Alboussière

A Morison
email: adrien.morison@ens-lyon.fr

S Labrosse

R Deguen

T Alboussière

Timescale of overturn in a magma ocean cumulate

Keywords: magma ocean, overturn, mantle dynamics, linear stability

Timescale of overturn in a magma ocean cumulate

Introduction

A common scenario considered for the formation of terrestrial planets is the crystallization of a global magma ocean from the bottom-up, because the liq-uidus of silicate magmas increases with pressure more steeply than the isentropic temperature, at least at low to moderate mantle pressure [START_REF] Andrault | Solidus and liquidus profiles of chondritic mantle: Implication for melting of the Earth across its history[END_REF][START_REF] Fiquet | Melting of peridotite to 140 gigapascals[END_REF]Thomas and Asimow, 2013;[START_REF] Boukaré | Thermodynamics of the MgO-FeO-SiO2 system up to 140 GPa: Application to the crystallization of Earth's magma ocean[END_REF]. The crystallization of the surface magma ocean is expected to be rapid, around 1 Myr (e.g. [START_REF] Abe | Thermal and chemical evolution of the terrestrial magma ocean[END_REF][START_REF] Lebrun | Thermal evolution of an early magma ocean in interaction with the atmosphere[END_REF]. This has led several authors to assume convection in the solid part of the crystallizing mantle does not start until the mantle is entirely crystallized (e.g. [START_REF] Hess | A model for the thermal and chemical evolution of the Moon's interior: implications for the onset of mare volcanism[END_REF][START_REF] Abe | Thermal and chemical evolution of the terrestrial magma ocean[END_REF][START_REF] Parmentier | Gravitational differentiation due to initial chemical stratification: origin of lunar asymmetry by the creep of dense KREEP?[END_REF][START_REF] Elkins-Tanton | Magma ocean fractional crystallization and cumulate overturn in terrestrial planets: Implications for Mars[END_REF][START_REF] Elkins-Tanton | Early magnetic field and magmatic activity on Mars from magma ocean cumulate overturn[END_REF]Zhang et al., 2013).

However, this assumption deserves scrutiny since the compositional and thermal structure of the mantle after complete crystallization could be widely different if solid-state convection does set in during its crystallization.

Two processes might lead to the destabilization of the solid mantle during its crystallization. First, assumming fractional crystallization, the surface magma ocean gets enriched in incompatible elements. As a secondary result, the new solid formed at the solid/liquid boundary gets richer and richer in these elements as crystallization progresses. Iron is such an element and its abundance is such that it affects significantly the density of both the solid and the liquid. The solid formed at the end of the crystallization is richer in iron than the solid formed at the beginning of the crystallization, leading to an unstable setup with material denser at the top than at the bottom of the solid mantle.

The second process that can further destabilize the solid mantle is the temperature gradient in the solid. The solidus temperature increases with pressure, and is steeper than the isentropic temperature profile. Assuming the temperature in the solid stays close to the solidus, the resulting profile is hence unstable. This effect is enhanced by fractional crystallization and the associated enrichment of the solid in incompatible elements: their presence further decreases the solidus temperature and the compositional gradient discussed above induces an even steeper solidus.

Numerical simulations including these processes suggest it is possible for solid-state convection to set in prior to the entire crystallization of the surface magma ocean (e.g. [START_REF] Maurice | Onset of solid-state mantle convection and mixing during magma ocean solidification[END_REF][START_REF] Boukaré | Timing of mantle overturn during magma ocean solidification[END_REF]. Whether convection in the mantle starts during or after the crystallization of the surface magma ocean is found to have profound implications on the preservation of compositional heterogeneities as well as the dynamics of the mantle [START_REF] Ballmer | Reconciling magma-ocean crystallization models with the present-day structure of the Earth's mantle[END_REF]Tosi et al., 2013). These results further confirm the need to assess the parameters controlling the onset of convection in the primitive mantle.

A dynamical feature of the solid cumulate in contact with a magma ocean that has not been accounted for in the past studies is the possibility of exchange of matter at the boundary between the solid and the ocean via melting and freezing. We use a boundary condition developed for the inner core boundary [START_REF] Deguen | Thermal convection in Earth's inner core with phase change at its boundary[END_REF] to take this effect into account. This boundary condition is expected to have important effects on the convection pattern and heat flux as well as the timescale with which convection sets in [START_REF] Deguen | Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries[END_REF][START_REF] Labrosse | Rayleigh -Bénard convection in a creeping solid with melting and freezing at either or both its horizontal boundaries[END_REF].

Our aim is to assess how the timescale at which convection starts in the solid cumulate compares with the time needed to crystallize a surface magma ocean. Different scenarios are explored to determine the parameters controlling the onset of convection in the magma ocean cumulate. We consider the case where fractional crystallization happens during the entire cooling history of the magma oceans as well as the case where no compositional fractionation occurs.

We explore the classical case for which no matter crosses the boundary between the magma ocean and the solid cumulate, and also the case with a boundary that allows matter transfer accross it. The study is applied to the Earth, Mars, and the Moon.

Methods

We consider a mantle that is initially fully molten and crystallizes from the bottom or some intermediate depth upward. The goal of the present study is to determine the timescale for convection to start in the solid part of the mantle as the magma ocean crystallizes.

For the sake of simplicity, we assume the compaction length to be small and neglect the thickness of a mush layer at the phase change interface. Matter on one side of the boundary is entirely liquid while matter on the other side is entirely solid. We nonetheless allow for compositional fractionation to occur as the mantle crystallizes. The temperature at the solid/liquid boundary is denoted T m and referred to as the melting temperature.

Depending on how the temperature profile in the magma ocean compares with the profile of the melting temperature, two situations can occur. Either the solidification of the ocean progresses from the bottom up, or the solidification starts from an intermediate depth leading to a setup in which the solid part of the mantle is surrounded by two magma oceans. In this second scenario, the crystallization of the surface magma ocean (SMO) is thought to be a lot faster than the crystallization of the basal magma ocean (BMO) [START_REF] Labrosse | A crystallizing dense magma ocean at the base of the Earth's mantle[END_REF].

We assume the solid mantle is a spherical shell of internal radius R -and external radius R + . Since the crystallization of the BMO is much slower than the crystallization of the SMO, we assume R -to be constant even for the case where the solid shell is surrounded by two magma oceans. The presence or absence of a BMO however affects the boundary condition applied at the bottom boundary of the solid mantle (see section 2.4).

As the magma ocean cools down, R + increases to reach the total radius of the planetary body, denoted by R T . The temperature at the top boundary of the solid follows the melting temperature. The composition of the solid changes as well with the radius if we assume fractional crystallization occurs. For the sake of simplicity, we only consider fractionation of iron. The mass fraction of FeO, denoted by C, varies between 0 (e.g. Forsterite) and 1 (e.g. Fayalite). Although simplistic, such a model allows us to study the effect of the density gradient due to fractional crystallization on the dynamics of the solid. Figure 1 shows the composition and temperature profiles at two different times. We assume the velocity of the freezing front Ṙ+ does not vary laterally and that the SMO is well mixed, the temperature and compositional fields in the resulting solid

� � � � � � � � ����� ��� � � � � � � � � � � � � � � � � � � ������ � � � � ������� ����� ��� � � � � � � � � � � � � � � � �
Figure 1: Temperature and composition reference profiles. Solid lines are the profiles at time t, dashed lines the profiles at time t + δt. The green area is the solid mantle at time t. The yellow area represents the part of the surface magma ocean (SMO, in red) that has crystallized during δt. All the annotations on the axes are written at time t (see table 1 for the meaning of symbols). Notice how the melting temperature decreases between the two instants owing to the enrichment in iron of the surface magma ocean. The slopes of the curves are exaggerated for readability purpose.

hence only vary with the radial position (as long as no solid-state convection operates).

In this section, we introduce the simple phase diagram we use to compute the resulting temperature and composition profiles in the solid under the assumption that no convection occurs in the solid (section 2.1). This serves as base state which stability against overturning motion is studied. We don't treat the full dynamics of the overturn but compute, using a linear stability analysis, the growth rate of an overturning instability to compare it to the crystallization rate of the magma ocean. The latter is computed using a magma ocean cooling model which gives R + as a function of time, as described in section 2.2.

Composition and temperature reference profiles

Under the assumption that no convection occurs during crystallization, one can determine the resulting temperature and compositional profiles in the cumulate. These profiles are used as reference profiles in order to perform the linear stability analysis (section 2.5).

L Ṙ+ /κ Γ(t)
Thickness of the solid part

L/L M Γ S (t) Thickness of the SMO (R T -R + )/L M Φ ±
Phase change number § 10 -2 ; ∞ We consider a magma ocean crystallizing from some depth R -up to its top radius R T . The mass fraction of the heavy component (FeO) is C(r) in the solid and C l (t) in the liquid, assuming that no diffusion (nor convection) operates in the solid (therefore C does not depend on time) and convection mixes the liquid efficiently (therefore C l depends only on time). At the freezing front, the phase relation is

C(R + (t)) = DC l (t) (2.1)
with D the partition coefficient (considered constant) and R + (t) the timeevolving radius of the freezing interface.

Assuming the magma ocean undergoes fractional crystallization, the composition profile in the cumulate is exponential. At the radial position r it is

C(r) =    C 0 � RT 3 -(R -) 3 RT 3 -r 3 � 1-D if r < R s 1 if r > R s , (2.2) with R s = � (R -) 3 C 1 1-D 0 + R T 3 � 1 -C 1 1-D 0 �� 1/3
(2.3) the value of R + at which C l reaches 1 (see appendix A for more details).

Since the diffusion timescale is much larger than the other time scales considered here, we assume the temperature profile in the cumulate stays close to the melting temperature. We take into account variations of the melting temperature T m due to both the pressure and the composition. A higher concentration in iron leading to a lower melting temperature, the resulting temperature profile in the solid is steeper than a constant-concentration solidus when fractional crystallization is accounted for (Figure 1). The melting temperature T m verifies:

dT m dr = ∂ T m ∂P ∂ P ∂r + ∂ T m ∂C ∂ C ∂r .
(2.4)

With ∂ P ∂r =
-ρg and eq. ( 2.2), one obtains

dT m dr = -ρg ∂ T m ∂P + 3C(1 -D) r 2 R T 3 -r 3 ∂ T m ∂C . (2.5)
For the sake of simplicity, we assume ∂ T m ∂P and ∂ T m ∂C to be constant (see table 1 for values).

We denote T = T -T isen the superisentropic temperature in the solid, with

T isen = T -exp � αg(R --r) C p � (2.6)
the isentropic temperature profile in the solid, with α the coefficient of thermal expansion, g the acceleration of gravity and C p the heat capacity. We assume the variations of α, C p and g with depth to be negligible. The reference superisentropic temperature (denoted T ) gradient is then:

d T dr = -ρg ∂ T m ∂P +3C(1-D) r 2 R T 3 -r 3 ∂ T m ∂C + αg C p T -exp � αg(R --r) C p � . (2.7)

Crystallization time scale

Assuming the temperature profile in the SMO to be isentropic and neglecting variations of α, g and C p with depth, the potential temperature at the surface is:

T p = T + exp � -αg(R T -R + ) C p � .
(2.8)

Note that we are neglecting the temperature drop across the boundary layer at the bottom of the magma ocean. This is justified by the very small viscosity of the magma and the main buoyancy force coming from cooling to the atmosphere at the top surface. [START_REF] King | Heat transfer by rapidly rotating Rayleigh-Bénard convection[END_REF] showed that the scaling law for the heat flux in a rotating fluid (such as the surface magma ocean) depends on how the quantity

Ra S E 3/2 S = αgΔT ν 1/2 κ(2Ω) 3/2
compares to 1, with E S the Ekman number and Ra S the Rayleigh number in the SMO. A conservative lower bound with the thermal expansivity α ∼ 10 -5 K -1 , the gravity g ∼ 10 m/s 2 , the super-isentropic temperature difference ΔT ∼ 1 K, the kinematic viscosity ν ∼ 10 -5 m 2 /s, the thermal diffusivity κ ∼ 10 -6 m 2 /s and the rotation rate Ω ∼ 10 -4 s -1 is

Ra S E 3/2 S ∼ 10 5 � 1.
We then consider the heat flux is not controlled by rotation and scales as Nu = 0.16Ra

2/7 S Γ 6/7 S with Γ S = (R T -R + )/L the dimension-
less thickness of the SMO [START_REF] King | Heat transfer by rapidly rotating Rayleigh-Bénard convection[END_REF]. Note that this scaling does not depend on the Prandtl number in the range of values explored by [START_REF] King | Heat transfer by rapidly rotating Rayleigh-Bénard convection[END_REF], i.e. 1 � Pr � 100. Since Pr ∼ 10 is a reasonable value for a magma ocean, we assume this scaling holds for our study. We neglect variations of Ra S with time and assume the magma ocean behaves like a gray body at its upper surface. Heat flow conservation at the surface gives the following equation for the surface temperature T s :

k(T p -T s ) L M 0.16Ra 2/7 S Γ -1/7 S = εσ(T 4 s -T 4 ∞ ) (2.9)
where T ∞ is the black body equilibrium temperature, σ is the Stefan-Boltzmann constant and ε the emissivity. The emissivity should depend on the atmosphere dynamics and composition (particularly its water content) and vary with time.

Taking this effet into account would require an atmosphere model (e.g. [START_REF] Abe | Thermal and chemical evolution of the terrestrial magma ocean[END_REF][START_REF] Lebrun | Thermal evolution of an early magma ocean in interaction with the atmosphere[END_REF]. For the sake of simplicity, we assume the emissivity to be constant, tuning its value to obtain a crystallization timescale that matches the ones of [START_REF] Lebrun | Thermal evolution of an early magma ocean in interaction with the atmosphere[END_REF]

(see table 1 for values).

As the SMO crystallizes (i.e. R + increases with time), we assume the temperature at the top of the solid mantle T + follows the solidus (eq. (2.5)), and the temperature profile in the SMO follows an isentropic profile. As R + grows, two phenomena produce heat that should be evacuated: the crystallization itself with an associated latent heat L h , and the cooling of the magma ocean.

Assuming this heat is entirely evacuated through radiation in the atmosphere modeled as a gray body, one obtains the following equation:

εσR 2 T (T 4 s -T 4 ∞ ) = ρL h R + 2 dR + dt -ρC p d dt � � RT R + T + exp � αg(R + -r) C p � r 2 dr � .
(2.10) The last term of this equation can be developed (keeping in mind that the lower bound of the integral R + depends on time). This yields the time derivative of T + , which is written as a derivative with respect to R + using the chain rule. One obtains an ordinary differential equation on R + (t) whose numerical integration gives the position of the interface between the solid and the surface magma ocean as a function of time.

Set of dimensionless equations

L = R + -R -, L 2 M /κ, κ/L, ηL 3 /κ, ΔT = T --T +
are used as scales for length, time, velocity, mass and temperature respectively. Note that R + and T + vary with time as the surface magma ocean crystallizes.

L M = R T -R -is the
thickness of the solid mantle once the SMO is entirely crystallized. Note that all scales depend on time except the one for time itself, which is why Γ = L/L M appears in the following equations. The dimensionless radial position is built as 1 + (r -R -)/L so that it is between 1 and 2 at all times. Similarly, the dimensionless temperature is chosen as (T -T + )/ΔT so that it is between 0 and 1 at all times.

Using the same symbols for dimensionless quantities, dimensionless conservation equations of mass, momentum, heat and iron fraction are written as:

∇ • u = 0 (2.11) 0 = -∇p + ∇ 2 u + Ra (Θ -�Θ�) r + Rc (c -�c�) r (2.12) Γ 2 ∂ Θ ∂t + u • ∇(Θ + T ) -∇ 2 Θ = W � (r -1) ∂ Θ ∂r + � ∂ T ∂r � + Θ � (2.13) Γ 2 ∂ c ∂t + u • ∇(c + C) = W (r -1) ∂ c ∂r .
(2.14) Note that these equations are written under the assumption that Ṙ-= 0.

u
Other terms would appear on the right hand side of eqs. (2.13) and (2.14) in the general case involving the crystallization of a basal magma ocean. For Earth's case, we assume the basal ocean crystallizes much slower than the surface ocean, and as such we neglect Ṙ- [START_REF] Labrosse | A crystallizing dense magma ocean at the base of the Earth's mantle[END_REF]. We assume the diffusion of the compositional field is negligible since the diffusion coefficient of composition is much smaller than that of heat. Moreover, diffusion of T is neglected while that of Θ is retained in order to ease the linear stability analysis. This is justified a posteriori by the fact that the diffusion timescale is much longer than the other timescales considered in this study.

Phase change boundary condition

In the classical Rayleigh-Bénard setup, convecting matter arriving near an horizontal boundary forms a topography whose height is limited by the weight viscous forces can sustain. This topography is often neglected and a nonpenetrative boundary condition is assumed at the interface (u r (R + ) = 0). However, in the system studied here, the boundary between the magma ocean and the cumulate is a phase change interface. A topography of the solid with respect to the equilibrium position of the solid/liquid interface can then be eroded by melting or freezing. Provided that the melting/freezing time is short compared to the time needed to build the topography by viscous forces, it is thus possible to have a non-zero normal velocity accross the interface. This is taken into account with the help of the boundary condition introduced for the inner core by [START_REF] Deguen | Thermal convection in Earth's inner core with phase change at its boundary[END_REF]. This boundary condition, which translates the continuity of normal stress across the interface, is written in dimensional form as:

Δρ m gτ φ u r + 2η ∂ u r ∂r -p = 0. (2.15)
where Δρ m is the density difference between the solid and liquid phases and τ φ is the phase change timescale. Note that our definition of the dynamic pressure (defined here as p = P -�P �) differs from that of p used by [START_REF] Deguen | Thermal convection in Earth's inner core with phase change at its boundary[END_REF].

The laterally constant term Δρ m gτ φ Ṙ is thus included in p instead of explicitly appearing in the boundary condition. The dimensionless form of the boundary condition is

±Φ ± u r + 2 ∂ u r ∂r -p = 0 (2.16)
where Φ is the phase change number defined as:

Φ ± = |Δρ m | ± gLτ φ η (2.17)
(the superscript + denotes the interface between the SMO and the solid at R + while the superscript -denotes the interface between the BMO and the solid at R -). Moreover, the continuity of tangential stress is simply written as a classic free-slip boundary condition.

The phase change timescale τ φ is related to the time needed to transport latent heat in the magma ocean from the areas that freeze to the areas that melt [START_REF] Deguen | Thermal convection in Earth's inner core with phase change at its boundary[END_REF]:

τ φ = ρL h (ρ -Δρ m ) 2 C p (∂ P T m -∂ P T isen )gu � (2.18)
where u � is the velocity scale in the magma ocean. A reasonable value for the latter is u � ∼ 1 m s -1 [START_REF] Lebrun | Thermal evolution of an early magma ocean in interaction with the atmosphere[END_REF]. Using nominal values for the other parameters, we find that τ φ ∼ 10 4 s. Plugging this in the expression of the phase change parameter eq. ( 2.17 which leads to a flow-through boundary (we use this value rather than 10 -5

because the resolution of radial modes is more computationally demanding as Φ decreases, while the overturn timescale is not affected as shown in the results).

For the Earth, these two values are also considered at the bottom of the solid, accounting for the possible presence of a basal magma ocean (BMO, [START_REF] Labrosse | A crystallizing dense magma ocean at the base of the Earth's mantle[END_REF]. For Mars and the Moon, we do not consider the possibility of a BMO and the bottom interface is hence non-penetrative, u r (R -) = 0. Rather than being realistic, these extreme constant values are used to study how the possibility of melting and freezing at the interface affects the stability of the solid, both in terms of onset time of overturn and preferred mode of motion.

The estimated nominal value being Φ ∼ 10 -5 , we expect the real system should be closer to the flow-through case than to the classical non-penetrative case.

Determination of overturn timescale

We start from a completely molten primitive mantle (R + = R -and T + = T -). We numerically integrate eq. (2.10) to obtain R + as a function of time (the potential surface temperature T p and the surface temperature T s are computed using eq. (2.8) and eq. (2.9)).

At each timestep of this integration, we compute the reference temperature and composition profiles in the solid as shown in section 2.1 as well as the dimensionless numbers Ra(t), Rc(t), W (t) and Γ(t). Using a Chebyshevcollocation approach (e.g. [START_REF] Guo | The Application of the Chebyshev-Spectral Method in Transport Phenomena[END_REF][START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF], the set of linearized equations around the reference state is written as an eigenvalue problem (see appendix B). Solving numerically this problem yields the growth rate and shape of the most unstable mode of overturn. The inverse of that growth rate is the timescale for convection to set in in the solid shell. We compute this timescale at each timestep of the evolution of the SMO. By comparing this timescale with the corresponding time in the evolution of the SMO, we can assess whether convection is able to take place before the entire magma ocean is crystallized. Three different models are considered for the bulk of the solid:

1. full model: compositional, thermal, and moving frame terms are taken into account;

2. thermal model: compositional terms are left out, modeling the ideal case where no fractional crystallization occurs and the sources of instability are purely thermal (eq. (2.14) and the corresponding buoyancy term in eq. ( 2.12) are ignored);

3. frozen-time model: moving frame terms (right-hand-side of eqs. (2.13) and (2.14)) are left out, resulting in a frozen-time approach where all long term evolution terms are ignored when studying the stability of the system at a given instant.

We also compare the timescale obtained by linear stability analysis with the Stokes time τ Stokes = ηL 2 /(ΔρgL 3 M ) computed at each time to check whether this time is a relevant proxy of the stability of the solid mantle.

Results

The The simplest cases are the one neglecting the compositional effects on density. For such cases, the destabilization timescale tends to infinity for a given non-zero thickness of crystallized mantle. This thickness corresponds to the one needed for instabilities to overcome diffusion of perturbations of the reference state. In other words, it corresponds to the thickness above which the Rayleigh number in the solid part is above the critical Rayleigh number. For the Moon, this thickness is never reached and the Moon's mantle stays stable with respect to purely thermal convection. For the Earth and Mars, this thickness is reached rather early, after ∼ 500 km of solid mantle is formed. As crystallization progresses, the thickness and the temperature contrast between the top and the bottom of the solid mantle increase. The available buoyancy in the system therefore increases. This leads to a strong decrease of the destabilization timescale, which becomes much shorter than the time needed to crystallize the remaining surface magma ocean (up to 6 orders of magnitude, depending on which boundary conditions are considered). This suggests that even in the purely thermal case, solid-state convection sets in before the mantle is completely crystallized for planets larger than Mars.

The cases taking compositional effects on density into account are always unstable. This contrasts with the purely thermal cases and is due to the fact

� ��� ���� ���� ���� ���� ������������������������������ ����� ������ ����� ����� ����� � � � � ����� ���������������������������������������� ������������������������� �������������������������� ����������������� � � � � ��������������������� � � �� � � � � �������������������������������� � � � �� � � �������������������� ����������� � ��� ��� ��� ���� ���� ������������������������������ ���� � ��� ��� ��� ��� ���� ������������������������������ ���� Figure 2:
Growth time of the most unstable mode as a function of the crystallized mantle thickness for the Earth, Mars, and the Moon. The solid black line is the time necessary to crystallize the remaining surface magma ocean. Colors represent different boundary conditions: both horizontal boundaries non-penetrative (blue); flow-through boundary condition between the solid and the surface magma ocean to model the possibility of melting and freezing (see section 2.4 for details) (green); and flow-through boundary conditions for both horizontal boundaries assuming the presence of a basal magma ocean (red). Linestyles represent different approximations regarding compositional effects (fractional crystallization and effect on density) and moving frame contributions: both are taken into account (solid lines), compositional effects are neglected (dash-dotted lines), or moving frame terms are neglected (dotted lines). The black dashed line is the Stokes time for each thickness, given for comparison.

that diffusion of the composition field is neglected. There is no mechanism to damp perturbations around the reference state, the latter is hence always unstable. Similarly to what is observed for the thermal cases, the destabilization timescale drops dramatically as the solid mantle thickens. For the Earth and Mars, the destabilization timescale ends up being shorter than the crystallization time of the remaining SMO by several orders of magnitude. The case where moving frame terms are neglected exhibits a shorter destabilization time scale at small thickness. The moving frame terms play a stabilising role only at the begining of mantle crystallization for the Earth and Mars but are significant through the entire Moon's mantle crystallization. The stabilising effect of the moving terms can be understood from the energy conservation eq. (2.13). Taking a temperature perturbation θ > 0 and the associated velocity perturbation u r > 0, one can notice there is a competition between the advection term u r ∂ r T < 0 and the moving frame term W (r -1)∂ r θ whose average is negative.

The same reasoning can be made with a negative perturbation and on the iron conservation eq. (2.14).

For the Moon, the destabilization timescale is always greater than the time needed to crystallize the SMO. However, it should be noted that in this study the time to crystallize the SMO is computed assuming a well-mixed SMO with a surface behaving like a black body. The formation of a light solid crust enriched in plagioclase when around 80% of the SMO is crystallized is expected to slow down the solidification of the SMO by a few million years (e.g. [START_REF] Elkins-Tanton | The lunar magma ocean: Reconciling the solidification process with lunar petrology and geochronology[END_REF]. This would leave enough time for convection to set in in the solid since the destabilization timescale we find is much shorter than that.

Φ + = ∞ Φ -= ∞ (a) classical case Φ - = ∞ Φ + = 10 -2 (
The three boundary conditions exhibits different destabilization timescales.

The case where both boundaries are non-penetrative (which is the case classically considered) needs more time to destabilize than the case where the boundary between the surface magma ocean and the solid allows melting and freezing.

Convective patterns obtained with a flow-through boundary are substantially different than the classical ones (Figure 3). Aspect-ratio-1 rolls are obtained with classical boundary conditions. However, when the top boundary allows phase change, a spherical-harmonic-degree-1 near-translation mode develops.

Matter freezes on one side of the spherical shell, goes around the core or basal magma ocean, and melts on the other side. In the case with a basal magma ocean and its boundary with the solid of flow-through type, matter also crosses the inner boundary of the spherical shell, resulting in a true translation mode.

These two translation modes involve very little or no deformation of the solid compared to the classical case, and therefore less viscous forces acting against convection. This explains the smaller destabilization timescale associated with these modes as well as the lower critical thickness in the purely thermal case. The most unstable mode is the one with the shortest destabilization timescale. One can notice that in the non-penetrative case (Φ + → ∞), the most unstable mode corresponds to aspect-ratio-1 rolls. The typical roll size of the most unstable mode (rπ/l) is roughly 900 km for the top case (l = 15) and 1850 km for the bottom case (l = 8). However, with a flow-through boundary (Φ + → 0), the most unstable mode is the near-translation mode for both cases.

� ��� ���� ���� ���� ���� ������������������������������ ��

� �� � �� � � � � ����� ����������������� � � � � ��������������������� � � �� � � � � �������������������������������� � � � �� � � � ��� ��� ��� ���� ���� ������������������������������

����

Figure 5: Ratio between the destabilization timescale obtained for the purely compositional case τ C (thermal terms are left out) and the timescale obtained for the purely thermal case τ T (compositional terms are left out). When this ratio is above one, it means the thermal reference profile is more unstable than the compositional reference profile. The Moon is not shown here since the purely thermal case is never unstable (τ T → ∞). The colors are the same as in Figure 2, blue: non-penetrative condition for both horizontal boundaries (Φ ± = ∞); green: flow-through condition at the boundary between the solid and the surface magma ocean; and red: flow-through condition at both horizontal boundaries.

timescales at the end of the crystallization when compositional terms are taken into account. That decrease is due to the strong (i.e. very unstable) compositional gradient appearing at the end of the crystallization. It does not affect the destabilization timescale obtained with non-penetrative boundary conditions; this can be explained by the fact that the strong compositional gradient is in a very thin layer at the top of the domain where vertical velocities vanish, and therefore does not contribute to the driving of the down-and up-welling currents.

A comparison between the purely thermal and purely compositional cases for the Earth and Mars is shown on Figure 5. The ratio between the destabilization timescales for theses two cases is 0 before the critical thickness for the purely thermal case is reached. For Mars, the compositional profile is always more unstable than the thermal profile and controls the destabilization timescale of the system. For the Earth, however, the ratio between the two cases is fairly close to 1 for a large part of the crystallization history: neither the thermal nor the compositional profile dominates the destabilization timescale of the system. The colors are the same as in Figure 2, blue: non-penetrative condition for both horizontal boundaries (Φ ± = ∞); green: flow-through condition at the boundary between the solid and the surface magma ocean; and red: flow-through condition at both horizontal boundaries.

�� � �� � �� � �� � ����������������� �� � �� � �� � �� � �� � �� � � � � � � � � � � � � � � � � � � � � ����� ����������������� � � � � ��������������������� � � �� � � � � �������������������������������� � � � �� � � �� � �� � �� � �� � �� � �� �� ����������������� ���� �� � �� � �� � �� � �� � �� �� ����������������� ����
on the boundary condition. It should be noted that the ratio τ LSA /τ Stokes depends on the body and the boundary conditions considered. Notably, permeable boundary conditions lead to a decrease of τ LSA .

Discussion

We showed for the Earth and Mars that the growth timescale of convective instabilities in a crystallizing mantle from the bottom up is several orders of magnitude smaller than the time needed to fully crystallize that mantle. This holds even without taking into account fractional crystallization and the unstable density gradient it induces. This contrasts with the assumptions made in several studies [START_REF] Hess | A model for the thermal and chemical evolution of the Moon's interior: implications for the onset of mare volcanism[END_REF][START_REF] Elkins-Tanton | Magma ocean fractional crystallization and cumulate overturn in terrestrial planets: Implications for Mars[END_REF]Tosi et al., 2013) where the overturn is assumed to take place because of the compositionally induced density gradient after the entire mantle is crystallized. The numerical simulations performed by [START_REF] Ballmer | Reconciling magma-ocean crystallization models with the present-day structure of the Earth's mantle[END_REF] for Earth-like objects lead to a destabilization of the solid after a few Myr, and those performed by [START_REF] Maurice | Onset of solid-state mantle convection and mixing during magma ocean solidification[END_REF] for Mars-like objects lead to a destabilization after roughly 1 Myr. These times are not easily comparable to the timescales we compute via linear stability analysis since the physical problems are different in non-trivial ways: the simulations of [START_REF] Ballmer | Reconciling magma-ocean crystallization models with the present-day structure of the Earth's mantle[END_REF] are in a 2D aspect-ratio-1 cartesian box, those of [START_REF] Maurice | Onset of solid-state mantle convection and mixing during magma ocean solidification[END_REF] are in cylindrical geometry with a variable viscosity, a melt extraction mechanism and a solidus temperature that depends only on pressure. However, despite these differences, the destabilization time uncovered by these simulations are rather similar to the one we predict for the non-penetrative cases: of the order of 1 Myr for the Earth and 0.5 Myr for Mars. This confirms the linear growth rate of instabilities is a relevant proxy for the timescale at which convection sets in.

Moreover, allowing transfer of matter via melting and freezing at the interface between the solid and the surface magma ocean reduces dramatically the timescale with which solid-state convection can set in. It also changes the shape and harmonic degree of the most unstable mode: a degree-one translation mode is preferred. Therefore, the possibility of melting and freezing at the interface should be accounted for when studying the overturn of the primitive mantle of planetary bodies. For example, the case of the Moon is an interesting potential application. This body has a strong dichotomy: the near-side presents more mare basalts, more KREEP material, and a thinner crust than the far-side.

Wasson and Warren (1980) already proposed that such features could be due to a slower cooling of the lunar magma ocean on the near side than on the far-side.

A permeable boundary would allow the solid mantle to overturn with a dominant degree-one before the entire crystallization of the mantle (keeping in mind that the end of the crystallization is much slower than what we predict with our simple model, see Elkins-Tanton et al. ( 2011)). The mechanisms proposed to build a degree 1 at the scale of the Moon involve the dynamics of an entirely crystallized lunar mantle (e.g. [START_REF] Parmentier | Gravitational differentiation due to initial chemical stratification: origin of lunar asymmetry by the creep of dense KREEP?[END_REF]Zhong et al., 2000). The possibility to form a degree one while the crystallization of the magma ocean is still ongoing is therefore worth exploring with more complete models to test whether this dominant degree-one can be conserved after crystallization of the magma ocean and/or helps the development of degree-one instabilities such as the ones predicted in the aforementioned studies. It is also tempting to associate the degree-one feature of Mars (the Marsian dichotomy) to the same process but, as explained above, the first degree-one overturn of the solid mantle is expected to happen long before its complete crystallization. Secondary overturning instabilities are possible after the first one that we cannot investigate with the tools presented above. A more complete study investigating the finite amplitude dynamics is necessary to understand the implications of this work to planets larger than the Moon.

It should be noted that several parameters involved in the problem are badly constrained. The viscosity of the solid mantle and even its rheology is such a parameter. It is highly dependent on how close the temperature in the solid is from the solidus and could easily vary by a few orders of magnitude (e.g. Solomatov, 2015). Since the destabilization timescale scales as the Stokes time (Figure 6), the growth time of instabilities (Figure 6). This validates the general approach proposed by [START_REF] Boukaré | Timing of mantle overturn during magma ocean solidification[END_REF] to assess whether solid-state convection sets in before the magma ocean is entirely crystallized: they compare the Stokes time with the time needed to crystallize the magma ocean and their numeri- cal experiments yields that syn-crystallization convection is possible when the ratio between these two times exceeds ∼ 5 × 10 4 . This value however was determined with non-penetrative boundary conditions, the actual threshold should be a few orders of magnitude higher (meaning syn-crystallization convection is allowed for shorter solidification timescales) since the flow-through boundary condition leads to a faster destabilization of the cumulate for the same Stokes time. Another aspect that deserves care is that for Earth-sized bodies, the Stokes time should incorporate both the thermal and compositional density constrasts. [START_REF] Boukaré | Timing of mantle overturn during magma ocean solidification[END_REF] compare the "compositional" Stokes time with the solidification timescale; while this is perfectly valid for the Moon and Mars for which the thermal density contrast is much smaller than the compositional one, this does not hold for the Earth where both terms have similar magnitudes (Figure 5). The tremendous importance of the viscosity is why a viscosity of 10 18 Pa s is assumed in this study since it is a higher bound for the near-solidus viscosity (see Solomatov, 2015, and references therein) and gives the most conservative estimate for the destabilizing time. The viscosity could be significantly lower if the melt fraction is important in the cumulate, Solomatov (2015) suggests 10 14 Pa s as a lower bound at 40% melt fraction (roughly the rheological transition). Another potential effect of viscosity that is neglected in the study is dynamical: since solid state convection occurs during the crystallization of the magma ocean, the temperature in the solid departs from the solidus temperature profile and as a result the viscosity increases. Moreover, the compositional profile becomes gravitationally stable with iron-enriched heavy material being transported from the top to the bottom of the solid. These two effects combined may lead to the stopping of the solid state convection (Solomatov, 2015).

��� ��� ��� ��� ��� � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � ����� ��� ��� ��� ��� ��� ��� ����������������������� ��� ��� ��� ��� ��� ��� � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � ���� ��������������� � � � � � ��������������������� � � �� � � � � �������������������������������� � � � �� � �
Depending on the size of the magma ocean considered, it could then be possible either that the magma ocean crystallizes completely before convection may start again in the solid, or that convection sets in again in the solid before it is entirely crystallized. Studying this scenario requires a more complex method that a simple linear stability analysis since it involves a non-linear feedback between the dynamics of the solid part and its viscosity, temperature, and compositional fields.

Another unconstrained parameter is the partition coefficient of iron between the solid and liquid. An exploration of this parameter shows that the effect of the partition coefficient is rather limited for the Earth, and slightly more important for Mars (Figure 7). This is in agreement with Figures 2 and5 showing the difference between the purely thermal case (corresponding to the extreme value D = 1) and the purely compositional case is rather small for the Earth but more important for Mars.

Finally, our choice of a constant emissivity results in a roughly constant solidification rate, whereas more sophisticated cooling models including an atmosphere predict most of the mantle crystallizes quickly, and the solidification slows down when only a shallow magma ocean remains. Although such effects are important to build realistic solidification models, they should not affect dramatically our results. Indeed, a faster crystallization at the beginning would lead to a destabilization of the solid mantle at a larger thickness, but we expect this difference to be small since the destabilization timescale is rapidly much lower than the solidification timescale.

Conclusions

Upward crystallization of the silicate mantle of planets within a magma ocean is expected to produce a unstably stratified situation, because of both temperature and composition. In this study, we have addressed the question of whether the overturning instability develops faster than the time it takes to crystallize the magma ocean. To that end, we have developed a linear stability analysis tool to compute the growth rate of the fastest overturning mode and studied systematically the effect of the most important parameters: the planet's size (Moon to Earth size), the partition coefficient and the type of boundary condition between the solid and the liquid. In particular, we have introduced a boundary condition that accounts for the possibility of melting and freezing at the interface between the solid mantle and the magma ocean.

This study shows convection is likely to start in the solid mantle of the Earth, Mars and the Moon before the entire crystallization of the surface magma ocean. Evolution models of the primitive mantle of planetary bodies should therefore account for convection and the associated mixing in the solid part of the crystallizing mantle.

This result holds for the Earth and Mars even without fractional crystallization and the unstable compositional gradient it creates in the cumulate. The value of the partition coefficient is found to have little impact on the timing of mantle overturn.

The timescale at which convection sets in scales as the Stokes time. Specifically, it is proportional to the viscosity of the solid. However, it should be kept in mind that these results are obtained assuming a newtonian rheology and a constant viscosity in the solid mantle. Given the central role of viscosity in this problem, better knowledge of the viscosity and rheology of the primitive solid mantle is of primary importance to study its dynamics.

Finally, the possibility of exchange of matter between the solid mantle and the magma ocean(s) should be accounted for in dynamical models of the primitive mantle since it greatly alters the pattern of convection as well as the destabilization timescale. It could even be a way of producing degree-one structures such as the ones observed on the Moon and Mars.

A. Composition profile resulting from the fractional crystallization of the surface magma ocean

Conservation of the heavy component implies that

d dt � � R + R - C(r)r 2 dr + 1 3 � R T 3 -R + 3 � C l � = 0, (A.1)
where no compressibility effect on density is considered, allowing the bulk density to drop out of the equation. Using eq. (2.1) and Ċl = Ṙ+ dC l dR + , assuming Ṙ+ > 0 at all time and R -constant:

1 C l dC l dR + = 3(1 -D) R + 2 R T 3 -R + 3 . (A.2)
Using eq. (2.1), eq. (A.2) can be written for C:

1 C dC dR + - 1 D dD dR + = 3(1 -D) R + 2 R T 3 -R + 3 . (A.3)
Since C(r) does not depend on time, this equation holds for any r � R + (t) (i.e.

everywhere in the solid) and can be written as: However, it is useful to consider the limiting case of a constant partition coefficient D. In that case, a solution to this equation is .5) with C 0 = DC l0 the mass fraction of FeO in the first solid formed.

C = C 0 � R T 3 -R -3 R T 3 -r 3 � 1-D , ( A
Note that eq. (A.5) diverges when r → R T but is in fact only valid as long as C < 1 and C l < 1. When C l reaches 1, the solid formed has the same composition as the liquid. The complete solution therefore is

C =    C 0 � RT 3 -R -3 RT 3 -r 3 � 1-D if r < R s 1 if r > R s , (A.6) with R s = � (R -) 3 C 1 1-D 0 + R T 3 � 1 -C 1 1-D 0 �� 1/3 (A.7)
the value of R + (t) such that C l (t) = 1.

B. Linear Stability

Since the solid is considered isoviscous and no source of toroidal flow is imposed at the boundaries, the velocity field can be expressed in terms of the scalar poloidal potential P: u = ∇×∇×(Pr) (e.g. [START_REF] Ricard | Mantle dynamics with induced plate tectonics[END_REF][START_REF] Ribe | Analytical approaches to mantle dynamics[END_REF]. Linearizing eqs. (2.11) to (2.14) around the reference state (u = 0; T ; C) gives:

Q = ∇ 2 P (B.1) ∇ 2 Q = Ra Θ -�Θ� r + λ + Rc c -�c� r + λ (B.2) Γ 2 ∂ Θ ∂t + ∂ T ∂r L 2 P r + λ -∇ 2 Θ = W � (r -1) ∂ Θ ∂r + � ∂ T ∂r � + Θ � (B.3) Γ 2 ∂ c ∂t + ∂ C ∂r L 2 P r + λ = W (r -1) ∂ c ∂r .
(B.4)

The boundary conditions on the temperature and composition perturbations are trivial: where l and m are the spherical harmonics degree and order and σ l is the growth rate associated to the harmonic degree l. The system is laterally degenerated and m does not affect the growth rate of the perturbation nor the shape of the radial modes P l (r), Q l (r), Θ l (r) and c l (r). These radial modes are discretized using a Chebyshev collocation approach (e.g. [START_REF] Guo | The Application of the Chebyshev-Spectral Method in Transport Phenomena[END_REF][START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF]. Each radial mode is expressed as a vector whose components are the

Θ ± =
L = 0 : N 0 : N 1 : N -1 1 : N -1                       d 2 + (l2 -2)r -2 λ 0 0 0 0 D 2 -1 0 0 1 : N -1 d 2 + (l2 -2)r -2 λ 0 0 0 N l2(Φ + r -1 λ -2r -2 λ + 2r -1 λ d) -(1 + r λ d) 0 0 0 0 D 2 -Rar -1 λ -Rcr -1 λ 1 : N -1 l2(-Φ -r -1 λ -2r -2 λ + 2r -1 λ d) -(1 + r λ d) 0 0 N -(∂r T )l2r -1 λ 0 D 2 + W + � (r -1)d + (∂r T ) + 1 � 0 1 : N -1 -(∂r C)l2r -1 λ 0 0 W + (r -1)d 1 : N -1 (B.12) R = 0 : N 0 : N 1 : N -1 1 : N -1       0 0 0 0 0 : N 0 0 0 0 0 : N 0 0 Γ 2 1 0 1 : N -1 0 0 0 Γ 2 1 1 : N -1 (B.13)
where 1 is the identity matrix, r ij = r i 1 ij , r λ = r + λ1, l 2 = l(l + 1) and

D 2 = d 2 + 2r -1 λ d -l 2 r -2 λ
The extra row and column on top and right of the matrices are respectively the column and row indices of each of the submatrices.

For example, the top left submatrix of the matrix L is only the first row (hence the 0 on the extra column) of the matrix d 2 + (l 2 -2)r -2 λ .

At a given instant during the crystallization, all the dimensionless numbers W , λ, Γ, Ra and Rc appearing in the matrices L and R are known. For any harmonic degree l of the perturbation, finding its growth rate σ l and associated vertical mode X is an eigenvalue problem. The largest eigenvalue is the growth rate, and the associated eigenvector represent the vertical modes. At a given instant, we look for the harmonic degree l with the highest growth rate σ l , which is then used to compute the dimensional destabilization time scale L 2 M /(κσ).

  is the velocity field, p the dynamic pressure, Θ the temperature perturbation with respect to the reference profile T and c the composition perturbation with respect to the reference profile C. �x� denotes the lateral average of the quantity x. Ra is the thermal Rayleigh number, Rc is the compositional Rayleigh number. The terms on the right hand side of eqs. (2.13) and (2.14) are due to the time dependence of the scales L and ΔT , which brings new advection terms associated with the change of frame, with W = L Ṙ+ /κ the dimensionless velocity of the freezing front. See table 1 for the definition and values of the various symbols.

  ) yields Φ ∼ 10 -5 . The phase change number Φ compares the phase change timescale τ φ (i.e. the time needed to erode topography via melting and freezing) to the viscous timescale (i.e. the time needed to build topography with viscous forces). The value of Φ allows to tune continuously the boundary condition between a nonpenetrative classical condition (Φ → ∞) and a fully permeable condition (Φ → 0). Although this number should depend on time since L depends on time and τ φ depends also on time but in a non trivial way, it is kept constant in this study. Two extreme values are tested for the SMO/solid interface: Φ = ∞ which leads to the classical non-penetrative boundary condition and Φ = 10 -2

  destabilization timescales for the Earth, Mars, and the Moon with various boundary conditions along with the time needed to crystallize the remaining SMO are shown on Figure 2. Comparison of the destabilization timescales obtained for various bulk setups and boundary conditions yields information regarding their contribution to the destabilization of the solid.

Figure 3 :

 3 Figure 3: Most unstable convection modes for the Earth when a 1700 km thick mantle has crystallized, for different boundary conditions represented by the values of the Φ parameters at the top and the bottom, as indicated. The dark zones represent negative temperature anomalies while the bright zones represent positive temperature anomalies. The streamlines are superimposed. Note that the linear stability analysis offers no constraint on the orientation and amplitude of these modes, only their harmonic degree and radial shape. (a): both boundaries non-penetrative, the convection rolls have an aspect ratio approximatively equal to 1; (b): flow-through top boundary, the flow pattern is of spherical harmonic degree one, the streamlines go through the top boundary but go around the central part; (c): flowthrough conditions at both boundaries, the flow pattern is of spherical harmonic degree one, the streamlines go through both boundaries, resulting in a translation mode of convection. Similar behavior is obtained for the other bodies.

Figure 4 Figure 4 :

 44 Figure4shows the transition between the non-penetrative and the flowthrough regime occurs over a rather short range of values of the phase change number. Φ + � 1 leads to near-translation while Φ + � 100 leads to classical aspect-ratio-one rolls.A notable feature on Figure2is the steep decrease of the destabilization

Figure 6

 6 Figure 6 shows that the destabilization timescale τ LSA is proportional to the Stokes time τ Stokes = ηL 2 /(ΔρgL 3 M ). Two effects alter this relation: moving frame terms whose effects are not included in the Stokes time, and the strong compositional gradient at the end of the crystallization whose effects depend

Figure 6 :

 6 Figure 6: Growth time of the most unstable mode versus the Stokes time for the Earth, Mars, and the Moon. The solid line is the destabilization timescale obtained with the linear stability analysis τ LSA (case with all terms accounted for). The dashed lines correspond to τ LSA ∝ τ Stokes . Composition, temperature and moving frame terms are all taken into account.The colors are the same as in Figure2, blue: non-penetrative condition for both horizontal boundaries (Φ ± = ∞); green: flow-through condition at the boundary between the solid and the surface magma ocean; and red: flow-through condition at both horizontal boundaries.

  it is directly proportional to the viscosity and could therefore vary by a few orders of magnitude. The strong relation between the viscosity of the cumulate and the overturn scaling has been investigated by Ballmer et al. (2017): their numerical experiments confirm the overturn onset scales as the Stokes time. It should be noted that our flow-through boundary conditions does not affect this result, it only reduces the proportionality factor between the Stokes time and

Figure 7 :

 7 Figure 7: Thickness of the solid cumulate at which the destabilization timescale equals the time needed to crystallize the rest of the SMO for several values of the partition coefficient, D ∈ [0.01, 0.99].The Moon is not shown here since the destabilization timescale is greater than the time needed to crystallize the SMO. The colors are the same as in Figure2, blue: non-penetrative condition for both horizontal boundaries (Φ ± = ∞); green: flow-through condition at the boundary between the solid and the surface magma ocean; and red: flowthrough condition at both horizontal boundaries.

  is general and allows to take into account variations of D.

  eq. (2.16) and the free-slip boundary condition are written in term of the poloidal potential as: -/L -1 is a curvature term due to the definition of the dimensionless radius. L 2 is the horizontal laplacian:L 2 • = ∂ r ((r +λ) 2 ∂ r •)-(r +λ) 2 ∇ 2 •.The quantity Q is introduced to ease the formulation of this system as an eigenvalue problem involving square matrices.The perturbations P, Q, Θ and c are developed using spherical harmonics, P l (r)Y m l (θ, φ)e σ l t (B.9)

values at the N + 1

 1 Chebyshev nodal points (respectively denoted P, Q, T and C). Radial derivatives evaluated at the nodal points r i = expressed with a differentiation matrix d, e.g. ∂ r P(r i ) = (dP) i . We formulate the system of linearized equations along with the associated

Table 1

 1

	: Symbols used in this paper. All quantities with a + superscript are evaluated at the top boundary (R + ), while quantities with a -superscript are evaluated at the bottom boundary (R -). ¶ The emissivity values for the Earth and Mars are chosen so that the crystallization time scale of the SMO is of the order of 1 Myr (Lebrun et al., 2013). For the Moon, we neglect the effects of the atmosphere and assume a black body cooling. * This choice assumes a 400 km thick basal magma ocean. Using R -= 3471 km does not change significantly the results. † From Andrault et al. (2011), ‡ from Andrault et al. (2012). § 10 -2 : flow-through, ∞: non-penetrative. For the Moon and Mars, the possibility of a BMO is not considered and Φ -= ∞ (see section 2.4 for details).

Acknowledgements

We thank the editor Bruce Buffet and an anonymous reviewer for their useful remarks that helped to improve the clarity of our paper. This study is funded by the French Agence Nationale de la Recherche (grant number ANR-15-CE31-0018-01, MaCoMaOc).


Thierry Alboussière

Bound of dissipation on a plane Couette dynamo

Keywords: numbers: 47.27.N-, 52.30.Cv, 91.25.Cw dynamo, energy dissipation, variational turbulence

Variational turbulence is among the few approaches providing rigorous results in turbulence. In addition, it addresses a question of direct practical interest, namely the rate of energy dissipation. Unfortunately, only an upper bound is obtained as a larger functional space than the space of solutions to the Navier-Stokes equations is searched. Yet, in some cases, this upper bound is in good agreement with experimental results in terms of order of magnitude and power law of the imposed Reynolds number. In this paper, the variational approach to turbulence is extended to the case of dynamo action and an upper bound is obtained for the global dissipation rate (viscous and Ohmic). A simple plane Couette flow is investigated. For low magnetic Prandtl number Pm fluids, the upper bound of energy dissipation is that of classical turbulence (i.e. proportional to the cubic power of the shear velocity) for magnetic Reynolds numbers below P -1 m and follows a steeper evolution for magnetic Reynolds numbers above P -1 m (i.e. proportional to the shear velocity to the power four) in the case of electrically insulating walls. However, the effect of wall conductance is crucial: for a given value of wall conductance, there is a value for the magnetic Reynolds number above which energy dissipation cannot be bounded. This limiting magnetic Reynolds number is inversely proportional to the square root of the conductance of the wall. Implications in terms of energy dissipation in experimental and natural dynamos are discussed.

I. INTRODUCTION

Natural dynamos exist whenever the conditions of their existence is possible, i.e. when a sufficiently large magnetic Reynolds is reached. One can imagine that this is an additional route for mechanical energy dissipation and that it is more likely for natural systems to take it than not to take it. It has been argued sometimes that a state of maximal dissipation rate should be reached and this idea has been used as a closure assumption for turbulence (e.g. Malkus [1]). It is not necessary to make such an assumption, as explained clearly by Howard [2], and yet one can draw useful information from the determination of rigorous lower and upper bounds on energy dissipation in turbulent flows. These bounds are obtained in a larger functional space than the solutions to the Navier-Stokes equation, hence they are not necessarily attained. However, in a number of cases, turbulent flows lead to a dissipation rate of the same magnitude as this upper bound [3].

The variational approach to turbulence was introduced by Malkus [4], Busse [5] and Howard [6]. More recently, this approach was reformulated by Doering and Constantin [7] and expressed in a simpler way, using the concept of a background function (not necessarily the mean flow of turbulence) following Hopf [8]. The objective is to bound energy dissipation under the constraint of horizontally averaged energy balance (for statistically plane * Electronic address: thierry.alboussiere@ujf-grenoble.fr invariant configurations). In a series of papers, the approach was improved by optimizing the spectral Lagrange parameter [9] and the background function [3]. The final bound is better than the nearly rigorous bound of Busse. Our objective here is not to exhaust these possibilities of optimization but rather to apply the general principle to a new configuration, i.e. a dynamo problem.

The idea of applying a variational approach to a magnetohydrodynamic flow has been applied already to a Couette and Poiseuille flow subjected to an applied transverse magnetic field [10]. This paper has been a source of motivation for the present work with two significant variations. First, there is now no imposed magnetic field and secondly, magnetic boundary conditions are different. In the paper by Alexakis et al., magnetic disturbances are constrained to vanish at the boundaries, which does not correspond to a physically plausible situation. It is assumed here that the fluid domain is bounded by an infinite domain of electrically insulating medium (or with a conducting solid wall of finite thickness in between).

A plane Couette flow configuration is considered and our objective is to find an upper bound to the total energy dissipation when a prescribed velocity is applied. There is no applied magnetic field of external origin and the problem may look like a purely hydrodynamical one. However, the flow may support dynamo action (see [11,12]) and the amount of dissipated energy must then take into account Joule dissipation. Upper bounds will be obtained as a function of two dimensionless parameters, the Reynolds number and the magnetic Reynolds number. In addition, the effect of electrically conducting walls of finite thickness will be investigated.

Section II provides details on the flow configuration, notation, dimensionless variables and equations. In section III, the principle of decomposition with background flow is presented. The horizontally averaged energy balance is obtained in section IV and the expression for the total energy dissipation (Joule plus viscous) is given in section V. Energy dissipation bounds are obtained in section VI and VII respectively, when velocity fluctuations and magnetic fluctuations are considered respectively in addition to the background flow. An improved bound is determined numerically in section VIII. Section IX is devoted to a discussion of the bounds obtained and their relevance to experimental and geophysical configurations and section X to the directions in which the variational approach could be extended to deal with more relevant models of dynamo.

It may be useful to provide some guidance on how to read this paper. The Hopf-Doering-Constantin method is explicitly introduced and subsequent calculations of upper bounds are also detailed explicitly. A reader with no prior knowledge of the method can check all results with pen and paper until the end of section VII. Section VIII does not provide any fundamentally new result and can be ignored in a first reading. This section has required some standard numerical calculations of eigenvalues related to the magnetic spectral constraint, and the method is only sketched. This section serves two purposes: first, it provides a better upper bound as it corresponds to a background function for which the spectral condition is only just satisfied (it is zero for a particular disturbance) and secondly, it provides a confirmation that the upper bound derived analytically in section VII is relevant as it obeys the same scaling law as the numerical bound at large magnetic Reynolds numbers.

II. PLANE COUETTE FLOW CONFIGURATION

The dimensionless Navier-Stokes and induction equations can be written

∂u ∂t + u • ∇u = -∇p + j × B + Re -1 ∇ 2 u, (1)
∂B ∂t + u • ∇B = B • ∇u + R -1 m ∇ 2 B. (2)
In the equations above, the dimensional length and velocity scales are chosen to be half the distance, H, and differential velocity, U , between the plates (see Fig. 1) while the quantities H/U , ρU 2 , U √ ρµ, ρ/µU/H are taken as dimensional scales for time, pressure, magnetic field and electric current density respectively. The dimensionless parameters are the Reynolds number Re = U H/ν and the magnetic Reynolds number R m = µσU H. The symbols ρ, µ, ν and σ denote density, vacuum magnetic permeability, kinematic viscosity and electrical conductivity respectively. In this dimensionless formulation, the As shown on Fig. 1, the thickness of the fluid layer is 2H -2E where E is the thickness of the electrically conducting layer of the wall. That thickness may be zero E = 0, in which case the wall is electrically insulating. In all cases, the overall thickness of the electrically conducting domain is 2H. For simplicity, the electrical conductivity of the solid conducting layer is identical to that of the fluid. Hence, at z = H -E and z = -H + E, there is continuity of all three components of the magnetic field. Moreover, the induction equation is identical in the fluid and in the conducting layer. Therefore, there is no boundary condition to consider at z = H -E or z = -H + E. The magnetic boundary condition, at z = ±H, is that the magnetic field in the electrically conducting domain is matched continuously to a curl-free and divergence-free magnetic field outside (see appendix A).

III. BACKGROUND VELOCITY DECOMPOSITION

A frame of reference (x, y, z) is chosen with x aligned with the direction of the imposed shear velocity and z in the direction perpendicular to the plates. Its origin is half-way between the plates.

The velocity field u is written as the sum of a steady parallel profile Φ(z)e x satisfying the boundary conditions Φ(±(1 -e)) = ±1 and of a field v with homogeneous boundary conditions v(x, y, ±(1 -e), t) = 0, i.e. u(x, y, z, t) = Φ(z)e x + v(x, y, z, t)

It must be stressed here that Φ -the so-called background velocity profile -needs not be the average value of the full velocity field u, neither in the temporal nor in the horizontal space average sense.

Substituting this decomposition into the Navier-Stokes and induction equations leads to evolution equations for v and B for any particular choice of a background profile:

∂v ∂t + v • ∇v + Φ ∂v ∂x + v z Φ ′ e x = -∇p + j × B +Re -1 Φ ′′ e x + ∇ 2 v , (4) ∂B ∂t + v • ∇B + Φ ∂B ∂x = B • ∇v + B z Φ ′ e x +R -1 m ∇ 2 B, (5)
where Φ ′ and Φ ′′ denote the first and second derivative of the background velocity profile Φ respectively.

IV. ENERGY BALANCE OF FLUCTUATIONS

The dot product of equations ( 4) and ( 5) with v and B respectively are integrated over V a large rectangular volume -L < x < L, -L < y < L and -1 < z < 1. They are then integrated over a long period of time T . The following space and time averaged equations are obtained:

Φ ′ v x v z = (j × B) • v -Re -1 Φ ′ ∂v x ∂z -Re -1 ∇v : ∇v , (6) 0 = Φ ′ B x B z + j • (v × B) -R -1 m j 2 , ( 7
)
where the Poynting flux going out of the control volume has been assumed to be zero (there is no external source or sink, see appendix A) and where the following global space-time average is introduced for any quantity f :

< f >= lim T →∞ lim L→∞ 1 L 2 T T 0 1 -1 L 0 L 0 f dxdydzdt. ( 8
)
The work of Lorentz forces (j × B) • v and electromotive work j • (v × B) are equal and opposite, so that the sum of ( 6) and ( 7) leads to the following energy balance:

-Φ ′ v x v z + Φ ′ B x B z -Re -1 Φ ′ ∂v x ∂z -Re -1 ∇v : ∇v -R -1 m j 2 = 0. ( 9
)
V. TOTAL ENERGY DISSIPATION

Let us take ρνU 2 /H 2 as a dimensional scale for energy dissipation per cubic metre. The average dimensionless energy dissipation due to viscous effects is ∇u : ∇u and the average dimensionless Joule dissipation is P -1 m j 2 . When u is expressed using the decomposition (3), the total dissipation D = ∇u :

∇u + P -1 m j 2 becomes D = ∇v : ∇v + Φ ′2 + 2 Φ ′ ∂v x ∂z + P -1 m j 2 . ( 10
)
Combining the energy constraint (9), so as to remove the linear term Φ ′ ∂v x /∂z , gives the following expression for the dissipation:

D = Φ ′2 -[ ∇v : ∇v + 2Re Φ ′ v x v z ] -P -1 m j 2 -2R m Φ ′ B x B z . ( 11
)
Dissipation is bounded by the background dissipation Φ ′2 when both conditions ∇v : ∇v

+ 2Re Φ ′ v x v z > 0 and j 2 -2R m Φ ′ B x B z > 0 are
satisfied for all admissible vector fields v and B, i.e. divergence-free vector fields satisfying the appropriate boundary conditions. Those two conditions are called spectral conditions in the framework of the background method because they can be treated as an eigenvalue problem and this is preceisely the method that will be followed in section VIII. However, it is also possible to ensure that those conditions are satisfied using other methods and they will be dealt with using functional inequalities in the following sections VI and VII.

VI. VELOCITY SPECTRAL CONDITION

Equation (11) shows that the spectral conditions can be treated independently for velocity and magnetic disturbances. The velocity spectral condition consists in ensuring that ∇v : ∇v + 2Re Φ ′ v x v z remains positive for all velocity fields v.

In order to satisfy the spectral constraint at large Reynolds numbers, it is convenient to choose a background velocity profile which is uniform in most of the fluid, with a linear profile on each side so as to recover correct boundary conditions (see fig. II). The thickness δ of these linear parts (virtual boundary layers) is a free parameter. The 'production term' 2Re Φ ′ v x v z is confined to these regions where it can hopefully be balanced by viscous 'dissipation' ∇v : ∇v for a sufficiently small thickness δ.

This problem is solved as follows (Doering and Gibbons [START_REF] Doering | Applied Analysis of the Navier-Stokes Equations[END_REF], Alexakis et al. [10]). Using the condition v = 0 at the lower wall, velocity at another position z in the fluid can be bounded as follows using Schwartz relationship:

v x = z -1+e ∂v x ∂z ′ dz ′ ≤ z -1+e dz ′ z -1+e ∂v x ∂z ′ 2 dz ′ ≤ √ z + 1 -e 1-e -1+e ∂v x ∂z ′ 2 dz ′ . ( 12
)
Combining with the corresponding equation for v z and using Young relationship leads to:

v x v z ≤ 1 2 (z + 1 -e) 1-e -1+e ∂v x ∂z ′ 2 + ∂v z ∂z ′ 2 dz ′ ≤ 1 2 (z + 1 -e)
1-e -1+e

∇v : ∇vdz.

Integrating v x v z Φ ′ over the lower boundary layer yields:

-1+e+δ -1+e v x v z Φ ′ dz ≤ δ 4 1-e -1+e
∇v : ∇vdz

Averaging over the x and y direction and taking into account the upper boundary layer leads to:

Φ ′ v x v z ≤ δ 2 ∇v : ∇v . ( 15
)
The velocity spectral condition is then satisfied as soon as:

δ ≤ Re -1 . (16)
In terms of energy dissipation (11), assuming magnetic effects to be absent, this background profile with sufficiently small boundary layers ( 16) provides an upper bound:

D ≤ 2Re (17)

VII. MAGNETIC SPECTRAL CONDITION

In the previous section, it has been possible to find a condition on δ fulfilling the spectral requirement for velocity perturbations. The corresponding task for magnetic perturbations cannot be exactly similar as the property v = 0 on the walls has to be changed into physically sound magnetic boundary conditions. We are considering that the magnetic field inside the fluid and boundary domains should match a potential field outside decaying to zero at infinity (or possibly to a non-zero constant if a non-zero mean electric current density is allowed in the electrically conducting domain).

We cannot easily find a simple expression for the boundary condition for B x or B z in the physical space, however we obtain below a boundary condition for the (x, y) average product B x B z . This will then been used to bound Φ ′ B x B z in terms of j 2 .

Let us denote the horizontal average over the (x, y) plane with an overbar:

f (z) = lim L→∞ 1 L 2 L 0 L 0 f (x, y, z) dxdy. ( 18
)
Outside the domain, the curl of the magnetic field is zero. So is the product of its y component with B z averaged over x and y:

B z (∂ z B x -∂ x B z ) = 0. ( 19
)
Using integration by part, divergence-free condition for B, another integration by parts and the z component of ∇ × B = 0 leads to:

B z (∂ z B x -∂ x B z ) = ∂ z B x B z -B x ∂ z B z , = ∂ z B x B z + B x ∂ y B y , = ∂ z B x B z -B y ∂ y B x , = ∂ z B x B z + B y ∂ x B y , = ∂ z B x B z = 0, (20)
As the product B x B z is vanishing at infinite z, we conclude that B x B z is zero at the boundary of the conducting domain, i.e. when z = ±1. Strictly speaking, this is enough to see that it will be possible to bound B x B z in terms of the mean square of the gradient of B using Poincaré's theorem, and then in terms of the mean joule dissipation. A convenient way of obtaining such a bound is to repeat the calculations leading to (20) while retaining a non-zero electric current density for -1 ≤ z ≤ 1. This leads to:

∂ z B x B z = j y B z -j z B y = (j × B) x , (21)
which is the expression of the x component of Lorentz forces in terms of the magnetic stress tensor, averaged on constant z planes. The product B x B z can be evaluated inside the electrically conducting domain using the boundary condition B x B z = 0 at z = ±1.

B x B z (z) = z -1 ∂ z B x B z dz ′ , = z -1 (j × B) x , (22)
which can be bounded as follows:

B x B z (z) ≤ z -1 dz ′ z -1 j 2 B 2 dz ′ . ( 23
)
As shown in appendix B, the x and y averaged square magnetic field in the electrically conducting domain is bounded as follows:

B 2 (z) ≤ 2 j 2 . ( 24
)
Hence, equation (23) can be written:

B x B z (z) ≤ √ 2 √ z + 1 j 2 . ( 25
)
The magnetic boundary condition is rather loose compared to the velocity boundary condition which leads here to a bound proportional to the square-root of z + 1 for the product B x B z rather than a linear dependence for the velocity product v x v z . Using equation ( 25) it is now possible to bound

B x B z Φ ′ : | B x B z Φ ′ | ≤ δ -1+e+δ -1+e B x B z (z) dz + 1 δ 1-e 1-e-δ B x B z (z) dz, ≤ √ 2 δ j 2 e+δ e √ udu. ≤ √ 2 3δ j 2 (e + δ) 3/2 -e 3/2 (26)
The magnetic spectral condition j 2 -2R m Φ ′ B x B z ≥ 0 is satisfied as soon as:

8 √ 2 R m 3δ (e + δ) 3/2 -e 3/2 ≤ 1. ( 27
)
This relationship is then used to express δ is terms of viscous dissipation for the background flow, hence to determine an upper bound for global dissipation. This bound is plotted on Fig. 2 for various values of the wall thickness e. The limiting case of vanishing wall thickness (or rather of an electrically insulating wall) can be easily derived from (27) and corresponds to:

8 √ 2 R m δ 1/2 3 ≤ 1, (28)
and since the dissipation of the background profile is 2/δ, we have: 17) for the velocity spectral condition ('plus' signs), and (28) for the magnetic spectral condition, for different wall conductances e = 10 -3 , 10 -5 , 10 -7 and 0, plotted here for Pm = 10 -5

D ≤ 256 R 2 m 9 . ( 29
)
One must take the higher upper bound derived from the velocity and magnetic spectral conditions in Fig. 2 as both conditions must be satisfied. It is straightforward from equations ( 17) and ( 27) that the magnetic spectral condition will take over the velocity spectral condition when the magnetic Reynolds number is exceeding 9/(128 P m ) ≃ 0.07 P -1 m . It can also be derived from (27) that the critical magnetic Reynolds number at which our upper bound estimate diverges towards infinity is R m = 1/(4 √ 2e).

VIII. NUMERICAL DETERMINATION OF UPPER BOUNDS

The velocity and magnetic spectral conditions appearing in (11) are tackled numerically in this section. The reason for this is not to achieve a better (i.e. lower) upper bound, although this is the case. The main reason is to check that our analysis in sections VI and VII provides results that are of the correct magnitude and not gross overestimates of dissipation upper bounds. For a given background profile, the numerical analysis provides Reynolds and magnetic Reynolds numbers such that the spectral conditions are satisfied for all admissible velocity and magnetic fluctuations, and this cannot be improved as one particular velocity disturbance and one particular magnetic disturbance make the spectral conditions just zero.

The method is sketched here: the spectral condition can be treated as a problem of energy stability for the background flow [START_REF] Doering | Applied Analysis of the Navier-Stokes Equations[END_REF]. This is changed into an eigenvalue problem and the game consists in finding a Reynold or magnetic Reynolds number such that the maximum eigenvalue is zero. The velocity and magnetic fluctuations are Fourier transformed in the x and y directions, while they are expanded using Chebyshev collocation polynomials in the z direction. The optimal magnetic disturbance is always found in the limit of zero wave numbers k x and k y , whereas the optimal velocity disturbance is found for k x = 0 and k y finite and increasing with the Reynolds number (see Fig. ( 4)).

For this numerical analysis, the background function is no longer piecewise linear. In order to avoid numerical difficulties due to discontinuities in the velocity gradient, a hyperbolic sine function class is chosen Φ = sinh(zδ). In the limit of large dissipation, the background velocity is zero everywhere except in thin boundary layers where the velocity profile is an exponential function. Similarly to what has been done in sections VI and VII, the background functions contain a single free parameter, i.e, the typical thickness δ of the boundary layers. Due to the new profile of the background function, the associated dissipation takes a different form that can be derived analytically:

1 -e + sinh 2 1-e δ 2δ sinh 1-e δ 2 .
(30)

For a given choice of δ, hence dissipation, a lower bound of Reynolds and magnetic Reynolds number is sought. This is repeated for various values of δ and the resulting curve can be read as an upper bound of dissipation for each value of Reynolds or magnetic Reynolds number (see Fig. 4 and Fig. 3 respectively). Both spectral conditions can be plotted on the same figure once a magnetic Prandtl number is specified.

It is clear that this numerical procedure leads to improved upper bounds -by a factor 10 approximatelycompared to those obtained analytically, while similar trends are obtained. For the simple Couette flow configuration under investigation, and for electrically insulating boundaries, it has been found that a classical 'hydrodynamical' upper bound for dissipation holds when R m is smaller than around P -1 m : the dimensional upper bound of dissipation per unit mass is of order U 3 /H. For larger magnetic Reynolds numbers though, another scaling is followed and the bound on dissipation becomes proportional to µ 2 σ 2 ν U 4 . This new bound is higher than the 'hydrodynamical' bound and is independent of the half-distance H between the plates.

When the boundary plates have the same electrical conductivity as the fluid and a finite thickness e, the resulting upper bound is changed dramatically. The upper bound corresponding to the magnetic part of the spectral condition increases suddenly and diverges to infinity for a value of the magnetic Reynold number depending on the thickness of the walls. This critical magnetic Reynolds number scales as e -1/2 . This value can be reached much before the change of regime between 'hydrodynamical' and (insulating) 'magnetic' branches of upper bound estimates (R m ∼ P -1 m ). It should not be too surprising that dissipation can be unbounded at some finite magnetic Reynolds number. Such a situation can be easily simulated with a Bullard dynamo model. This is a solid-rotation dynamo with no possible feedback of the magnetic field on the structure of the flow. Hence magnetic energy (and dissipation) grows without limit above the threshold of linear instability. As soon as electrical conductivity is restricted to the fluid (e = 0), this behaviour is no longer possible. For any finite magnetic Reynolds number, energy dissipation is certainly finite. Yet the scaling law for this dissipation is distinctly above the 'hydrodynamical' a priori estimate.

Whether one is inclined to believe that upper bounds of dissipation are generally relevant to turbulence or not, there is one firm and strong result from the present work which is the extended validity of the Kolmogorov scaling up to R m ∼ P -1 m . It might have been a reasonable guess to predict a distinct increase of the dissipation as soon as R m exceeds a critical threshold of order unity. This analysis shows that this is not possible with electrically insulating boundaries. This certainly has implications on the design of future Couette-type dynamo experimental setups, depending on whether a small or large dissipation is sought. Although there is no mathematical statement on how dissipation is divided between mechanical and electrical phenomena, it might be a better choice to have electrically conducting boundaries, should one want to increase the ratio of Joule to viscous dissipation.

It is appropriate to discuss here how energy dissipation was affected by dynamo action in those three experimental setups where a self-sustained liquid dynamo was observed. They are the Riga (Latvia) [START_REF] Gailitis | [END_REF], Karlsruhe (Germany) [15] and Cadarache (France) [16] dynamos. From these references, a table can be shown (table I), We have concentrated mainly on the low magnetic Prandtl number case, because our primary interest lies in the Earth's core dynamics and in the relevant liquid metal experiments. Nevertheless the results we have obtained are equally valid for large magnetic Prandtl numbers. It is however anticipated that our dissipation bound will be a gross overestimate at low values of the hydrodynamic Reynolds number. In that case, the flow is the simple uniform shear flow between the moving walls and such a flow is not going to sustain dynamo action easily. This perspective should not be ruled out as strong Lorentz forces might change the flow in such a way that it could drive a dynamo, however this transition -if it exists -will be severely subcritical.

Saturation can somehow be addressed within this work. The total energy dissipation is an upper bound for Joule dissipation, which can be used through equation (B1) to obtain an upper bound of the magnetic energy in a saturated regime of dynamo action. This bound cannot give realistic estimates near dynamo threshold, as dissipation is still dominated by the viscous contribution. At large magnetic Reynolds numbers, this bound is much larger than equipartition between kinetic and magnetic energy.

In the derivation of velocity and magnetic spectral conditions, there are strong similarities. The only difference is due to boundary conditions for velocity disturbances and magnetic field disturbances. This mere difference in the boundary conditions leads to fundamentally different upper bounds for dissipation. Magnetic boundary conditions are sometimes not well treated in general textbooks (see, for instance, the energy stability of hydromagnetic flows in the otherwise excellent book by Joseph [START_REF] Joseph | Stability of fluid motions I[END_REF], volume II). It is confirmed again in the present work that one should pay great attention to the relevant physical boundary conditions.

X. PERSPECTIVES

It is certainly useful to derive such upper bounds for dissipation in the field of magnetofluid dynamics, and more particularly regarding dynamos. When complex situations are considered (magnetic field, global rotation, ...) there is no reliable heuristic approach that can provide good estimates for dissipation. Ordinary hydrodynamic turbulence can be roughly tackled by Kolmogorov's theory but it is not safe to extend it to other types of turbulence. When the 'Doering-Constantin-Hopf' background profile method can provided an upper bound, this is a solid reliable result. The bound can be rather constraining for instance in the case of a Couette flow with electrically insulating boundaries, at low P m : it is not an obvious result that dissipation must remain similar to hydrodynamic dissipation up until R m ∼ P -1 m . There are actually a lot of other configurations for which the method can be applied.

Upper bounds for dissipation can be calculated for experimental setups involving liquid metal flows. In Grenoble , LGIT laboratory, we have the DTS setup involving 40 liters of sodium in an imposed magnetic field [START_REF] Nataf | [END_REF][19][20]. Upper bounds for dissipation would provide some information about the importance of size scaling of the experimental setup. This is crucial when planning a larger setup that might sustain dynamo action.

Thermal convection can also be taken into account to derive dissipation bounds. This is the driving force for the flow of liquid iron in the Earth's core. Upper bounds of dissipation may prove useful in terms of the thermal budget of Earth throughout its history. Although the configuration of the core of the Earth is not similar to the simple Couette flow considered here, it could be very relevant to examine the role of a thin electrically conducting layer at the bottom of the mantle, as it could make a large change in the upper bound for dissipation. end of this appendix. Poloïdal and toroïdal contributions lead to

B = -∂ y q -∂ 2 xz p ∂ x q -∂ 2 yz p ∇ 2 S p (B4) j = ∇ × B = ∂ y ∇ 2 p -∂ 2 xz q -∂ x ∇ 2 p -∂ 2 yz q ∇ 2 S q (B5)
Hence, x and y averaged magnetic energy and ohmic dissipation can be written

B 2 (z) = k 2 2 |Q| 2 + k 2 2 |∂ z P | 2 + k 4 2 |P | 2 (B6) j 2 (z) = k 2 2 -k 2 P + ∂ 2 z P 2 + k 2 2 |∂ z Q| 2 + k 4 2 |Q| 2 (B7)
Integration by parts leads to the following expression for the global dissipation using the condition that the poloïdal and toroïdal components must decay far away from the fluid layer

j 2 = 1 -1 j 2 dz = ∞ -∞ j 2 dz = ∞ -∞ k 2 2 ∂ 2 zz P 2 + k 4 |∂ z P | 2 + k 6 2 |P | 2 dz + ∞ -∞ k 2 2 |∂ z Q| 2 + k 4 2 |Q| 2 dz (B8)
The toroïdal scalar function Q must vanish at z = ±1. Hence:

Q(z) = z -1 ∂ z Qdz (B9)
Hence by Cauchy-Schwartz

|Q| 2 ≤ z -1 dz z -1 |∂ z Q| 2 dz ≤ 2 1 -1 |∂ z Q| 2 dz (B10)
The treatment of the poloïdal is slightly more involved as the boundary conditions available are less straightforward: ∂ z P ± P = 0 at z = ±1. They can be used as follows leading to a similar result

|∂ z P | 2 ≤ k 2 1 -1 |∂ z P | 2 dz ′ + 1 -1 ∂ 2 zz P 2 dz ′ (B15)
All three terms in (B6) can be bounded using (B10), (B13) and (B15) and then compared to the contributions of the ohmic dissipation in (B8), leading to the expected result:

B 2 (z) ≤ 2 j 2 (B16)
Importantly, this result is independent of k, so that it will equally apply to a sum of contributions of different wavenumbers k x and k y , hence proving the expected result for the fields B and j.

The horizontally independent contributions b 0x (z)e x + b 0y (z)e y must now be considered. Their associated ohmic dissipation can be written A similar inequality holds for b 0y , so comparing (B17) and (B18) is also compatible with the inequality to prove.

FIG. 1 :

 1 FIG. 1: Plane Couette flow configuration

FIG. 2 :

 2 FIG.2: Bound on dissipation obtained from equations[START_REF] Joseph | Stability of fluid motions I[END_REF] for the velocity spectral condition ('plus' signs), and (28) for the magnetic spectral condition, for different wall conductances e = 10 -3 , 10 -5 , 10 -7 and 0, plotted here for Pm = 10 -5

5 FIG. 3 :

 53 FIG. 3: Numerical magnetic upper bound for energy dissipation for different wall thicknesses

FIG. 4 :

 4 FIG.4: Numerical upper bound for energy dissipation from the velocity spectral condition and associated spanwise wavenumber ky of marginal perturbation

1 z 1 - 1 k 2 1 - 1 k 2

 1112112 2kP = (kP -∂ z P ) + (kP + ∂ z P ) = z -1k∂ z P -∂ 2 zz P dz ′ -k∂ z P + ∂ 2 zz P dz ′ (B11) from which the modulus of P can be bounded using the triangle and Cauchy-Schwartz inequalities 2k |P | ≤ |∂ z P | + ∂ 2 zz P dz ′ ≤ |∂ z P | 2 dz ′ + 2 is made on ∂ z P 2∂ z P = -(kP -∂ z P ) + (kP + ∂ z P ) (B14)

  electrical current flux in the x and y direction is not zero, then this is an unphysical situation as magnetic energy in the space above and below the fluid would be infinite, and such a situation would take an infinitely long time to be established. If it is zero, then b 0x and b 0y vanish on z = ±1, and one can derive a bound for |B| 2 using Cauchy-Schwartz inequality|b 0x | 2 (z) =

TABLE I :

 I Extra dissipation due to dynamo action slightly above threshold in experimental dynamos . the Karlsruhe dynamo experiment can be clearly distinguished by the large amount of extra dissipation compared to the other two. This experiment is also one where differential velocity is forced deep inside the setup. For the Riga and Cadarache dynamos, a turbine or propellers are set in rotation near the outer boundary of the setup. There can only be a gross comparison to the pure Couette flow configuration presented here, but the case of thick electrically conducting boundaries is also putting shear velocities well inside the electrically conducting domain and results in a huge increase of the dissipation bound. Conversely, it cannot be ruled out that Riga and Cadarache dynamos are still obeying the same classical 'hydrodynamical' law of dissipation, similarly to the Couette flow with electrically insulating walls. The authors of the dynamo experiments themselves specify that the extra dissipation is obtained by difference of the actual dissipation and the extrapolated curve of dissipation below threshold. Even though it has been used to estimate Ohmic dissipation by some authors, it should be stressed here that there is no theoretical reason why this extra dissipation should be the Ohmic contribution to dissipation. For instance, in the present work, we use the fact that there is an exact cancellation between the work of Lorentz forces and the electromotive work, but we have no access to either of these terms.

		Riga [14] Karlsruhe [15] Cadarache [16]
	% above critical Rm % extra dissipation	6 -7 10	6.7 42	30 15 -20
	where the estimated extra dissipation when dynamo ac-
	tion is present is given (in percentage) for a certain mag-
	netic Reynolds number (expressed in percentage above
	the threshold)

Acknowledgments

Thanks are due to Alexia Gorecki and Maylis Landeau who have been working on student projects on closely related topics and provided stimulating discussions. I have benefited from a discussion with Christophe Gissinger on the role of the Poynting flux and with Alexandros Alexakis who pointed out to me a mistake in the original proof regarding the analytical bound for the magnetic spectral constraint. I am also grateful to the members of the Geodynamo team (LGIT lab in Grenoble) for their comments on this work.

Appendix A: Magnetic boundary conditions and associated Poynting flux

Outside our electrically conducting domain -H < z < H (i.e. the fluid layer and electrially conducting part of the plates), it is assumed that there exist no other electrically conducting domains, and no regions with magnetic properties: let us just think of it as empty space. Coherent with the classical magnetostatic approximation used in the electrically conducting domain to derive the induction equation ( 2), the displacement current is also neglected outside, so that the magnetic field satisfies the following Maxwell equations:

from which it follows that the magnetic field is harmonic.

The magnetic field can be written as a Fourier integral:

and the harmonicity of B implies that each component (k x and k y fixed) must satisfy:

where

As a consequence B must decay exponentially away from the fluid layer on a typical distance 1/k. There is one exception when there is a significant magnetic contribution near k = 0. However, we shall exclude that possibility on physical grounds: in that case, the exponential decay length increases without limit and there would be an infinite amount of magnetic energy stored oustside the fluid. That would take an infinitely long time to settle and is thus incoherent with our assumption of the existence of a stationary flow solution.

Faraday's equation (another equation from the set of Maxwell's equations) is used to determine the electric field:

showing that it must also decay exponentially. Hence, the Poynting flux density E × B/µ decays exponentially and is thus as close to zero as one wishes some distance away from the fluid layer. From the assumption of statistical stationarity, there can be no accumulation of magnetic energy, hence this Poynting flux E × B/µ must vanish when averaged along x and y, at any position z outside the electrically conducting domain. This reasoning is valid for each wavenumber and the total Poynting flux is the sum of individual wavenumber components, hence this proves that the averaged Poynting flux must be zero at any position z and in particular at the boundary z = H and z = -H (or z = ±1 in dimensionless coordinates).

Let us insist that this result is true when the Poynting flux is averaged on a plane of constant z (not pointwise) and when it is also averaged in time, assuming stationary turbulence.

It should be noted that the above reasoning fails in two circumstances. One is when there are other electrically conducting domains. In that case, magnetic energy can be dissipated or generated in each domain and there can be a net exchange of energy between the fluid layer and the other domains. The other possibility is that energy is radiated away from the fluid layer and goes away infinitely far. This is only possible when the displacement current is taken into account in Ampere's equation (A1). Then a significant fraction of energy can be radiated away only when the timescale τ and lengthscale l of the magnetic field are such that l/τ is comparable to the speed of light [START_REF] Alboussière | Fundamentals of MHD, a chapter of[END_REF]. This is a necessary condition for our system to work as an antena, and we shall not consider this limiting case.

Appendix B: Bound on magnetic energy for a given amount of Joule dissipation

We show here that magnetic energy in the conducting domain is bounded pointwise by the integral Joule dissipation. In our dimensionless terms, the following relationship hods for any value of z:

The magnetic field B is decomposed in poloïdal-toroïdal contributions:

B = ∇q ×e z +∇×(∇p × e z )+b 0x (z)e x +b 0y (z)e y (B2)

The assumption of spatial statistical invariance along x and y directions allows us to write the poloïdal and toroïdal contributions of B in terms of a Fourier integral [p, q] = [P, Q] (k x , k y , z, t) e (ikxx+ikyy) dk x dk y (B3) Each Fourier component can be considered individually, as magnetic energy and ohmic dissipation will be just the sum of energy and dissipation of each contribution. It is enough to prove (B1) for all contributions. The contribution b 0x (z)e x + b 0y (z)e y will be dealt with at the