Marie Běhounková
email: marie.behounkova@univ-nantes.fr

Gaël Choblet
email: gael.choblet@univ-nantes.fr

Onset of convection in a basally heated spherical shell, 1 application to planets 2

Keywords: onset of convection, scaling laws, temperature-dependent 7 viscosity, shell geometry, sluggish rid regime 8

Convective instabilities related to the early dynamics of planetary mantles just after core formation play an important role in the subsequent evolution. Although these early stages of planetary dynamics are likely to imply more complex phenomena such as global melting and fractional solidification, and although density variations of compositional origin are likely to play an important role, little is known about the onset of solid-state convection in a fluid with temperature-dependent viscosity heated from below. Here, we investigate onset times of convection in order to obtain scaling relationships for the influences of Rayleigh number, viscosity parameter describing the dependency on the temperature and geometry of spherical shell (measured by f , ratio between the

inner and outer radii). We performed three dimensional numerical experiments and we concentrate on the dynamical regime described by global viscosity contrast smaller than 10 4 . Onset times and wavelengths of the first instabilities using both dynamical (free-slip) and kinematical (no-slip) boundary conditions are investigated. For both boundary conditions, the scaling may be written in the form t ′ ∝ (Ra * ) a , where a is approximately -2 /3 and Ra * = Ra(µ(θ * )) is a Rayleigh number specifically associated with a relevant temperature (viscosity) value (θ * ≈ 0.25). In addition, the dimensionless onset times (using the shell thickness as a characteristic length scale) are almost independent on the geometry of the shell for large range of the geometrical factor (f 0.2). In order to better understand these processes, 3D results are compared with two simple methods: the linear stability (LS) analysis and the growth of Rayleigh-A c c e p t e d M a n u s c r i p t

Taylor (R-T) instabilities. The LS analysis values of the onset times are much smaller due to the "frozen time" approach (i.e. the conductive propagation of the hot front is not taken into account). The dependency of the onset time on the Rayleigh number is overestimated, especially for the free-slip conditions, where the "frozen time" effect is even more significant. For the R-T instability analysis, although the onset times are also underestimated, the agreement with 3D simulations is good in terms of efficient scaling relationships. When applied to the dimensions and plausible initial state of terrestrial planets (Mars, the Earth and Venus), the scaling relationships provide an idealized framework to investigate early dynamics. Due to uncertainties associated with the "initial" temperature field and viscosity parameters, the computed onset times vary by several orders of magnitude (between 0.1 Myr and 500 Myr). These are likely to be smaller than the ones obtained for the onset of convection at the base of the lithosphere. For the investigated range of parameters, the minimal preferred degree for the onset instabilities is estimated to be approximately 10 so that, other ingredients or a different dynamical regime, have to be considered to promote the very low degree convective instabilities suggested for the early evolution of Mars.

Introduction 9

The internal distribution of temperature and of chemical composition dur-10 ing the early stages of planetary evolution are shaped by a series of primordial 11 energetic events such as possible large scale impacts and core differentiation.

12 Although in a purely thermal framework, the primordial nature of the temper-13 ature field within the planet will be forgotten due to the strong temperature-14 dependence of the viscosity (cf. [START_REF] Tozer | Heat transfer and convection currents[END_REF], especially in the case of a hot start 15 (e.g. [START_REF] Schubert | Whole planet cooling and the ra-1082 diogenic heat source contents of the Earth and Moon[END_REF], variations in composition may strongly affect this 16 simple scheme. One possible origin of chemical stratification may result from A c c e p t e d M a n u s c r i p t the solidification of a global magma ocean that likely resulted from the energetic processes mentioned above. Despite the fact that a precise timing may be controversial [START_REF] Wood | Cooling of the Earth and core formation after giant impact[END_REF][START_REF] Allègre | The major differentiation of the Earth at 4.45 Ga[END_REF], the existence of such an early magma ocean is supported by geochemical measurements on certain isotopes of tungsten and neodymium (cf. [START_REF] Kleine | Rapid accretion and early core formation on asteroids and the terrestrial planets from Hf-W chronometry[END_REF][START_REF] Blichert-Toft | Short-lived chemical heterogeneities in the archean mantle with implications for mantle convection[END_REF]. Although the solidification of the magma ocean is a complex phenomenon involving many processes (e.g. [START_REF] Solomatov | Fluid dynamics of a terrestrial magma ocean[END_REF] some of which still lack a full understanding, it has been suggested that the density distribution subsequent to the crystallization may be highly unstable and would result in a large scale overturn (see for Mars, [START_REF] Elkins-Tanton | Magma ocean fractional crystallization and cumulate overturn in terrestrial planets: Implications for Mars[END_REF], 2005).

The new stratification could then be stable enough to prevent the mixing of this deep denser layer possibly strongly enriched in heat producing elements (cf. [START_REF] Elkins-Tanton | Early magnetic field and magmatic activity on Mars from magma ocean cumulate overturn[END_REF]. A stratification of possibly different origin, is also suggested in some models of the present day dynamics of the Earth's mantle in order to reconcile geochemical and geophysical observations. The dynamical implications are studied both in the laboratory (e.g. Le [START_REF] Bars | Whole layer convection in a heterogeneous planetary mantle[END_REF] and numerically (e.g. [START_REF] Tackley | Strong heterogeneity caused by deep mantle layering[END_REF]. Since the above dynamical processes are complicated and because they are controlled by poorly constrained factors, simple idealized models lead to a useful insight into the onset of solidstate convection. Indeed, while these models may lack fundamental aspects of the early dynamics of planetary interiors, they allow a precise understanding of the phenomena and permit to consider large ranges for the (unconstrained) parameters.

The onset of convection of a homogeneous fluid with a temperature-and/or stress-dependent viscosity has been investigated in several contexts. The first instabilities induced by the sudden surface cooling of an initially hot fluid and its subsequent evolution was extensively studied (e.g., [START_REF] Davaille | Transient high-Rayleigh number thermal convection with large viscosity variations[END_REF]Choblet and Sotin, 2000;[START_REF] Korenaga | Physics of multiscale convection in earth's mantle: Onset of sublithospheric convection[END_REF][START_REF] Huang | Controls on sublithospheric smallscale convection[END_REF][START_REF] Zaranek | The onset of convection in fluids with strongly temperature-dependent viscosity cooled from above with implications for planetary lithospheres[END_REF][START_REF] Dumoulin | Onset of small-scale instabilities at the base of the lithosphere: scaling laws and role of pre-existing lithospheric structures[END_REF]. These studies mostly focus on the case of a fluid with strongly dependent viscosity (viscosity contrast > 10 6 ) in set-ups where the asymptotic stagnant lid regime occurs (see e.g. Solo-

A c c e p t e d M a n u s c r i p t matov, 1995). Applications to the evolution of the oceanic lithosphere on the Earth and the onset of small scale convection at its base have been proposed, for example, in order to explain the apparent heating from below away from hot spot tracks [START_REF] Davaille | Onset of thermal convection in fluids with temperature-dependent viscosity -applications to the oceanic mantle[END_REF], and many subsequent studies). A few studies also applied these results to the early dynamics of planetary interiors (see Choblet and Sotin, 2001, in the case of Mars). In the symmetrical case of heating from below, for the regime corresponding to strongly temperature dependent viscosity, hot instabilities develop that are initially not strong enough to penetrate the stiff cold material: small-scale convection is restricted to a hot sublayer whose boundary gradually extends upwards. This layer becomes eventually unstable and large scale convection develops. This regime was investigated by [START_REF] Thompson | Generation of mega-plumes from the coremantle boundary in a compressible mantle with temperature-dependent viscosity[END_REF] in a study concerning the formation of a superplume and by [START_REF] Solomatov | Small-scale convection in the D" layer[END_REF] as a possible explanation of the origin of D" layer within the Earth. [START_REF] Ke | Early transient superplumes and the origin of the Martian crustal dichotomy[END_REF] employed numerical simulations for this regime as a possible mechanism leading to the formation of the crustal dichotomy on Mars. In the case of icy-satellites, the non-Newtonian (power-law) rheology of ice is supposed to be a key parameter. The onset of convection is then often considered for an initially conductive steady-state, i.e., the critical value of the Rayleigh number defined classically for the whole layer is investigated rather than the transient process caused by the instantaneous heating/cooling of the fluid layer. [START_REF] Barr | Onset of convection in the icy Galilean satellites: Influence of rheology[END_REF] study the case of temperature and strain-rate dependent viscosity in this context. A conclusion is that the non-Newtonian rheology is important mainly in the case of ice with large grain sizes. Based on a different approach (an initially convective state is subjected to a decrease, step by step, of the Rayleigh number until convection vanishes), Solomatov andBarr (2006, 2007) focus on the influence of the rheology, demonstrating the differences between Newtonian and power-law viscosities and the dependence on the initial perturbation.

Here, we propose to investigate the onset of convection after the solidification of a magma ocean in the terrestrial planets, using three dimensional numerical experiments based on numerical tool OEDIPUS (Choblet, 2005;Choblet et al., A c c e p t e d M a n u s c r i p t 2007) in a basally heated spherical shell with either isoviscous or temperature dependent viscosity. This set-up is described in section 2. Intermediate values of the viscosity contrast (∆µ) across the hot boundary layer are considered (< 10 4 ).

This corresponds to the transitional regime between the isoviscous case and the asymptotic regime investigated in previous studies. Due to poorly constrained viscosity parameters at the appropriate conditions for deep planetary mantles (see e.g. [START_REF] Yamazaki | Some mineral physics constraints on the rheology and geothermal structure of Earth's lower mantle[END_REF][START_REF] Hirth | Rheology of the upper mantle and the mantle wedge: A view from the experimentalists[END_REF][START_REF] Korenaga | A new analysis of experimental data on olivine rheology[END_REF] and initial temperature distribution, both regimes (transitional or asymptotic regime) are indeed plausible. The methods are detailed in section 3 and appendices A and B. For a better understanding of the onset of convection, the full 3D numerical solutions (3.1) are systematically compared with the results of two simplified methods. The first approach is the numerical solution of "frozen-time" linear stability analysis (e.g. Chandrasekhar, 1961;[START_REF] Yang | The onset of thermmal convection in a horizontal fluid layer heated from below with time-dependent heat flux[END_REF], investigating the critical Rayleigh number for temperature dependent viscosity (3.2 and App. A). The second approach uses linearized Rayleigh-Taylor analysis in a simple two layers model (see e.g. [START_REF] Zaranek | The onset of convection in fluids with strongly temperature-dependent viscosity cooled from above with implications for planetary lithospheres[END_REF]Ke andSolomatov, 2004, 2006) where the onset time is determined as the time when the growth of the R-T instability exceeds the growth of the thermally induced boundary layer in the convection framework (3.3 and App. B). The differences between dynamical (in our case free-slip) and kinematical (no-slip) boundary conditions are investigated. The free-slip boundary condition approximates the interface between solid and liquid material, i.e. this corresponds to the boundary between the solid mantle and the liquid outer core of terrestrial planets or to icy satellites with an icy crust overlying deep internal oceans. The no-slip condition reflects interfaces between two solid layers, this may represent the boundary between an icy mantle and silicate-rich core, for example. Section 4 describes the numerical results: we focus especially on the influence of the vigor of convection (4.1) and of the temperature dependency of viscosity (4.2). The geometry of the shell is addressed as an additional parameter (4.3). Finally, a generic scaling relationship describing these effects is introduced (4.4). All the aspects concerning the scaling relationship are further A c c e p t e d M a n u s c r i p t discussed in section 5. Such laws are then used with the appropriate scaling factor in order to assess the onset of convection within terrestrial planets (section 6).

Governing equations

In the following, the Boussinesq approximation for infinite Prandtl number is taken into account. The dimensionless equations (conservation of mass, momentum and energy with neglected viscous and internal heating) are

0 = ∇ ′ • v ′ , (1) 0
= -∇ ′ p ′ + ∇ ′ • µ ′ (θ ′ ) ∇ ′ v ′ + ∇ ′ T v ′ -Raθ ′ e r , (2)
∂θ ′ ∂t ′ = -v • ∇ ′ θ ′ + ∇ ′2 θ ′ . ( 3
)
where Ra is Rayleigh number. • ′ denotes the dimensionless variables and the following scaling is used:

x = dx ′ , time t = d 2 κ t ′ , velocity v = κ d v ′ , pressure p = µ0κ d 2 p ′ , temperature T = T 0 + ∆T θ ′ , µ ′ (θ ′ ) = µ(T ) µ0 = µ(T ) µ(T0) , where d = r t -r b
is the thickness of the fluid layer (r b and r t are associated with inner and outer boundary, respectively), T 0 is the surface temperature, ∆T is the temperature difference across the shell and κ is the thermal diffusivity.

The viscosity is supposed to depend exponentially on temperature

µ(T ) = µ 0 exp -a vis T -T 0 ∆T , µ ′ (θ ′ ) = exp (-a vis θ ′ ) , (4)
where a vis is the variable viscosity parameter. This dependency is an approximation of the viscosity described by the Arrhenius law:

µ(T ) = A exp Q * RT , (5)
where Q * is the activation enthalpy and R the gas constant. In the framework of a basally heated layer, the two laws are similar for a vis = Q * ∆T RT0(T0+∆T ) .

6

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

The initial and boundary conditions for temperature are θ ′ (r, t = 0) = 0, θ ′ 0 (r = r b , t) = 1 and θ ′ (r = r t , t) = 0. No internal heat sources are introduced.

Besides the initial and boundary conditions, the convective behavior described by conservation laws (1-3) and by relationship (4) is controlled by three nondimensional parameters: the Rayleigh number Ra, the viscosity parameter a vis and the spherical shell geometrical factor f defined by the ratio of the inner and outer radius of the shell f = r b rt = rt-d rt .

The Rayleigh number Ra(θ ′ ) for a given viscosity (temperature) is

Ra(θ ′ ) = ρ 0 g∆T αd 3 κµ(θ ′ ) (6)
with ρ 0 reference density, α thermal expansion, g the gravity acceleration. The surface Rayleigh number Ra 0 = Ra(θ ′ = 0) and the bottom Rayleigh number Ra b = Ra(θ ′ = 1) are considered in the following.

Method

3D convection

The numerical method described in Choblet (2005) and Choblet et al. (2007) is used to obtain a three dimensional solution of the system (1-3) in the spherical shell. The composite mesh based on the "cubed sphere" [START_REF] Ronchi | The "cubed sphere": a new method 1079 for the solution of partial differential equations in spherical geometry[END_REF] transformation is employed, the resulting grid consisting in six identical blocks.

Due to the time demands, most of the convection simulations are carried out in one block. Additional vertical boundaries are hence introduced where free-slip and no-heat-flux conditions are prescribed. In order to test the influence of these artificial boundaries, several tests in the whole spherical shell are also performed.

The computational grid in one block consists of 32 × 64 × 64, 64 × 64 × 64 or 128×64×64 discrete cells depending on the Rayleigh number, viscosity contrast and the geometry (the unstable hot boundary layer is described at least 5-6 points in the vertical direction at the onset of convection).

At the beginning of the simulation, the temperature field is T 0 within the shell and it is stochastically perturbed with a maximum amplitude of 10 -3 ∆T .

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

The onset time for convection is defined as the time when the maximum deviation between the horizontally averaged temperature and the conductive heating profile reaches locally a value of 0.1% (see Fig. 1).

The preferred degree (wavelength) of the instabilities at the onset time are also estimated. The "degree" in one block can be determined by two dimensional fast Fourier transform [START_REF] Press | Numerical recipes 1065 in FORTRAN: the art of scientific computing[END_REF]. Assuming the preferred wavelength remains identical, the coresponding estimate of the preferred degree for the whole sphere is then obtained by the multiplication of this value by factor 4 (one block corresponds to one fourth of the sphere, i.e. to π /2 of the sphere in both equatorial and meridional directions). This provides a guidance for the relative behavior (i.e. increase or decrease) of the preferred degree as a function of the studied parameters. For experiments in the whole shell, the spherical harmonic expansion coefficients are computed by integration and lead to a precise evaluation of the preferred degree.

Linear stability analysis

The linear stability analysis belongs to a traditional approach for computing the "classical" critical Rayleigh number, i.e. for computing critical Rayleigh number for steady-state conduction solution. For constant viscosity and spherical geometry, this approach is described in detail (including the influence of the shell geometry and preferred degree) in Chandrasekhar (1961). The effect of the temperature dependent viscosity as well as the influence of the shell geometry on "classical" critical Rayleigh number and preferred degree is analyzed in [START_REF] Ratcliff | Steady tetrahedral and cubic patterns 1068 of spherical shell convection with temperature-dependent viscosity[END_REF]. Here, we look for the critical Rayleigh number and its corresponding degree (i.e. associated to the smallest time at which convection occurs) in a model with depth-dependent viscosity. The temperature profile first evolves conductively in the basally heated shell. A critical Rayleigh number Ra c (t ′ ) is computed for a given time t ′ and thus a given temperature profile. These results may be related to the ones obtained for 3D numerical experiment supposing that the onset of convection corresponds to the first time when the Rayleigh number of the 3D convection experiment exceeds Ra c (t ′ ). This method, however, does

A c c e p t e d M a n u s c r i p t not take into account the propagation of the wave front, i.e. the "frozen time"

problem (e.g. [START_REF] Yang | The onset of thermmal convection in a horizontal fluid layer heated from below with time-dependent heat flux[END_REF] see differences between the "frozen time" model with "propagation theory") is solved. This aspect is discussed further in the next section. Technical details concerning this method are described in appendix A.

Rayleigh-Taylor instability analysis

Similarly to linear stability analysis, the Rayleigh-Taylor instability also belongs to classical (semi)analytical approaches with many possible applications to mantle dynamics (see e.g. [START_REF] Ribe | Analytical approaches to mantle dynamics[END_REF] for a summary). For example, both

Cartesian configurations with layers of different viscosities (e.g. Canright and Morris, 1993) and spherical layers (e.g. [START_REF] Ribe | The global hotstpot distribution and instability 1074 of D[END_REF]) have been

considered. Some models also include more complex geometry of the buoyant structure such as the cylindrical anomaly studied by [START_REF] Lister | The effect of geometry on the gravitational instability 1059 of a buoyant region of viscous fluid[END_REF] in the context of diapirism beneath mid-ocean ridges. Here the onset time for "convection" is defined as the time when the growth of any R-T instability exceeds the conductive propagation of the heat front [START_REF] Zaranek | The onset of convection in fluids with strongly temperature-dependent viscosity cooled from above with implications for planetary lithospheres[END_REF][START_REF] Ke | Plume formation in strongly temperaturedependent viscosity fluids over a very hot surface[END_REF]. A two layers model is considered so that the viscosity profile of the 3D convection experiment is modeled by a step-like function characterized by the viscosity contrast ∆µ. Several definitions of ∆µ are possible based on either a constant temperature fraction ∆θ (∆µ = exp(∆θa vis )) or on the temperature average in each layer (∆μ = exp(∆ θa vis )). These various definitions are discussed in the following. This approach is described in details in appendix B.

Results

Here we systematically compare 3D numerical results with the two simplified approaches. We study the onset times depending on the Rayleigh number Ra, the shell geometrical factor f , viscosity parameter a vis and either free-slip or no-slip.

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

Tab. 1 lists results of 3D numerical experiment for free-slip conditions (Tab. 2 for no-slip conditions). The estimate of the wavelength (degree) of the instabilities and the results for the computations in six blocks are also shown in these tables. The results obtained for one block (i.e. with additional vertical boundaries) do not differ significantly from those where the whole shell is considered.

The largest difference of the onset time between two similar experiments in one and six blocks is ≈ 5%. This discrepancy can also be partly explained, however, by the different initial perturbations (for stochastic nature of the onset time see e.g. [START_REF] Korenaga | Physics of multiscale convection in Earth's mantle: Evolution of sublithospheric convection[END_REF]. The preferred degrees also agree rather well despite the uncertainty inherent to the case of runs performed on one block (the value is forced to be a multiple of 4).

Furthermore, the numerical experiments were performed for different initial perturbations in the case of constant viscosity and Rayleigh number Ra = 10 6 .

As expected, the onset time is decreasing with increasing amplitude of the initial perturbation. Decreasing the initial perturbation from 10 -3 to 10 -5 for free-slip conditions, the onset time increases by factor ≈ 1.3. This corresponds relatively well to results previously obtained by [START_REF] Korenaga | Physics of multiscale convection in Earth's mantle: Evolution of sublithospheric convection[END_REF] who found the approximately constant factor ≈ 1.5.

Influence of Rayleigh number

As an example, we discuss the results obtained for a given geometry (f = 0.55) and an isoviscous fluid (a vis = 0). We focus on the dependency of the onset time on the value of the Rayleigh number. The onset times and corresponding wavelength for all three approaches are summarized in Fig. 2 for both freeslip (Fig. 2a) and no-slip (Fig. 2b) boundary conditions. All three approaches induce a linear character in the log-log scale, hence the interpolation by power law t ′ = ARa a is used as an analytical tool: inverted values for A and a are listed in Tab. 3. The dependency of the onset time is rather similar for both mechanical boundary conditions (a = -0.67 for free-slip and a = -0.69 for no-slip). The obtained slopes agree well with the typical onset time scaling t ′ ∝ Ra -2 /3 (e.g. [START_REF] Blair | The onset of cellular convection in the fluid layer with time-dependent density gradients[END_REF][START_REF] Jhaveri | Randomly forced Rayleigh Bérnard convection[END_REF]; Choblet and A c c e p t e d M a n u s c r i p t [START_REF] Solomatov | Fluid dynamics of a terrestrial magma ocean[END_REF][START_REF] Huang | Controls on sublithospheric smallscale convection[END_REF].

Sotin

Results of the linear stability analysis for both boundary conditions also exhibit a linear character for a log-log scale when t ′ < 0.01 (Fig. 2, solid thick line). For larger times, the temperature on the upper boundary using formula (A12) is non-zero and the power-law scaling is not valid anymore (this effect of the heat front reaching the upper boundary is less pronounced for 3D convection results). For this reason, onset times larger than 0.01 are not considered when computing the slope. The influence of Rayleigh number is significantly stronger for LS than for 3D results (a = -1.02 for free-slip and a = -0.72 for no-slip). Moreover, the values of the onset time for LS are more than one order of magnitude smaller. Both effects should probably be attributed to the use of the "frozen time" approach-in the 3D convection results, the upward advective propagation of the instabilities needs to be faster than the conductive propagation of the heat front which is neglected in this approach (for influence of the "non-frozen" time approach see e.g. [START_REF] Yang | The onset of thermmal convection in a horizontal fluid layer heated from below with time-dependent heat flux[END_REF]. If the boundary layer is defined according to a thickness proportional to the square root of the time (h ′ ∼ √ t ′ ) then the propagation velocity of the heat front evolves according to the relationship ḣ′ ∼ 1 / √ t ′ . The "frozen time" effect is thus higher for lower onset times than for higher ones. This behavior may explain the larger negative slope in the case of free-slip conditions. It may also partly cause the discrepancy between the slopes observed for no-slip and free-slip results: the onset times associated to no-slip boundary conditions are higher and the effect of the "frozen time" approach is thus less significant.

"Onset times" determined by the Rayleigh-Taylor instability analysis also underestimate the 3D values (however, significantly less than the LS analysis).

Again, the curve is almost linear for the studied range. The power-law scaling leads to a = -0.70 for free-slip and a = -0.70 for no slip and agrees relatively well with slopes obtained from 3D simulations.

The degree of the fastest growing anomalies increases with increasing Ra for all methods (i.e. thinner boundary layers lead to smaller preferred wavelengths of the first instabilities). The preferred degrees are higher for no-slip than A c c e p t e d M a n u s c r i p t

for free-slip condition. However, they vary strongly among the methods. The lowest degree is obtained for the linear stability analysis. For 3D numerical experiments, the estimate of the degree is limited by the common multiple factor 4 when only one block is used. The increase of the degree with increasing Rayleigh number is however obvious. The preferred degree is highest for the R-T analysis. This may be due to the fact that, in the vicinity of the minimum, the dependency of the onset time on the degree is rather low (see App. B; Fig. B1a). Considering that the reciprocal value of the fastest growing degree is proportional to the thickness of the boundary layer and assuming h ′ ∼ √ t ′ , the power law scaling for Rayleigh number Ra leads to

1 l ∼ h ∼ √ t ∼ Ra -1 /3 ⇒ l ∼ Ra 1 /3 . (7)
The values obtained for preferred degree by the R-T instability analysis follow well this prediction l ∼ Ra 0.33 .

Influence of the viscosity variations

Until now, we have discussed only results for constant viscosity. In this section, we focus on the temperature dependent viscosity. In our computations, we consider values of a vis smaller than 10. This corresponds to the transitional regime observed between the isoviscous regime and the asymptotic regime associated to very large viscosity contrasts investigated earlier by some authors [START_REF] Thompson | Generation of mega-plumes from the coremantle boundary in a compressible mantle with temperature-dependent viscosity[END_REF][START_REF] Solomatov | Small-scale convection in the D" layer[END_REF][START_REF] Ke | Plume formation in strongly temperaturedependent viscosity fluids over a very hot surface[END_REF].

Fig. 3 summarizes results for temperature dependent viscosity and free-slip conditions. Fig. 3a shows the dependency of the onset time t ′ on the viscosity parameter a vis keeping a constant value for the surface Rayleigh number Ra 0 .

In all the results presented in this paragraph, the initial viscosity within the spherical shell is thus identical. The viscosity variations with temperature within the region where the hot front propagated, increase with a vis . For all methods, a nearly linear character in the log scale is observed. We thus use a relationship in the form t ′ = B exp(ba vis ). The interpretation of this scaling is discussed in convection results and both simplified methods when viscosity variations are introduced. For R-T analysis, different definitions of the viscosity contrast ∆µ between the two prescribed layers were tested. First, a series of fractions of the global temperature difference are used: ∆µ = exp(a vis ∆θ) with ∆θ equal to 0.2, 0.4, 0.6, 0.8 and 1 (cf. Eq. (B14), dashed lines in Fig. 3). Second, the difference ∆ θ between the average temperature in each layer is used (dashed dotted line in Fig. 3). The slope of the curve is influenced significantly by this definition.

The strongest dependency is obtained for ∆θ = 1 (b = -0.25). The case where ∆ θ is used to define the viscosity contrast lies between ∆θ = 0.4 and 0.6 and leads to a smaller slope (b = -0.15). This value is close to the one obtained for the 3D results (b = -0.16). For the linear stability analysis, the slope is almost zero (Fig. 3a), reflecting the fact that the onset time is not strongly influenced by the increasing viscosity parameter a vis . In this case with free-slip boundaries, the onset time is predominantly controlled by the Rayleigh number on the surface. While the difference between 3D and LS results was found to increase with decreasing values of the onset time in the isoviscous case (Fig. 1a), this difference diminishes here (Fig. 3a). In fact, while the conduction of the heat front is identical whether the fluid is isoviscous or not, the growth rate of instabilities is larger when a vis increases. And since the onset time decreases more slowly in this case than for constant a vis and increasing Ra, the frozentime effect is less important. No-slip results (not presented in Fig. 3 onset of convection. Consequently, the scaling t ′ ∝ Ra -a is not valid for these large values of the onset time. In the case of the LS analysis, this deviation of large values of the onset times from the general trend is observed even in the isoviscous case (Fig. 2). However, in the case of the R-T analysis, the influence of the upper boundary is detected for smaller values of the onset time due to the definition of the viscosity contrast between the layers: while results up to t ′ = 5 • 10 -2 seem to follow the same unique slope in the isoviscous reference case (Fig. 2), results above t ′ > 10 -2 deviate when viscosity variations are introduced (Fig. 3b). For LS stability and R-T instability analysis, we therefore take here into account only onset times t ′ < 0.01. The corresponding segments of the curves are then also linear so that a fit based on a scaling similar to the calculations with constant Ra 0 is adopted; parameters are noted (B ′ , b ′ ) and the inverted values are presented in Tab. 3. For this restricted range, the agreement is good between the values obtained for the 3D convection results (b ′ = 0.56) and the R-T instability when the viscosity contrast is based on ∆ θ (b ′ = 0.54).

Again the results for viscosity based on ∆ θ lie between the results obtained for ∆θ = 0.4 and 0.6.

For all methods, the degree of the fastest growing anomaly does not vary monotonically with increasing a vis . This more complex behavior compared with the isoviscous case, may be caused by the existence of several regimes as discussed by [START_REF] Ke | Early transient superplumes and the origin of the Martian crustal dichotomy[END_REF].

Influence of the shell geometry

We now focus on the influence of the geometrical factor f of the shell. Results

for constant viscosity and Rayleigh number Ra = 10 6 are shown in Fig. 4a-b.

Values of the onset times for the LS analysis (solid line) and the R-T instability analysis (dashed line) are again systematically lower than those of 3D numerical experiment (solid circles). For the interpolation, we propose a purely mathe-

matical function t ′ = C • f c • exp(-cf ).
This scaling satisfies the expected limits based on the energetic expectations for the onset times for f approaching 0 or 1 when the shell thickness is used as a characteristic length scale. For f → 1, the A c c e p t e d M a n u s c r i p t onset time is expected to tend smoothly to the Cartesian limit: lim f →1

∂t ′ ∂f = 0.
For f → 0, the onset time should be infinity (which is the case of the proposed function for c < 0). Due to the asymptotic behavior of the onset time corresponding to low values of f , the inverted value for c depends on the investigated range and on the sampling of the geometrical factor f .

When results of the linear stability analysis in the case of no-slip boundaries are considered, the best fit returns a positive value of c (Tab. 3) in contradiction with the assumption inherent to the above scaling. This is due to the complex behavior of the onset time caused by increasing preferred wavelengths with decreasing f (Fig. 4b). This can also be observed to a lesser extent for free-slip LS and for R-T. A key result is that the maximal curvature encompassed by the c value in the proposed scaling is significantly different between the LS results

for free-slip (c = -0.29) and the R-T results (|c| < 10 -2 for both free-slip and no-slip). This reflects the fact that onset times are almost independent of f above a given value and this value is larger for LS than for R-T.

In the case of 3D results, the lack of numerical experiments in the region f ∈ (0, 0.2) prevents a precise assessment of the c value. For this reason, we do not compare it with the results obtained for LS and R-T. A slight increase of the onset times with f can be observed for the studied range both for free-slip and no-slip. This gradual increase is also detectable for LS (and, to a lesser extent, R-T), especially in the no-slip case. A candidate effect for this range of f , where the geometrical effect on the heat budget gets less and less significant, might be that since the preferred wavelength is forced to increase when f gets close to 1, the onset time is delayed. The fact that the R-T results seem to be less affected would then be related to the small dependency of the growth rate as a function of preferred degree, as already mentioned (see Fig. B1a).

We also investigate the scaling in terms of Rayleigh number for different values of the geometrical factor f using R-T instability analysis. The functions a(f ) and b(f ) are introduced to quantify this effect (t ′ ∝ Ra A general conclusion is that, although small values of f induce that a smaller amount of energy is supplied in order to prescribe the hot temperature on the inner boundary, the values of the dimensionless onset time are globally independent of f . Thus, the hot front always reaches a similar fraction of the shell thickness d before it becomes unstable. As could be expected, the preferred wavelength of the first instabilities grows with decreasing geometrical factor f .

Scaling relationships

Considering the scaling of the onset time as a function of both Rayleigh number Ra and viscosity parameter a vis and neglecting the influence of the geometry on dimensionless onset time (a and b are not functions of f , thus assuming f 0.2), we propose a combination of the relationship obtained by varying the two parameters independently:

t ′ ∝ Ra a 0 exp (ba vis ) . (8)
The values obtained for 3D numerical experiments for free-slip conditions are a ≈ -0.67 and b ≈ -0.16 (see Tab. 3). Using Eq. ( 8

Since the runs with constant Ra b provide t ′ ∝ exp (0.56a vis ) (see Tab. 3) and Eq. ( 9) gives b -a ≈ 0.51, the assumptions in Eq. ( 8) seem to be reasonably valid as a general description of the two dimensional parametrical space. If the same procedure is applied for 3D numerical experiments with no-slip boundaries and a ≈ -0.69 and b ≈ -0.20 (see Tab. 3), we get b -a ≈ 0.49. Again, this agrees rather well with results obtained for constant Ra b (t ′ ∝ exp (0.52a vis )).

The equation ( 8) can also be interpreted as

t ′ ∝ (Ra * ) a (10)
where

Ra * = ρ0g∆T αd 3 κµ *
is associated with a viscosity value µ * = µ (θ * ) that "controls" the onset time of the convection, replacing the explicit scaling in terms of viscosity parameter. Since Ra * = Ra 0 exp (θ * a vis ), the temperature defining the viscosity µ * can be expressed as θ * = b a (cf. Eqs. 8 and 10). Using the values listed in Tab. 3, we obtain θ * ≈ 0.24 for free-slip and θ * ≈ 0.29 for no-slip. This is summarized for all methods in Fig. 5 where the scaled onset time t ′ Ra(θ * ) -a is shown as a function of a vis .

The validity of this scaling is demonstrated in Fig. 5 for 3D calculations.

The onset time t ′ is normalized using the dependency predicted by Eq. 10 ( t′ = times show an approximately linear dependence on viscosity parameter a vis in the log-scale. In the case of R-T results and no-slip LS results, the value of the slope is lower than 0.05. However, due to significantly different scaling for freeslip LS results, the slope is higher (β ≈ 0.15). The normalized onset times for a constant bottom Rayleigh number Ra b show more complex variations. These follow an almost linear trend up to a given (transitional) value of the viscosity parameter a vis . For higher values, the onset of convection is affected by the presence of the upper boundary. The transitional value of a vis is approximately 5 in the case of R-T instability analysis. For LS results, it is shifted toward slightly higher values (≈ 7) as a result of lower onset times.

t ′ • (Ra * ) -a )

Discussion

For the isoviscous 3D runs, the scaling of the onset time approximately follows the typical scaling t ′ ∝ Ra -2 /3 for both free-slip and no-slip conditions.

The obtained slopes agree well with previous Cartesian numerical and laboratory studies for a fluid cooled from above (e.g. [START_REF] Blair | The onset of cellular convection in the fluid layer with time-dependent density gradients[END_REF][START_REF] Jhaveri | Randomly forced Rayleigh Bérnard convection[END_REF]Choblet and Sotin, 2000;[START_REF] Huang | Controls on sublithospheric smallscale convection[END_REF]. For an isoviscous fluid described by Eqs. ( 1)-( 3), the two cases (heating from below, cooling from above) are symmetric and the scaling should be identical. This symmetry also holds more-or-less for the spherical shell (f 0.2, depending on Rayleigh number). Moreover, while no tests were performed in 3D, it was shown for the R-T results that this slope is also rather constant for a broad range of values for the geometrical factor f (in our experiments, the slopes obtained for R-T instability analysis are systematically comparable to 3D numerical experiments).

However, when the viscosity contrast is increased the symmetry of the flow disappears (e.g. [START_REF] Solomatov | Scaling of temperature-and stress-dependent viscosity convection[END_REF]. For cooling from above and large viscosity contrast, a cold thermal boundary layer first develops. The onset of convection occurs when the thin sublayer at the base of the cold and viscous boundary layer reaches the critical Rayleigh number (see e.g. [START_REF] Dumoulin | Onset of small-scale instabilities at the base of the lithosphere: scaling laws and role of pre-existing lithospheric structures[END_REF].

The subsequent convective motions develop in the whole domain except for the A c c e p t e d M a n u s c r i p t stagnant lid. The following asymptotic scaling for the viscosity parameter is used

t ′ ∝ Ra a • a β vis , (11)
where 5 6 < β < 7 6 is proposed [START_REF] Huang | Controls on sublithospheric smallscale convection[END_REF][START_REF] Zaranek | The onset of convection in fluids with strongly temperature-dependent viscosity cooled from above with implications for planetary lithospheres[END_REF][START_REF] Dumoulin | Onset of small-scale instabilities at the base of the lithosphere: scaling laws and role of pre-existing lithospheric structures[END_REF]. In the case of a fluid layer heated from below and again considering large limit of ∆µ(> 10 6 ), a "symmetrical" regime has been described [START_REF] Thompson | Generation of mega-plumes from the coremantle boundary in a compressible mantle with temperature-dependent viscosity[END_REF][START_REF] Solomatov | Small-scale convection in the D" layer[END_REF]Ke andSolomatov, 2004, 2006). Convective instabilities develop in the low viscosity boundary layer. These are not strong enough at first to penetrate into the stiff cold material above and small-scale convection is restricted to the low-viscosity sublayer. With further heating, this sublayer grows and becomes eventually unstable: a large scale, low degree pattern develop.

Our study address the intermediate range between the isoviscous case and this asymptotic regime (contrast ∆µ < 10 4 ). At steady-state, this regime would correspond to the "sluggish lid" regime where cold and more viscous instabilities participate to convection (e.g. [START_REF] Solomatov | Scaling of temperature-and stress-dependent viscosity convection[END_REF]. This transitional regime has been less studied than the classical asymptotic stagnant-lid regime: [START_REF] Korenaga | Physics of multiscale convection in earth's mantle: Onset of sublithospheric convection[END_REF] propose a "unifying" scaling of the various regimes, introducing a functional based on the concept of available buoyancy (c.f. Conrad and Molnar, 1999). In their laboratory experiments mainly focused on the stagnant lid regime, [START_REF] Davaille | Transient high-Rayleigh number thermal convection with large viscosity variations[END_REF] also investigate the range between the isoviscous case and the asymptotic regime (∆µ < 10 6 ). and θ * denotes the value of temperature corresponding to the viscosity value that "controls" the instability.

Although the dynamics of the onset of convection are complex when viscosity gradients of various amplitudes are considered, using this specific value µ(θ * ) (a constant fraction of the viscosity contrast) as a characteristic viscosity in the definition of the Rayleigh number is sufficient to describe the influence of a vis . Our results indicate that a ≈ -2 /3 and θ * ≈ 0.25 for both free-slip and no-slip boundary conditions. Furthermore, the scaling obtained independently for constant surface and bottom Rayleigh numbers show that the assumptions inherent to this scaling are reasonable (see section 4.4).

Implication for planets

Applicability to planets Although the onset of solid-state flow in planetary interiors is a complicated process influenced by many factors, the values of the onset times reported here and, most importantly, the observed slopes in the simple three parameters space provide interesting guidelines with regard to the early dynamics of terrestrial planets. We propose a set-up where a homogeneous layer is heated from below. This provides an idealized model for the "initial" thermal state of a newly formed planet, once core formation is completed and subsequent to the solidification of a large fraction of a plausible magma ocean.

In fact, a global magma ocean resulted likely from the rapid occurrence of a series of primordial energetic events. Whether hafnium-tungsten chronometry based on rocks of presumably planetary origin-Mars, the Moon, the Earth, Vesta-e.g. [START_REF] Yin | A short timescale for terrestrial planet formation from hf-w chronometry of meteorites[END_REF]; [START_REF] Kleine | Rapid accretion and early core formation on asteroids and the terrestrial planets from Hf-W chronometry[END_REF] provides a precise timing for the metal-silicate separation is controversial, e.g. [START_REF] Wood | Cooling of the Earth and core formation after giant impact[END_REF]; [START_REF] Allègre | The major differentiation of the Earth at 4.45 Ga[END_REF]. In all cases, geochemical data indicate that the completion of this process is likely to occur before 100 Myr after the formation of the solar system. In this context, the few models proposing colder scenarios where core formation is significantly delayed (as for Mars, [START_REF] Senshu | Thermal evolution of a growing Mars[END_REF] are problematic. The precise partition of the heat associated with core formation between 20

A c c e p t e d M a n u s c r i p t the silicate mantle and the iron-rich core is, however, still largely uncertain.

These aspects are further detailed below.

Simplifications Neither radioactive heating nor dissipation are considered in the models presented here. Although the value of the viscous dissipation may depend on the initial amplitude of the perturbation analogically to what has been proposed in the case of a power-law viscosity [START_REF] Solomatov | Onset of convection in fluids with strongly temperature-dependent, power-law viscosity 2. dependence on the initial perturbation[END_REF], the influence on the onset of convection is probably relatively small. Short-lived radiogenic isotopes such as 26 Al (and to a lesser extent 60 Fe) are likely to be present in the primordial bricks building the silicate part of planets. If accretion is achieved shortly after the formation of CAIs (< 10 Myr), these may contribute efficiently to the heat budget. In the case where heating from below remains the prominent energy supply of the mantle, these additional internal sources would reduce the onset time.

Furthermore, possible chemical variations are not considered. If the mantle is heterogeneous, a two layers system could develop (where a stable deep layer may be preserved depending on compositional buoyancy). The specific dynamical regimes associated with such a set-up are described by Le [START_REF] Bars | Whole layer convection in a heterogeneous planetary mantle[END_REF], for example. In the case of the Earth, if such a stratification exists at present-day, it might be a consequence of plate dynamics; the deep reservoir could be created by slab remnants or delaminated continental material. Therefore, the presence of this stratification throughout the evolution is not required. In the case of Mars, one key issue is the possibility of a magma ocean cumulate overturn as proposed by [START_REF] Elkins-Tanton | Magma ocean fractional crystallization and cumulate overturn in terrestrial planets: Implications for Mars[END_REF], 2005). This would result in a stable stratification of the deep mantle due to composition and would strongly inhibit whole mantle convection. In this case, radiogenic elements would be highly concentrated in the deep layer. The problem of the onset of purely thermal convection would then mostly concern the above mantle layer in a configuration similar to that described by Le [START_REF] Bars | Whole layer convection in a heterogeneous planetary mantle[END_REF] in the case of a very large value of the buoyancy number (measuring the ratio between the density variations of compositional origin to the ones caused We also omit the effect of the phase transitions which may affect the onset times and the wavelength of the preferred instabilities especially if these are located close to the core-mantle boundary (CMB). In the case of the Earth, the exothermic post-perovskite phase transition near CMB is a candidate for such a deep transition, however, its presence is rather unlikely for the early Earth due to possibly higher temperatures (e.g. [START_REF] Oganov | Theoretical and experimental evidence for a post-1061 perovskite phase of MgSiO 3 in Earth's D" layer[END_REF]. On Mars, the presence of the endothermic transition to perovskite near core-mantle boundary has also been mentioned as a possibility (e.g. Breuer et al., 1997;[START_REF] Harder | Phase transition and three-dimensional planform of thermal convection in the Martian mantle[END_REF][START_REF] Roberts | Degree-1 convection in the Martian mantle and 1076 the origin of the hemispheric dichotomy[END_REF] and will be discussed later.

Parameter values for terrestrial planets In addition, the parameter values for the simple convection model we investigate, are rather uncertain especially under the condition of early planets. Two of the parameters considered in the present study (namely, the Rayleigh number Ra and the viscosity parameter a vis ) strongly depend on (i) the internal temperature (T 0 in our models), (ii) the temperature contrast through the boundary layer (∆T in our models) and

(iii) the values of viscosity parameters, especially activation parameters associated with planetary materials. As already mentioned, both (i) and (ii) directly result from the energy budget of planetary accretion and core formation and from the heat partitioning between mantle and core. The energy associated with impacts by bodies of varying size [START_REF] Tonks | Magma ocean formation due to giant impacts[END_REF], including the probable Moon-forming impact on the Earth at the end of accretion (Canup, 2004) but also, maybe, the large impact responsible for the Martian hemispheric dichotomy (e.g. [START_REF] Andrews-Hanna | The Borealis basin and the origin of the martian crustal dichotomy[END_REF], is most probably large enough to [START_REF] Fei | The interior of Mars[END_REF]. The upper limit corresponds to the solidus temperature 607 (≈ 2300 K) [START_REF] Takahashi | Speculations on the Archean mantle: missing link between komatiite and depleted garnet peridotite[END_REF]. For the Earth, the current temperature at depth [START_REF] Peltier | Mantle viscosity and ice-age ice sheet tomography[END_REF][START_REF] Lambeck | Sea-level change, glacial rebound and mantle viscosity for northern europe[END_REF].

646

The relationship between parameter a vis used in the 3D experiments and actual 647 values of activation parameters E * and V * is a vis = (E * +pV * )∆T RT0(T0+∆T ) . However, both Mars. We use the scaling relationship obtained for the 3D numerical experi-667 ments (Eq. 8) for free-slip boundaries in order to compute the onset time t onset .

668

We assume that our scaling is valid up to at least a viscosity parameter value 669 of a vis = 10. Above this value, as mentioned in section 5, a different regime 670 may occur where small-scale convection first develops with a convective thick-671 ening of a low-viscosity sublayer before large scale low-degree instabilities affect 672 the whole layer (e.g. [START_REF] Ke | Early transient superplumes and the origin of the Martian crustal dichotomy[END_REF]. We thus only report results

673

corresponding to values of the viscosity parameter smaller than 10. As noted 674 above a vis depends on T 0 and Q * . In Fig. 6a-b, the region in the parameter 675 space where a vis > 10 is indicated as an empty area. In the case of Mars, this 676 occurs for a temperature contrast above the CMB of more than 1700 K. In the 677 Earth, the viscosity parameter is lower than 10 in the whole range considered the parameter space. The lower limit is obtained for the lowest viscosity considered and the highest temperature contrast. In contrary, the upper limit is reached for the highest viscosity and lowest temperature contrast. Note that, as indicated by the isocontours, there is a systematic trade-off between the two parameters proposed in Fig. 6a-b, η 0 and ∆T . The main effect is caused by viscosity: as can be seen from the scaling relationship, for a constant ∆T , an increase of one order of magnitude is obtained when η 0 is decreased by a factor of 30. For a constant η 0 , an increase of 1000 K for ∆T results in 4 times smaller onset time.

As mentioned above, the viscosity parameter a vis and thus the onset times depends besides the free parameter ∆T also on two rather uncertain parameters Q * and T 0 . In the case of the Earth and Venus, the influence of the internal temperature T 0 is demonstrated in Fig. 6, plates c (T 0 = 2500 K), plate b (T 0 = 3000 K) and plate d (T 0 = 3500 K). As we already discussed, the parameter a vis increases with decreasing T 0 . The empty area corresponding to a vis > 10 thus increases with decreasing T 0 (cf. Fig. 6 plates b and c). Nevertheless, the onset times are not influenced significantly for the investigated range of the internal temperatures T 0 : a slight shift towards the higher onset times is observed for constant temperature contrast ∆T if T 0 increases and the upper and lower boundaries of the estimated onset times are rather similar. In the case of variable activation enthalpy Q * , a similar effect (the increase of geometrical factor with increasing activation enthalpy) is observed.

Although highly variable, these values of the onset times are systematically smaller than the values obtained for the onset of convection beneath a stagnant lid forming due to cooling from above (see Choblet and Sotin, 2001, for a simple scaling) or [START_REF] Zaranek | The onset of convection in fluids with strongly temperature-dependent viscosity cooled from above with implications for planetary lithospheres[END_REF], where a more sophisticated model where compositional stratification is taken into account, thus delaying the onset of convection). A consequence is thus that, once the core is formed and the solidification of the magma ocean is achieved, whether an overturn occurs as is possible for Mars, or not, the bulk mantle is likely to be destabilized first by upwelling instabilities first and downwellings from the cold lithosphere will ing possibly due to a change of creep mechanism at depth. In these models, a viscosity increase higher than 8 at a depth of ≈ 1000 km systematically produces a degree one instability that seems stable independently of the Rayleigh number. In a recent paper [START_REF] Zhong | Migration of Tharsis volcanism on Mars caused by differential rotation of the litosphere[END_REF], the idea of a weak asthenosphere is further developed and, together with lithospheric variations, is proposed as a mechanism to explain not only the formation of the crustal dichotomy but also the spatial and time evolution of volcanism on the Tharsis Rise (namely the migration from south to north): differential rotation of the lithosphere excited by degree one convection is obtained in such models. This configuration cannot be investigated by purely basally heated convection and is not addressed in this study. Whether the viscosity increases in the mid-mantle of Mars remains an open question. A better assessment of the present-day relationship between topography and areoid could help to resolve this issue.

Conclusions

In order to investigate the onset of solid-state convection within the interior of terrestrial planets, 3D numerical experiments are presented for the idealized model of a spherical shell heated from below. These results are systematically compared with two simplified methods: linear stability (LS) analysis and Rayleigh-Taylor (R-T) instability analysis. We focus especially on the dynamical regime obtained for a viscosity contrast ∆µ ranging between 1 and 10 4 . The results of 3D numerical experiments are well approximated by a scaling relationship t ′ ∝ (Ra * ) a where a ≈ -2 /3 and Ra * is a specific Rayleigh number corresponding to temperature θ * ≈ 0.25. This latter parameter includes the effect of various viscosity contrasts in a simple framework. Although this scaling is valid for both mechanical boundary conditions, the onset times for free-slip are systematically lower than for no-slip, as expected. The aspect ratio of the shell measured by parameter f (ratio between the inner and outer radii) does not influence significantly the dimensionless onset times as long as the inner sphere is large enough (f 0.2).

Page 29 of 60

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

For linear stability analysis, values of the onset time are much smaller and the obtained slope a is lower especially for free-slip boundary conditions due to the "frozen time" approach. In the case of Rayleigh-Taylor instability, even though the onset times are also underestimated (however significantly less than for LS analysis), the above scaling based on the 3D results remains valid, to a lesser extent. Two views are classically considered for the onset of convection:

(i) the need to reach a critical thickness for the boundary layer as supposed in the LS analysis, (ii) the need for convective instabilities to reach a specific growth rate (larger than the conductive propagation rate) as assumed in the R-T. An important conclusion is that the latter process provides a better description of the first convective instabilities for a variable-viscosity fluid heated from below.

The preferred degree monotonically increases with increasing Rayleigh number Ra and geometrical factor f in the isoviscous case. In the case of temperaturedependent viscosity, the phenomenon is more complex. However, in general, the low degrees are preferred for high viscosity contrasts and low Rayleigh numbers.

Application of these scaling relationships to Venus, the Earth, and Mars

shows that values of the onset times vary by more than three orders of magnitude for reasonable parameter ranges.

∇ ′ • ṽ′ = 0, (A1) Ra g(r ′ ) g 0 ∇ ′ θ′ × e r + µ ′ ∇ ′2 (∇ ′ × ṽ′ ) + ∂µ ′ ∂r ′ e r × ∇ ′2 ṽ′ + ∇ ′ × e r • ∇ ′ ṽ′ + (∇ ′ ṽ′ ) T + + ∂ 2 µ ′ ∂r ′2 e r × e r • ∇ ′ ṽ′ + (∇ ′ ṽ′ ) T = 0 (A2) ∇ ′2 θ′ -ṽr ∂θ ′ 0 ∂r ′ = ∂ θ′ ∂t ′ , (A3)
where Eq. ( A2) is a curl of momentum Eq. ( 2) with depth-dependent viscosity The equations (A1-A4) lead to the following set of the equations for poloidal 822 scalar P and the temperature field θ′ :

823 -Ra g(r) g 0 θ′ r ′ = µ∇ ′ 4 P + ∂µ ∂r ′ 2 ∂ 3 ∂r ′ 3 P + 4 r ′ ∂ 2 ∂r ′ 2 P - - 2 r ′ 2 ∂ ∂r ′ P + 2 r ′ 3 P + 2 r ′ 2 B ∂ ∂r P - 2 r ′ 3 BP + + ∂ 2 µ ∂r ′ 2 ∂ 2 ∂r ′ 2 P - 2 r ′ 2 P - 1 r ′ 2 BP (A5) ∂ θ′ ∂t ′ = ∇ ′2 θ′ - 1 r ∂θ ′ 0 ∂r ′ BP, ( A6
)
where

824 1 r 2 B• = 1 r 2 sin ϑ ∂ ∂ϑ sin ϑ ∂ ∂ϑ • + 1 r 2 sin 2 ϑ ∂ 2 ∂ 2 ϕ • = = ∇ 2 • - ∂ 2 ∂r 2 • - 2 r ∂ ∂r • . (A7)
In order to solve the linear stability problem, the spectral decomposition of

-Ra g(r ′ ) g 0 Θ lm r ′ = = µ ′ d 4 dr ′ 4 + 4 r ′ d 3 dr ′ 3 - 2l(l + 1) r ′ 2 d 2 dr ′ 2 + l 2 (l + 1) 2 -2l(l + 1) r ′ 4 P lm + + dµ ′ dr ′ 2 d 3 dr ′ 3 + 4 r ′ d 2 dr ′ 2 - 2l(l + 1) + 2 r ′ 2 d dr ′ + 2l(l + 1) + 2 r ′ 3 P lm + A c c e p t e d M a n u s c r i p t + d 2 µ ′ dr ′ 2 d 2 dr ′ 2 + l(l + 1) -2 r ′ 2 P lm , (A8) d 2 dr ′ 2 + 2 r ′ d dr ′ - 1 r ′ 2 l(l + 1) Θ lm + 1 r ′ l(l + 1)P lm dθ ′ 0 dr ′ = 0. (A9)
These equations are solved in a spherical shell bounded by spheres with radii

834 r ′ = f 1-f = r b d (inner boundary) and r ′ = 1 1-f = rt d (
P ′ lm r ′ = 1 1 -f = P lm r ′ = f 1 -f = 0 and d dr P ′ lm r ′ = 1 1 -f = d dr P ′ lm r ′ = f 1 -f = 0 for no-slip or d 2 d 2 r P ′ lm r ′ = 1 1 -f = d 2 d 2 r P ′ lm r ′ = f 1 -f = 0 for free-slip. (A10)
The critical Rayleigh number is then obtained solving the condition for in-

842

In our case, the basally heated problem with θ ′ (r ′ , t

′ = 0) = 0, θ ′ 0 r ′ = f 1-f , t ′ = 843 1 and θ ′ r ′ = 1 1-f , t ′ =
846 θ ′ 0 (r ′ , t ′ ) = f 1 -f 1 r ′ - f 1 -f r ′ -f 1-f r ′ - - 2 r ′ π ∞ n=1 f 1 -f 1 n sin nπ r ′ - f 1 -f exp -n 2 πt ′ . (A11)
For small times, the relationship (A11) can be approximated by the solution for 847 an infinite region bounded internally by a sphere (Carslaw and Jaeger, 1959):

848 θ ′ 0 (r ′ , t ′ ) = f 1 -f 1 r ′ erfc r ′ -f 1-f 2 √ t ′ . (A12)
A c c e p t e d M a n u s c r i p t

Hence, for a given time, we obtain the temperature profile and then compute the critical Rayleigh number and its corresponding degree.

The example of the onset times obtained for a constant viscosity a vis = 0, f = 0.55 and free-slip conditions is shown in Fig. A1 as a function of Rayleigh number. For times lower than ≈ 0.01 (Fig. A1, solid line), the onset time varies almost linearly with Rayleigh number in the log-log scale and its slope is

t ′ ∝ Ra -1.02
. The degree of the first instability is growing from l = 3 to l = 4

with increasing Rayleigh number. For t ′ > 0.01, a sudden change of the slope is observed. This effect can be explained by the presence of the upper boundary.

The temperature profile can be described by the time-dependent solution of pure conduction in the infinite region bounded internaly by sphere (Eq. A12) for low t ′ . The critical Rayleigh number thus varies with t ′ . If the hot front reaches the upper boundary (the temperature based on Eq. ( A12) is non-zero there), the temperature profile should be described by Eq. (A11). For high t ′ , the temperature profile becomes very close to the steady-state conduction solution.

The critical Rayleigh number for the onset of convection thus reaches (for t ′ > 0.5) the asymptotic value corresponding to the classical defition (i.e. the value needed to destabilize a steady-state conductive profile): in this example, the value of the classical critical Rayleigh number is 712 and the preferred degree is 3 (see e.g. [START_REF] Ratcliff | Steady tetrahedral and cubic patterns 1068 of spherical shell convection with temperature-dependent viscosity[END_REF]. As expected, the dependency of onset time for higher degrees (8, 12, 16 and 20) also illustrated in Fig. A1 (thin solid lines) also display this asymptotic behavior.

B. Rayleigh-Taylor instability analysis

In this approach, we define the onset time as the time when the maximum growth velocity ξ of the Rayleigh-Taylor instability is equal to the propagation velocity of boundary layer ( ḣ) (e.g. [START_REF] Zaranek | The onset of convection in fluids with strongly temperature-dependent viscosity cooled from above with implications for planetary lithospheres[END_REF][START_REF] Ke | Plume formation in strongly temperaturedependent viscosity fluids over a very hot surface[END_REF]:

ξ = ḣ. (B1)
Page 34 of 60

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

Therefore, the growth of the R-T instability in the spherical shell needs be determined. The shell is bounded by r b (inner radius) and r t (outer radius)

and a simple two layers model is considered. In this model, the less-denser fluid layer (lower layer) is described by thickness h, viscosity µ 1 and density ρ 1 , the denser fluid (upper layer) is characterized by the viscosity µ 2 and the density ρ 2 = ρ 1 + ∆ρ. Thus, the internal interface between the layers is located at r i = r b + h. The interface between the layers is presumed to be initially infinitesimally perturbed with amplitude ξ 0 . This initial perturbation ξ 0 of the internal boundary is supposed to grow exponentially with time

ξ = ξ 0 exp st, ( B2
)
where s is the growth rate of the instability.

The growth rate is computed as follows [START_REF] Ribe | Analytical approaches to mantle dynamics[END_REF]: In each layer, the flow in the i-th layer follows the equation of motion for the constant viscosity and without body force

∇ 4 P i = 0, (B3)
where P i is a poloidal scalar (see Eq. A4) representation of the velocity v i = (v i r , v i ϑ , v i ϕ ) within the i-th layer. On the boundaries, either no-slip or free-slip is prescribed. On the internal interface, the continuity of the normal and tangential velocity and of the tangential stresses and discontinuity of the normal stress is prescribed. Supposing spectral decomposition P = ∞ l=0 l m=-l exp (st)P lm (r)Y lm (ϑ, ϕ)

and the orthonormality of the spherical harmonic functions, the boundary conditions on the outer boundaries give

P 2 lm (r t ) = P 1 lm (r b ) = 0 and d dr P 2 lm (r t ) = d dr P 1 lm (r b ) = 0 for no-slip or d 2 d 2 r P 2 lm (r t ) = d 2 d 2 r P 1 lm (r b ) = 0 for free-slip. (B4)
A c c e p t e d M a n u s c r i p t

On the internal interface (r i ), the linearized conditions for the continuity of 896 the normal and tangential velocities are prescribed

897 [v r ] + -= P 2 lm (r i ) -P 1 lm (r i ) = 0 898 [v ϑ ] + -= [v ϕ ] + -= d dr P 2 lm (r i ) - d dr P 1 lm (r i ) = 0. (B5)
The linearized continuity of the tangential stresses gives:

899 [σ rϑ ] + -= [σ rϕ ] + -= 0 µ 2 d 2 dr 2 - l(l + 1) -2 r 2 P 2 lm (r i ) -µ 1 d 2 dr 2 - l(l + 1) -2 r 2 P 1 lm (r i ) = 0 . (B6)
The linearized discontinuity of the normal stress is computed as

900 [-p + σ rr ] + -= -∆ρgξ µ 2 r d 3 dr 3 + 3 d 2 dr 2 -3 l(l + 1) r d dr + 3 l(l + 1) r 2 P 2 lm (r i ) - -µ 1 r d 3 dr 3 + 3 d 2 dr 2 -3 l(l + 1) r d dr + 3 l(l + 1) r 2 P 1 lm (r i ) = = -∆ρgξ. (B7)
Using linearized kinematic condition and Eq. (B7), we can write

901 ξ = v r (r i ) ξ = - s[-p + σ rr ] + - g∆ρ s[-p + σ rr ] + - g∆ρ + v r = 0, s[-p + σ rr ] + - g∆ρ - 1 r l(l + 1)P = 0. (B8)
Considering the Green function of the of the poloidal scalar in the i-th layer

902 P i l (r, r i ) = A i n r l + B i n r -l-1 + C i n r l+2 + D i n r -l+1 , (B9)
A c c e p t e d M a n u s c r i p t

we get 8 unknowns of 8 equations (B4, B5, B6 and B8). And the growth rate s 903 for the given degree l is computed then from the condition of the existence of 904 the solution. The results for the thin layer limit (ǫl ≪ 1) agrees well with the 905 results obtained by [START_REF] Ribe | The global hotstpot distribution and instability 1074 of D[END_REF]. A more detailed description 906 of the growth of R-T instability in the Cartesian coordinates can be found 907 e.g. in [START_REF] Ribe | Analytical approaches to mantle dynamics[END_REF], chap. 7.04.9.1.

908

In order to compare the results of the Rayleigh-Taylor instability with the 3D 909 convection results, the dimensionless growth rate s ′ may be related to Rayleigh 910 number as follows [START_REF] Ke | Early transient superplumes and the origin of the Martian crustal dichotomy[END_REF] 911

s ′ = d 2 κ s = d 2 κ ∆ρgh µ 2 s = f ǫ∆µ 1 -f Ra 0 s, ( B10
)
where ǫ = h /r b , ∆µ = µ2 /µ1 is a viscosity contrast and s is a function dependent 912 on l, ∆µ, ǫ and f .

913

In order to compare these results with the 3D convection experiments, the We restrict ourselves only to small times (t ′ < 0.01 for f ≈ 0.5) using 919 the relationship (A12)

920 t ′ = h ′ 2 4 1 erfc -1 θ ′ L 1 + h ′ 1-f f 2 , (B11)
leadint to an implicit equation for h ′ solved numerically. The velocity of the 921 propagation of the boundary ḣ′ can be computed for given t ′ and h ′ as follows

922 ḣ′ = 2 √ π exp h ′2 4t ′ h ′ 4t ′ 3 /2 θ ′ L 1-f f + 2 √ πt ′ exp h ′2 4t ′ . (B12)
Apparently, both the thickness h ′ and its time derivative depend on the defini- As already mentioned, we are looking for the time when the condition (B1)

926 ξ 0 s ′ exp(s ′ t ′ ) = ḣ′ (B13)
is satisfied for the fastest growing degree l and for given Rayleigh number Ra 0 , 927 geometrical factor f and viscosity contrast between the layers defined by ∆µ.

928

In order to compute the time, the bisection method in time is used (Press et al., 929 1992).

930

Another key issue for this method is the estimation of ∆µ. Supposing the 931 relationship 4, the viscosity contrast is defined as follows

932 ∆µ = exp (-θ 2 a vis ) exp (-θ 1 a vis ) = exp ((θ 1 -θ 2 )a vis ) = exp (∆θa vis ) , (B14)
where θ 1 and θ 2 are values of temperature representing lower and upper layer.

933

Several solutions are proposed here: either with a prescribed (constant) value 934 of ∆θ or using the actual averages of both layers to define θ 1 and θ 2 . The results corresponding to ξ 0 = 10 -3 , 10 -4 and 10 -5 , respectively. Hence, the slope is influenced only weakly by the initial perturbation.

In order to compare with the 3D convection results, the value of the initial perturbation can be roughly estimated as follows: the grid resolution in the vertical direction is typically dr ≈ 0.02 for 64 cells. Supposing a maximal amplitude of 10 -3 for the stochastic temperature perturbations, the undulation of the isothermic surface could be approximately 2 • 10 -5 at the beginning of the simulation. This value is the one we select for the prescribed initial amplitude ξ 0 for the R-T analysis. Note that the nonlinear effect of the temperature diffusion on the estimate of ξ 0 should also be taken into account-for decreasing Rayleigh number the onset time is increasing, hence the diffusion effect grows and the initial value of the amplitude should decrease. This effect was however neglected.

Nevertheless, as we show here, the value of ξ 0 influences the slope of the curve only weakly.

A c c e p t e d M a n u s c r i p t and R-T instability analysis (dashed line) for ξ 0 = 10 -5 and its corresponding degrees. A c c e p t e d M a n u s c r i p t

List of Tables Table 1 List of onset times for the numerical experiments corresponding to free-slip boundary conditions; a.i.p-amplitude of the initial perturbations; nb-number of blocks; l-estimated preferred degree at onset time.

Table 2 List of onset times for the numerical experiments corresponding to no-slip boundary conditions; a.i.p-amplitude of the initial perturbations; nb-number of blocks; l-estimated preferred degree at onset time.

Table 3 Summary of the proposed scaling (using least squares fit)

for both free-slip and no-slip boundary conditions and all methods; 3D-3D numerical solution, R-T-Rayleigh-Taylor instability analysis, LS-linear stability analysis.

Table 4 Reference parameters for Venus, the Earth and Mars.

Page 49 of 60 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 1 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 2 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 3 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 4 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 5 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 6 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 7 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 8 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 9 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 10 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 15 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 20 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] 23 -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

f [1] -5 -4 -3 -2 t' [1]
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 f [1]

  . Values of (B, b) are listed in Tab. 3. Note that linear stability analysis consider only depth-dependent viscosity, and R-T instability analysis uses an even simpler description of the viscosity variations based on a two layered viscosity stratification. For this reason, specific differences arise between 3D

A

  but whose best scaling fit are reported in Tab. 3) show a significantly different behavior, closer to the observed slopes for 3D results and R-T. Results obtained for a constant value of the bottom Rayleigh number Ra b are presented in Fig. 3b. The surface Rayleigh number Ra 0 decreases with increasing a vis and hence the onset time increases. For relatively high a vis , the onset time is delayed and the heat front reaches the upper boundary before the 13

  a(f ) and t ′ ∝ exp (b(f )a vis )). The absolute value of function a(f ) is expected to increase with f : by definition of the time analogous to the onset of convection in the R-A c c e p t e d M a n u s c r i p t T framework, smaller values of Ra imply larger thicknesses for the lower layer. Since the growth of the thermal boundary layer decreases more rapidly with time for low f than for high f , critical thicknesses of the boundary layer take more time to develop (see Eq. A12). Our results were obtained for Rayleigh numbers larger than 10 3 (corresponding roughly to the critical value for convection to occur within the shell) and smaller than 10 9 . The following values of the power a(f ) are obtained for free-slip conditions and isoviscous fluid: a(0.01) = -0.74, a(0.1) = -0.72, a(0.2) = -0.71, a(f > 0.3) = -0.70. In the case of the function b(f ), the dependency for free-slip conditions has been found for Ra = 10 6 and a vis ranging between 0 and 10: b(0.01) = -0.05, b(0.1) = -0.13, b(0.2) = -0.15, b(f > 0.3) = -0.15.

  ), the results obtained for a constant bottom Rayleigh number Ra b = Ra 0 exp (a vis ) may be interpreted as follows Ra a b exp ((b -a)a vis ) .

  and the variations of this normalized onset time with a vis are reported for all the calculations performed in this study. Note that the value of θ * used in this normalization is based solely on results obtained for a prescribed value of the surface Rayleigh number. Fig.5shows that the proposed scaling also describes the other 3D calculations (i.e. with other values of Ra 0 ): a leastsquares fit with the function t′ = α exp (βa vis ) returns values of β smaller than 10 -2 for both free-slip and no-slip and the scattering of 3D convection runs around the interpolated curve is low, see Fig. 5. As expected, the R-T and LS results for constant surface Rayleigh number Ra 0 show a similar behavior even though the normalized onset times are underestimated. The normalized onset Page 18 of 60 A c c e p t e d M a n u s c r i p t

  Fig. 12 in their study indicates that an exponential scaling such as the one proposed here (Eq. 8) describes correctly the experiments: onset times for convection increase linearly with ∆µ in a log-log scale. In the transient set-up of the present study, onset times behave similarly. The hot (less viscous) thermal boundary layer above the inner interface first thickens. Convective instabilities occur so that the whole thermal boundary layer takes part in the large-scale motion. For this reason, we propose an alternative dependency on a vis (cf. Eq. 8) equivalent to formulation A c c e p t e d M a n u s c r i p t in Eq. 10: t ′ ∝ (Ra * ) a , where Ra(θ * ) = ρ0g∆T αd 3 κµ(θ * )

  s c r i p t by temperature gradients), i.e. a stratified regime with a "flat" interface. Indeed while this layer should remain stable for a significant fraction of the age of the planet, if radiogenic heat sources are present with concentrations ten times larger or more than chondritic[START_REF] Elkins-Tanton | Early magnetic field and magmatic activity on Mars from magma ocean cumulate overturn[END_REF], the increase of temperature due solely to these internal heat sources may easily reach several hundreds of Kelvins in a few tens of million years and lead to hot convective instabilities in the less dense layer above.

  part of the planet's mantle[START_REF] Elkins-Tanton | Early magnetic field and magmatic activity on Mars from magma ocean cumulate overturn[END_REF]. Fur-594 thermore, core formation on any terrestrial planet raises its temperature due to 595 the conversion of potential energy into heat. The presence of a magma ocean is 596 thus highly probable and it is supported by the geochemical data. The solidifica-597 tion of such a magma ocean is described by[START_REF] Solomatov | Magma oceans and primordial mantle differentiation[END_REF], for example. The 598 temperature profile subsequent to the solidification of the magma ocean could 599 then be close to the solidus value. Again, it has been noted by Elkins-Tanton 600 et al. (2005) that the density stratification resulting from this solidification may 601 not be stable and could promote a gravitational overturn. 602 (i) Internal temperature. For all reasons mentioned above, large ranges for 603 the internal temperature value T 0 need to be considered. For Mars, a temper-604 ature range for the mantle near CMB of (1900, 2300)K is supposed. The lower 605 limit is based on the estimate of the current temperature near CMB in Mars 606

608

  660 km is supposed to be approximately 1900 K (e.g.[START_REF] Ito | Postspinel transformations in the system Mg 2 SiO 4 -Fe 2 SiO 4 and some geophysical implications[END_REF] 609Boehler, 2000). Considering an adiabatic profile, the current temperature near 610 CMB is approximately 2500 -2900 K providing an estimate for the lowest pos-611 sible value for temperature before the onset of convection. The upper bound is 612 supposed to lay around the solidus of the lower mantle, i.e. 4300 K(Holland and 613 Ahrens, 1997;[START_REF] Zerr | Solidus of Earth's deep mantle[END_REF]. For Venus, we assume the same temperature 614 range: Venus is slightly smaller and (at least present-day) surface conditions 615 are different than on the Earth so that the energy budget may differ. However, 616 we suppose that these differences are not substantial especially when compared 617 to the uncertainties of the studied parameters. The thermal evolution within 618 the Earth and Venus is discussed in[START_REF] Stevenson | Magnetism and thermal evolution of the terrestrial planets[END_REF]. The main effect of 619 varying T 0 results in large variations of the viscosity value which we take into 620 account in the following by considering viscosity η 0 as a free parameter. T 0 also 621 influences the value of a vis . We first focus on two specific values of T 0 indicated 622 in Tab. 4 and then consider results obtained for other values of T 0 in the case 623 of the Earth and Venus (Fig.6b-d).

  Temperature contrast. The minimal temperature contrast between the 625 mantle and the core ∆T is limited by the liquidus temperature of the core 626 material. The liquidus of iron near CMB in Mars is estimated to be approxi-627 mately 2100 K(Boehler, 1996). In the Earth, the iron liquidus is around 3000 K 628(Boehler, 1996). The presence of light elements within the core would signifi-629 cantly reduce these values. Supposing the values mentioned above for the man-630 tle temperature before the onset of convection, i.e. 2100 K (Mars) and 3000 K 631 (Earth/Venus), a high temperature increase near CMB is not required. We thus 632 investigate the following range for ∆T : 500-2000 K consistent with the current 633 estimate within the Earth[START_REF] Lay | Core-mantle boundary heat flow[END_REF]. We will see below that larger 634 values of ∆T are likely to be associated with a different dynamical regime than 635 the one investigated here where the associated viscosity contrast is moderate. 636 (iii) Viscosity parameters. The parameters controlling the viscosity value 637 and its variations near the bottom of the mantle are another key issue. Due to 638 the lack of constraints, we basically suppose them to be free parameters (com-639 bined together with the internal temperature T 0 ). Indeed, besides temperature, 640 viscosity also depends on the activation energy and volume, grain size and water 641 content. Values of these parameters are subjected to large uncertainties espe-642 cially for the lower mantle material. We thus consider the following viscosity 643 range (10 18 , 10 22 )Pa • s. The lower limit is a viscosity near solidus (Solomatov, 644 2007) in the Earth. The upper viscosity bound is based on a estimate of the vis-645 cosity in the lower mantle in the Earth (e.g.

648

  the activation parameters Q * = E * + pV * and the temperature T 0 are rather 649 uncertain. In order to compute the activation enthalpy Q * = E * + pV * , we use 650 the parameters listed in Tab. 4 and consider a pressure value corresponding to 651 mantle close to the CMB. For Mars, the values are based on what is proposed 652 for the Earth's upper mantle material (Korenaga and Karato, 2008). For the 653 Earth and Venus, the activation energy and volume are rather uncertain for the 654 lower mantle material under high pressure; we use values based on the results of 655 A c c e p t e d M a n u s c r i p t Yamazaki and Karato (2001). This leads to Q * = E * + pV * ≈ 380 kJ • mol -1 for 656 Mars and ≈ 540 kJ • mol -1 for the Earth and Venus. However, the value of the 657 activation volume is rather poorly constrained. Korenaga and Karato (2008) 658 report a value of the activation volume of 6 ± 5 cm 3 • mol -1 , hence, values in the 659 range (1, 11) cm 3 • mol -1 are plausible. Supposing these values, the activation 660 enthalpy near CMB may vary between 280 and 480 kJ • mol -1 for Mars and 661 between 340 and 1700 kJ • mol -1 for the Earth and Venus. Similarly to high 662 temperature contrasts ∆T , the upper bounds of the estimated enthalpy range 663 lead to a different dynamical regime which is not addressed here. 664 Onset of convection Fig. 6a-b shows the onset times of the first convec-665 tive instability for parameter values corresponding to Venus or the Earth and 666

  678here. However, if we had considered smaller values for T 0 or higher values of Q * , 679 the empty region may have represented a larger fraction of the studied range of 680 values for ∆T . 681 The onset times for the studied range of parameters vary between 0.1 Myr 682 and 700 Myr Fig. 6a-b, the values for the Earth and Venus (< 400 Myr) being 683 only slightly smaller than the ones obtained for Mars (> 0.3 Myr). The extreme 684 values, although they are very different, are reached only for small areas within

  is the determination of a preferred wavelength of the anomalies associated with the onset of the convection. Even though our estimate of the degree (wavelength) for 3D runs in one single block is limited by the common factor 4, it was demonstrated that the low degree instabilities are favored by rather low viscosity and/or high viscosity contrast. For the investigated range of the parameters (internal viscosity µ 0 and temperature contrast ∆T ), i.e. supposing Rayleigh number Ra 0 10 7 (for the Earth and Venus)and Ra 0 5 • 10 5 (for Mars) and addressing intermediate viscosity contrast a vis 10, it is highly improbable to achieve low degree-anomalies. Based on the 3D runs in six block, we roughly estimate that the lowest degree in our case is approximately 10.The onset of convection in the martian mantle as a the fluid heated from below was investigated by[START_REF] Ke | Early transient superplumes and the origin of the Martian crustal dichotomy[END_REF] as a possible cause for the hemispheric crustal dichotomy. The authors report the possibility of a transient low-degree superplume resulting from the destabilization of a deep layer where small-scale convection occurs. Large viscosity contrasts (∆µ = 10 6 ) and relatively low viscosity values need to be considered to promote this regime corresponding to the empty area in Fig. 6a. The present study does not rule out the possibility of such a large viscosity contrast regime, especially for early Mars (this regime does not appear in the parameter space we considered for the Earth and Venus) since the values of the activation parameters are subjected to large uncertainties.Following a first study based on the dynamics of a Rayleigh-Taylor instability(Zhong and Zuber, 2001),[START_REF] Roberts | Degree-1 convection in the Martian mantle and 1076 the origin of the hemispheric dichotomy[END_REF], investigate other mechanisms possibly leading to the formation of low-degree instabilities early in the martian mantle. The first one is related to the presence of an endothermic phase transition deep within the mantle and leads to low-degree convection only for a restricted region in the parameter space. Furthermore, this pattern takes a long time to develop which makes it unrealistic. The other mechanism investigated by[START_REF] Roberts | Degree-1 convection in the Martian mantle and 1076 the origin of the hemispheric dichotomy[END_REF] is the presence of a viscosity layer-

815µ′

  = µ ′ (r) and equation (A3) is a linearized energy equation for the depth-816 dependent reference temperature profile θ ′ 0 = θ ′ 0 (r). This set of equations is 817 solved for the solenoideal (Eq. A1) velocity field represented by a poloidal scalar 818 (toroidal part of the field is identically equal to zero see e.g. Schubert et al. 819 (2001)) and the temperature perturbations θ′ . The poloidal scalar P is defined 820 by 821 ṽ′ = ∇ ′ × (r ′ e r × ∇ ′ P ) .

  exp (st)P lm (r)Y lm (ϑ, ϕ) and of temperaexp (st)Θ lm (r)Y lm (ϑ, ϕ) are used, Y lm are 827 fully normalized spherical harmonics (e.g.[START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF]. Further-828 more, we restrict ourselves to the eigenvalue s = 0. Hence, the "frozen time" 829 approximation is considered and the propagation of the temperature front is 830 not taken into account. Under these conditions and due to the orthogonality 831 of the spherical harmonics functions, considering Bf (r)Y lm = -l(l + 1)f (r)Y lm , 832 the equations (A5-A6) can be rewritten into 833

  837stability (at least one non-trivial solution of the system (A8-A9) for boundary 838 condition (A10) exists) with a method similar to Chandrasekhar (1961). The re-839 sults for the conductive profile reproduce the results obtained by Chandrasekhar 840 (1961) for constant viscosity and Ratcliff et al. (1996) for temperature dependent 841 viscosity.

914

  time dependency of the thickness of the lower layer h = r i -r b (corresponding 915 to the hot thermal boundary layer in the numerical experiments) and its time 916 derivative ḣ need to be prescribed. The interface between the layers is defined 917 by a specific value of the temperature θ ′ L taking into account purely conductive 918 solution.

  h ′ = √ πt ′ for f → 1.

935A

  Fig. B1b. As expected, the onset times for given Rayleigh number decreases 948

  Figure 1Example of the determination of the onset time for Ra = 10 6 , f = 0.55, a vis = 0 and free-slip conditions; the deviation of the horizontally averaged temperature (solid line) from the conduction temperature profile (dashed line) for the six discrete layers closest to core-mantle boundary is reported. The onset time corresponding to a 0.1% deviation is indicated by an arrow.

Figure 2

 2 Figure 2 Onset time as a function of Rayleigh number Ra for constant f = 0.55, a vis = 0 for a) free-slip and b) no-slip. Results of 3D numerical simulations (solid circles) and associated fit (thin solid line); linear stability analysis (solid and dash-dotted line)and R-T instability analysis (dashed line) for ξ 0 = 10 -5 and corresponding degrees.

Figure 3

 3 Figure 3Onset time as a function of a vis for f = 0.55 and free-slip a) Ra 0 = 10 6 and b) Ra b = 10 6 . Results of 3D convection experiment (solid circles) and associated fit (thin solid line), linear stability analysis results (solid line), R-T instability analysis results(dashed line for ∆θ=1, the lowest onset times, 0.8, 0.6, 0.4 and 0.2, dash-dotted line for temperature average (∆ θ), see Eq. (B14)and text in App. B for further detail) for ξ 0 = 10 -5 and corresponding degrees. In case of constant bottom Rayleigh number (panel b), solid parts of the lines denote results for t ′ < 0.01.

Figure 4

 4 Figure 4Onset time depending on spherical shell geometry f for constant Ra = 10 6 , a vis = 0 for a) free-slip and b) no-slip conditions, the results of 3D convection experiment (solid circles), the linear stability analysis results (solid and dash-dotted line)

Figure 5

 5 Figure 5 Scaled onset time t ′ • Ra -a (θ * ) as a function of viscosity parameter a vis . Results of 3D numerical solution (solid circles-Ra = const, solid triangles-Ra 0 = 10 6 , solid squares-Ra b = 10 6 ) and the associated regression (solid line), results for R-T instability analysis (thick lines, ∆µ based on temperature averaging) and linear stability analysis (thin lines) for Ra 0 = 10 6 (dashed line) and Ra b = 10 6 (dash-dotted line), results for a) free-slip and b) no-slip boundary conditions.

Figure 6

 6 Figure 6 The estimate of the onset time based on 3D convection experiments and reference temperatures T 0 for a) Mars and b) the Earth and Venus; c-d) the dependency of onset times on T 0 in the Earth and Venus.

Figure A1

 A1 Figure A1Linear stability analysis results for f = 0.55, a vis = 0 and free-slip conditions. Minimal onset time (solid line-the solution corresponding to the approximation (A12), dashed line-the solution corresponding to (A11)) and corresponding degree l onset times for degrees 8, 12, 16 and 20 (thin solid lines) are also presented.

Figure

  Figure B1 Rayleigh-Taylor instability analysis for f = 0.55, a vis = 0 and free-slip, a) the dependency of the onset time on the degree l for different Rayleigh numbers Ra and ξ 0 = 10 -5 , b) the dependency of onset time on the Rayleigh number for different initial amplitude ξ 0 and its corresponding degrees.

  Figure 1:

  Fig. 1

  Fig. 2

  Fig. 3

  Figure 4:

  Fig. 4

  Fig. 5

  Figure 6:

  Fig. 6

  Fig. A1

44 •

 44 Fig. B1

  outer boundary). On the

	835
	836	impermeable boundaries, either no-slip or free-slip condition are prescribed

  Boehler, R., 1996. Melting temperature of the Earth mantle and core: Earth's Choblet, G., Sotin, C., 2001. Early transient cooling of Mars. Geophys. Res.

	978
	979	thermal structure. Annu. Rev.Earth Planet. Sci 24, 15-40.
	980	Boehler, R., 2000. High-pressure experiments and the phase diagram of lover
	981	mantle and core materials. Reviews of Geophysics 38, 221-245.
	982	Breuer, D., Yuen, D., Spohn, T., 1997. Phase transitions in the Martian mantle:
	983	Implications for partially layered convection. Earth Plan. Sci. Let. 148, 457-
	984	469.
	985	Canright, D., Morris, S., 1993. Buoyant instability of a viscous film over a
	986	passive fluid. J. Fluid Mech. 255, 349-372.
	987	Canup, R., 2004. Simulations of a late lunar-forming impact. Icarus 168, 433-
	988	456.
	989	Carslaw, H., Jaeger, J., 1959. Conduction of heat in solids, 2nd Edition. Oxford
	990	University Press.
	991	Chandrasekhar, S., 1961. Hydrodynamic and hydromagnetic stability. Oxford
	992	University Press.
	993	Choblet, G., 2005. Modelling thermal convection with large viscosity gradients
	994	in one block of the 'cubed sphere'. J. Comput. Phys. 205, 269-291.
	995	Choblet, G., Čadek, O., Couturier, F., Dumoulin, C., 2007. OEDIPUS: A new
	996	tool to study the dynamics of planetary interiors. Geophys. J. Int. 170, 9-30.
	997	Choblet, G., Sotin, C., 2000. 3D thermal convection with variable viscosity:
	998	can transient cooling be described by a quasi-static scaling law? Phys. Earth
		Planet. Inter. 119, 321-336.
	1000
	1001	Let. 28(15), 3035-3038.
	1002	Conrad, C., Molnar, P., 1999. Convective instability of a boundary layer with
	1003	temperature-and strain-rate-dependent viscosity in terms of available buoy-
		ancy. Geophys. J. Int. 139, 51-68.

999

Table 3 :

 3 Tab. 3

thickness of the mantle d 2900 km 1700 km geometrical factor f 0.55 0.50 activation parameters for dry diffusion creep in olivine [START_REF] Korenaga | A new analysis of experimental data on olivine rheology[END_REF] activation energy E * 261 kJ • mol -1 activation volume V * 6 cm 3 • mol -1 activation parameters for diffusion creep in perovskite, based on [START_REF] Yamazaki | Some mineral physics constraints on the rheology and geothermal structure of Earth's lower mantle[END_REF] activation energy E * 200 kJ • mol -1 activation volume V * 2.5 cm 3 • mol -1

Table 4:

Tab 4.