Yanick Ricard

J Matas

F Chambat

Seismic attenuation in a phase change coexistence loop 1

Keywords:

Most phase transformations in the mantle occur across regions of multi-phase coexistence. Inside these regions, the long term incompressibility becomes very low because the density can increase both by compression and by changing phase. This difference between long term and elastic incompressibilities is a typical situation where seismic attenuation may happen. In this paper, we discuss the various differences between the classical theory of sound attenuation in a reacting fluid and the case of seismic propagation in a two-phase loop. We derive a simple analytical model of a two-phase loop to show that the phase change should affect both the bulk and the shear attenuation and in rather similar proportion. We show that attenuation occurs over two different frequency ranges. For the olivine-wadsleyite phase change, the low frequency attenuation occurs for periods larger than hundreds of years but the high frequency band occurs between 1 mn and 1 h (from 16 to 0.27 mHz) in the domain of surface waves and seismic modes.

We predict both bulk and shear quality factors between 1-10 in the middle of the 410 km phase loop.

A c c e p t e d M a n u s c r i p t uation. In spite of observational difficulties, several global models of 1D attenuation have however been published (e.g. [START_REF] Dziewonski | Preliminary reference earth model[END_REF]Widmer et al., 1991;[START_REF] Durek | A radial model of anelasticity consistent with long-period surface wave data[END_REF]. The disagreement between them is however large and often larger than the uncertainties suggested by each individual model (see e.g., Romanowicz and Mitchell, 2007, for a discussion). Three dimensional models are also available but are still a challenge and only the structures of the largest wavelengths have been mapped [START_REF] Gung | Q tomography of the upper mantle using threecomponent long-period waveforms[END_REF]. A better knowledge of attenuation is necessary to interpret the tomographic images and would however significantly improve our knowledge of mantle temperature [START_REF] Anderson | Absorption-band q model for the earth[END_REF][START_REF] Karato | Origin of lateral variation of seismic wave velocities 544 and density in the deep mantle[END_REF][START_REF] Matas | On the anelastic contribution to the temperature 556 dependence of lower mantle seismic velocities[END_REF][START_REF] Brodholt | Chemical versus thermal heterogeneity in the lower mantle: The most likely role of anelasticity[END_REF][START_REF] Lekic | Frequency dependence of 550 shear wave attenuation[END_REF].

In the last 40 years, (e.g., [START_REF] Jackson | Physical mechanisms of seismic-wave attenuation[END_REF][START_REF] Anderson | The earth as a seismic absorption band[END_REF][START_REF] Karato | Defect microdynamics in minerals and solid-state 541 mechanisms of seismic-wave attenuation and velocity dispersion in the mantle[END_REF], various attenuation mechanisms have been discussed including those due to phase changes. The attenuation of sound in a media undergoing a phase change is indeed a classical example of irreversible process that leads to attenuation [START_REF] De Groot | Non-Equilibrium Thermodynamics[END_REF]. Recently [START_REF] Li | Effect of phase transitions on compressional-wave veloci-552 ties in the earth's mantle[END_REF] have succeeded in the very difficult laboratory measurement of attenuation that takes place across the mantle transition zone, due to the presence of phase changes. Their paper assumes that the phase change attenuation is only related to compression, i.e., to what seismologists call the "bulk attenuation" (the quantity of energy lost during an oscillation of a pure isotropic compression). This attenuation is accounted for by the quality factor Q κ on which seismologists have very little resolution. Seismologists tend to ascribe most of the attenuation to the "shear attenuation" accounted for by the quality factor Q µ .

In the laboratory, the pressure changes used to drive the phase change are at the gigapascal level, while those due to seismic wave propagation are much smaller, typically of order 10 -7 GPa (e.g., [START_REF] Aki | Quantitative seismology[END_REF]. To rescale their observations, [START_REF] Li | Effect of phase transitions on compressional-wave veloci-552 ties in the earth's mantle[END_REF] propose a qualitative model of attenuation where the pressure perturbation δP associated with the seismic wave would drive the phase change at the interface between two grains by a distance d ∝ δP and this length should be compared to the time t necessary for cation diffusion, with t ∝ d 2 . Their model suggests therefore that the attenuation and the relaxation times are related to the amplitude of the seismic perturbation. This non-linearity would invalidate various assumptions of seismology, A c c e p t e d M a n u s c r i p t like the principle of linear superposition or the ability to describe the wavefield observation of an instrument as a series of convolutions. It would imply that seismic waves from large earthquakes see a more attenuating mantle than those from small ones.

The model of [START_REF] Li | Effect of phase transitions on compressional-wave veloci-552 ties in the earth's mantle[END_REF] is based on general considerations that are explained in more details in Jackson (2007). However, Jackson (2007) warns us that "no attempt has been made to model the time-dependent stress at the phase boundary or the transformation kinetics, potentially strongly influenced by the rheology of the surrounding medium". This is what we do in this paper where we propose a micromechanical model of a coexistence loop. We show that dissipation occurs in two different time periods and that the resulting attenuations are independent of the seismic wave amplitudes. We confirm that phase change loops may be the zones of large attenuations. We show that they should affect rather similarly the compressibility and the shear modulus.

Reaction rates of phase changes

The mechanism of attenuation due to a phase change in the mantle (or in a fluid) is easy to understand (de [START_REF] De Groot | Non-Equilibrium Thermodynamics[END_REF]. The changes of pressure due to a propagating sound wave affect differently the chemical potentials of the various coexisting phases and thus modify locally the thermodynamic equilibrium. This drives a mineralogical phase change, a possible source of dissipation. However, the theory of seismic attenuation in the mantle cannot be directly derived from that of sound attenuation in fluids because the physics differs by at least four aspects.

• First, the propagation of elastic waves is related to the rigidity µ (entirely for the S waves, and partly for the P waves) which is not considered for sound wave attenuation in fluids.

• Second, contrary to gases or fluids that are usually used in textbooks to illustrate thermodynamics, the thermodynamic equilibrium in complex solid aggregates is related to stresses, not to pressure. The pressure is not a continuous quantity across the grain interfaces. According to [START_REF] Shimizu | The non-equilibrium thermodynamics of intracrystalline diffusion under non-hydrostatic stress[END_REF], a chemical potential A c c e p t e d M a n u s c r i p t tensor should be defined on interfaces and the reaction rate should also depend on the crystal orientation. We assume here that the equilibrium on an interface only depends on the stresses normal to this interface σ n which is a continuous variable [START_REF] Paterson | Non hydrostatic thermodynamics and its geological applica-562 tions[END_REF]. In the absence of any viscoelastic stresses, the normal stress and the pressure can be identified and the usual thermodynamic rules are recovered.

• Third, the rheology of the mantle is not only simply elastic but viscoelastic.

Deviatoric stresses can relax for times larger than the Maxwell time of the viscoelastic mantle, the ratio of viscosity to rigidity (see e.g., [START_REF] Ricard | Physics of mantle convection[END_REF].

• Fourth, the phase transformations in the mantle are not univariant. As mantle materials are solid solutions and involves various cations, the phase changes occur across phase loops where two of more phases of various compositions coexist.

For example, around 410 km, an olivine with a Fe/Mg ratio typically of 1/10 [START_REF] Ringwood | Phase transformations and differentiation in subducted lithosphere: implications for mantle dynamics, basalt petrogenesis, and crustal evolution[END_REF], enters a phase loop where wadsleyite, with a larger Fe/Mg ratio, nucleates and then grows. This larger ratio is balanced by a symmetrical decrease of the Fe/Mg ratio in the remaining olivine. Across the phase loop the percentage of wadsleyite increases with depth and this wadsleyite has a decreasing Fe/Mg ratio until the ratio of 1/10 which corresponds to the disappearance of the last grains of olivine. Notice that in a phase loop, the two phases are already present and the nucleation of new grains should not control the kinetics of transformation, contrary to what may happen when a single phase moves through a phase transition [START_REF] Rubie | Kinetics of the olivine-spinel transformation in subducting lithosphere: experimental constraints and implications for deep slab processes[END_REF].

The sound propagation theory shows that the attenuation is ultimately related to the difference between the elastic parameters at very high frequency (the unrelaxed parameters) and those at very low frequency (the relaxed parameters) if the relaxation occurs within the period of the sound wave. The time dependent pressure variations δP (t) due a a high frequency seismic body wave, and the associated density variations δρ(t) are related by

δP (t) = κ ∞ δρ(t) ρ , (1)
A c c e p t e d M a n u s c r i p t

where κ ∞ is the elastic incompressibility (or elastic bulk modulus). The subscript ∞ indicates that this corresponds to the limit of infinite frequency, ω = +∞. More precisely, κ ∞ should be the isentropic incompressibility κ S but we will not distinguish in this paper between the isothermal and isentropic elastic incompressibilities, κ T and κ S , that are at any rate, equal within 1%.

At thermodynamic equilibrium, inside a multi-phase loop where the density jumps by ∆ρ over a depth range ∆P , depth dependent pressure variations and depth dependent density variations are roughly proportional and related by

dP dr = κ 0 ρ dρ dr with κ 0 = ρ ∆P ∆ρ . (2)
The equilibrium relation (2) defines the relaxed bulk modulus in the limit of zero frequency.

If we take the example of the phase change around 410 km depth, between olivine and wadsleyite, the unrelaxed incompressibility (elastic bulk modulus) is around κ ∞ = 180 GPa. With an average density of 3630 kg m -3 and a density jump of 180 kg m -3 over a thickness of 10 km (this value is reasonable although estimates ranging from 5 to 30 km have been proposed (see e.g., [START_REF] Shearer | Earth's Deep Interior: mineral physics and tomography from the atomic to the global scale[END_REF][START_REF] Van Der Meijde | Seismic evidence for water deep in earth's upper mantle[END_REF][START_REF] Ricard | Mineral physics in thermod-chemical mantle models[END_REF]), the relaxed incompressibility is κ 0 = 7 GPa. Outside a coexistence loop, the elastic incompressibility κ ∞ that can be measured by a time dependent phenomenon (the propagation of elastic waves) and the incompressibility measured along a radial profile κ 0 are usually considered as equal (or at least very close, see e.g., [START_REF] Bullen | The problem of the earth's density variation[END_REF].

The numerical expression of κ 0 in eq. ( 2) can be expressed in a more physical way that demonstrates that κ 0 is bounded by κ ∞ [START_REF] Li | Effect of phase transitions on compressional-wave veloci-552 ties in the earth's mantle[END_REF]. The density jump across a phase change ∆ρ is due both to an intrinsic density jump ∆ρ χ (the density difference between the two phases at a given pressure and temperature) and to the compression of the material across the coexistence loop. In other term, the relaxed compressibility in the phase loop, 1/κ 0 , is due both to the elastic compressibility 1/κ ∞ and to an apparent compressibility to the density jump ∆ρ χ existing between the two A c c e p t e d M a n u s c r i p t phases so that

1 κ 0 = 1 κ ∞ + ∆ρ χ ρ 1 ∆P > 1 κ ∞ . ( 3
)
131

When the transition thickness becomes very large, κ ∞ and κ 0 become therefore equal.

132

The evolution of an interface interacting with an elastic wave can be computed from

137 ρ α (v α -V).n = ρ β (v β -V).n = -∆Γ, (4)
138 where ρ i and v i are the densities and the velocities of each phase, n the normal to the 139 interface of the two media, directed from α to β, V the interface velocity and ∆Γ the 140 reaction rate of the α → β reaction (in kg m -2 s -1 ). The velocity jump across the

141 interface is therefore 142 (v α -v β ).n = - ∆ρ ρ α ρ β ∆Γ, (5) 143
where ∆ρ stands for ρ β -ρ α .

144

Although the expression of the reaction rate ∆Γ might be very complex, it must been chosen in previous modelling to be (e.g. [START_REF] Morris | Coupling of interface kinetics and transformation-induced strain 559 during pressure-induced solid-solid phase changes[END_REF][START_REF] Krien | The accomodation of volume changes in phase transition 546 zones: implications for mantle dynamics and metamorphism[END_REF])

∆Γ ∝ -δσ n , (6)
(the minus sign comes from the convention sign for normal stress, opposite to pressure, in fluid mechanics). For a very slow perturbation, however, the equation ( 6) cannot hold. On a time scale for which inter atomic diffusion occurs, the Fe/Mg content of each phase evolves, and as the chemical potentials are also functions of composition, a new equilibrium is found. The reaction occurs until the pressure change and the density change are related by the condition (2). This implies to choose

∆Γ ∝ -δσ n -κ 0 δρ ρ . (7)
This relation could be rigourously obtained by following the more formal derivation of de [START_REF] De Groot | Non-Equilibrium Thermodynamics[END_REF] provided the pressure replaces the normal stress. As κ 0 is very small, we will see that the difference in the attenuation predictions between using ( 6) or ( 7) is however only sensible at very long periods, outside the seismic frequency band.

We now need to define the geometrical distribution of the phases inside the loop and use ( 4) with ( 6) or ( 7) to be able to predict the effect of a seismic wave on the interface and thus on the attenuation. We assume that in each half of a phase loop, the minor phase is made of spherical grains surrounded by the major phase [START_REF] Morris | Coupling of interface kinetics and transformation-induced strain 559 during pressure-induced solid-solid phase changes[END_REF].

For example, in the shallower half of the 410 km depth transition, we consider the βphase as surrounded by the α matrix (see Figure 1). The outer radius R e represents the average distance between the grains with radii r ≤ R e of the minor phase. In the deepest part of the loop, the major β-phase is supposed to surround the last grains of α phase. This model will be used to describe the whole loop, although it is obvious that none of the two phases surrounds the other one in the middle of the loop.

When the normal stress on an interface is increased, a new film of the high pressure β-phase grows at the expense of the low pressure α-phase. This reaction, associated with minor changes of volume, facilitates the deformation, decreases the effective GPa, the thickness of the reacting film is only nanometric [START_REF] Li | Effect of phase transitions on compressional-wave veloci-552 ties in the earth's mantle[END_REF].

187

This very small perturbation is however a significant source of dissipation.

188

Radial deformation and complex incompressibility 189

We assume for simplicity that the two phases have the same elastic properties with 190 bulk modulus κ and rigidity µ. We neglect the difference of these parameters for the 191 two phases which is of the order of the jumps in incompressibility and rigidity, about 192 10% in PREM, at 410 km depth [START_REF] Dziewonski | Preliminary reference earth model[END_REF]. We introduce the an elastic shell of compressibility κ ∞ and rigidity µ can be written

198 u i r = a i s + b i s 2 , ( 8
)
199 with i = α in the outer shell and i = β, inside (b β is obviously 0 to insure that u β r is 200 finite at s = 0). The change of density is

201 δρ = -3 ρ u r (R e ) R e , (9)
202 and the effective compressibility of the medium is κ, 6) or ( 7), we can write

203 κ = ρ δP δρ = - R e 3u r (R e ) δP, (10)
v α -v β = C δσ n (S) -3κ 0 u(R e ) R e , ( 11
)
209
where C is a kinetic factor in m s -1 Pa -1 (the factor C includes the ∆ρ/ρ α ρ β term 210 of ( 5)). By choosing a vanishingly small κ 0 or the κ 0 deduced from the observed 211 thickness of the transition we will be in agreement with (6) or with (7).

219

We get

220 κ HF = κ ∞ + κ i -κ ∞ 1 + iωτ 1 , (12) 221
(HF stands for "high frequency" as it will be explained below) where

222 κ i = κ ∞ -S 3 (κ ∞ -κ 0 ) κ ∞ + 4µ ∞ /3 κ ∞ S 3 + 4µ ∞ /3 , ( 13
)
223
and where the relaxation constant τ 1 is

224 τ 1 = R e S 3Cκ ∞ κ ∞ + 4µ ∞ /3 κ ∞ S 3 + 4µ ∞ /3 . ( 14
)
225

The use of the subscripts i for "intermediate" will soon be explained.

236 µ * = µ ∞ iωτ M 1 + iωτ M , ( 15
)
237 where τ M is the Maxwell time, ratio of viscosity to rigidity (the real elastic rigidity).

238

If we redo the same modeling using µ * instead of µ ∞ we get the same expressions 239 as ( 12) and ( 13) but where µ * replaces µ ∞ . The new expression of the viscoelastic 240 incompressibility can be simplified and becomes after some algebra

241 κ = κ ∞ + κ i -κ ∞ 1 + iωτ 1 + κ 0 -κ i (1 + iωτ 1 )(1 + iωτ 2 ) , ( 16
)
242
with the long relaxation time τ 2

243 τ 2 = τ M κS 3 + 4µ/3 κS 3 . ( 17
)
244

We can safely assume, at least for the olivine-wadsleyite transformation (see e.g., Li

251 κ ≃ κ HF = κ ∞ + κ i -κ ∞ 1 + iωτ 1 , ( 18
)
252
where HF stands for high frequency and a second one for ωτ 2 ∼ 1 (and ωτ 1 << 1) in 253 which we get

254 κ ≃ κ LF = κ i + κ 0 -κ i 1 + iωτ 2 . ( 19
)
255 This is the low frequency approximation of κ. The low frequency limit of the high 256 frequency incompressibility, ω = 0 in (18), is κ i which of course is also the high

Pure shear deformation and complex rigidity

The existence of a phase change has also an effect on the rigidity. When at uniform pressure, the stresses are not uniform, the high pressure phase grows in the direction of the maximum stress and the reverse reaction occurs in the perpendicular direction.

This eases the deformation and therefore reduces the effective rigidity. Notice that it is only because the chemical potential in solids is related to the normal stress not to the pressure (uniform in a pure shear deformation) than reactions occur.

An analytical expression can be obtained, although the derivation is more cumbersome and less rigourous than in the spherical case. Let us consider the deformation of a nucleus of β-phase surrounded by a shell of α-phase when a pure shear deformation is applied to the external boundary. The pure shear deformation far from the central nucleus is in cartesian coordinates u z = -γz, u x = γx/2 and u y = γy/2 (γ is the strain)

and can be written in spherical coordinates after a standard change of coordinates

u r = - 1 2 γr(3 cos 2 θ -1), and u θ = 3 2 γr cos θ sin θ, (20)
where θ is the colatitude.

In a pure shear experiment performed in a laboratory, the rigidity would be the ratio between the vertical stress applied on the surface of a core sample, at position z, and the vertical strain measured at the same position -zσ zz (z)/2u z (z). In our analytical model, instead of imposing the deformation on surfaces of constant cartesian coordinates, they are imposed on the sphere or radius R e . We consider that the effective rigidity can however be estimated by

µ = - σ rr (R e , θ))R e 2u r (R e , θ) . (21)
We can solve for the deformation inside the two-phase aggregate by assuming that it keeps the same degree 2 geometry. As the radial stress and the radial deformation have the same geometry, µ is independent of θ. In this case, the general solution of the A c c e p t e d M a n u s c r i p t momentum balance yields

u i r = a i s + b i s 3 + c i s 2 + d i s 4 (3 cos 2 θ -1) (22)
314 and 315

u i θ = -3a i s -b i 7µ + 5λ ∞ λ ∞ s 3 -6 µ ∞ 3λ ∞ + 5µ c i s 2 + 2 d i s 4 cos θ sin θ ( 23
)
316
where s is again the normalized radius

s = r/R e , λ ∞ = K ∞ -2µ ∞ /3 is an elastic 317
Lamé parameter, and i stands for α or β.

318

Using ( 22) and ( 23), the stress tensor can be computed in the inner sphere and in

327 µ HF = µ ∞ + µ i -µ ∞ 1 + iωτ 3 . ( 24
)
328 where the intermediate rigidity µ i is

329 µ i = µ ∞ -µ ∞ F (S, κ ∞ µ ∞ ), ( 25
)
330
where F is a cumbersome function of S and κ ∞ /µ ∞ . We can choose for simplicity and the relaxation times τ 3 is

331 κ ∞ = 5µ ∞ /3
F (S, κ ∞ µ ∞ = 5 3 ) = 105S 3 8 + S 2 604 + 280S 3 + 56S 5 + 5S 7 , ( 26
)
τ 3 = 315 SR e Cµ ∞ 1 604 + 280S 3 + 56S 5 + 5S 7 . ( 27
)
The constant τ 3 of shear modulus relaxation is also depicted in Figure 2 and its value is comparable to that appearing in the high frequency incompressibility τ 2 .

Like for the incompressibility, the intermediate rigidity µ i that corresponds to the relaxed limit of the rigidity with respect to the short time scale of phase change can in turn be seen as an unrelaxed rigidity for periods larger than τ 3 but much smaller than the Maxwell time. Replacing directly in ( 24) and ( 25), µ ∞ by µ * , see ( 15), is straightforward but leads to a very complex expression of µ. We checked numerically that the behavior of µ * F (κ ∞ /µ * ) can be qualitatively approximated by µ * F (κ ∞ /µ ∞ ). For periods much larger than τ 3 , the rigidity varies therefore as

µ = µ ∞ + µ i -µ ∞ 1 + iωτ 3 iωτ M 1 + iωτ M . ( 28
)
As τ 3 << τ M , the expression of the rigidity at high frequency, when ωτ M >> 1 is µ ≃ µ HF . At low frequency, when ωτ 3 << 1 we get

µ ≃ µ LF = µ i - µ i 1 + iωτ M , ( 29
)
which reaches zero after complete relaxation when ω = 0, i.e., when the medium behaves viscously rather than elastically.

The real and imaginary parts of the rigidity are plotted in Figure 3. Like for the incompressibility, two transitions are predicted. The low frequency transition only occurs when the elastic stresses in the surrounding shell relax and stop screening the inner nucleus from the outside stresses. At zero frequency (infinite period), the rigidity vanishes contrary to the incompressibility that remains finite.

A c c e p t e d M a n u s c r i p t

Bulk and shear attenuations

In the case of complex elastic parameters, the seismic waves propagates with a 359 frequency-dependent attenuation usually defined by the quality factors Q ξ where ξ 360 stands for κ or µ. These quality factors are defined by

361 Q ξ = - Re(ξ) Im(ξ) = ξ R + ξ U ω 2 T 2 R (ξ U -ξ R )ωT R , ( 30
)
362
where T R is the appropriate relaxation time, Re(ξ) and Im(ξ) the real and imaginary 363 parts of ξ that varies between the relaxed and unrelaxed limits, ξ R (for ω = 0) and ξ U

364

(for ω = +∞). Attenuation is maximum at the frequency

365 ω 0 = 1 T R ξ R ξ U , ( 31
)
366
where Q ξ reaches its minimum 367

Q 0 ξ = 2 ξ R ξ U ξ U -ξ R . ( 32
)
368

Notice that in our model, each elastic parameter has low and high frequency modes.

369

The bulk modulus can relax from ξ U = κ to ξ R = κ i with the time constant τ 1 , then 370 from ξ U = κ i to ξ R = κ 0 with the time constant τ 2 . Similarly, the rigidity relaxes 371 from ξ U = µ to ξ R = µ i with the time constant τ 3 , then from ξ U = µ i to ξ R = 0 with viscously rather than elastically. The general behavior for the shear attenuation is that 394 of a linear solid called a Burger body (see e.g. [START_REF] Karato | Defect microdynamics in minerals and solid-state 541 mechanisms of seismic-wave attenuation and velocity dispersion in the mantle[END_REF]. The bulk 395 attenuation behaves differently as contrary to the Burger body, the quality factor does 396 vanish at infinite periods.

397

The fact that the bulk quality factor is lower than the shear quality factor is not a 398 general result of our model. The difference κ i -κ ∞ that controls the bulk attenuation is 399 proportional to κ 0 -κ ∞ (see ( 12) and ( 13)). This difference tend to decrease when the To avoid unnecessary complexities, we assume λ ∞ = µ ∞ or κ ∞ = 5µ ∞ /3. This time, the boundary conditions on the outer shell imply

a α + b α + c α + d α = - γ 2 , (B-3) 1 4 (-12a α -48b α -3c α + 8d α ) = 3γ 2 , (B-4)
the continuity of radial and shear stress at the phase-change interface implies

(8a α -4b α s 2 -19 c α s 3 -32 d α s 5 ) = (8a β -4b β s 2 ), (B-5) -(24a α + 120b α s 2 + 15 c α s 3 + 64 d α s 5 ) = -(24a β + 120b β s 2 ), (B-6)
and the continuity of θ-displacement

1 4 (-12a α s -48b α s 3 -3 c α s 2 + 8 d α s 4 ) = 1 4 (-12a β s -48b β s 3 ). (B-7)
As there is no global volume change in this pure shear experiment the jump of radial displacement is given by ( 6) Figure 3: Incompressibility and rigidity as functions of the period. The real part is depicted by lines, the imaginary part by shadows. The thin line corresponds to the incompressibility computed using (6), the thick line using (7). In the former case the incompressibility goes to zero at long period while in the latter, it reaches the relaxed compressibility κ 0 . The high frequency relaxation occurs for seismic and tidal periods.

iω (a α s + b α s 3 + c α s 2 + d α s 4 ) -(a β s + b β s 3 ) = C 1 4 µ ∞ R e ( 8a
Page 27 of 30 The elastic incompressibility κ, seen by high frequency body waves, the totally relaxed compressibility κ 0 that can be deduced from the density jump and the thickness of the transition, and an intermediate incompressibility κ i , see equation ( 13). The transition between κ and κ i and an associated attenuation should occur for periods corresponding to surface waves. The transition between κ i and κ 0 should take a much longer time (a Maxwell time larger than 100 yrs).

133

  the mechanical properties of the two phases and the boundary conditions on the inter-134 face. At the interface between grains, the total normal stress (pressure plus deviatoric 135 stress) and shear stress are continuous. The boundary condition for the velocity across 136 the reacting interface is expressed by

145

  cancel when the two phases are at thermodynamic equilibrium. The definition of the 146 thermodynamic equilibrium is however subtle in the two-phase loop, and we discuss 147 here two possible equilibrium conditions. We show in the following that the choice of 148 one or the other expression does not change our conclusions on the seismic attenuation 149 within the two-phase loop.150When the material inside a phase loop is at equilibrium, the reaction rate is zero. 151 When the system is perturbed, for example by a change of the far-field pressure, 152 the normal stress on interfaces changes. The rules of irreversible thermodynamics 153 (de Groot and Mazur, 1984) suggests that the reaction rate close to equilibrium is pro-154 portional to the distance to equilibrium, i.e. to the change of affinity of the reaction. 155 The associated changes of chemical potentials are initially related to the normal 156 stress perturbations on the grain interfaces δσ n . Therefore the reaction rate has often

  the elastic energy and therefore leads to attenuation. Notice that 185 with the low pressure changes involved during the propagation of a seismic front, ∼ 186 10 -7

  193normalized radius s = r/R e where R e is the external radius and we use S for the value 194 of s on the two-phase interface.195In the case of isotropic radial compression of the two medium of Figure1, it is 196 straightforward to show that with spherical symmetry the radial deformation u r (s) in 197

  204where u r (R e ) and δP are the radial displacement and the pressure perturbation at the 205 outer radius R e .206 By using the general expression of the jump of normal velocity (5) and one of the 207 A c c e p t e d M a n u s c r i p t kinetic laws (

212

  If a sinusoidal pressure perturbation of frequency ω, δP exp(iωt), is applied on 213 the external rim, the deformations and therefore the coefficients a i and b i , will also 214 vary at the same frequency. From the general expression of the deformation, (8), the 215 normal stress can be expressed (see Appendix A). By matching the normal stresses 216 across the phase boundary at s and using the normal velocity jump condition (11) with 217 v α = iωu α r and v β = iωu β r , the three constants a α , b α and a β can be found from which 218 the effective incompressibility (10) is readily obtained (see Appendix A).

245

  and[START_REF] Li | Effect of phase transitions on compressional-wave veloci-552 ties in the earth's mantle[END_REF], that τ 1 << τ 2 ∼ τ M , i.e., that the reaction occurs in a time 246 shorter than the Maxwell time of a few hundred years. As is physically expected, the 247 high frequency limit, when 1 << ωτ 1 << ωτ 2 is the elastic value κ = κ ∞ and the 248 low frequency limit ωτ 1 << ωτ 2 << 1 is the incompressibility deduced from the 249 thickness of the phase loop, κ = κ 0 . The incompressibility variations occur within two 250 frequency bands, one for ωτ 1 ∼ 1 (and thus ωτ 2 >> 1) in which we recover (12),

Figure 3

 3 Figure 3 depicts the evolution of the real and imaginary part of the incompressibil-

  319 the outer shell (see Appendix B). The final resolution involves the determination of 6 320 parameters a α , b α , c α , d α , a β and b β . They can be obtained by matching 6 boundary 321 conditions; the continuity of shear stress, normal stresses and tangential deformation 322 on the interface, the jump condition for the normal velocity on the interface, and the 2 323 external boundary conditions (20) (see Appendix B). Notice that no density variation 324 occurs in the assemblage submitted to a pure shear and the kinetic law (6) is therefore 325 appropriate. Similarly to what we obtain for the incompressibility, the effective rigidity 326 deduced from the model can be written as

372 the time constant τ M . 373 Figure 5

 3735 Figure 5 depicts the minimum quality factor Q 0 κ and Q 0 µ for the various relaxation

Figure 6

 6 Figure 6 depicts the attenuation quality factor Q κ and Q µ computed at the center

400

  loop thickness increases (see (3)). On the contrary the difference µ i -µ ∞ controlling 401 the shear attenuation is independent of the loop thickness (see(24)). The ratio Q κ /Qµ 402 decreases therefore with the loop thickness. Increasing the loop thickness over 50 km 403 for the 410 km transition (keeping the other parameters unchanged) would lead to Q µ 404 lower than Q κ . In other words, a thin loop is mostly attenuating because of its Q κ , a 405 thick one because of its Q µ . attenuation in coexistence loop is certainly simplified in particular in 408 the description of the geometry of the two phases. However we believe that various 409 aspects of our model are very robust.410We confirm that the phase loops, in agreement with Li and Weidner (2008of a significant attenuation for periods belonging to the low frequency 412 seismic band (surface waves and seismic modes). However in their paper, only the bulk 413 attenuation was taken into account whereas we have shown that the shear attenuation 414 should be affected as well. Moreover, we found that the two attenuations and their 415 relaxation times are independent of the amplitude of the seismic wave in agreement 416 with the usual assumption of seismology. The attenuation bands are narrow because 417 of our assumption of a unique grainsize across the coexistence loop. A distribution of 418 grainsizes would probably broaden the frequency bands of attenuation.419Another large difference between[START_REF] Li | Effect of phase transitions on compressional-wave veloci-552 ties in the earth's mantle[END_REF] and the present paper is 420 that the transition between the elastic incompressibility, κ ∞ and that deduced from 421 PREM κ 0 (very low incompressibility in coexistence loops) does not occur in a single 422 step but in two steps. This is summarized in Figure7where we plot as a function of 423 depth the instantaneous κ ∞ , the intermediate κ i and the κ 0 bulk moduli. We derived 424 a κ ∞ (r) that follows the PREM values outside the loop and varies linearly within the 425 loop. We introduce the same κ ∞ (r) in the computation of κ i , (13). These depth-426 dependent incompressibilities look more realistic, but truly our model assumes the κ ∞ 427 is uniform within the loop. We predict that the transformation is first limited by the 428 shielding of the external stresses by the elastic rigidity of the matrix. This behavior 429 is in agreement with[START_REF] Morris | Coupling of interface kinetics and transformation-induced strain 559 during pressure-induced solid-solid phase changes[END_REF]. This first "fast" transition between κ ∞ and κ i 430 (see Figure7) occurs within minutes or hours. Therefore the incompressibility seen 431 by seismology is somewhere between these two values. The final transformation (slow 432 kinetics) corresponding to the transition between κ i and κ 0 occurs in a second step 433 when the elastic stresses are released, and a new two-phase equilibrium is found after 434 diffusion and reequilibration of the composition. A very similar figure could be drawn 435 for the rigidity, but the totally relaxed rigidity would simply be zero.436The minimum values of Q κ and Q µ predicted for the olivine-wadsleyite loop are 437 significantly lower than the typical range of attenuation of low frequency seismologic 438 models. This is true for the shear attenuation and even more for the bulk attenuation. 439 In a recent review, Resovsky et al. (2005) propose for depths around 400 km a Q µ 440 roughly between 150 and 200 and Q κ larger than 2000. These radial models have used 441 a very simple parametrization of the attenuation in a limited range of layers (typically 442 A c c e p t e d M a n u s c r i p t larger than 200 km) and in layers than coincide with the seismic discontinuities. This 443 parametrization is certainly an inappropriate choice to detect a narrow depth range of 444 attenuation astride a velocity discontinuity. 445 The Q κ and Q µ quality factors seen by given seismic or tidal perturbations are 446 not necessary the minimum values shown in Figures 5 or 6, but are functions of their 447 frequencies, of the kinetic constant C and of the grain size R e . The exact frequency de-448 pendent expressions for incompressibility, rigidity and attenuation can be easily com-449 puted from (12), (24) and (30). The values of C and R e are still uncertain for the 450 olivine-wadsleyite phase change and mostly unknown for other mantle phase changes.451The narrow depth range of the coexistence loop, ∼ 10 km, may make this localized 452 attenuation difficult to detect but may significantly bias the average determined over 453 a large zone of sensibility. The presence of bulk attenuation in the upper mantle has 454 indeed be proposed by various studies (Resovsky et al., 2005; Durek and Ekstrom, a strong attenuation associated with seismic velocity jumps should 457 also affect the reflection and transmission factors of short period seismic waves hitting 458 these interfaces. Notice however that body waves seem to have a too high frequency to 459 be right on the attenuation maximum, but according to Figure 6 might see an attenua-460 tion of order 10-100 (for frequencies between 1 Hz to 20 mHz). This could possibly be 461 seen from S waves multiply reflected under the Earth surface and with a turning point 462 below or above the 410 km discontinuity. At any rate, a better estimate of attenuation 463 in phase loops would constrain the kinetic behavior of these phase changes (Chambat 464 et al., 2009). 465 We have discussed the attenuation due to the olivine-wadsleyite transition because 466 this transition is simple and only involves two phases of similar composition. Other 467 mantle transitions should behave similarly like the wadsleyite-ringwoodite loop around 468 520 km deep (but, of course, with their appropriate kinetic laws that may not lead 469 to characteristic times for the phase changes in tune with the seismic periods). The 470 mantle transitions involving three phase like the ringwoodite to ferropericlase+perov-471 skite, should also lead to attenuation because of a large difference between the relaxed 472 and unrelaxed properties. However because of the necessary long distance motion of

Figure 1 :

 1 Figure1: The growing β-phase is surrounded by the low pressure α-phase. The interface at radius r moves up or down depending of the applied external pressure that controls the reaction rate ∆Γ. The external diameter 2Re will be interpreted as the average distance between grains.

Figure 2 :

 2 Figure2: Relaxation times of incompressibility (τ 1 and τ 2 ) and shear modulus (τ 3 and Maxwell time τ M ) across a phase loop. The two short relaxation times are in the seismic band of surface waves, for the bulk attenuation τ 1 and the shear attenuation τ 3 . The long bulk and shear relaxation times τ 2 and τ M correspond to viscoelastic behaviours occurring after at least one Maxwell time. The horizontal lines correspond to periods of 1 h and 1 yr (frequencies of 0.27 mHz and 0.03 µHz).

Figure 4 :Figure 6 :Figure 7 :

 467 Figure4: To compute the shear attenuation, we submit the two-phase to a pure shear experiment. The pressure remains constant, but the high pressure β-phase starts growing in the direction of the maximum stress, the α-phase in the direction of the lowest stress as the reaction rate is related to the normal stress on interfaces.

  elastic matrix even at infinite time. To clarify this point we have to remember that at 233 large times, the Earth mantle does not behave as elastic but as viscous. The rigidity in 234 the definition of κ i , is in fact the high frequency limit of a viscoelastic rheology (see

	235
		e.g., Ricard, 2007)
		M a n u s c r i p t
	230 A c c e p t e d
	231	correspond to the incompressibility κ 0 (see (2)) obtained from a radial seismological
	232	model. This is because in a purely elastic model, the deviatoric stresses remain in the

The incompress-226 ibility at infinite frequency which is usually called the unrelaxed incompressibility is 227 simply the elastic incompressibility κ ∞ , while κ i is the relaxed incompressibility ob-228 tained for ω = 0. 229 It might be surprising that the relaxed incompressibility κ i in this model does not

  which corresponds to a Poisson ratio of 1/4 (or to the equality of the

332

two Lamé Parameters λ ∞ and µ ∞ , (e.g.

[START_REF] Malvern | Introduction to the mechanics of a continuum medium[END_REF]

) which is a common rule of 333 thumb for elasticity of silicates. In this case, the function F is 334

Acknowledgments

This work has been supported by the CNRS-INSU Dyeti program. We thank Tine Thomas for interesting discussions on seismic attenuation.

Appendix A. Radial deformation

The deformation given by (8) corresponds to the radial stress

The continuity of normal stress across the phase transition boundary satisfies

while the application of an external pressure corresponds to

The condition of jump of normal displacement where the phase change occurs is

The resolution of equations (A-2)-(A-4) leads to the unknowns a α , b α and a β and therefore to the expression of the effective incompressibility (see ( 10))

δP.

(A-5)

Appendix B. Pure shear deformation

From expressions ( 22)-( 23), stresses can then be readily obtained,

where i stands for α or β. The high frequency quality factor for both Q 0 κ and Q 0 µ are below 10 in the center of the loop. The exact value of κ 0 (7 10 9 GPa, thick solid line, or 0, thin solid line) is not very important. The very low quality factor of Q 0 κ at low frequency does not belong to the seismic or tidal domain.


M Wolstencroft
email: wolstencroftm@cardiff.ac.uk

J H Davies

D R Davies

Nusselt-Rayleigh number scaling for spherical shell Earth mantle simulation up to a Rayleigh number of 10 9

Keywords: Planetary Mantle, Rayleigh Number, Nusselt Number, Thermal Evolution, Spherical Geometry

An investigation of the power law relationship between Nusselt number (Nu) and Rayleigh number (Ra) for Earth's convecting mantle is presented. The Nu(Ra) relationship was calculated from the results of a model with three dimensional spherical geometry and free slip boundary conditions. Both basally and internally heated convection has been examined.

For Nu(Ra) = aRa β , β was found to be 0.294 ± 0.004 for basally heated systems, which is lower than the value of 1 3 suggested by conventional boundary layer theory. The exponent β = 0.337±0.009 for internally heated systems, when the internally heated Ra is converted to a basally heated equivalent for comparison. The influence of the method used to calculate β was also considered, with particular attention paid to high Ra. As an example of the significance of β = 0.29 rather than 1 3 , a Ra of 10 9 results in a surface heat flux which is ≈ 32% lower. Within the range of Ra used in this study, there is no evidence that β changes at high Ra. Therefore, that mechanism cannot be used to moderate mantle temperature when projecting back to early Earth conditions. The differing planform of basally and internally heated models was shown to result in different scaling relationships between root mean square surface velocity and Ra for the two modes of heating, in particular, a much lower surface velocity for internally heated cases relative to equivalent Ra basally heated cases.

Introduction

The aim of geophysical research is to understand Earth and its processes over the entire history of the Planet on all scales. The endeavour is hindered by incomplete knowledge of some of the basic mechanisms of planetary development. The thermal evolution of Earth is one such area. Earth is, to a large extent, a heat engine. Assumptions regarding the evolution of temperature over time appear in most, if not all, conceptual models of plate tectonic theory, geochemistry and geophysics.

A favoured method for investigating Earth's thermal evolution has been one dimensional (1D) parameterised modelling of mantle convection (for example: [START_REF] Turcotte | Finite amplitude convective cells and continental drift[END_REF], the mantle comprising by far the largest portion of the Earth by volume (≈85%). A number of valuable insights have been gained into convective processes [START_REF] Sharpe | Parameterized mantle convection and the Earth's thermal history[END_REF][START_REF] Schubert | Subsolidus convective cooling histories of terrestrial planets[END_REF][START_REF] Davies | Thermal histories of convective Earth models and constraints on radiogenic heat production in the Earth[END_REF][START_REF] Turcotte | On the thermal evolution of the Earth[END_REF][START_REF] Stevenson | Magnetism and thermal evolution of the terrestrial planets[END_REF][START_REF] Honda | A simple parameterized model of Earth's thermal history with the transition from layered to whole mantle convection[END_REF][START_REF] Choblet | 3D thermal convection with variable viscosity: can transient cooling be described by a quasi-static scaling law?[END_REF][START_REF] Korenaga | Onset of convection with temperature and depth-dependent viscosity[END_REF][START_REF] Korenaga | Energetics of mantle convection and the fate of fossil heat[END_REF][START_REF] Nimmo | The influence of potassium on core and geodynamo evolution[END_REF]. Even when much more complex mantle convection simulations can be undertaken in 3D spherical geometry [START_REF] Tackley | Effects of multiple phase transitions A c c e p t e d M a n u s c r i p t in three-dimensional spherical model of convection in Earth's mantle[END_REF][START_REF] Bunge | A sensitivity study of three-dimensional spherical mantle convection at 10 8 Rayleigh number: Effects of depth-dependent viscosity, heating mode and an endothermic phase change[END_REF][START_REF] Zhong | Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection[END_REF][START_REF] Davies | Steady plumes produced by downwellings in Earth-like vigor spherical whole mantle convection models[END_REF], there is still value in parameterised models. They (i) provide a convenient way to run a large number of simulations, (ii) allow very high Rayleigh numbers and (iii) permit simulations of thermal evolution over almost all of Earth history.

A critical aspect of these parameterised models is the assumed relationship between the ratio of convective to conductive heat transfer, characterised by the non-dimensional Nusselt number (Nu):

Nu = q q K (1)
Page 3 of 40

A c c e p t e d M a n u s c r i p t

where: q is the heat transferred by convection, while q K = k∆T /D, is the amount of heat that would be conducted through a layer of thickness D with a temperature difference ∆T across it, k being the thermal conductivity.

The other important parameter is the convective vigour, characterised by the Rayleigh Number (Ra):

Ra = gρα∆T D 3 κµ (2)
where: g is the acceleration due to gravity, ρ is density, α is the coefficient of thermal expansion, ∆T is superadiabatic temperature drop across the shell, D is domain thickness (2900 km), κ is thermal diffusivity and µ is dynamic viscosity. Both equations are presented in their most common forms after [START_REF] Davies | Dynamic Earth : plates, plumes, and mantle convection[END_REF].

A series of models with sub Earth-like parameters may be used to derive the Nusselt-Rayleigh number relationship:

Nu(Ra) = aRa β ( 3
)
where a is a constant and β, the index of the power law relation, is the main value of interest.

Nu(Ra) can then be extrapolated to present day and early Earth conditions (higher Ra).

These values can be compared with heat flux and temperature measurements to constrain the thermal budget of Earth.

The value of Ra (Equation 2), is often quoted to demonstrate how Earth-like a particular model is in terms of its convective vigour. Earth's mantle is considered to have a basally

Page 4 of 40

A c c e p t e d M a n u s c r i p t heated Ra in excess of 3 × 10 6 [START_REF] Davies | Dynamic Earth : plates, plumes, and mantle convection[END_REF] and probably in the range 10 8 [START_REF] Bunge | A sensitivity study of three-dimensional spherical mantle convection at 10 8 Rayleigh number: Effects of depth-dependent viscosity, heating mode and an endothermic phase change[END_REF][START_REF] Weeraratne | Transitions in the style of mantle convection at high Rayleigh numbers[END_REF]. The assumed Ra of Earth is based on the Ra obtained when Earth-like parameters are used in Equation 2. The accuracy of these values is open to interpretation, leading to the uncertainty as to the exact Ra of Earth's mantle.

Historically, β was derived from boundary layer theory in a unit square box (aspect ratio 1) and found to be 1 3 [START_REF] Turcotte | Finite amplitude convective cells and continental drift[END_REF]. [START_REF] Christensen | Heat transport by variable viscosity convection and implications for the Earth's thermal evolution[END_REF] suggested that β may be closer to 0.1 for systems with a rigid lid at the upper boundary. However, [START_REF] Gurnis | A reassessment of the heat transport by variable viscosity convection with plates and lids[END_REF] argued that the presence of plate tectonics on the Earth allowed much greater cooling efficiency, similar to that of a free-slip upper boundary. In his more Earth-like models, β was generally around 0.3. More recently, Solomatov and colleagues [START_REF] Solomatov | Scaling of temperature-and stress-dependent viscosity convection[END_REF][START_REF] Reese | Stagnant lid convection in a spherical shell[END_REF] have derived the Nu(Ra) relationship for the stagnant-lid and sluggish-lid regimes. This paper follows [START_REF] Gurnis | A reassessment of the heat transport by variable viscosity convection with plates and lids[END_REF], by considering mobile surface models as being most representative of Earth. Stagnant and sluggish lid regimes might be more meaningful for other terrestrial bodies and therefore any relationships derived for Earth should not be applied to other bodies without first considering which regime best describes their mode of convection. [START_REF] Korenaga | Onset of convection with temperature and depth-dependent viscosity[END_REF] and [START_REF] Korenaga | Energetics of mantle convection and the fate of fossil heat[END_REF] confirmed the value of around 0.3 for models of the onset of convection with temperature dependent viscosity, implying that although simpler models lack some of the more subtle features of the real mantle, the Nu(Ra) scaling relationship is robust. [START_REF] Stevenson | Styles of mantle convection and their influence on planetary evolution[END_REF] also noted that β is rarely observed to be exactly 1 3 . Experimental work, such as that of [START_REF] Castaing | Scaling of hard thermal turbulance in Rayleigh-Benard convection[END_REF] and some low vigour modelling work [START_REF] Bercovici | Three-dimensional thermal convection in a spherical shell[END_REF][START_REF] Bercovici | Three-dimensional convection of an infinite-Prandtlnumber compressible fluid in a basally heated spherical shell[END_REF] have produced results where β = 2 7 or ≈ 0.28. This has led to the suggestion that there is a transition in the behaviour of a thermally convecting system at Ra beyond 10 7 , with increasing convective vigour becoming less efficient at increasing the thermal throughput of the system. While such a change in β is an intriguing possibility, neither of the studies operated entirely in the physical domain of mantle convection on Earth. [START_REF] Bercovici | Three-dimensional thermal convection in a spherical shell[END_REF][START_REF] Bercovici | Three-dimensional convection of an infinite-Prandtlnumber compressible fluid in a basally heated spherical shell[END_REF] modelled significantly sub-Earth Ra, while [START_REF] Castaing | Scaling of hard thermal turbulance in Rayleigh-Benard convection[END_REF] Previous modelling work to evaluate Nu(Ra) has concentrated on 1D, 2D and 3D Cartesian box domains, with only a small number of studies employing Earth-scale 3D spherical geometry [START_REF] Bercovici | Three-dimensional thermal convection in a spherical shell[END_REF][START_REF] Bercovici | Three-dimensional convection of an infinite-Prandtlnumber compressible fluid in a basally heated spherical shell[END_REF][START_REF] Iwase | An interpretation of the Nusselt-Rayleigh number relationship for convection in a spherical shell[END_REF]. Although models that do not employ 3D spherical geometry can provide valuable insights into mantle mechanisms, they are often limited by imposed 'side' boundary conditions and non Earth-like aspect ratios.

Many studies carried out in the 1980's and early 1990's lacked sufficient computing resources to model high degrees of structural complexity, with many using significantly lower than Earth-like Ra and limiting convection to low spherical harmonic degree or prescribed initial conditions to improve computational efficiency (e.g. [START_REF] Bercovici | Three-dimensional thermal convection in a spherical shell[END_REF][START_REF] Bercovici | Three-dimensional convection of an infinite-Prandtlnumber compressible fluid in a basally heated spherical shell[END_REF].

This study advances the work of constraining the relationship between Nu and Ra by using 3D spherical geometry and utilising computationally efficient multi resolution solution algorithms, which allow the examination of more Earth-like Ra [START_REF] Davies | Applying multi-resolution numerical methods to geodynamics[END_REF]. The relationship can then be characterised to higher Ra, giving greater confidence in the resulting Nu(Ra) relationship and limiting the degree of extrapolation when considering Earth's mantle. Very long model runs have been undertaken, allowing the use of short-scale random initial conditions. This removes the need for specific initial conditions. Since it is known that more than one large scale planform can be stable at lower Ra [START_REF] Bercovici | Three-dimensional thermal convection in a spherical shell[END_REF], our approach avoids the need to make a choice a priori regarding the convective pattern.

Simulation Method

The most important controlling parameter for this study is the Ra of the modelled man- . Previous work has used the approach of simulating a range of low Ra and using the derived scaling to extrapolate the results to higher Ra. This study also utilises this approach, however, recent innovations applied to the TERRA code [START_REF] Davies | Applying multi-resolution numerical methods to geodynamics[END_REF] and the availability of large computing resources allows modelling up to Ra bh = 10 8 (basal heating)

and Ra ih = 10 9 (internal heating). Consequently, the degree of extrapolation is marginal.

The TERRA code utilised here [START_REF] Baumgardner | Three dimensional treatment of convective flow in the Earths mantle[END_REF][START_REF] Bunge | Mantle convection modeling on parallel virtual machines[END_REF][START_REF] Yang | Matrix-dependent transfer multigrid method for strongly variable viscosity infinite Prandtl number thermal convection[END_REF] is a benchmarked, 3-D spherical finite element mantle dynamics code. In order to simulate very high Ra cases a novel extension of TERRA was used. [START_REF] Davies | Applying multi-resolution numerical methods to geodynamics[END_REF] applied a multigrid refinement approach, which allows a priori selective radial refinement of model resolution. Thus, it was possible to model the mantle with a resolution of ≈ 14 km at upper and lower boundary layers, thereby resolving thermal boundary layers at a basally heated Ra bh = 10 8 and an internally heated Ra ih = 1.44 × 10 9 . Lower resolution can be used for the sake of efficiency as lower Ra result in thicker thermal boundary layers, which are resolvable at lower spatial resolution.

The dynamical problem solved by TERRA can be defined in terms of conservation of mass

∇ • v = 0 (4) momentum 1 ρ ∇P = ν∇ 2 v + αg∆T a r (5)
and energy

Page 7 of 40

A c c e p t e d M a n u s c r i p t

∂T ∂t + v • ∇T = κ∇ 2 T + J ρC p (6)
where

129 κ = k ρC p (7)
v is velocity, P is dynamic pressure, ν is kinematic viscosity, t is time, T is temperature, ∆T a is adiabatic temperature drop, J is rate of internal heat generation per unit volume, C p is specific heat at constant pressure and r is the inward directed unit radial vector.

A Newtonian viscosity constitutive equation was assumed [START_REF] Bunge | Mantle convection modeling on parallel virtual machines[END_REF].

In this work, we assume the Boussinesq approximation, where density differences are neglected except in the buoyancy term of the momentum equation. The Prandtl number is assumed to be infinite.

The equations are solved dimensionally in TERRA but can be non-dimensionalised using the following relations (after [START_REF] Davies | Physical model of source region of subduction zone volcanics[END_REF]:

x ′ = x D , t ′ = κt D 2 , T ′ = T ∆T and P ′ = D 2 P µκ (8)
where 'primes' indicate non-dimensional terms and x is distance.

From these one obtains:

v ′ = vD κ and Q ′ = J ρC p D 2 κ∆T ( 9
)
Page 8 of 40 where Q ′ is non-dimensional internal heating and v ′ is non-dimensional velocity. Dropping the primes, Equations 4, 5 and 6 can then be written as:

∇ • v = 0 (10) ∇P = ∇ 2 v + RaT r (11) ∂T ∂t + v • ∇T = ∇ 2 T + Q (12)
Thus it becomes evident that the Ra is the key controlling non-dimensional parameter.

The technical details of TERRA and its solvers are extensively documented [START_REF] Baumgardner | Three dimensional treatment of convective flow in the Earths mantle[END_REF][START_REF] Bunge | Mantle convection modeling on parallel virtual machines[END_REF][START_REF] Yang | Matrix-dependent transfer multigrid method for strongly variable viscosity infinite Prandtl number thermal convection[END_REF] only details directly relevant to this work are described here. The computational domain (i.e. the spherical shell) is discretised by means of a regular icosahedron, where the 20 triangular sides (paired as 10 diamonds) are repeatedly divided into 4 sub-diamonds until the desired resolution is achieved. The degree of refinement is referred to as 'mt' and increases in powers of 2 i.e. 16, [START_REF]Table 5: Resulting Nu for Ra=10 9 with different values of a for β = 1 3 . Note that for β = 0.29, a remains 0[END_REF]64,128,256 and so on. The grid is extended radially by placing several of these spherical shells above one another, generating a mesh of triangular prisms (layers) with spherical ends.

In this study, nr, the number of radial layers, is set to:

nr = mt 2 (13)
The following assumptions were used: isochemical convection, free slip upper and lower boundaries, no mineral phase changes, incompressible rheology, whole mantle convection, and constant viscosity. For basally heated cases, both boundaries were isothermal, whilst for internally heated cases the upper boundary was isothermal and the lower boundary insulating. The manner in which Ra and Nu are calculated depends on the style of heating.

For a model where all heat enters through the base of the shell, the Ra is calculated as, and the Nu (after [START_REF] Bunge | A sensitivity study of three-dimensional spherical mantle convection at 10 8 Rayleigh number: Effects of depth-dependent viscosity, heating mode and an endothermic phase change[END_REF] is

Nu bh = h s D 4πk(r s r b )∆T (15)
As internally heated cases have no fixed lower boundary temperature, different methods for calculating Ra and Nu must be adopted. Internally heated Ra is

Ra ih = αρ 2 gHD 5 µκk ( 16
)
while internally heated Nu (after [START_REF] Reese | Scaling laws for time-dependent stagnant lid convection in a spherical shell[END_REF] is

Nu ih = Hρr 2 s (1 -η 2 (3 -2η)) 6k(T av -T s ) (17)
where: η = inner to outer shell radius ratio (0.546... for the Earth), H is the rate of internal heat generation (W/kg), r s is the outer radius of the shell (6370 km), r b is the inner radius of the shell (3480 km), T av is the average temperature over all radial layers, T s is the surface temperature, and h s is the surface heat flux.

While the concept of Nu and Ra remains the same, previously one was able to use ∆T , now the rate of internal heating combined with the shell thickness and thermal conductivity

Page 10 of 40

A c c e p t e d M a n u s c r i p t is used to produce a ∆T equivalent, based on heat flux across the domain [START_REF] Davies | Dynamic Earth : plates, plumes, and mantle convection[END_REF].

Comparing equations 14 and 16 it is evident that the internally heated Ra will be 1-2 orders of magnitude larger than a basally heated Ra from an otherwise similar model. An internally heated Ra can be converted into its basally heated equivalent [START_REF] Davies | Dynamic Earth : plates, plumes, and mantle convection[END_REF]:

Ra bh i = Ra ih Nu ih (18)
In this work, Ra bh and Ra bh i are used unless otherwise stated.

Input Parameters

While the vigour of convection is controlled by the Ra, TERRA solves the equations dimensionally. For completeness, the actual input values that remain unchanged from case to case are presented in Table 1. In this study, Ra was set by altering the dynamic viscosity of the model. This produced Ra bh in the range 10 3 to 10 8 and Ra ih from 10 4 to 10 9 . The initial conditions for each case were identical and based on low amplitude, short wavelength random temperature variations. Such an initial condition inevitably extends the time taken for a given model to evolve but has the advantage of not imposing an initial large-scale planform on the model.

Table 1 here.

The cases investigated in this study are outlined in Table 2. Eight cases were basally heated and seven cases were internally heated. Minimum viscosity was ≈ 10 20 P a•s (near Earth-like).

Table 2 here.

Model runs were allowed to evolve until a statistically stable condition was achieved. This

Page 11 of 40

A c c e p t e d M a n u s c r i p t differs from true steady state in that the models are intrinsically unsteady at high Rayleigh numbers, a feature which is likely to be shared by the real mantle. In this study, a simulation is considered to have reached 'steady state' if there is no strong long-term trend in the surface heat flux as a function of time.

An important question is whether a case has been sufficiently well resolved. As an approximate measure, the number of radial model layers within the thermal boundary layer should be ≥ 5 (after [START_REF] Lowman | Steady plumes in viscously stratified, vigorously convecting, three-dimensional numerical mantle convection models with mobile plates[END_REF]. It is this requirement that dictates which resolution a case is most efficiently modelled at.

Results

Degree of Model Evolution

Model cases were allowed to evolve such that they were no longer influenced by their initial condition. Figs. 1 and 2 demonstrate the evolution of each case.

Figures 1 and2 here.

The most apparent difference between basally and internally heated cases is that internally heated models converge towards the same heat flux, whereas basally heated models each find their own level. This is because heat output must equal heat generation at 'steady state' in the internally heated cases, otherwise the model is in a transient heating or cooling state. Basally heated models, by contrast, need to balance their heat input at the base and loss at the surface. With the viscosity controlling the thickness of the thermal boundary layer, each basally heated case has a different heat flux through the mantle.

Relationship of Nusselt Number to Rayleigh Number

A potential source of error is that the Nu will be time dependent in cases that are not truly steady state, a feature of most realistic mantle models. It can be seen from Table 3 here.

The Nu(Ra) results are presented in Table 3 and Figs. 3 and4. The mean Nu values were used to produce a power law fit by simple linear regression (Fig. 3 and4). This fit, while good, can alternatively be generated in what might be considered a more rigorous manner, using a weighted least squares method. In the weighted least squares method the weights were taken to be the variance in the mean of the Nu, which were greater for the higher Ra runs due to their temporal variability. It is important to note that the value for β produced by these two fits is very similar for basally heated cases but significantly different for internally heated cases. The data appear to be generally well fit by a single slope within range of Ra investigated, with the possible exception of the lowest Ra cases, which appear to lie on a slightly lower gradient for basally heated cases. This is discussed further in Section 4.

Figures 3 and4 here.

Scaling laws relating Nu and Ra were found to be:

Nu = 0.284Ra 0.294±0.004 bh ( 19
)
for basally heated cases (Fig. 3) and:

Nu = 0.164Ra 0.337±0.009 bh i (

Surface velocity -Ra scaling

The Nu(Ra) relationship for basally and internally heated cases have a similar value of a. In contrast, when the root mean square (V RM S ) surface velocity relationship (in ms -1 ) to the Ra is examined (Fig. 5) the resulting values of a differ more markedly, with a given Ra producing a significantly lower velocity for the internally heated cases than for the basally heated cases.

Figure 5 here.

The respective relationships are:

V RM S bh = 10 -13 Ra 0.538 (22)

V RM S ih = 7 × 10 -14 Ra 0.500 (23)
In general, the V RM S results display a much greater scatter than Nu(Ra). The cause of this is the greater time dependence of the V RM S when compared to the surface heat flux, the key variable for the Nu. Therefore, at high Ra, V RM S is more unsteady than surface heat flux. As a result, interpretations of the V RM S scaling relationship should be considered more approximate than Nu(Ra).

Spectral Analysis

In order to better understand how the convective wavelength of the models varies with Ra, spectral plots of spherical harmonic power were generated for a subset of cases (Figure 6). As would be expected, the power shifts progressively to higher harmonic degree (l) as

Page 14 of 40

A c c e p t e d M a n u s c r i p t

Ra increases and convective length scale reduces. For basally heated cases, the CMB and surface are generally symmetrical in terms of where the power peaks are located. At the highest basally heated Ra there appear to be two peaks at l = 6 and l = 10, indicating that different scales of flow may be co-existing. Internally heated cases demonstrate higher l and show an increasing asymmetry between surface and CMB as Ra increases. The surface generally demonstrates a higher l than the CMB. In comparison with basally heated cases, internally heated cases have power over a wide range of l and depth with many smaller power peaks suggesting a broad range of convective scale. Internally heated cases with high Ra show more power at the surface over a range of l from ≈ 9 -28.

Figure 6 here.

Discussion

There are a number of factors, particularly relating to how the Ra and Nu numbers are calculated, that need to be considered so this work can be seen in the correct context. Ra and Nu are calculated from average values. For example, surface heat flux will not be the same at all points globally. Therefore, as the Nu is calculated from model output, the nature of the model, be it 1D, 2D or 3D will have a bearing on the final value. One must also consider that not all Rayleigh numbers quoted in the literature represent the same measure. Ra for internally and basally heated models are, by necessity, calculated differently (Equations 14 and 16 and cannot be directly compared (unless the internal heating Ra is divided by the Nu, as is done here).

As a further example, the work of [START_REF] Christensen | Thermal evolution models for the Earth[END_REF], involving temperature dependent viscosity, defined Ra using the viscosity at average temperature, which produces somewhat lower values for β when using basal heating and higher values when using internal heating.

The advantage of Christensen's method for calculating the Ra is that it allows the relationship to be defined for temperature dependent viscosity convection. The disadvantage is that comparisons with other work are not straightforward. Christensen obtained a value

Page 15 of 40 A c c e p t e d M a n u s c r i p t of β=0.2 when using his definition of viscosity as opposed to the value of 0.28 when using the more commonly found definition of viscosity, the average viscosity over the entire domain.

The models presented here assumed free slip upper and lower boundaries (after [START_REF] Gurnis | A reassessment of the heat transport by variable viscosity convection with plates and lids[END_REF]. While free slip is a good approximation of the very large viscosity contrast at the CMB, it must be remembered that Earth's surface boundary layer is composed of rigid plates. [START_REF] Gurnis | A reassessment of the heat transport by variable viscosity convection with plates and lids[END_REF] argues that a stress free surface boundary condition is a better approximation to Earth than a rigid lid. In addition to the plate boundary condition, Earth also has phase transitions at 410 km, 660 km and D ′′ [START_REF] Murakami | Post-perovskite phase transition in MgSiO3[END_REF][START_REF] Oganov | Theoretical and experimental evidence for a post-perovskite phase of MgSiO3 in Earths D" layer[END_REF]. These features combined with other aspects such as: temperature and depth dependent viscosity, non uniform internal heating and heterogeneous composition may impose some limits on the applicability of the values for β presented here. It is also possible that grain size dependent viscosity in the diffusion creep regime is important [START_REF] Solomatov | Grain size-dependent viscosity convection and the thermal evolution of the Earth[END_REF]. If mantle minerals behave in a similar manner to laboratory analogues, the viscosity of the lower mantle could be influenced by factors affecting grain growth rate. As this study confines itself to thermal effects, the influence of such a mechanism was not investigated.

However, if such a mechanism proved important, it could exert a similar influence on thermal evolution as varying β.

Comparable previous studies (e.g. [START_REF] Korenaga | Energetics of mantle convection and the fate of fossil heat[END_REF] range of values for β, depending on their modelling and measuring methods, but have either concluded that 1 3 is a point which their models 'will' trend to at higher Ra (e.g. [START_REF] Iwase | An interpretation of the Nusselt-Rayleigh number relationship for convection in a spherical shell[END_REF], or that the value of β is entirely different due to particular boundary conditions [START_REF] Christensen | Heat transport by variable viscosity convection and implications for the Earth's thermal evolution[END_REF]. [START_REF] Iwase | An interpretation of the Nusselt-Rayleigh number relationship for convection in a spherical shell[END_REF] used 3D spherical harmonic models to obtain a value closer to 1 4 , very close to the β ≈ 0.24 of the uncorrected internally heated cases presented here. They argued that at higher Rayleigh numbers their β would tend to The power law scaling of RMS surface velocity and Ra shows notable differences between the basally and internally heated cases (Fig. 5), with internally heated cases producing lower velocities for any given Ra. Fig. 7 demonstrates how the planform of the two heating modes is different. In particular, the characteristic lateral length of individual convection 'cells' for internally heated cases is approximately half that of basally heated cases. In a basally heated case, the organising presence of the upper and lower thermal boundary layers produces larger scale features with correspondingly larger buoyancy anomalies. For an internally heated case, thermal anomalies are smaller and have lower buoyancy anomalies associated with them. Therefore, in an internally heated case, convection 'cells' are smaller than in an equivalent basally heated case (Fig. 7). This results in basally heated cases displaying higher velocities than internally heated cases for a given Ra.

The link between the above observation and the Nu(Ra) scaling differences (β = 0.29 for basal heating, 0.33 for internal heating) lies in the way in which heat is input and passes through the system. A basally heated system must transfer a given 'parcel' of heat from the bottom to the top of the mantle, whereas with internal heating, heat is input throughout the volume and the majority of heat 'parcels' have less distance to travel before reaching the surface. This effect appears to scale better with increasing Ra than it does for basally heated cases, leading to a higher value for β for internally heated cases. That this β is close to the 1 3 predicted by boundary layer theory, may be the result of the aspect ratio of the convection cells, as opposed to that of the domain, which remains constant. The generally lower aspect ratios of the internally heated convection cells appear to be closer to the 1:1 aspect ratio assumed by boundary layer theory and shown to give β ≈ 0.33 in 2D modelling [START_REF] Olsen | A comparison of heat transfer laws for mantle convection at very high Rayleigh numbers[END_REF].

Analysis is an important and often under emphasised-aspect of modelling. Many mod- A weighted fit which discards the lowest 3 Ra cases (Fig. 8) provides a much closer fit to the data and agrees more closely to the unweighted value for β (0.290 unweighted to 0.309 weighted). The cause of this improvement lies in the fact that the 3 cases with the lowest Ra for basal heating appear to lie on a slightly lower trend than the rest of the data points (see Fig. 3). This is a consequence of the large scale cubic/tetrahedral convective patterns that low Ra cases tend to adopt (Fig 9). These patterns appear to be less efficient (in a heat transfer sense) than the more complex patterns found at higher Ra. This has been observed previously for low Ra spherical geometry cases [START_REF] Bercovici | Three-dimensional thermal convection in a spherical shell[END_REF][START_REF] Bercovici | Three-dimensional convection of an infinite-Prandtlnumber compressible fluid in a basally heated spherical shell[END_REF].

Table 4 shows a selection of previously published values for β, calculated from a variety of model types.

Table 4 here. [START_REF] Bercovici | Three-dimensional thermal convection in a spherical shell[END_REF][START_REF] Bercovici | Three-dimensional convection of an infinite-Prandtlnumber compressible fluid in a basally heated spherical shell[END_REF] vigour, the convective pattern has a significant influence on β. The models used in this study explored much larger Ra, however the work of [START_REF] Bercovici | Three-dimensional thermal convection in a spherical shell[END_REF][START_REF] Bercovici | Three-dimensional convection of an infinite-Prandtlnumber compressible fluid in a basally heated spherical shell[END_REF] remains valid for low vigour convection.

Work by [START_REF] Giannandrea | Variable viscosity convection experiments with a stress-free upper boundary and implications for the heat transport in the Earth's mantle[END_REF] produced a value of 0.28 in large aspect ratio tank experiments. The aspect ratio of convection cells (λ) seems to have an impact in some geometries [START_REF] Olsen | A comparison of heat transfer laws for mantle convection at very high Rayleigh numbers[END_REF]. At larger values of λ, β seems to tend to 0.25, while at lower λ, β is closer to 0.33. It is likely that this effect was responsible for the value of 0.28 found by [START_REF] Giannandrea | Variable viscosity convection experiments with a stress-free upper boundary and implications for the heat transport in the Earth's mantle[END_REF], with Rayleigh numbers up to 10 5 . An arbitrary choice of aspect ratio is not an issue with 3D spherical geometry models, at least where whole mantle convection is considered.

The work of [START_REF] Korenaga | Energetics of mantle convection and the fate of fossil heat[END_REF][START_REF] Korenaga | Archean geodynamics and the thermal evolution of Earth[END_REF] examines the effect of chemically controlled viscosity in the mantle. Korenaga proposes that the dehydration of the upper mantle caused by the melting generating ocean crust causes the upper mantle to become more viscous. Logically, higher mantle temperatures in the past would be expected to result in more melting and a greater stiffening effect. Thus, Korenaga proposes that a hotter Archean mantle would in fact have a lower effective Ra than has been assumed. This effect would have implications for the value of β as it would have to be adjusted to accommodate this stiffening effect. [START_REF] Korenaga | Energetics of mantle convection and the fate of fossil heat[END_REF] suggests that β may be very low or even negative when this effect is taken into account. Without this effect included the value for β produced by Korenaga's models is 0.3, generally in agreement with other recent convection models and the values presented here.

It is worth noting that, while [START_REF] Korenaga | Energetics of mantle convection and the fate of fossil heat[END_REF][START_REF] Korenaga | Archean geodynamics and the thermal evolution of Earth[END_REF] uses a value for β derived from a uniform viscosity model, he then uses the scaling relationship in parameterised models where viscosity is temperature dependent. This is a slight clash of concept as the β being used was not explicitly defined for such a case. In practice, this does not appear to introduce significant error as β appears to be robust in such a situation [START_REF] Korenaga | Energetics of mantle convection and the fate of fossil heat[END_REF]Korenaga Page 20 of 40 , 2002). This again highlights the need for β to be very well defined. [START_REF] Zhong | Dynamics of thermal plumes in three-dimensional isoviscous thermal convection[END_REF] used CITCOM, a 3D Cartesian box model with 100% basal heating and the same definition of basally heated Ra as this study found that plume number (the global number of plumes transporting a significant amount of heat, see [START_REF] Zhong | Dynamics of thermal plumes in three-dimensional isoviscous thermal convection[END_REF] for details)

scales at Ra 0.31 . The work presented here agrees broadly with [START_REF] Zhong | Dynamics of thermal plumes in three-dimensional isoviscous thermal convection[END_REF] as large scale plumes, originating at the lower thermal boundary layer, are likely to be the primary transporters of heat, thus controlling surface heat flux in a 100% basally heated system. Internal heating by contrast, tends to produce a similar heat flux but over a broader area. Plume number would not be expected to scale in the same manner for internally heated cases.

An additional factor that may affect β is the type of convection. [START_REF] Castaing | Scaling of hard thermal turbulance in Rayleigh-Benard convection[END_REF] divide convection into soft and hard turbulent regimes, with the transition from soft to hard at Ra = 4 × 10 7 . Their experiments to extremely high vigour (Ra = 10 12 ) using liquid Helium, give β = 2 7 (0.2857...). While this type of experiment is obviously not intended to exactly simulate mantle convection, it suggests that convecting systems are capable of displaying behaviour at odds with boundary layer theory. The value of β derived by this study (Equation 19) is closer to 0.29 than 1 3 , indicating that a certain amount of non-classic boundary layer behaviour may be occurring. The broad range of Ra, reaching Earth-like vigour, lend weight to the argument that this is a valid result. The tendency of any change in behaviour at high Ra would be to reduce the value of β for Earth. The mechanism suggested to cause this reduction is the entrainment of parts of the lower boundary layer, possibly as fingers of hot material [START_REF] Castaing | Scaling of hard thermal turbulance in Rayleigh-Benard convection[END_REF]. This study finds no evidence of a change in β at high Ra, as illustrated in Figs. [START_REF] Moore | Heat transport in a convecting layer heated from within and below[END_REF] proposes that boundary layer theory applies to 'ambient' convection but does not account for plumes. The 'overshoot' of plumes thins the boundary layer, operating in a similar manner to the additional advection described by [START_REF] Jarvis | Convection Models and Geophysical Observations: Mantle Convection, Plate Tectonics and Global Dynamics[END_REF]. This effect is greater with larger plumes and cases in the low Ra range would have a higher surface heat flux than boundary layer theory would predict. For higher Ra cases, where plumes are smaller, this effect is reduced.

In the context of either Jarvis and Peltier or Moore, Nu(Ra) derived from a wide range of Ra will have a value of β < 1 3 . In effect, this additional heat flux represents the interference of the lower boundary layer on the upper, something which is not included in boundary layer theory assumptions. A hint of such flow may be present in Fig. 6 at higher Ra, where there is more than one peak in the power spectra. Care must be taken when comparing this study directly to the work of Moore as his conclusions are based on mixed mode heating, where the upper and lower boundary layers are not of equal thickness. General observations on boundary layer theory should however be comparable.

Other ideas also involve the disruption of the upper thermal boundary layer by plumes, a dramatic example of which is provided by the degree-1 convection found by [START_REF] Zhong | Degree-1 mantle convection and the crustal dichotomy on Mars[END_REF] A feature of Fig. 3 is that if the trend is extended to Nu = 1 (the onset of convection), our predicted critical Ra would be of order 10 2 . The critical Ra has been extensively investigated and is accepted to be ≈ 712 [START_REF] Bercovici | Three-dimensional thermal convection in a spherical shell[END_REF]. This apparent inconsistency in our results is explained by the fact that at Nu = 1 (critical Ra) there is effectively no convection and therefore no boundary layers (or the whole layer can be considered a single boundary layer), boundary layer theory of convection is no longer valid and conduction is dominant. [START_REF] Moore | Heat transport in a convecting layer heated from within and below[END_REF] suggests that at near critical Ra a single slope does not fit because conduction interferes significantly with convection.

There is also the question of whether the data are indeed best fit by a single slope (β).

At very low Ra, basally heated cases seem to lie on a slightly flatter trend (lower β). This is caused by the presence of near steady state convective planforms similar to those observed by [START_REF] Bercovici | Three-dimensional thermal convection in a spherical shell[END_REF]. Apart from these very low Ra cases, results from the range of Ra modelled in this study do appear to be fit by a single slope (e.g. Fig. 3). Based on 2D Cartesian models, [START_REF] Moore | Heat transport in a convecting layer heated from within and below[END_REF] suggests that at Ra above those investigated in this work, the data are best fit by several slopes. In particular, Moore cites [START_REF] Lenardic | Thermal convection below a conducting lid of variable extent: Heat flow scalings and two-dimensional, infinite Prandtl number numerical simulations[END_REF], who suggested that β = 1 3 is recovered above Ra ≈ 10 9 . Unfortunately, their use of 2-D modelling with a rigid lid surface boundary condition limits the applicability of the result to the data presented here. No studies have been carried out in 3D spherical geometry at Ra > 2 × 10 9 , leaving open the question of whether β = 1 3 can be recovered a extremely high Ra. Some laboratory convection experiments to high Ra in plane layers such as those of [START_REF] Castaing | Scaling of hard thermal turbulance in Rayleigh-Benard convection[END_REF] have been able to fit a single slope to their data over a wide range of Ra. It is interesting to note that the basally heated Ra of 10 9 -10 11 , suggested by Moore as the point at which β = 1 3 , is probably much higher than the likely Ra over most of Earth's history.

The implication of β = 0.29 as opposed to 1 3 is that the extrapolated surface heat flux

22

A c c e p t e d M a n u s c r i p t in the past (at higher Ra) would be lower than otherwise estimated. This lower heat flux would arise from the efficiency of the convective system, before other factors, such as rheology, are taken into account.

To demonstrate this, consider the relationship Nu = aRa β where a is 0.284 (from Equation 19). When using an upper bound present-day Earth Ra bh of 10 9 , a value of β = 0.29 results in Nu = 116 (Surface Heat Flux (q) of 114 TW). The value of q was calculated from the Nu using the following equation re-arranged and simplified from equation 15:

q = Nu × 9.832 × 10 11 W (24)
To compare this to β = 1 3 another value for a must be chosen to prevent the two trends from being pinned at the origin, leading to an artificially large Nu at high Ra. In order to obtain reasonable projection, a was calculated using Nu and Ra values for the mid-point of the data on Fig. 3 and also for the lowest value Ra in the data. These two values pin the projected values to points within the data. The three possibilities for projecting the power law are shown in sketch form in Fig. 10.

Figure 10 here.

A comparison of Nu values for β = 1 3 and β = 0.29 at Ra = 10 9 is shown in Table 5.

The favoured values are highlighted bold and show a 32% reduction in the expected Nu when β = 0.29 rather than 1 3 . The difference in q would be 55 TW. These results are particularly significant when one considers that this is a fairly conservative choice of a when generating the β = 1 3 Nu(Ra) relationship. Clearly, when extrapolating to higher Ra, the difference between the outcomes predicted by the two values of β would be greater still.

Page 24 of 40 It can be seen that the Earth's ancient heat flux could be lower than conventional boundary layer theory would suggest for a relatively small change in β. This perhaps offers a mechanism that could partially assist in resolving the thermal catastrophe paradox which occurs when projecting the Earth's temperature into the distant geological past as described by [START_REF] Korenaga | Energetics of mantle convection and the fate of fossil heat[END_REF]. However, the reduction of β is not sufficiently large to solve the problem by itself. Further mechanisms, possibly only significant in early Earth, are required to reconcile the present day mantle heat flux.

Conclusion

The Nu(Ra) scaling relationship has been investigated in 3D spherical geometry over a large range of convective vigour, using both internally and basally heated models. Basally and internally heated systems were directly compared by converting the internally heated Ra to a basally heated equivalent. This was achieved by dividing the internally heated Ra by the Nu (Equation 18). It was found that β is closer to 0.29 than

  used liquid Helium, a fluid with a low Prandtl number.Thus, investigations at infinite Prandtl number and high Ra are required to test the possi-

  Figs. 1 Page 12 of 40 A c c e p t e d M a n u s c r i p t and 2 that the surface heat flux of a statistically stable model tends to fluctuate about its average value. For this reason the Nusselt numbers were calculated by averaging over the last 1000 timesteps of each model at 10 timestep intervals. The last 1000 timesteps of each run are higlighted in Figs. 1 and 2.

  20) for internally heated cases (Fig. 4 using Ra bh i ) when simple linear regression was used to fit the data. The value of β obtained by weighted least squares fit was not used further in this Page 13 of 40 A c c e p t e d M a n u s c r i p t study as the variance of residuals was larger than for the unweighed fit. For completeness, the uncorrected internally heated Ra gives the following relationship: Nu = 0.297Ra 0.24 ih (21)

  the weighted line fit is notably different. Previous studies have found a

1 3 ,

 3 an assumption somewhat at odds with high Ra experimental work such as[START_REF] Castaing | Scaling of hard thermal turbulance in Rayleigh-Benard convection[END_REF] and indeed with this study. A significant difference in the work of[START_REF] Iwase | An interpretation of the Nusselt-Rayleigh number relationship for convection in a spherical shell[END_REF] is their use of local Nu and Ra (calculated across the thermal boundary layer) as opposed to global Nu and Ra as is traditionally the case. Therefore, their value for β is based on somewhat different assumptions to the classical case and may be more sensitive to variations in planform, possibly similar to the differences in Rayleigh-RMS velocity scaling of different heating modes observed in this work.

Figure 7

 7 Figure 7 here.

  Figures 8 and 9 here.

  have proposed a value for β of around 0.28. Their work, in 3D spherical geometry, demonstrates tetrahedral and cubic patterns of convection at relatively low Ra. They note differing values of β for the different planforms, with β for their cubic cases given as 0.28, while tetrahedral solutions give 0.26. This indicates that at low

  1 and 2. The mechanism reducing β relative to the classical theory operates consistently over the range of Ra studied. The work presented here does not support β reducing with Ra as a solution to the thermal catastrophe when modelling thermal evolution backwards from the present. The question of why β derived from experiment and simulation is not completely consistent Page 21 of 40 A c c e p t e d M a n u s c r i p t with boundary layer theory has been addressed by a number of studies. What many have in common, is the idea that the simplifications underlying boundary layer theory are more suited to some Ra than others. Jarvis and Peltier (1989) demonstrated that the assumption in boundary layer theory of zero vertical velocity along the top of the unit box is more valid for higher Ra than lower Ra. Lower Ra cases have a greater degree of advection of heat in the upwelling corner of the box, which enhances heat flux out of the surface. Higher Ra cases result in less vertical heat advection, more closely resembling boundary layer theory.

  when considering Martian convection. Degree-1 convection was not observed Page 22 of 40 A c c e p t e d M a n u s c r i p t in any of the cases presented here.

1 3 (Figure 1 .

 31 Figure captions

Figure 2 .

 2 Figure 2. Summary of model evolution for internally heated cases, boxes indicate the time interval over which the average Nu was calculated. Logarithmic axes.

Figure 3 .

 3 Figure 3. Two possible fits to the results from basally heated cases. Trends of 1 3 and 2 7 are also shown. A value of 1 3 being the value predicted by boundary layer theory, and 2 7 being the value suggested by Castaing et al. (1989) from their high Ra, but low Prandtl number experiments.

Figure 4 . 2 PageFigure 5 .

 425 Figure 4. Two possible fits for output from the internally heated cases. Trends of 1 3 and 2

Figure 6 .

 6 Figure 6. Spectral plots showing power at a given spherical harmonic degree (l) for the temperature field up to l = 32. l increases up each plot, the CMB is on the left and surface on the right. Row (a) consists of basally heated cases and row (b) internally heated cases. Ra increases across each row from left to right.

Figure 7 .

 7 Figure 7. Comparison of temperature field planform at 22 km depth for cases with a dynamic viscosity of 1 × 10 22 Pas. (a) basally heated model and (b) internally heated model. Colour figure can be found in the electronic version of this paper.

Figure 8 .

 8 Figure 8. An improved fit may be achieved when the lowest 3 Ra cases are excluded from the analysis. Basally heated cases only.

Figure 9 .

 9 Figure 9. Case B7, 900 km depth. Temperature plot in Kelvin. Note tetrahedral planform. Colour figure can be found in the electronic version of this paper.

Figure 10 .

 10 Figure10. Sketch representation of how the chosen value for a effects the Nu when extrapolating to higher Ra (see text for details).

  Figure 1:

  Figure 2:

  Figure 3:

  Figure 5:

  Figure 7:

  have shown that simplified models are capable of capturing the essence of the more complex real Earth. This study has charac-

	terised Nu(Ra) using a model that is more Earth-like, in terms of convective vigour and
	geometry, than the majority of previous work.
	Nu(Ra) derived by boundary layer theory of thermal convection gives a β value of 1 3 (Tur-
	cotte and Oxburgh, 1967). It is therefore important to consider why our model derived

value for β is closer to 0.29. The value of β=0.337 for internally heated cases appears to agree, though we note that classical boundary layer theory assumes no internal heating. The error for β is however somewhat higher for the internally heated value than for the basally Page 16 of 40

Table 5 here .

 5.

Table 1 :

 1 Common input values used for all cases. Note: the governing equations can be expressed nondimensionally, these values are presented for completeness.

	Page 33 of 40

are shown for comparison.

Acknowledgements

The authors would like to thank John Baumgardner for help and support with TERRA. Models were run on supercomputing clusters at Cardiff (Helix), Liverpool (NESSC) Universities and HECToR (the UK National Supercomputer). Special thanks go to John Brodholt for enabling access to HECToR as part of the NERC Mineral Physics Consortium. The authors would like to thank Chris Reese and Mark Jellinek for their thoughtful and detailed reviews. MW was funded by NERC e-Science studentship: NER/S/A/2005/13131.