Thibaut T Clarté

Nathanaël Schaeffer

Stéphane Labrosse

Jérémie Vidal

The Effects of Robin Boundary Condition on Thermal Convection in a Rotating Spherical Shell

Convection in a spherical shell is widely used to model fluid layers of planets and stars. The choice of thermal boundary conditions in such models is not always straightforward. To understand the implications of this choice, we report on the effects of the thermal boundary condition on thermal convection, in terms of instability onset, fully developed transport properties and flow structure. We use the Boussinesq approximation, and impose a Robin boundary condition at the top. This enforces the temperature anomaly and its radial derivative to be linearly coupled with a proportionality factor β. Using the height H of the fluid layer, we introduce the non-dimensional Biot number Bi * = βH. Varying Bi * allows us to transition from fixed temperature for Bi * = +∞, to fixed thermal flux for Bi * = 0. The bottom boundary of the shell is kept isothermal. We find that the onset of convection is only affected by Bi * in the non-rotating case. Far from onset, considering an effective Rayleigh number and a generalized Nusselt number, we show that the Nusselt and Péclet numbers follow standard universal scaling laws, independent of Bi * in all cases considered. However, the large-scale flow structure keeps the signature of the boundary condition with more vigorous large scales for smaller Bi * , even though the global heat transfer and kinetic energy are the same. Finally, for all practical purposes, the Robin condition can be safely replaced by a fixed flux when Bi * 0.03 and by a fixed temperature for Bi * 30.

Introduction

Thermal convection is the main heat transfer mechanism in the interiors of planets and many stars [START_REF] Phillips | The physics of stars[END_REF] and is responsible for most of their dynamics. Rayleigh-Bénard convection, which is the flow developing in a plane layer subject to a potentially destabilising temperature difference, is commonly used as model for the dynamics of such astrophysical objects (Busse & Carrigan 1976;Busse 1989). It is also one of the most studied example of pattern-forming system [START_REF] Cross | Pattern formation outside of equilibrium[END_REF].

A horizontal layer of fluid heated from below and cooled from above can be unstably stratified -assuming, as is most usually the case, that the thermal expansion coefficient is positive. Indeed, in that case, the upper cold fluid is denser than the hot one at the bottom. It has been demonstrated [START_REF] Rayleigh | On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side[END_REF][START_REF] Chandrasekhar | Hydrodynamic and Hydrodynamic stability[END_REF]) that for a large enough temperature gradient, the fluid adopts a convective state in which heat transport is partially realised by fluid motion. Convection between two isothermal planes (i.e. with fixed temperature boundary condition, FT) modelled with the Boussinesq approximation constitutes the classical Rayleigh-Bénard thermal convection setup. Lord Rayleigh offered a model for such a phenomenon in 1916 Rayleigh as an explanation for the experiments of [START_REF] Bénard | Les tourbillons cellulaires dans une nappe liquide propageant de la chaleur par convection en régime permanent[END_REF]. Other setups of the Rayleigh-Bénard system have been explored, for example by choosing a fixed flux boundary condition (FF), which is closer to usual experimental conditions [START_REF] Busse | Nonlinear convection in a layer with nearly insulating boundaries[END_REF]. However, for geo-and astrophysical applications it can be difficult to choose between those two mathematical usual conditions: fixed field or fixed derivative, also known as Diriclet and Neumann boundary conditions, respectively. Indeed, in some physical situations, the boundary conditions is more accurately described by a thermal flux proportional to the temperature itself, leading to

β T + ∂T ∂r = 0 (1.1)
with β a proportionality factor homogeneous to an inverse length-scale. Using a reference length-scale H, one can construct the non-dimensional Biot number Bi = βH.

(1.

2)

The boundary condition (1.1) is intermediate between the Neumann (FF) and Dirichlet (FT) conditions which are obtained for respectively Bi = 0 and Bi → ∞, and is called the Robin boundary condition. Robin boundary conditions appear naturally in the presence of weakly conductive boundaries [START_REF] Busse | Nonlinear convection in a layer with nearly insulating boundaries[END_REF] or conductive-radiative empirical Newton's law [START_REF] O'sullivan | Newton's law of cooling-a critical assessment[END_REF]. Note that linearised black body radiative equilibrium are often interpreted as Newton's law (for example in meteorological models, [START_REF] Savijärvi | Meteorology of the terrestrial planets[END_REF]).

Thermal contact between two layers of distinct length and thermal conductivity can also be modelled by a Robin boundary condition as done by [START_REF] Guillou | On the effect of continents on mantle convection[END_REF].

Modelling the radiative equilibrium between a liquid shell and an atmosphere is our initial motivation. More precisely, we are interested in understanding the behaviour of surface global magma oceans crystallizing from the bottom upward [START_REF] Labrosse | Fractional melting and freezing in the deep mantle and implications for the formation of a basal magma ocean. The Early Earth: Accretion and Differentiation[END_REF]. Following the current scientific consensus, Earth may have undergone giant impacts during its formation [START_REF] Solomatov | 9.04 -magma oceans and primordial mantle differentiation[END_REF]. Such impacts were able to melt a large part of the nowadays solid mantle, producing what is called a magma ocean. This physical object is not known with extreme precision but was characterised, at its formation, by large depth (∼ 1 × 10 6 m), low viscosity (∼ 0.1 Pa s) and rapid flows (∼ 10 m s -1 ) [START_REF] Solomatov | 9.04 -magma oceans and primordial mantle differentiation[END_REF]. Modelling the cooling and the crystallization of this terrestrial global magma ocean is of large interest in order to get information about the early evolution of the Earth and other terrestrial planets. In that case, the radiative equilibrium between an atmosphere and the surface of the magma ocean may be modelled using a Robin boundary condition. Indeed, approximating this radiative equilibrium by an interaction between two black bodies, leads to set the heat flux at the upper boundary of the magma ocean to the difference between the radiative flux leaving the surface and the one received from the atmosphere (assumed isothermal). Using Stefan-Boltzman's law (Boltzmann 1884) (which can be derived from Planck's law) and linearizing the flux variation with temperature leads to the Robin boundary condition.

Up to now, the influence of Robin boundary conditions has received scant attention. Linear stability analyses were conducted in a Cartesian system without rotation by [START_REF] Sparrow | Thermal instability in a horizontal fluid layer: effect of boundary conditions and non-linear temperature profile[END_REF], who determined that increasing the Biot number enables transitioning monotonously from FF to FT cases (in terms of critical Rayleigh numbers). This transition was found to occur in the range Bi =0.1 to 100.

Several studies consider only the extremal cases of FT and FF, the latter being a model for the poorly conductive wall used in some experiments. In 2D Cartesian geometry, [START_REF] Chapman | Nonlinear rayleigh-bénard convection between poorly conducting boundaries[END_REF] showed, thanks to stability analysis calculations, the generation of much longer convective cells with symmetric FF, compared to the ones observed in symmetric FT configurations. A similar 3D setup was also considered by [START_REF] Busse | Nonlinear convection in a layer with nearly insulating boundaries[END_REF] who showed that square-pattern convection planforms were preferred to two-dimensional rolls produced in the case of FT boundary conditions. Later, [START_REF] Ishiwatari | The effects of thermal conditions on the cell sizes of two-dimensional convection[END_REF] explored, by numerical integration in a Cartesian twodimensional system, the effects of the different combinations of boundary conditions and internal heating on thermal convection. They observed similar effects on the wavelength of convection cells of the FF boundary conditions (independently of the presence of internal heating) as previously mentioned by Chapman and Proctor. All these Cartesian geometry studies show that FF produce convection patterns with larger wavelength than FT.

In geophysical fluid dynamics, planetary rotation is often significant and rotational effects (Coriolis force) must be included in the physical model. In fact, rotation is able to impact dramatically Rayleigh-Bénard convection, altering convection patterns (elongated in Taylor columns) and increasing the intensity of the thermal gradient needed to start convection [START_REF] Chandrasekhar | Hydrodynamic and Hydrodynamic stability[END_REF].

Adding the influence of global rotation, [START_REF] Takehiro | Linear stability of thermal convection in rotating systems with fixed heat flux boundaries[END_REF] studied horizontal layers of fluid with lateral boundaries and rotation axis inclined. They showed that, for a vertical axis of rotation, the behaviour of FF convection at the threshold becomes similar to the FT case when rotation is increased. In the case of a horizontal rotation axis, the two setups remain different, the FF one being characterised by a zero critical wavenumber.

More recently, [START_REF] Falsaperla | Stability in the rotating bénard problem with newton-robin and fixed heat flux boundary conditions[END_REF] conducted studies on a Cartesian system with Robin boundary conditions influenced by weak rotation and also observed a monotonuous transition from FF to FT when varying the Biot number.

Influence of Robin boundary conditions has been mainly studied in Cartesian geometry. On the other hand, both rotating and non-rotating convection in spherical shells have been widely studied for its geophysical and planetary relevance but mostly using traditional boundary conditions (FF or FT).

The crystallization of the terrestrial magma ocean is a complex question which cannot be reduced to the influence of Robin boundary condition. Nevertheless, we have chosen to first focus on one of the fundamental aspect of our system of interest. The goal of this work is to study Rayleigh-Bénard convection in a spherical shell subjected to a Robin boundary condition at its upper boundary, with and without rotation. We first determine numerically the characteristics of the onset of convection with a linear stability analysis, and then we explore the non-linear properties of the same system in terms of convection efficiency. In all the following we focus first on the two end-member situations, the FF and FT configurations, before studying intermediate cases by varying the value of Bi . The robustness of these observations is tested by varying the value of the Prandtl number.

Methods

Physical System

We consider a spherical shell of aspect ratio η = r i /r e = 2/3 (all the parameters and fields used are listed in the appendix A), with r e and r i the external and internal radii, respectively, subjected to rotation of vector Ωe z and filled with a Newtonian fluid of constant kinematic viscosity ν, thermal conductivity k, reference density ρ 0 , specific heat capacity c p and thermal expansivity α. The system is subjected to a radial gravity field g = -ge r , g being uniform. Note that our choice of a uniform value is relevant for thin shells and generally speaking for mantle studies (e.g. Bercovici et al. (1989)), but other radial dependencies of the gravity field can be used in other contexts (e.g. g ∝ r for constant density models, or g ∝ 1 r 2 for centrally condensed mass, or to conduct exact computations as in [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF]).

The physical behaviour of the system is basically monitored by three equations of conservation (conservation of momentum -the Navier-Stokes equation -, mass and energy) associated with boundary conditions. We assume the Boussinesq approximation, that is variations of density are neglected except in their contributions to the buoyancy force in which the density varies as function of

ρ ρ 0 = 1 -α(T -T 0 ), (2.1)
ρ being the density, while ρ 0 and T 0 are reference values of density and temperature and T the temperature anomaly around T 0 . The reference state is motionless (uniform zero velocity) and follows a radial conductive temperature profile T ref between the two boundaries with given temperatures (FT)

T ref (r) = ∆T r i r e -r i r e r -1 + T 0 , (2.2)
where ∆T is the total temperature difference of the reference profile across the shell. In the Boussinesq approximation, the value of T 0 plays no role and any value can be used such that the reference temperature profile satisfies the Robin boundary condition. In the following, temperatures are defined with respect to this reference. The problem is rendered dimensionless using r e as length scale, ν re as velocity scale, ∆T as temperature scale and ρ 0 νΩ as pressure scale. The system is characterised by the following dimensionless numbers: the Biot number Bi , the Prandtl number Pr = νρ 0 c p /k, the Rayleigh number Ra =

∂Θ ∂t + u • ∇ Θ + T ref ∆T = 1 Pr ∇ 2 Θ (2.3)
for the energy, (2.4) for the momentum and ∇ • u = 0, (2.5)

Ek ∂u ∂t + u • ∇u + 2e z × u + ∇Π = Ek ∇ 2 u + RaPr Ek Θe r ,
for the mass, with u, Θ and Π the fluid velocity, the temperature anomaly and the reduced pressure, respectively. Mechanical boundary conditions imposed on u are non-penetrative (u r = 0), stress-free (i.e. tangential stress is null, that is with [σ] the tensor of viscosity stresses of the fluid, t the tangential and n normal vector to the interface, ([σ] • n) • t = 0 at the top (r = r e ) and no-slip (u = 0) at the bottom (r = r i ). For the thermal boundary conditions, we impose a fixed temperature at the bottom (FT, i.e. Θ(r i ) = 0) and a Robin boundary condition at the top (i.e. Bi Θ(r e ) = ∂Θ ∂r (r e ), where Bi is the Biot number). These conditions can be relevant for modelling a fluid undergoing a phase change at its bottom interface and subject to a radiative equilibrium at the top surface. Indeed, pressure at the bottom is set by the fixed depth of the interface, which imposes the value of the temperature thanks to the liquid-solid equilibrium. Concerning the top interface, radiative equilibrium between a layer of fluid, with surface temperature T (r e ) and surface heat flux -k ∂ T ∂r (r e ), and an isothermal atmosphere at temperature T a leads to an energy balance,

σ(T 4 a -T 4 (r e )) = k ∂ T ∂r (r e ), (2.6)
where σ 5.67 × 10 -8 W.m -2 .K -4 is the coefficient of the Stefan-Boltzmann law giving the amount of heat flux produced by a black-body at a given temperature. We further introduce the reference depth-dependent temperature T ref such that its value at the surface is set by (2.8)

4σ(T 4 a -T 4 ref (r e )) = k ∂T ref ∂r (r e ), ( 2
Generally speaking, a Robin boundary condition applied on a field Y and its spatial derivative with respect to a space variable l at a boundary r 0 can be expressed as Bi Y + ∂ Y ∂l = 0. With that definition Bi = 0 imposes a FF boundary condition (zero spatial derivative i.e. flux) while Bi = +∞ imposes a FT boundary condition (temperature anomaly is zero at the interface).

The Biot number in the usual thermal context is written as Bi = hL/k, with h the convective heat transfer coefficient having the dimension of a W m -2 K -1 , k the thermal conductivity of the affected body and L its typical length scale. Bi can be estimated by different means depending on the type of heat transfer modelled. For material pieces in contact with a conductive-convective atmosphere, the thermal exchanges can be modelled by a thermal Newton's law with h a coefficient of the order of 10 W m -2 K -1 . As an illustration a meter-size piece of copper in contact with the atmosphere is characterised by Bi 0.1. Conversely, for a same size polystyrene piece, Bi 1000. Last example, a 5 cm glass of water cooled by contact with air in usual conditions of pressure and temperature is characterised by Bi 5. All these estimations can be computed from generic or specific handbooks (e.g. [START_REF] Kreith | The CRC Handbook of Thermal Engineering[END_REF]).

In the case described in eq. (2.8), h = 4σT 3 a so that

Bi = 4σT 3 a L k .
(2.9)

Even if our modelling remains unspecific to magma oceans (we basically model a spherical rotating shell of fluid) we should mention the estimated values of the dimensionless numbers characterizing a terrestrial magma ocean modelled by our method. Parameters given are those of a terrestrial fully liquid magma ocean.

Values are computed in both cases and reported in Tab. 1.

Name (symbol) Magma Ocean

Range in this paper Ekman number (Ek ) 3 × 10 -13 to 3 × 10 -9 10 -4 to ∞ Prandtl number (Pr ) 5 × 10 1 to 4 × 10 8 0.03 to 30 Rayleigh number (Ra) 2 × 10 18 to 9 × 10 25 10 4 to 10 9 Biot number (Bi) 3 × 10 5 to 3 × 10 11 0 to ∞ Table 1. Dimensionless numbers characterizing a terrestrial global magma ocean, maximal and minimal values are given, computed from [START_REF] Solomatov | 9.04 -magma oceans and primordial mantle differentiation[END_REF], [START_REF] Maas | Effects of earth's rotation on the early differentiation of a terrestrial magma ocean[END_REF] and [START_REF] Snyder | Experimental determination of the thermal conductivity of molten CaMgSi2O6 and the transport of heat through magmas[END_REF] (for typical values of the thermal conductivity). The extremal values for each parameter cover the full range between a fully liquid to a quasi solid (because of a high proportion of crystals in suspension) magma ocean. Finally, the range of parameters used in this study is precised.

These parameters constitute extreme conditions for numerical simulations (for example in terms of Ek ) which cannot be considered with current computational ressources. Since we focus here on a specific aspect of such a system (modelled by a Robin boundary condition) we restricted our study to a more convenient parameter space (see Tab. 1).

In addition to radiative equilibrium, Robin boundary condition can also model thermal contact between two layers of distinct conductivity (λ 1 and λ 2 ) and thickness (L 1 and L 2 ). In that case, following [START_REF] Guillou | On the effect of continents on mantle convection[END_REF], a Robin boundary condition can be written at the frontier. It links linearly dimensionless temperature and flux by a factor Bi = λ2L1

L2λ1 . Such a model could be used for internal oceans of icy satellites of Jupiter [START_REF] Sotin | Internal structure and dynamics of the large icy satellites[END_REF] or Saturn [START_REF] Collins | Enceladus' south polar sea[END_REF] with a layer of liquid water lies below an ice crust. In that case, the ratio of the conductivities would be typically around 3 [START_REF] Kreith | The CRC Handbook of Thermal Engineering[END_REF]). This means that a thick ice layer (typically 300 km) above a thin ocean (typically 10km) produces a small Bi 0.1 while a deep ocean below a thin icy crust (reverse thickness ratio) produces a large Bi 10 configuration.

Numerical Tools

We use the eigenvalue solver SINGE (freely available at https://bitbucket.org/ vidalje/singe; [START_REF] Vidal | Quasi-geostrophic modes in the Earth's fluid core with an outer stably stratified layer[END_REF][START_REF] Monville | Rotating double-diffusive convection in stably stratified planetary cores[END_REF] in order to determine the critical Rayleigh number of the system, Ra C , and, in terms of spherical harmonics, the order of azimutal symmetry m of the most unstable mode at the onset of convection in our system.

SINGE is a Python code conceived to determine the eigenvalues-eigenmodes of incompressible and stratified fluids in spherical cavities. SINGE uses the SLEPc solver [START_REF] Vidal | Quasi-geostrophic modes in the Earth's fluid core with an outer stably stratified layer[END_REF], finite difference (2 nd order scheme with N r points) in the radial direction and spherical harmonic decomposition, with degree l and azimuthal wave number m truncated at (l max , m max ). Equatorial symmetry or anti-symmetry is imposed.

The non-linear problem is solved by the XSHELLS code (freely available at https: //bitbucket.org/nschaeff/xshells). This geodynamo code has been used to solve convection in spherical rotating shell (see for example [START_REF] Kaplan | Subcritical thermal convection of liquid metals in a rapidly rotating sphere[END_REF] or [START_REF] Guervilly | Turbulent convective length scale in planetary cores[END_REF]). XSHELLS solves the Navier-Stokes equation in spherical setup thanks to the toroidal-poloidal decomposition and a pseudo-spectral approach using the fast SHTns library [START_REF] Schaeffer | Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations[END_REF]. Poloidal and toroidal scalars as well as temperature are expanded over spherical harmonics truncated in a similar way to SINGE. In the radial direction, XSHELLS uses a second order finite difference differentiation. Each simulation is run long enough to reach a statistically steady state. We then dispose of a 3D-distribution of the anomalies of the main fields (temperature, velocity etc., see fig. 1 for example). We extract from these latter 2D maps of anomalies or global quantities at a given time or compute time-average data.

Both numerical codes have been modified for this study in order to take into account the Robin boundary conditions.

Linear Stability Analysis

Differences between FF and FT cases with respect to Ekman number

We determine the critical Rayleigh number and azimuthal wave number (when it is relevant) at the linear onset of convection for different rotation rates (Ek = +∞ to 10 -6 ), in both FF and FT end-member cases. We checked that the symmetric modes are always the preferred unstable modes at the linear onset of convection (i.e.their Ra C is the smallest).

In fig. 2 we plot Ra C as a function of Ek , for different values of Pr , in order to illustrate the existence of two regimes (expected from previous studies, e.g. [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF]): in the first one (Ek → +∞) rotation does not affect Ra C while this number diminishes with Ek in the second regime. For small enough Ek , the expected power law Ra C ∝ Ek -4/3 is observed [START_REF] Chandrasekhar | The instability of a layer of fluid heated below and subject to Coriolis forces[END_REF]. Note that, although Ra C is a relevant criterion to discriminate FF and FT setups at high values of the Ekman number, it is no longer the case when rotation impacts notably the critical Rayleigh number. The effect of the Prandtl number is contrasted: the general behaviour of Ra C as a function of Ek is not affected by the value of Pr but its influence on the exact value of Ra C depends on the Ekman number. For Ek 10 -1 , Pr has no visible effect since cases are differentiated only by their boundary conditions. This is in line with the classical behaviour of nonrotating Rayleigh-Bénard convection for which the onset of convection is known to be independent of the value of the Prandtl number (e.g. [START_REF] Chandrasekhar | Hydrodynamic and Hydrodynamic stability[END_REF]). On the other hand, for Ek 10 -3 , the Pr = 0.3 cases become distinct from the cases with larger values of Pr (Pr = 3 or 30).

Let us now turn to the azimuthal structure of the flow at the onset of convection. Fig. 3 shows the critical value of the azimuthal number m for the most unstable mode for both boundary conditions as a function of the Ekman number. We observe no difference between FF and FT configurations for high rotation rates: in the limit of small values of Ek , m scales as Ek -1/3 , as expected (Busse 1970). For large values of Ek , the critical value of m tends toward a constant. In the case of Bi = 0, the transition to a constant value of m happens sharply at around Ek = 2 × 10 -3 . This particularity may be compared to the transition described by [START_REF] Falsaperla | Stability in the rotating bénard problem with newton-robin and fixed heat flux boundary conditions[END_REF] at Ek 10 -2 .

We should also mention that in the non-rotating case, Ek = ∞, the m number is not relevant. Indeed in a non rotating setup, the most unstable mode is degenerated in terms of m number. Nevertheless, in that case, the critical degree l can be determined and, with our chosen setup, is equal to 5 for a FF boundary condition and 6 for a FT one.

To observe more precisely the distinction between FF and FT cases, we plot in fig.

4

Ξ = Ra C (Bi = ∞, Ek ) -Ra C (Bi = 0, Ek ) Ra C (Bi = ∞, Ek = ∞) -Ra C (Bi = 0, Ek = ∞) , (3.1)
which quantifies the normalised difference between the FF and FT cases as a function of Ek . By definition, we have Ξ → 1 when Ek → +∞.

For Ek < 10 -3 , as shown on fig. 4, decreasing the Ekman number decreases the difference between FF and FT cases, independently from Pr . This growing similarity of Ra C (F F ) and Ra C (F T ) when rotational effects are increased is in good accordance with the conclusions of [START_REF] Takehiro | Linear stability of thermal convection in rotating systems with fixed heat flux boundaries[END_REF]. Indeed it has been shown, in the case of top and bottom FF boundary conditions and for Cartesian models in which Ω is vertical, that when Ek goes to zero Ra

C (F F ) → Ra C (F T ).
Even though high similarities of Ra C and azymuthal periodicity of the flow are observed at high rotation rates between the FF and FT cases, the structure of the solution, in particular its radial shape, depends strongly on the choice of boundary condition. In the FT case, the temperature is uniform at the boundary and lateral heterogeneity is restricted to the bulk of the domain whereas, in the FF case, temperature is variable at the boundary. This is shown later, on figure 7.

Rac

No rotation Ek = 10 3 Ek = 3.10 4 Ek = 10 5 Infinite Bi value limit

Transition from FF to FT owing to Robin boundary condition

We determine the critical Rayleigh number for different values of the Biot number with various Ekman numbers in the same range as formerly. Increasing Bi from 0 to +∞ tunes the boundary condition from a fixed flux upper boundary condition to a fixed temperature one (FF to FT transition). We can reasonably expect that the amplitude of this transition is limited by the extremal cases described in section 3.1. Since the convective setup can evolve continuously from FF to FT configuration thanks to the Robin boundary condition, we expect (as the simplest hypothesis) significant influence of this condition on the onset only in the first regime (high Ekman numbers), in which a large Ra C difference is observed.

In order to determine the evolution of the critical Rayleigh number with respect to the Biot number, we plot on fig. 5, for a given Ekman number, its normalised values,

Ra C = Ra C (Bi ) Ra C (Bi = ∞) . (3.2)
A continuous and monotonous transition can be observed in terms of the critical Rayleigh number. Without rotation, the effect is the largest, with a minimum ratio of 0.5 (Ra C (F F ) = Ra C (T F )/2). As expected from the previous section, the amplitude of the transition is reduced when the Ekman number is decreased. Finally, most of the transition occurs for Bi ∈ [1, 100].

In order to further characterise the transition of the solution as a function of the Biot number, we define

Γ = Ra C (Bi ) -Ra C (Bi = 0) Ra C (Bi = ∞) -Ra C (Bi = 0) (3.
(Bi t = 9) (Bi t = 7) (Bi t = 13) (Bi t = 30) Figure 6. Γ = Ra C (Bi )-Ra C (Bi =0) Ra C (Bi =∞)-Ra C (Bi =0
) as a function of Bi for various velocity boundary conditions at the top of the sphere and values of Ek , Pr = 3. It illustrates the transition of RaC from FF values (Γ = 0) to FT (Γ = 1) ones. If mechanical boundary condition is mentioned, this one is imposed at both boundaries, if not, the usual setup (free-slip at the top, no-slip at the bottom) is used, see 2.1 set uniformly to 1. It is worth noting that the transition shows a common behaviour for different combinations of mechanical boundary conditions (symmetric free-slip or noslip boundary conditions) and rotation rates. The only noticeable difference from one configuration to the other is the location of the transition. Indeed, we are able to model every transition by

γ = Bi Bi t + Bi , (3.4)
where Bi t is an unknown function of the Ekman number. This model, in which we estimated the suitable Bi t , shows good agreement with data (fig. 6). We observe that, for moderate rotation rates (Ek > 10 -3 ), this threshold is almost constant with Bi t ≈ 10. This sets the intermediate regime roughly between Bi = 1 and Bi = 100. The evolution of the Bi t parameter with the Ekman number is not monotonous: it first decreases with decreasing Ek then increases and becomes larger than its value for Ek = ∞ when Ek < 10 -3 . Nevertheless, as mentioned previously in section 3.1, the amplitude of the transition -in term of Ra C -becomes insignificant with such small values of the Ekman number.

We further discuss the evolution of the solution structure at the onset of convection as a function of the Biot number. Following the results obtained for the two end-member cases of FF and FT situations, discussed in section 3.1, we first study whether the azimuthal m number of the selected mode depends on the Biot and/or Ekman numbers. This enables us to build a schematic phase diagram of the flow structure at the onset, shown in fig. 7.

In the far left part (high Ek ), solutions are weakly affected by rotation and can be described as FF like or FT like structures. When the Ekman number decreases, the solution becomes affected by rotation showing smaller azimuthal number, i.e. a Nevertheless, FF and FT cases remain distinct since the radial extension is still a discriminating diagnostic. Indeed, FF configurations can be characterised by a smaller radial extension of the temperature perturbation compared to FT setups. Also, in the case of FF (and for Bi 100), the temperature perturbation extends to the boundary whereas it is zero at the boundary for FT cases.

Possible dependence on the Prandtl number

We have shown in fig. 4, that the relative difference of Ra C measured by the quantity Ξ, decreases when Ek decreases. This general behaviour is not substantially affected by the variation of Pr in the range 0.3 Pr 30 that we investigated.

In order to discuss the robustness of our observations concerning the effect of Bi against variations of the Prandtl number, we studied the FF to FT transition in a non-rotating 0 10 1 10 0 10 1 10 2 10 3 10 4 Bi

0.0 0.2 0.4 0.6 0.8 1.0 Pr=0.3 Pr=3 Pr=30 Bi T = 9 Figure 8. Γ = Ra C (Bi )-Ra C (Bi =0) Ra C (Bi =∞)-Ra C (Bi =0
) as a function of Bi with variable Pr and Ek = +∞. Γ goes to 0 (respectively 1) when the boundary condition goes to FF, i.e, Bi = 0 (respectively FT i.e. Bi → ∞).

case for Pr = 0.3 and 30. A representation of the transition in terms of critical Rayleigh number is plotted in fig. 8. We can note that the transitional regime is almost independent from Pr . In particular, the range of Biot number over which the transition from FF to FT occurs is 1 Bi 100.

Finally we consider the evolution of the flow structure in a meridional plane when Pr is varied in the presence of rotation (Ek = 10 -4 ). The radial cylindrical component u s is plotted for FF and FT setups in fig. 9 II for Pr = 0.3 or 0.03. Note that for larger value (Pr = 3 and 30), patterns similar to the Pr = 0.3 case were observed (not shown). For such large values (Pr 0.3), we observe that the flow at the onset is totally contained outside the tangent cylinder -the cylinder with axis Ω tangent to the inner boundary at the equator, see fig. 9 I. For lower values of Pr , we observe that the flow at the onset mainly develops within the tangent cylinder. This behaviour seems unaffected by the boundary condition imposed at the top interface.

In other words, when Pr is decreased, we observe a transition from a symmetric equatorial mode (localised near the equator) to a polar mode (localised near the poles, inside the tangent cylinder). We can compare our results with those for a thinner shell obtained by Garcia et al.. They determined a phase diagram of the selected modes at the onset, as a function of the rotation rate and Pr . For the same (Pr , Ek ) couples they observed another transition than ours: from an equatorial attached mode to a spiralling columnar one. This difference can possibly be attributed to the different aspect ratio. Nevertheless in a previous study conducted with a similar shell as ours, Garcia et al. observed a transition from symmetric equatorial mode to anti-symmetric polar one for Ek < 10 -6 , without mention of polar modes for higher Ek . Our results suggest that such polar mode can be found up to Ek = 10 -4 .

Non Linear Simulations

In order to investigate the behaviour of the convective flow in the non-linear regime we conducted various non-linear numerical simulations for different values of Ra, Ek , Pr and Bi (see section 2.1). We aim to determine if thermal convection subjected to a Robin boundary condition shows specific behaviours for global parameters (convection efficiency or kinetic energy) as well as local ones (structure of the flow). We first consider end-member configurations (FF and FT), before comparing them with solutions for intermediate Biot number.

In the following, we use a global diagnostic, the Nusselt number, which we define as

Nu = 1 -η η Q ∆T ef f = 1 -η η Qr e ∆T ef f k (4.1)
where Q is the total heat flux leaving the system and ∆T ef f is the observed temperature difference across the system. The tilde symbolises the dimensionless character of the physical quantities. This definition of the Nusselt number is valid in the whole range of boundary conditions (from FF to FT) as imposed by the value of Bi . In order to facilitate the understanding, we note Q ref the heat flux, computed in the motionless reference state (see eq. 2.2), leaving the system at the upper interface. When Bi → 0, in a FF situation,

Q = Q ref , Q(1-η) reη = ∆T and Nu = ∆T T (ri)-T (re) . When Bi → ∞, in a FT situation, T (r i ) -T (r e ) = ∆T , η∆T re(1-η) = Q ref and Nu = Q Q ref .
In both endmember cases, the usual definitions of Nu for FF (see for example [START_REF] Long | Scaling behaviour in spherical shell rotating convection with fixed-flux thermal boundary conditions[END_REF] or FT (see for example [START_REF] Plumley | Scaling laws in rayleigh-benard convection[END_REF]) setups are obtained and Nu provides a dimensionless measure of the efficiency of heat transfer in all cases. According to this definition, Nu = 1 in the absence of convection. Consequently Nu -1 will be preferred on some plots where Ra ≈ Ra C . Nu is used to quantify the efficiency of convection in the non-linear system in the statistically stationary regime, consequently its average over time will be presented in this regime. First of all, we study end-member configurations in terms of Bi values, that is Bi = 0 or Bi → +∞, for rotating (Ek = 10 -4 ) or non-rotating cases, for Pr = 0.3 or 3. The Rayleigh number has been set between Ra ≈ Ra C and Ra/Ra C ≈ 10 3 . As expected from numerous previous studies (e.g. [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF][START_REF] Long | Scaling behaviour in spherical shell rotating convection with fixed-flux thermal boundary conditions[END_REF], several successive convection regimes are observable when Ra increases above Ra C .

We first consider Nu = f (Ra) as represented on fig. 10. We observe the well known convective behaviour (e.g. [START_REF] King | Heat transfer by rapidly rotating Rayleigh-Bénard convection[END_REF] or [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF]): an increase of Nu with Ra starting from Nu = 1 for Ra Ra C . Rotating cases show a higher Ra C and a smaller Nu at a given Ra compared to the non-rotating cases. Since the slope of Nu = f (Ra, Ek = ∞) is also steeper (it can approach 1/2 for Ra ≈ 10 7 ) rotating and non-rotating cases show increasingly similar Nu when Ra is increased. We do not reach large enough values of Ra to see the convergence of moderately rotating and non-rotating cases that has been documented [START_REF] King | Heat transfer by rapidly rotating Rayleigh-Bénard convection[END_REF][START_REF] Gastine | Turbulent rayleigh-bénard convection in spherical shells[END_REF][START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF]. Indeed the convective Rossby number, R oc = ( RaEk 2

Pr ) 1/2 , used to compare a priori the Coriolis force to the buoyancy force [START_REF] Gilman | Nonlinear dynamics of Boussinesq convection in a deep rotating spherical shell-i[END_REF]) of Ek = 10 -4 cases is always below 10 (see appendix B). Nevertheless, for weakly rotating cases, this convergence is well observed for Ra > 10 6 . The comparison of FT and FF cases on the basis of Nu values at given Ra shows that the two configurations almost always show distinct behaviour. Generally speaking Nu(F F ) < Nu(F T ). The only exception is noted for the Ek = 10 -4 cases when Ra ≈ Ra C which shows no real difference between FF and FT configurations. This agrees with the fact that the Ra C of these cases were found with the linear stability analysis to be very similar. Indeed, it is expected (Busse & Or 1986;[START_REF] Gillet | The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder[END_REF]) that, near the threshold, the behaviour of the system is controlled by Ra Ra C -1. In the non-rotating case, for large enough Ra, it is known that the Nusselt number scales as Ra α where α is expected to be similar to 1/3 (1/3 in [START_REF] Plumley | Scaling laws in rayleigh-benard convection[END_REF] for example or 2/7 are both frequently used, see for example [START_REF] Gastine | Turbulent rayleigh-bénard convection in spherical shells[END_REF] or [START_REF] Iyer | Classical 1/3 scaling of convection holds up to ra = 1015[END_REF] for FT configuration and [START_REF] Long | Scaling behaviour in spherical shell rotating convection with fixed-flux thermal boundary conditions[END_REF], [START_REF] King | Heat transfer by rapidly rotating Rayleigh-Bénard convection[END_REF] or [START_REF] Johnston | Comparison of turbulent thermal convection between conditions of constant temperature and constant flux[END_REF] for symmetric FF one). Here, we observe α = 2/7. This behaviour is well observed in our data set for the FT configuration. Nevertheless, a different exponent (2/9 here, as justified below) seems to be necessary to model the FF data. All these general observations can be applied to both Pr = 0.3 and Pr = 3 cases.

The differences observed in fig. 10 between FF and FT cases (Nu(F T ) > Nu(F F ) at a given (Ra, Pr , Ek )) can be explained by fundamental differences of the setups. In the first case, the stationary state is characterised by an effective difference of temperature ∆T ef f smaller than the fixed ∆T imposed in the FT case. To compare the values of Nu produced by FF and FT cases on an equal footing, we measure the effective Rayleigh number Ra ef f established in the stationary state defined as (see [START_REF] Calkins | The asymptotic equivalence of fixed heat flux and fixed temperature thermal boundary conditions for rapidly rotating convection[END_REF]; [START_REF] Johnston | Comparison of turbulent thermal convection between conditions of constant temperature and constant flux[END_REF]; [START_REF] Verzicco | A comparison of turbulent thermal convection between conditions of constant temperature and constant heat flux[END_REF]):

Ra ef f = Ra ∆T ef f ∆T . (4.2)
One can note that Ra ef f = Ra when FT boundary condition is imposed at the upper interface. In order to precisely compare FF and FT configurations at the same Ra ef f we first run a FF case then determine its Ra ef f in order to run a FT case at the same Ra = Ra ef f . We plot, in fig. 11, Nu as a function of Ra ef f . One can notice that, for a given (Pr , Ek ), Nu(Ra ef f ) follows the same behaviour independently of the top boundary condition 10 4 10 5 10 6 10 7 10 8 10 9 10 10 Ra 10 0 except for cases near the threshold for non-rotating convection. This is reminiscent of the behaviour observed with the linear stability approach ( §3) that showed that the threshold is highly dependent on the boundary conditions only in the absence of rotation.

10 1 Nu Bi = ∞, Ek = 10 -4 , Pr = 3 Bi = 0, Ek = 10 -4 , Pr = 3 Bi = ∞, Ek = ∞, Pr = 3 Bi = 0, Ek = ∞, Pr = 3 Bi = ∞, Ek = 10 -2 , Pr = 3 Bi = 0, Ek = 10 -2 , Pr = 3 Bi = ∞, Ek = 10 -4 , Pr = 0. 3 Bi = 0, Ek = 10 -4 , Pr = 0. 3 Bi = ∞, Ek = ∞, Pr = 0. 3 Bi = 0, Ek = ∞, Pr = 0. 3 Nu = Ra 2/7 Nu = Ra 2/9
As expected, a fixed-flux boundary condition implies a lower critical Rayleigh-number (which has similar implication on the critical Ra ef f since Ra Ra ef f for Ra Ra C ). Finally, as shown in previous studies (e.g. [START_REF] Calkins | The asymptotic equivalence of fixed heat flux and fixed temperature thermal boundary conditions for rapidly rotating convection[END_REF] FF and FT cases become indistinguishable for high enough values of Ra ef f (here typically Ra ef f > 3 × 10 5 ), in terms of Nu and are in good agreement with the 2/7 power law mentioned above.

Since one can write, in the FF case, Ra ef f,F F ∝ Ra Nu F F , the new scaling Nu F F ∝ Ra , where u is the rms flow velocity. As one can see on fig. 12, where Pe is plotted as a function of Ra ef f , for a given (Ek , Pr ) the general behaviour far from the onset of convection is independent from the top boundary condition. For the rotating case, Pe increases sharply with Ra ef f near the onset, making it difficult to observe low values of Pe. Pr = 0.3 cases can be distinguished by a systematically smaller value of Pe for Ek = ∞ while they follow the same general behaviour as for the Pr = 3 points for rotating cases. We also observe the decrease of the gap between rotating and non-rotating configurations for high values of Ra ef f . They almost reach a unified regime characterised by a Pe ∝ Ra 1/2 scaling as expected from [START_REF] Gastine | Turbulent rayleigh-bénard convection in spherical shells[END_REF] or [START_REF] Long | Scaling behaviour in spherical shell rotating convection with fixed-flux thermal boundary conditions[END_REF] (since Pe can be deduced from the Reynolds number by the multiplication by a factor Pr independent of Ra). We plot two models Pe ∝ Ra 1/2 , one for each Pr value. Note that the difference of the proportionality factor between the models is in good accordance with a Pe ∝ Pr 1/2 dependency -that is a Re ∝ Pr -1/2 relation found with a somehow different setup in [START_REF] Calzavarini | Rayleigh and Prandtl number scaling in the bulk of Rayleigh-Bénard turbulence[END_REF], for example.

We should comment on the short range of Rayleigh numbers (around 10 8 ) for which Ek = 10 -4 cases show Pe values higher than their non-rotating counterparts, which is visible for both values of Pr . This overshoot can be explained by our use of the total velocity field to compute Pe. Indeed in a configuration of convection in a rotating spherical shell with upper free slip boundary condition, we expect the development of strong zonal flows [START_REF] Yadav | Effect of shear and magnetic field on the heat-transfer efficiency of convection in rotating spherical shells[END_REF], responsible for up to 90% of the total kinetic energy, in a certain range of Ra. We then plot on fig. 13 similar quantities as in fig. 12 but using only the non-zonal (i.e. m = 0 components) velocity field to compute Pe. With this choice, the mentioned overshoot disappears. Moreover, the two branches attributed to low and high Pr values for Ek = 10 -4 which are distinct on fig. 12 collapse in fig. 13. Consequently, the distinction between these two branches is probably solely due to the presence of a strong zonal component at low Pr that is absent or weaker in the Pr = 3 cases with Ra ef f < 10 7 . Some more comments are given about the contribution of the zonal-flows to the total kinetic energy in appendix B.

As a conclusion from the global diagnostics we just studied (Nu and Pe), for a given (Ek , Pr , Ra ef f ), the behaviour of thermal convection in our system is independent of the top boundary condition (FT or FF) if Ra ef f Ra C . Therefore, the behaviour of the system can be predicted based on results with classical boundary conditions, provided Ra ef f is used instead of Ra. Moreover, this Ra ef f can be determined a priori in the non-rotating case from Ra value. This computation can likely be extended to finite Ek values as soon as Ra is large enough to be in the weakly rotating regime.

Temperature profiles

We also compare the two extreme configurations with respect to their temperature profile. We plot in fig. 14 radial profiles of mean temperature of FF and FT cases characterised by the same Ra ef f . More precisely we compute the mean radial temperature profile for several temperature fields saved at different times in the stationary state. Then, these mean radial profiles are used to produce a time averaged radial profile T (r). We 10 4 10 5 10 6 10 7 10 8 Ra eff

10 0 10 1 Nu Bi = ∞, Ek = 10 -4 , Pr = 3 Bi = 0, Ek = 10 -4 , Pr = 3 Bi = ∞, Ek = ∞, Pr = 3 Bi = 0, Ek = ∞, Pr = 3 Bi = ∞, Ek = 10 -2 , Pr = 3 Bi = 0, Ek = 10 -2 , Pr = 3 Bi = ∞, Ek = 10 -4 , Pr = 0. 3 Bi = 0, Ek = 10 -4 , Pr = 0. 3 Bi = ∞, Ek = ∞, Pr = 0. 3 Bi = 0, Ek = ∞, Pr = 0. 3 Figure 11.
Nu as a function of Ra ef f (4.2) for the same data as in fig. 10. For a given Pr , FF and FT points follow the same behaviour except near the threshold for non-rotating convection. Nu ∝ Ra 2/7 is plotted with a blue dotted line.

finally plot T -Tmin

Tmax-Tmin in order to compare the shape of the profiles since, by construction, ∆T ef f is the same in these two configurations. We can see in fig. 14 that temperature profiles of FF and FT configurations are almost identical when Ra ef f is large enough. In the case of a weakly surcritical simulation (Ra Ra C ), temperature profiles are more similar to diffusive ones since convection is not well developed.

Flow structure

Convection can also be studied on the basis of the spatial structure of the flow. We consider first the structures in the meridional plane and then in equatorial ones.

The meridional slices enable us to investigate the disposition of the flow with respect to the tangent cylinder in the rotating cases (without rotation this concept is not relevant). As shown in fig. 15(a) moderate supercriticality implies that flow is not active in the tangent cylinder. When Ra is increased enough (b), this spatial organization is lost and the flow occupies the whole shell. Nevertheless in both cases, the influence of rotation can be noticed thanks to the alignment of the flow along columns of axis Ω. RaC , FF and FT configurations follow the same behaviour at a given Ek , Pr . For large enough Ra ef f , Pe becomes independent of Ek . Solid and broken lines are 1/2 slope scalings for Pr = 3 and Pr = 0.3.

As a remark, we noticed in some cases (e.g. Bi = ∞, Ek = 10 -4 , Pr = 3, Ra = 3 × 10 7 i.e. in rotating cases at moderate Ra/Ra C ) the development of an anti-symmetric component of the velocity field with respect to the equatorial plane. This component was never stronger than the symmetric one. The anti-symmetric component of the flow is mostly seen within the tangent cylinder where the Taylor-Proudman constraint does not link both hemispheres.

Then, we choose to compare the respective properties of structures at equal Ra/Ra C in order to compare slightly supercritical configurations. For higher Ra, equal Ra ef f plots are compared since it has been shown that this number was an accurate criterion to fairly compare FF and FT configurations.

We can see on fig. 16 that it is usually possible to discriminate FF from FT cases without consideration of (Ek , Pr ) by checking if temperature anomalies extend to the upper boundary or not. This is illustrated clearly in the non-rotating configuration (see first line for slightly surcritical case and third one for higher Ra/Ra C ) and in the rotating case (see second line for slightly surcritical simulations and last one for moderate Ra/Ra C ). Except for the vertical extent of the thermal anomalies near the upper boundary, it is nearly impossible to visually distinguish between FF and FT equatorial temperature snapshots on the basis of the size of the anomalies or their spatial organization. We extend this comparison to 3D structures in the appendix D.

In the slightly supercritical rotating cases (second line) we observe the lateral periodicity (m = 17) predicted by the linear stability analysis (see fig. 7 and3), which is identical for both FF and FT boundary condition. It is worth noting to remark the high temperature equatorial anomaly which is characteristics of the fixed flux configuration.

FF to FT transition

We now focus on the transition from FF to FT configurations produced by an increase of the Biot Number from 0 to infinite limit values. We restrict our set of parameters to the Pr = 3 simulations. For various Ra and Ek we determined Nu when the Biot number evolves from one limit to the other.

Transition with respect to Nu

Intermediate Bi values cases can be compared with FF and FT ones on the basis of their Nu, as plotted on fig. 17 as a function of Ra ef f . They follow the same behaviour as for FF and FT cases, which shows that Ra ef f constitutes the relevant variable characterizing this different types of thermal convection subjected to different boundary conditions. We did similar observations on the basis of Pe rather than Nu.

As in the section 4.1.1 we compute Ra ef f as a function of Bi for (Ek = ∞, Pr = 3). We found that for fixed Ra, Ra ef f increases continuously and monotonically from FF to FT values following

Ra ef f (Bi ) = Ra ef f,F F + ∆Ra ef f Bi Bi + Bi t , (4.4)
with ∆Ra ef f the difference between FF and FT values of Ra ef f (which can be determined by eq. 4.3) and Bi t a threshold of approximate value 10 -we noticed that Bi t increased slightly with Ra. Our data show a good agreement with this model on three orders of magnitude of Ra (discrepancies are always smaller than 10%).

Transition for a given effective Rayleigh number

Even if equal Ra ef f implies similar Nu regardless of Bi , is it possible to detect differences in the flow structure at constant Ra ef f with respect to the Bi value? As an illustration we plot the radial profiles of temperature (fig. 18) of seven cases of approximately equal Ra ef f and variable Bi . Note that we could not obtain a collection of cases with exactly constant Ra ef f but variable (Ra, Bi ), since the value of Bi t used in eq. 4.4 evolves slightly with Ra, nevertheless differences are below 10%. These profiles have been time averaged over fields extracted from the statistically stationary regime.

Here we focused on non-rotating systems with Ra Ra C . As one can see, the general structure of the profiles is the same: the temperature difference is concentrated in narrow boundary layers located at both shell limits, while the interior of the shell is homogenised by convection and exhibits a slightly stabilizing temperature profile. While the profiles are nearly identical near the boundaries, the middle temperature is affected by Bi . This dependency is not monotonous: intermediate (i.e. between 1 and 30) Bi profiles have the lowest mean temperature, FF ones show medium temperature and FT cases the highest mean temperature.

In order to give a quantitative description of the spatial structure of the flow, we plotted the velocity spectrum (fig. 19) for each of theses cases. With the information provided by the spectra it is possible to identify the three regimes suggested by the radial profiles of the fig. 18. Indeed, the low, intermediate and high Bi value cases are characterised respectively by a dominance of l = 1, 2 and 3. Finally small-scale flows are favoured at higher Bi (that is FT-like top boundary condition) while the different spectra become indistinguishable as soon as l 7.

To conclude this section, we have shown that even if Ra ef f is the relevant parameter to compare convection efficiency of FF and FT cases, differences in the flow structure are observed at constant Ra ef f but variable Bi (without rotation). For rotating configurations (particularly at Ek = 10 -4 ), it has been impossible to bring out important differences of the flow structure between FF and FT cases. Usually these successive regimes are spotlit by scaling laws deduced from Nu = f (Ra) or Pe = f (Ra) graphs. Here we choose to combine these two in a Pe N Z = f (Nu -1) graph (more usual graphs have been presented in the beginning of this section 4.1, see fig. 12 and11). The non-zonal component of Pe is solely used since its exclude zonal winds that do not participate to the radial thermal transfer. We observe on fig. 20 different branches on which points are distributed mainly following their Ekman number (Ek = 10 -2 points are poorly distinguishable from non-rotating ones). For low Pe and Nu, two branches of specific Ek are visible. They are well modelled by a Pe ∝ (Nu -1) 1/2 relation. For higher Nu, these two branches join a new one with a steeper slope (Pe ∝ (Nu -1) 7/4 ). The non-rotating branch joins this steeper branch for Nu 1 while the rotating branch joins it for Nu 20. We observe no strong influence of Bi or Ek on these branches except in the non-rotating one for small Nu -1: FT points seem to show systematically slightly higher Pe compare to FF ones. Remarkably, there is almost no dependence on Pr in this representation -although only Pr = 0.3 and 3 are considered.

We further interpret these different branches relatively to convection regimes. Small Pe and Nu -1 are characteristic of a slightly surciritical convection that is a weakly non-linear (WNL) regime in which the convective flow is mainly monitored by Ra Ra C . From usual scaling laws (see Busse & Or (1986) andGillet &Jones (2006) 2.2 × 10 6 but variable (Ra, Bi), (Ek = ∞, Pr = 0.3). Similar profiles are obtained, nevertheless, three different behaviour can be distinguished: for large (FT), small (FF) and intermediate Bi values. Exact values of (Ra ef f , Ra) for each Bi are: Bi = 0, Ra ef f = 2.18 × 10 6 , Ra = 9 × 10 6 ; Bi = 0.1, Ra ef f = 2.24 × 10 6 , Ra = 9 × 10 6 ; Bi = 1, Ra ef f = 2.14 × 10 6 , Ra 6.6 × 10 6 ; Bi = 3, Ra ef f = 2 × 10 6 , Ra = 4.5 × 10 6 ; Bi = 10, Ra ef f = 2.2 × 10 6 , Ra = 3.3 × 10 6 ; Bi = 100, Ra ef f = 2.2 × 10 6 , Ra = 2.34 × 10 6 ; Bi = ∞, Ra ef f = 2.18 × 10 6 , Ra = 2.34 × 10 6 .

the WNL regime, we expect that on the one hand

Nu -1 ∝ Ra Ra C -1 , (4.5)
and on the other hand, when the mean thermal profile remains similar to the diffusive one (which is appropriate for Ra Ra C ), we expect following the equilibrium between viscous dissipation rate and mean buoyancy power (see eq. 3.3 in [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF] and eq. 2.16 in [START_REF] Gastine | Turbulent rayleigh-bénard convection in spherical shells[END_REF]) points between a rotating and a non-rotating branch can be understood. As shown in the linear stability analysis section 3, rotation impacts poorly the onset of convection for Ek > 10 -3 , which explains why results for Ek = 10 -2 fall on the non-rotating branch. Since both Ra C and L C can be expressed as Ek functions using appropriate scaling laws, a Pe ∝ Ek -1/2 dependency can be predicted. Nevertheless we are not able to test it since the scaling laws are established for smaller Ek than those we used here. Finally, Bi has generally a weak or no influence on Ra C except in the non-rotating configuration which explains the small difference between FF and FT points in the non-rotating branch.

Ra(Nu -1) Pr 2 ∝ νU 2 C L 2 C , ( 4
The highest values of Pe, Nu are characteristic of the non-or weakly rotating regime. We presented above ( § 4.1.1) the usual scaling laws used in the non-rotating limit for Pe and Nu (see fig. 12 and10). The combination of both these laws leads to Pe ∝ Nu 7/4 , (4.8) which is the slope observed in fig. 20, for high Nu -1 Nu. Note that our scaling law is very similar to the one found for this regime in another study of thermal convection with a different system and mechanical boundary conditions [START_REF] King | Turbulent convection in liquid metal with and without rotation[END_REF]. In the non-rotating branch, Ek is no longer a discriminating factor as expected. Remarkably, lower Pr points are not easily distinguished from the Pr = 3 group of points. This is worth noting since both Pe and Nu exhibit dependencies on Pr (see Figs. 11 and 12). Between these WNL and non-rotating regimes, a transitional regime is expected for rotating simulations. We observe it roughly between Nu 1.1 and Nu 20, it is characterised by an intermediate slope between the two end-members. Following a former study [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF] we expect a somehow higher (Nu 1.5) value of Nu as a lower limit of this regime. Note that since we never use Ek smaller than 10 -4 , we do not expect to significantly observe a rapidly rotating regime [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF][START_REF] Long | Scaling behaviour in spherical shell rotating convection with fixed-flux thermal boundary conditions[END_REF].

As a conclusion, we observe that variable Bi convective simulations follow the same general behaviour in term of convective regimes as the two better know end-members.

Discussion

Our objective was to study the influence on thermal convection of Robin boundary condition (characterised by the Biot number) on a rotating spherical shell of fluid heated from below by an isothermal lower interface. This boundary condition has been implemented on the anomalies of temperature and thermal flux at the top interface.

In order to compare with different situations, it appears more convenient to use a more universal Biot number Bi * = Bi/3 which is based on the fluid height H = r e /3 instead of the radius of our sphere r e . We used, first, linear stability analysis to determine the properties of this boundary condition through the transition from fixed flux to fixed temperature at the upper boundary (see 3). This transition is obtained by varying the Biot number which monitors the linear relationship between the flux and the temperature for a Robin boundary condition at the interface. We observe a continuous transition in terms of critical Rayleigh number of the system between the two end members. This transition is significant only at weak rotation for which the two end-member regimes (FF and FT) remain distinct. In the regime of high Ek (roughly Ek > 10 -4 ), a transition occurs approximately for Bi * ∈ [0.3, 30]. We give a simple model for the evolution of the critical Rayleigh number as a function of Bi * during this transition. Outside of this range, the general behaviour of the onset of convection can be described as the traditional FF or FT configuration. These results are in good accordance with previous work of [START_REF] Sparrow | Thermal instability in a horizontal fluid layer: effect of boundary conditions and non-linear temperature profile[END_REF]. In terms of solution structure at the onset of the convection, we offer a schematic phase diagram showing the existence of FF-like, FT-like or intermediate structures depending on Ek and Bi * . We have extended our results concerning the FF to FT transition, at least without rotation, to other velocity boundary conditions and other Pr numbers.

Non-linear simulations have been conducted in a second time (see 4) to study vigorous convection. We have shown that the relevant parameter to compare convective cases subjected to the whole range of boundary conditions authorised by a variable Bi * is a Rayleigh number based on the effective mean temperature difference across the shell. We provided a simple model for the evolution of Ra ef f as a function of Bi * for non rotating configurations.

In terms of global diagnostics quantifying convection (for example the Nusselt number or the Péclet number), the transition at a given Ra ef f shows valuable amplitude only for non-rotating and Ra ef f Ra C configurations. Indeed, in this case Ra C depends significantly on Bi * , while by construction Ra ef f ignores this effect. For other values of the parameters, the transition amplitude is considerably reduced. Concerning the organization of the flow, the transition is observable as long as the solution remains not too turbulent (that is Ra ef f Ra C ). For respectively small and large enough Bi * the system behaviour is very similar to the end-member setups. It is only when Bi * ∈ [3.10 -2 , 3.10 1 ] that intermediate specific structures can be noticed while global diagnostics remain unchanged quasi identical to those of FF and FT configuration. Since Ra ef f is not an input parameter, we also studied the transition at fixed Ra and gave a model for the evolution of Nu as a function of Bi (see appendix C). We observe a noticeable transition only for Bi * ∈ [3.10 -2 , 3.10 1 ], this one can be modelled as a simple function of Bi . The transition is in fact characterised by a unique parameter determined by (Ra, Pr , Ek ).

Consequently, the range of Biot number for which the results significantly differ from end-member FF or FT boundary conditions is quantified in this paper. Depending on the object of application, computing the Biot number allows one to decide whether it is worth considering the full Robin boundary condition or, otherwise, which end-member model to use. In practice we have shown that Robin boundary conditions can be safely replaced within the framework described in 2.1 -by fixed flux boundary condition for Bi * < 0.03 and by fixed temperature for Bi * > 30.

Considering the application to convection in magma ocean planets, we obtained (following the definition given in eq. 2.9 and the references given in 2.1) a Biot number of order 10 10 which implies that a fixed temperature at the upper surface is the proper boundary condition to apply.This remains true even if we should bear in mind that a geophysical object like a terrestrial magma ocean evolved with time by cooling and freezing. We must expect that the huge Bi * computed here decreases when the ocean evolves and becomes thinner. Indeed Bi * is directly proportional to the thickness of the system, inducing a reduction of a factor 100 to 1000 of the Biot number with the crystallization of the ocean.

Nevertheless it is only for very thin layers of high thermal conductivity below a cold atmosphere that FT modelling becomes irrelevant. For a caricatural example, a liquid mercury layer of a hundred of meters below an atmosphere at 300 K sets a Biot number of about 60 which places this system in the intermediate regime. Extreme cases like metric sized liquid layers authorise low Biot number, that is FF boundary conditions. We emphasize that this observation concerns only boundary conditions applied to anomalies and not to the determination of reference profiles which is a distinct work.

However, we also mentioned in the introduction other physical thermal systems which can be modelled using Robin boundary condition. One of them is the superposition of two thermally conductive layers. In that case it is easier to produce configurations for which Bi < 1, as thin water ocean (of thickness H) below thick ice layer (thickness H S ). Except top mechanical boundary condition, our physical model is appropriate to describe this kind of systems (namely internal oceans of icy satellites). Furthermore, in this situation, the boundary condition depends on the lateral length-scale of the temperature anomaly. For large length-scale, the Biot number is unchanged, but for small length-scales H S an effective Bi = Bi H/ (see [START_REF] Guillou | On the effect of continents on mantle convection[END_REF]. In this context, it may be appropriate to use a Robin boundary condition with a Biot dependent on the length-scale, or at least to use fixed-temperature for the smallest scales of the system.

We have not considered the influence of a symmetric top and bottom Robin boundary conditions or the effects of varying the aspect ratio. However, our results bear much generality, and a are very likely relevant to most convection setup, as soon as Nu and Ra ef f are used to characterise the system.

Global and deep oceans of hypothetical exo-planets could constitute another appropriate geophysical subject for such modelling. Robin boundary condition could be chosen to model various relations between a field and its derivative at an interface. We could, for example, contemplate modelling partial crystallization of a two component liquid using this method. Indeed the evolution of the liquid composition C at the solid/liquid interface during the crystallization (due to the fact that the produced solid does not have the same composition as the initial liquid) could be modelled by fixing a chemical flux j C determined by the local composition of the fluid. Such a relation could be made linear by a relation similar to K C + j C = 0 with K a coefficient similar to the Biot number. that we do not really observe the disappearance of the influence of the rotation on such simulations (see fig. 10 for example).

Appendix C. Transition at a given Ra

It was also possible to study the transition observed at constant Ra, when Bi varies from 0 to infinity. As shown previously this choice is not optimal to compare equivalent configurations but only uses input parameters.

Fig. 22 shows Nu -1 as function of Bi for various values of (Ra, Ek ). The amplitude of the transition by definition depends on Nu(Bi = 0) -Nu(Bi = ∞) which is a function of (Ra, Ek , Pr ). For all cases but one, Nu increases with Bi -at equal Ra, FT cases are more efficient in terms of heat transfer. The exception is explained below. Fig. 22 shows that in the most common cases the transition from FF to FT is monotonous and continuous. Moreover, the transition occurs mainly between Bi = 10 -1 and Bi = 10 2 .

In order to study more specifically the transition itself, we then plot for a given global diagnostics X(Ra, Ek , Bi ) (here X = Nu -1) the quantity

Σ = X -X(Bi = ∞) X(Bi = 0) -X(Bi = ∞) . (C 1)
By construction, Σ evolves from 0 to 1 when Bi goes from 0 to +∞ and it can be interpreted as a transition parameter. This quantity has been computed for the same subset that was used to produce fig. 22 and is plotted on fig. 23. First of all, we notice 10 -2 10 -1 10 0 10 1 10 2 10 3 ← FF Bi FT → 10 0 10 1

Nu

-1 Ra = 9. 00 × 10 4 , Ek = ∞ Ra = 9. 00 × 10 6 , Ek = 10 -4 Ra = 3. 00 × 10 7 , Ek = 10 -4 Ra = 9. 00 × 10 7 , Ek = ∞ that the large majority of the data subsets associated to a given (Ra, Ek ) couple follows the same continuous and monotonous transition. The quantity

Σ 0 = Bi Bi + Bi T (C 2)
approximates correctly the evolution of Σ, when the central value Bi T ≈ 3 to 10. Note that Σ is very similar to the quantity γ used to model Ra C = f (Bi ) in section 3.2. We observe a slight evolution of Bi T with Ra and Ek . The only notable exception to these observations is the slightly supercritical, Ra = 9 × 10 4 non-rotating case (blue dots), which is anomalous. For this case, the transition is not monotonous and is characterised by a minimum for Bi ≈ 1. In fact, this case is one of the few configurations in which the most efficient heat transfer is obtained for intermediate value of Bi . This specific behaviour is linked to the vicinity of the critical Rayleigh number in this simulation as well as to the large relative difference between Ra C for FF or FT boundary conditions in this case. As seen in the section 3.1, Ra C (FF, Pr , Ek ) < Ra C (FT, Pr , Ek ) and Σ Ra = 9. 00 × 10 6 , Ek = 10 -4 Ra = 3. 00 × 10 7 , Ek = 10 -4 Ra = 2. 34 × 10 6 , Ek = ∞ Ra = 3. 00 × 10 6 , Ek = 10 -4 Ra = 9. 00 × 10 4 , Ek = ∞ Ra = 9. 00 × 10 5 , Ek = 10 -2 Ra = 6. 00 × 10 6 , Ek = 10 -4 Σ 0 (Bi T = 3) Σ 0 (Bi T = 5) Σ 0 (Bi T = 10) 2006)). Thus, for this case characterised by (Ra = 9 × 10 4 ≈ Ra C , Ek = ∞, Pr = 3) we can reasonably expect that Nu(FF) > Nu(FT) which means that increasing Bi leads to a transition from higher to smaller values. Moreover, since Ek = ∞ the relative amplitude of the transition is large and implies high variations of Ra Ra C with Bi since Ra remains constant. This produces the non-monotonous variation of Nu with Bi .

We also consider other global diagnostics: Et which is the root mean square value of the temperature anomaly and P conv which is the root mean square of uΘ and found similar behaviour as for Nu with respect to Bi .

  .7) and the depth-dependent temperature anomaly θ = T -T ref . Assuming θ T 0 , equation (2.6) can be expanded to give -4σT 3 ref (r e )θ = k ∂ θ ∂r (r e ).

Figure 1 .

 1 Figure 1. 3D visualization of radial velocity (left) and temperature anomaly (right) of a nonlinear computation. Simulation run for Ra = 9 × 10 6 , Pr = 3, Ek = 10 -4 and Bi = 0. Internal and external surfaces have been defined slightly outside boundary layers.

Figure 2 .

 2 Figure 2. Critical Rayleigh number RaC of FF and FT cases as a function of Ekman number Ek for Prandtl number values Pr = 0.3, 3 and 30. For evanescent Ek RaC scales as Ek -4/3 (represented by the dotted line).

Figure 3 .

 3 Figure 3. Values of the critical azimuthal wave number m as a function of Ek for both FF (Bi = 0) and FT (Bi = ∞) cases. For small enough values of Ek , m scales as Ek -1/3 (represented by the dashed line).

Figure 4 .

 4 Figure 4. Ξ i.e.difference of RaC between FF and FT cases for a given Ek normalised by the non rotating difference (cf eq. 3.1) as a function of Ek . Ξ decreases with the same general behaviour independently from Pr .

Figure 5 .

 5 Figure 5. RaC (cf eq. 3.2) as a function of Bi for Pr = 3. The amplitude of the continuous transition from FF to FT is decreased when Ek is reduced. By definition, RaC → 1 when Bi → ∞ (i.e.when the boundary condition goes to FT).

  3)which is plotted as function of Bi on fig.6, for various values of the Ekman number. This quantity allows us to compare the behaviour of the transition regardless of its amplitude rotation, free slip BC Ek = 10 3 , free slip BC No rotation, no-slip BC Ek = 10 5 , no-slip BC Ek = 3.10 3 Ek = 3.10 4 Ek = 10 5

Figure 7 .

 7 Figure 7. Phase diagram of the solution structure at the linear onset of convection as a function of Ek and Bi, typical structures of temperatures anomalies and dependence on Bi or Ek is mentioned for each area. I: FT solution independent of rotation with constant value of m, II: FF solution independent of rotation with constant value of m, III: FT solution affected by rotation (m increases with rotation), IV: intermediate structure, V: FF solution affected by rotation, VI: common structure at high rotation rate.

Figure 9 .

 9 Figure 9. I: Schematic of the tangent cylinder (grey) as a part of a spherical shell rotating about an axis. II: Quarter of meridional slices representing radial cylindrical component of the velocity at the onset of convection for FF (c and d) and FT (a and b) top boundary conditions and Pr = 0.3 (b and d) or 0.03 (a and c), Ek = 10 -4 . We observe an effect of Pr on the flow, which appears outside of the tangent cylinder for Pr = 0.3 and inside it for Pr = 0.03. Since these maps are those of the flow at the onset, they are defined within any arbitrary multiplicative factor.

4. 1 .

 1 Behaviour of purely FT or FF setups 4.1.1. Global diagnostics

Figure 10 .

 10 Figure 10. Nusselt number, Nu, as a function of Ra for Pr = 0.3 (open symbols) or 3 (filled symbols) with (Ek = 10 -4 : squares and Ek = 10 -2 : circles) or without (Ek = ∞ circles) rotation and with fixed flux (red points) or fixed temperature (blue ones) at the upper boundary. The red dotted line represents the scaling Nu ∝ Ra 2/7 relevant to model FT configurations while the Nu ∝ Ra 2/9 -blue dotted line-models the FF ones.

  Ra ef f,F F ∝ Ra 7/9 F F and Nu F F ∝ Ra 2/9 which was observed on fig.10. Using our data to determine the coefficient of proportionality between Ra ef f,F F and Ra F F , we found for non rotating convection and Pr = 3 that, far enough from the onset of convection and along at least three orders of magnitude, Ra ef f,F F 8.6 × Ra below 5%. The similarity of the FF and FT configurations observed at equal Ra ef f on the basis of Nu prompted us to look at other quantities. Aside from determining Nu of each configuration, we computed their Péclet Number in statistically steady state, Pe = re uρ0cp k

Figure 12 .

 12 Figure 12. Péclet number Pe as a function of Ra ef f for the same data as plotted in fig. 11. Except for Ra ef fRaC , FF and FT configurations follow the same behaviour at a given Ek , Pr . For large enough Ra ef f , Pe becomes independent of Ek . Solid and broken lines are 1/2 slope scalings for Pr = 3 and Pr = 0.3.

Figure 13 .

 13 Figure13. Non-zonal component of Pe as a function of Ra ef f for Ek = 10 -4 cases. Similar behaviour is observed regardless of (Pr , Bi): Pe values approaches the non-rotating limit when Ra ef f is increased. The overshoot observed in fig.12is suppressed with the subtraction of the zonal component.

Figure 14 .

 14 Figure 14. Radial temperature profiles averaged in time and normalised by the mean temperature difference across the shell, T -T minTmax-T min . Reference, FF and FT profiles are compared, Ek = ∞, Pr = 3, Ra ef f ≈ 5.6 10 6 . For large enough Ra ef f , equal Ra ef f configurations show similar profiles.

Figure 15 .

 15 Figure 15. Meridional (a and b) and equatorial (c and d) maps of radial velocity. a and c: Bi = ∞, Ek = 10 -4 , Pr = 3, Ra = 9 10 6 , b and d: Bi = ∞, Ek = 10 -4 , Pr = 3, Ra = 9 10 7 . a: the flow is only active outside the tangent cylinder, b: the flow is active in the whole shell. c and d: flow loses its lateral organization when Ra is increased of a factor 10.

Figure 16 .

 16 Figure 16. Equatorial maps of temperature anomalies (l = 0 component is subtracted, Pr = 3 for all cases); left Bi = 0 (FF), right Bi = ∞ (FT); c,d,g,h: rotating cases (Ek = 10 -4 ), a,b,e,f: non rotating ones; a,b,c and d: Ra/RaC 1.1, e,f,g and h: Ra ef f = 3.10 8 . Except difference of thermal upper boundary layer, structures of FF and FT thermal anomalies are similar for each line.

Figure 17 .

 17 Figure 17. Nusselt number Nu as a function of Ra ef f for the same representation code as in fig. 10. Intermediate Bi cases have been added to the fig. 10 data. They are plotted with big symbols, coloured according to their Bi value. As in the fig. 10, points are distributed along two main branches determined by the value of Ek . Except for (Ek = ∞, Ra ≈ RaC ) points, intermediate Bi points are located on the same branches as extremal ones and then follow the same scalings as described in fig. 10.

Figure 18 .

 18 Figure 18. Normalised and time-averaged radial temperature profiles for quasi constant Ra ef f2.2 × 10 6 but variable (Ra, Bi), (Ek = ∞, Pr = 0.3). Similar profiles are obtained, nevertheless, three different behaviour can be distinguished: for large (FT), small (FF) and intermediate Bi values. Exact values of (Ra ef f , Ra) for each Bi are: Bi = 0, Ra ef f = 2.18 × 10 6 , Ra = 9 × 10 6 ; Bi = 0.1, Ra ef f = 2.24 × 10 6 , Ra = 9 × 10 6 ; Bi = 1, Ra ef f = 2.14 × 10 6 , Ra 6.6 × 10 6 ; Bi = 3, Ra ef f = 2 × 10 6 , Ra = 4.5 × 10 6 ; Bi = 10, Ra ef f = 2.2 × 10 6 , Ra = 3.3 × 10 6 ; Bi = 100, Ra ef f = 2.2 × 10 6 , Ra = 2.34 × 10 6 ; Bi = ∞, Ra ef f = 2.18 × 10 6 , Ra = 2.34 × 10 6 .

Figure 19 .

 19 Figure19. Time average energy spectrum of the velocity (sum of toroidal and poloidal components over the whole shell) for roughly constant Ra ef f 2.2 × 10 6 , Ek = ∞, Pr = 3 but variable Bi. For Bi 0.1, l = 1 is the dominant mode, intermediate values (1 Bi < 10) favour l = 2, and l = 3 dominates at high Bi like in FT configuration. The total kinetic energy in all these configurations varies by less than 10%.

W

  Figure20. Non-zonal component of Pe as a function of Nu -1 for various Ek , Pr , Bi coded by shape, filling and colour of the symbols respectively (as in fig.17 and 10). We distinguish two main regimes: the weakly non-linear regimes at low Pe and Nu 1 following the two broken-line models and the weakly or non-rotating regimes represented by the dotted line.

Figure 21 .

 21 Figure 21. Kinetic energy ratio, E E N Z as a function of the convective Rossby number Roc for various (Ek , Pr , Bi). Legend is similar to fig. 22. Zonal flows are mainly active for small Ek and Roc ∈ [0.1, 1].

Figure 22 .

 22 Figure 22. Nusselt Number as a function of the Biot Number for various (Ra, Ek , Pr = 3) simulations. Each data subset defined by (Ra, Ek ) is plotted by symbols of the same colour and shape. FF (Bi = 0) and FT (Bi → ∞) limits are represented respectively by solid and dotted lines of the same colour.

(

  Ra C FT,Ek )-Ra C (FF,Ek ) Ra C (FT,Ek =∞)-Ra C (FF,Ek =∞) decreases with Ek . Moreover, Nu -1 is proportional to Ra Ra C -1 near the threshold of convection (Gastine et al. (2016) and Gillet & Jones

Figure 23 .

 23 Figure 23. Σ transition parameter (cf. eq. C 1, computed with X = Nu -1) from FF (Bi = 0) to FT (Bi = ∞) as a function of the Biot Number for various (Ra, Ek , Pr = 3) simulations. Σ0 (cf. eq. C 2) are plotted for three values of BiT , BiT = 3, 5, 10

Figure 24 .

 24 Figure 24. 3D representation of the temperature anomaly (l=0 subtracted) inside a shell of radial limits 0.7 and 0.95 with Pr = 3. Rotational axis is z. Left column FF, right one FT. For each line (Ra ef f , Ek ) is kept approximately constant. First line (3 × 10 8 , ∞). Second line (3 × 10 8 , 10 -4 ). Third line (10 7 , 10 -4 ).

† Email address for correspondence: thibaut.clarte@univ-grenoble-alpes.fr

Acknowledgements

Computations were performed on the Froggy platform of CIMENT (https: //ciment.ujf-grenoble.fr), supported by the Rhône-Alpes region (CPER07 13 CIRA), OSUG2020 LabEx (ANR10 LABX56) and EquipMeso (ANR10 EQPX-29-01). ISTerre is part of Labex OSUG2020 (ANR10 LABX56). JV acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 847433). The opensource numerical codes xshells and singe are available at https://nschaeff. bitbucket.io/xshells/ and https://bitbucket.org/vidalje/singe/downloads/. Most figures were produced using matplotlib (http://matplotlib.org/) or paraview (http://www.paraview.org/).

Declaration of interest

The authors report no conflict of interest. Pr ) 1/2 is frequently used as a parameter to study the evolution of zonal flows. R oc estimates roughly the ratio of buoyancy to Coriolis forces [START_REF] Gilman | Nonlinear dynamics of Boussinesq convection in a deep rotating spherical shell-i[END_REF]. Small R oc implies strong influence of rotation while the infinite limit corresponds to a flow not affected by rotation. In order to study this phenomenon we plot on fig. 21, the ratio of the total kinetic energy E over the non-zonal component of this energy, E N Z , as a function of R oc . The larger is the ratio, the stronger is the zonal flow.

Appendix A. Notations

For the two values of Ek , we observed the zonal flow first grows with R oc then decreases. Moreover the zonal flow is almost negligible for the weakly rotating configuration (a flat maximum is observed for R oc 8). For Ek = 10 -4 we observe a maximum for R oc 0.5 at which the zonal flow represent more than 90% of the total kinetic energy of the system. This maximum is compatible with previous results realised with different boundary conditions (see [START_REF] Yadav | Effect of shear and magnetic field on the heat-transfer efficiency of convection in rotating spherical shells[END_REF] or [START_REF] Christensen | Zonal flow driven by strongly supercritical convection in rotating spherical shells[END_REF]). We also observe that Pr and Bi do not seem to alter significantly the general behaviour of the zonal flows, only small modification of the maximum being noticed. Finally, fig. 21 shows that the Ek = 10 -4 data do not reach very high R oc values (∼ 3 at most) which explains Appendix D. 3D representations Fig. 24 shows 3D maps of temperature anomalies (l=0 component subtracted) for subshells of our system -boundary layers excluded -for Pr = 3. FF and FT configurations are compared for quasi equal Ra ef f , Ek . For cases with Ra ef f = 3 10 8 , the effect of boundary condition is hardly visible. Conversely, cases with lower Ra ef f with a strong effect of rotation show differences in the location of the anomalies with respect to the tangent cylinder (see fig. 9 I).