A Table of results for the numerical models with different geometry and gravity profiles

# Turbulent Rayleigh-Bénard convection in spherical shells

## Abstract

We simulate numerically Boussinesq convection in non-rotating spherical shells for a fluid with a unity Prandtl number and Rayleigh numbers up to . In this geometry, curvature and radial variations of the gravitational acceleration yield asymmetric boundary layers. A systematic parameter study for various radius ratios (from to ) and gravity profiles allows us to explore the dependence of the asymmetry on these parameters. We find that the average plume spacing is comparable between the spherical inner and outer bounding surfaces. An estimate of the average plume separation allows us to accurately predict the boundary layer asymmetry for the various spherical shell configurations explored here. The mean temperature and horizontal velocity profiles are in good agreement with classical Prandtl-Blasius laminar boundary layer profiles, provided the boundary layers are analysed in a dynamical frame that fluctuates with the local and instantaneous boundary layer thicknesses. The scaling properties of the Nusselt and Reynolds numbers are investigated by separating the bulk and boundary layer contributions to the thermal and viscous dissipation rates using numerical models with and a gravity proportional to . We show that our spherical models are consistent with the predictions of Grossmann & Lohse’s (2000) theory and that and scalings are in good agreement with plane layer results.

B

Turbulent Rayleigh-Bénard convection in spherical shells]Turbulent Rayleigh-Bénard convection in spherical shells T. Gastine, J. Wicht and J. M. Aurnou]Thomas Gastine1,\nsJohannes Wicht,\nsJonathan M. Aurnou 2015 \volume650 \pagerange119–126

énard convection, boundary layers, geophysical and geological flows

## 1 Introduction

Thermal convection is ubiquitous in geophysical and astrophysical fluid dynamics and rules, for example, turbulent flows in the interiors of planets and stars. The so-called Rayleigh-Bénard (hereafter RB) convection is probably the simplest paradigm to study heat transport phenomena in these natural systems. In this configuration, convection is driven in a planar fluid layer cooled from above and heated from below (figure 1a). The fluid is confined between two rigid impenetrable walls maintained at constant temperatures. The key issue in RB convection is to understand the turbulent transport mechanisms of heat and momentum across the layer. In particular, how does the heat transport, characterised by the Nusselt number , and the flow amplitude, characterised by the Reynolds number , depend on the various control parameters of the system, namely the Rayleigh number , the Prandtl number and the cartesian aspect ratio ? In general, quantifies the fluid layer width over its height in classical planar or cylindrical RB cells. In spherical shells, we rather employ the ratio of the inner to the outer radius to characterise the geometry of the fluid layer.

Laboratory experiments of RB convection are classically performed in rectangular or in cylindrical tanks with planar upper and lower bounding surfaces where the temperature contrast is imposed (see figure 1b). In such a system, the global dynamics are strongly influenced by the flow properties in the thermal and kinematic boundary layers that form in the vicinity of the walls. The characterisation of the structure of these boundary layers is crucial for a better understanding of the transport processes. The marginal stability theory by Malkus (1954) is the earliest boundary layer model and relies on the assumption that the thermal boundary layers adapt their thicknesses to maintain a critical boundary layer Rayleigh number, which implies . Assuming that the boundary layers are sheared, Shraiman & Siggia (1990) later derived a theoretical model that yields scalings of the form and (see also Siggia, 1994). These asymptotic laws were generally consistent with most of the experimental results obtained in the 1990s up to . Within the typical experimental resolution of one percent, simple power laws of the form were found to provide an adequate representation with exponents ranging from 0.28 to 0.31, in relatively good agreement with the Shraiman & Siggia model (e.g. Castaing et al., 1989; Chavanne et al., 1997; Niemela et al., 2000). However, later high-precision experiments by Xu et al. (2000) revealed that the dependence of upon cannot be accurately described by such simple power laws. In particular, the local slope of the function has been found to increase slowly with . The effective exponent of roughly ranges from values close to near to when (e.g. Funfschilling et al., 2005; Cheng et al., 2015).

Grossmann & Lohse (2000, 2004) derived a competing theory capable of capturing this complex dynamics (hereafter GL). This scaling theory is built on the assumption of laminar boundary layers of Prandtl-Blasius (PB) type (Prandtl, 1905; Blasius, 1908). According to the GL theory, the flows are classified in four different regimes in the phase space according to the relative contribution of the bulk and boundary layer viscous and thermal dissipation rates. The theory predicts non-power-law behaviours for and in good agreement with the dependence and observed in recent experiments and numerical simulations of RB convection in planar or cylindrical geometry (see for recent reviews Ahlers et al., 2009; Chillà & Schumacher, 2012).

Benefiting from the interplay between experiments and direct numerical simulations (DNS), turbulent RB convection in planar and cylindrical cells has received a lot of interest in the past two decades. However, the actual geometry of several fundamental astrophysical and geophysical flows is essentially three-dimensional within concentric spherical upper and lower bounding surfaces under the influence of a radial buoyancy force that strongly depends on radius. The direct applicability of the results derived in the planar geometry to spherical shell convection is thus questionable.

As shown in figure 1(c), convection in spherical shells mainly differs from the traditional plane layer configuration because of the introduction of curvature and the absence of side walls. These specific features of thermal convection in spherical shells yield significant dynamical differences with plane layers. For instance, the heat flux conservation through spherical surfaces implies that the temperature gradient is larger at the lower boundary than at the upper one to compensate for the smaller area of the bottom surface. This yields a much larger temperature drop at the inner boundary than at the outer one. In addition, this pronounced asymmetry in the temperature profile is accompanied by a difference between the thicknesses of the inner and the outer thermal boundary layers. Following Malkus’s marginal stability arguments, Jarvis (1993) and Vangelov & Jarvis (1994) hypothesised that the thermal boundary layers in curvilinear geometries adjust their thickness to maintain the same critical boundary layer Rayleigh number at both boundaries. This criterion is however in poor agreement with the results from numerical models (e.g. Deschamps et al., 2010). The exact dependence of the boundary layer asymmetry on the radius ratio and the gravity distribution thus remains an open question in thermal convection in spherical shells (Bercovici et al., 1989; Jarvis et al., 1995; Sotin & Labrosse, 1999; Shahnas et al., 2008; O’Farrell et al., 2013). This open issue sheds some light on the possible dynamical influence of asymmetries between the hot and cold surfaces that originate due to both the boundary curvature and the radial dependence of buoyancy in spherical shells.

Ground-based laboratory experiments involving spherical geometry and a radial buoyancy forcing are limited by the fact that gravity is vertically downwards instead of radially inwards (Scanlan et al., 1970; Feldman & Colonius, 2013). A possible way to circumvent this limitation is to conduct experiments under microgravity to suppress the vertically downward buoyancy force. Such an experiment was realised by Hart et al. (1986) who designed a hemispherical shell that flew on board of the space shuttle Challenger in May 1985. The radial buoyancy force was modelled by imposing an electric field across the shell. The temperature dependence of the fluid’s dielectric properties then produced an effective radial gravity that decreases with the fifth power of the radius (i.e. ). More recently, a similar experiment named “GeoFlow” was run on the International Space Station, where much longer flight times are possible (Futterer et al., 2010, 2013). This later experiment was designed to mimic the physical conditions in the Earth mantle. It was therefore mainly dedicated to the observation of plume-like structures in a high Prandtl number regime () for . Unfortunately, this limitation to relatively small Rayleigh numbers makes the GeoFlow experiment quite restricted regarding asymptotic scaling behaviours in spherical shells.

To compensate the lack of laboratory experiments, three dimensional numerical models of convection in spherical shells have been developed since the 1980s (e.g. Zebib et al., 1980; Bercovici et al., 1989, 1992; Jarvis et al., 1995; Tilgner, 1996; Tilgner & Busse, 1997; King et al., 2010; Choblet, 2012). The vast majority of the numerical models of non-rotating convection in spherical shells has been developed with Earth’s mantle in mind. These models therefore assume an infinite Prandtl number and most of them further include a strong dependence of viscosity on temperature to mimic the complex rheology of the mantle. Several recent studies of isoviscous convection with infinite Prandtl number in spherical shells have nevertheless been dedicated to the analysis of the scaling properties of the Nusselt number. For instance, Deschamps et al. (2010) measured convective heat transfer in various radius ratios ranging from to and reported for , while Wolstencroft et al. (2009) computed numerical models with Earth’s mantle geometry () up to and found . These studies also checked the possible influence of internal heating and reported quite similar scalings.

Most of the numerical models of convection in spherical shells have thus focused on the very specific dynamical regime of the infinite Prandtl number. The most recent attempt to derive the scaling properties of and in non-rotating spherical shells with finite Prandtl numbers is the study of Tilgner (1996). He studied convection in self-graviting spherical shells (i.e. ) with spanning the range and . This study was thus restricted to low Rayleigh numbers, relatively close to the onset of convection, which prevents the derivation of asymptotic scalings for and in spherical shells.

The objectives of the present work are twofold: (i) to study the scaling properties of and in spherical shells with finite Prandtl number; (ii) to better characterise the inherent asymmetric boundary layers in thermal convection in spherical shells. We therefore conduct two systematic parameter studies of turbulent RB convection in spherical shells with by means of three dimensional DNS. In the first set of models, we vary both the radius ratio (from to ) and the radial gravity profile (considering ) in a moderate parameter regime (i.e. for the majority of the cases) to study the influence of these properties on the boundary layer asymmetry. We then consider a second set of models with and up to . These DNS are used to check the applicability of the GL theory to thermal convection in spherical shells. We therefore numerically test the different basic prerequisites of the GL theory: we first analyse the nature of the boundary layers before deriving the individual scaling properties for the different contributions to the viscous and thermal dissipation rates.

The paper is organised as follows. In § 2, we present the governing equations and the numerical models. We then focus on the asymmetry of the thermal boundary layers in § 3. In § 4, we analyse the nature of the boundary layers and show that the boundary layer profiles are in agreement with the Prandtl-Blasius theory (Prandtl, 1905; Blasius, 1908). In § 5, we investigate the scaling properties of the viscous and thermal dissipation rates before calculating the and scalings in § 6. We conclude with a summary of our findings in § 7.

## 2 Model formulation

### 2.1 Governing hydrodynamical equations

We consider RB convection of a Boussinesq fluid contained in a spherical shell of outer radius and inner radius . The boundaries are impermeable, no slip and at constant temperatures and . We adopt a dimensionless formulation using the shell gap as the reference lengthscale and the viscous dissipation time as the reference timescale. Temperature is given in units of , the imposed temperature contrast over the shell. Velocity and pressure are expressed in units of and , respectively. Gravity is non-dimensionalised using its reference value at the outer boundary . The dimensionless equations for the velocity , the pressure and the temperature are given by

 \boldmath∇⋅\boldmathu% \boldmathu=0, (1)
 ∂\boldmathu∂t+\boldmathu⋅\boldmath∇\boldmathu=−\boldmath∇p+RaPrgT\boldmather+\boldmathΔ\boldmathu, (2)
 ∂T∂t+\boldmathu⋅\boldmath∇T=1PrΔT, (3)

where is the unit vector in the radial direction and is the gravity. Several gravity profiles have been classically considered to model convection in spherical shells. For instance, self-graviting spherical shells with a constant density correspond to (e.g Tilgner, 1996), while RB convection models with infinite Prandtl number usually assume a constant gravity in the perspective of modelling Earth’s mantle (e.g. Bercovici et al., 1989). The assumption of a centrally-condensed mass has also been frequently assumed when modelling rotating convection (e.g. Gilman & Glatzmaier, 1981; Jones et al., 2011) and yields . Finally, the artificial central force field of the microgravity experiments takes effectively the form of (Hart et al., 1986; Feudel et al., 2011; Futterer et al., 2013). To explore the possible impact of these various radial distribution of buoyancy on RB convection in spherical shells, we consider different models with the four following gravity profiles: . Particular attention will be paid to the cases with , which is the only radial function compatible with an exact analysis of the dissipation rates (see below, § 2.3).

The dimensionless set of equations (1-3) is governed by the Rayleigh number , the Prandtl number and the radius ratio of the spherical shell defined by

 Ra=αgoΔTd3νκ,Pr=νκ,η=riro, (4)

where and are the viscous and thermal diffusivities and is the thermal expansivity.

### 2.2 Diagnostic parameters

To quantify the impact of the different control parameters on the transport of heat and momentum, we analyse several diagnostic properties. We adopt the following notations regarding different averaging procedures. Overbars correspond to a time average

 ¯¯¯f=1τ∫t0+τt0fdt,

where is the time averaging interval. Spatial averaging over the whole volume of the spherical shell are denoted by triangular brackets , while correspond to an average over a spherical surface:

 ⟨f⟩=1V∫Vf(r,θ,ϕ)dV;⟨f⟩s=14π∫π0∫2π0f(r,θ,ϕ)sinθdθdϕ,

where is the volume of the spherical shell, is the radius, the colatitude and the longitude.

The convective heat transport is characterised by the Nusselt number , the ratio of the total heat flux to the heat carried by conduction. In spherical shells with isothermal boundaries, the conductive temperature profile is the solution of

 ddr(r2dTcdr)=0,Tc(ri)=1,Tc(ro)=0,

which yields

 Tc(r)=η(1−η)21r−η1−η. (5)

For the sake of clarity, we adopt in the following the notation for the time-averaged and horizontally-averaged radial temperature profile:

 ϑ(r)=¯¯¯¯¯¯¯¯¯¯⟨T⟩s.

 Nu=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨urT⟩s−1Prdϑdr−1PrdTcdr=−ηdϑdr(r=ri)=−1ηdϑdr(r=ro). (6)

The typical rms flow velocity is given by the Reynolds number

 Re=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯√⟨u2⟩=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯√⟨u2r+u2θ+u2ϕ⟩, (7)

while the radial profile for the time and horizontally-averaged horizontal velocity is defined by

 Missing or unrecognized delimiter for \right (8)

### 2.3 Exact dissipation relationships in spherical shells

The mean buoyancy power averaged over the whole volume of a spherical shell is expressed by

 P=RaPr¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨gurT⟩=4πVRaPr∫rorigr2¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨urT⟩sdr,

Using the Nusselt number definition (6) and the conductive temperature profile (5) then leads to

 P=4πVRaPr2(∫rorigr2dϑdrdr−Nuη(1−η)2∫rorigdr).

The first term in the parentheses becomes identical to the imposed temperature drop for a gravity :

 ∫rorigr2dϑdrdr=r2o[ϑ(ro)−ϑ(ri)]=−r2o,

and thus yields an analytical relation between and . For any other gravity model, we have to consider the actual spherically-symmetric radial temperature profile . Christensen & Aubert (2006) solve this problem by approximating by the diffusive solution (5) and obtain an approximate relation between and . This motivates our particular focus on the cases which allows us to conduct an exact analysis of the dissipation rates and therefore check the applicability of the GL theory to convection in spherical shells.

Noting that for , one finally obtains the exact relation for the viscous dissipation rate :

 ϵU=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨(\boldmath∇×\boldmathu\boldmathu)2⟩=P=31+η+η2RaPr2(Nu−1). (9)

The thermal dissipation rate can be obtained by multiplying the temperature equation (3) by and integrate it over the whole volume of the spherical shell. This yields

 ϵT=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⟨(∇T)2⟩=3η1+η+η2Nu. (10)

These two exact relations (9-10) can be used to validate the spatial resolutions of the numerical models with . To do so, we introduce and , the ratios of the two sides of Eqs (9-10):

 χϵU =(1+η+η2)Pr23Ra(Nu−1)ϵU, (11) χϵT =(1+η+η2)3ηNuϵT.

### 2.4 Setting up a parameter study

#### Numerical technique

The numerical simulations have been carried out with the magnetohydrodynamics code MagIC (Wicht, 2002). MagIC has been validated via several benchmark tests for convection and dynamo action (Christensen et al., 2001; Jones et al., 2011). To solve the system of equations (1-3), the solenoidal velocity field is decomposed into a poloidal and a toroidal contribution

 \boldmathu\boldmathu=\boldmath∇×(\boldmath∇×W\boldmather)+\boldmath∇×Z%\boldmath$er$,

where and are the poloidal and toroidal potentials. , , and are then expanded in spherical harmonic functions up to degree in the angular variables and and in Chebyshev polynomials up to degree in the radial direction. The combined equations governing and are obtained by taking the radial component and the horizontal part of the divergence of (2). The equation for is obtained by taking the radial component of the curl of (2). The equations are time-stepped by advancing the nonlinear terms using an explicit second-order Adams-Bashforth scheme, while the linear terms are time-advanced using an implicit Crank-Nicolson algorithm. At each time step, all the nonlinear products are calculated in the physical space (, , ) and transformed back into the spectral space (, , ). For more detailed descriptions of the numerical method and the associated spectral transforms, the reader is referred to (Gilman & Glatzmaier, 1981; Tilgner & Busse, 1997; Christensen & Wicht, 2007).

#### Parameter choices

One of the main focuses of this study is to investigate the global scaling properties of RB convection in spherical shell geometries. This is achieved via measurements of the Nusselt and Reynolds numbers. In particular, we aim to test the applicability of the GL theory to spherical shells. As demonstrated before, only the particular choice of a gravity profile of the form allows the exactness of the relation (9). Our main set of simulations is thus built assuming . The radius ratio is kept to and the Prandtl number to to allow future comparisons with the rotating convection models by Gastine & Wicht (2012) and Gastine et al. (2013) who adopted the same configuration. We consider 35 numerical cases spanning the range . Table 1 summarises the main diagnostic quantities for this dataset of numerical simulations and shows that our models basically lie within the ranges and .

Another important issue in convection in spherical shells concerns the determination of the average bulk temperature and the possible boundary layer asymmetry between the inner and the outer boundaries (e.g. Jarvis, 1993; Tilgner, 1996). To better understand the influence of curvature and the radial distribution of buoyancy, we thus compute a second set of numerical models. This additional dataset consists of 113 additional simulations with various radius ratios and gravity profiles, spanning the range with . To limit the numerical cost of this second dataset, these cases have been run at moderate Rayleigh number and typically span the range for the majority of the cases. Table 2, given in the Appendix, summarises the main diagnostic quantities for this second dataset of numerical simulations.

#### Resolution checks

Attention must be paid to the numerical resolutions of the DNS of RB convection (e.g. Shishkina et al., 2010). Especially, underresolving the fine structure of the turbulent flow leads to an overestimate of the Nusselt number, which then falsifies the possible scaling analysis (Amati et al., 2005). One of the most reliable ways to validate the truncations employed in our numerical models consists of comparing the obtained viscous and thermal dissipation rates with the average Nusselt number (Stevens et al., 2010; Lakkaraju et al., 2012; King et al., 2012). The ratios and , defined in (11), are found to be very close to unity for all the cases of Table 1, which supports the adequacy of the employed numerical resolutions. To further highlight the possible impact of inadequate spatial resolutions, two underresolved numerical models for the two highest Rayleigh numbers have also been included in Table 1 (lines in italics). Because of the insufficient number of grid points in the boundary layers, the viscous dissipation rates are significantly higher than expected in the statistically stationary state. This leads to overestimated Nusselt numbers by similar percent differences ().

Table 1 shows that the typical resolutions span the range from () to (). The two highest Rayleigh numbers have been computed assuming a two-fold azimuthal symmetry to ease the numerical computations. A comparison of test runs with or without the two-fold azimuthal symmetry at lower Rayleigh numbers () showed no significant statistical differences. This enforced symmetry is thus not considered to be influential. The total computational time for the two datasets of numerical models represents roughly 5 million Intel Ivy Bridge CPU hours.

## 3 Asymmetric boundary layers in spherical shells

### 3.1 Definitions

Several different approaches are traditionally considered to define the thermal boundary layer thickness . They either rely on the horizontally-averaged mean radial temperature profile or on the temperature fluctuation defined as

 Missing or unrecognized delimiter for \right (12)

Among the possible estimates based on , the slope method (e.g. Verzicco & Camussi, 1999; Breuer et al., 2004; Liu & Ecke, 2011) defines as the depth where the linear fit to near the boundaries intersects the linear fit to the temperature profile at mid-depth. Alternatively, exhibits sharp local maxima close to the walls. The radial distance separating those peaks from the corresponding nearest boundary can be used to define the thermal boundary layer thicknesses (e.g. Tilgner, 1996; King et al., 2013). Figure 2(a) shows that both definitions of actually yield nearly indistinguishable boundary layer thicknesses. We therefore adopt the slope method to define the thermal boundary layers.

There are also several ways to define the viscous boundary layers. Figure 2(b) shows the vertical profile of the root-mean-square horizontal velocity . This profile exhibits strong increases close to the boundaries that are accompanied by well-defined peaks. Following Tilgner (1996) and Kerr & Herring (2000), the first way to define the kinematic boundary layer is thus to measure the distance between the walls and these local maxima. This commonly-used definition gives () for the inner (outer) spherical boundary. Another possible method to estimate the viscous boundary layer follows a similar strategy as the slope method that we adopted for the thermal boundary layers (Breuer et al., 2004). () is defined as the distance from the inner (outer) wall where the linear fit to near the inner (outer) boundary intersects the horizontal line passing through the maximum horizontal velocity.

Figure 2(b) reveals that these two definitions lead to very distinct viscous boundary layer thicknesses. In particular, the definition based on the position of the local maxima of yields much thicker boundary layers than the tangent intersection method, i.e. . The discrepancies of these two definitions are further discussed in § 4.

### 3.2 Asymmetric thermal boundary layers and mean bulk temperature

Figure 2 also reveals a pronounced asymmetry in the mean temperature profiles with a much larger temperature drop at the inner boundary than at the outer boundary. As a consequence, the mean temperature of the spherical shell is much below . Determining how the mean temperature depends on the radius ratio has been an ongoing open question in mantle convection studies with infinite Prandtl number (e.g. Bercovici et al., 1989; Jarvis, 1993; Vangelov & Jarvis, 1994; Jarvis et al., 1995; Sotin & Labrosse, 1999; Shahnas et al., 2008; Deschamps et al., 2010; O’Farrell et al., 2013). To analyse this issue in numerical models with , we have performed a systematic parameter study varying both the radius ratio of the spherical shell and the gravity profile (see Table 2). Figure 3 shows some selected radial profiles of the mean temperature for various radius ratios (panel a) and gravity profiles (panel b) for cases with similar . For small values of , the large difference between the inner and the outer surfaces lead to a strong asymmetry in the temperature distribution: nearly 90% of the total temperature drop occurs at the inner boundary when . In thinner spherical shells, the mean temperature gradually approaches a more symmetric distribution to finally reach when (no curvature). Figure 3(b) also reveals that a change in the gravity profile has a direct impact on the mean temperature profile. This shows that both the shell geometry and the radial distribution of buoyancy affect the temperature of the fluid bulk in RB convection in spherical shells.

To analytically access the asymmetries in thickness and temperature drop observed in figure 3, we first assume that the heat is purely transported by conduction in the thin thermal boundary layers. The heat flux conservation through spherical surfaces (6) then yields

 ΔToλoT=η2ΔTiλiT, (13)

where the thermal boundary layers are assumed to correspond to a linear conduction profile with a temperature drop () over a thickness (). As shown in Figs. 2-3, the fluid bulk is isothermal and forms the majority of the fluid by volume. We can thus further assume that the temperature drops occur only in the thin boundary layers, which leads to

 ΔTo+ΔTi=1. (14)

Equations (13) and (14) are nevertheless not sufficient to determine the three unknowns , , and an additional physical assumption is required.

A hypothesis frequently used in mantle convection models with infinite Prandtl number in spherical geometry (Jarvis, 1993; Vangelov & Jarvis, 1994) is to further assume that both thermal boundary layers are marginally stable such that the local boundary layer Rayleigh numbers and are equal:

 Raiλ=Raoλ→αgiΔTiλiT3νκ=αgoΔToλoT3νκ. (15)

This means that both thermal boundary layers adjust their thickness and temperature drop to yield (e.g., Malkus, 1954). The temperature drops at both boundaries and the ratio of the thermal boundary layer thicknesses can then be derived using Eqs. (13-14)

 ΔTi=11+η3/2χ1/4g,ΔTo≃Tm=η3/2χ1/4g1+η3/2χ1/4g,λoTλiT=χ1/4gη1/2, (16)

where

 χg=g(ri)g(ro), (17)

is the ratio of the gravitational acceleration between the inner and the outer boundaries. Figure 4(a) reveals that the marginal stability hypothesis is not fulfilled when different radius ratios and gravity profiles are considered. This is particularly obvious for small radius ratios where is more than 10 times larger than . This discrepancy tends to vanish when , when curvature and gravity variations become unimportant. As a consequence, there is a significant mismatch between the predicted mean bulk temperature from (16) and the actual values (figure 4b). Deschamps et al. (2010) also reported a similar deviation from (16) in their spherical shell models with infinite Prandtl number. They suggest that assuming instead might help to improve the agreement with the data. This however cannot account for the additional dependence on the gravity profile visible in figure 4. We finally note that for the database of numerical simulations explored here, which suggests that the thermal boundary layers are stable in all our simulations.

Alternatively Wu & Libchaber (1991) and Zhang et al. (1997) proposed that the thermal boundary layers adapt their thicknesses such that the mean hot and cold temperature fluctuations at mid-depth are equal. Their experiments with Helium indeed revealed that the statistical distribution of the temperature at mid-depth was symmetrical. They further assumed that the thermal fluctuations in the center can be identified with the boundary layer temperature scales and , which characterise the temperature scale of the thermal boundary layers in a different way than the relative temperature drops and . This second hypothesis yields

 θi=θo→νκαgiλiT3=νκαgoλoT3, (18)

and the corresponding temperature drops and boundary layer thicknesses ratio

 ΔTi=11+η2χ1/3g,ΔTo=Tm=η2χ1/3g1+η2χ1/3g,λoTλiT=χ1/3g. (19)

Figure 5(a) shows for different radius ratios and gravity profiles, while figure 5(b) shows a comparison between the predicted mean bulk temperature and the actual values. Besides the cases with which are in relatively good agreement with the predicted scalings, the identity of the boundary layer temperature scales is in general not fulfilled for the other gravity profiles. The actual mean bulk temperature is thus poorly described by (19). We note that previous findings by Ahlers et al. (2006) already reported that the theory by Wu & Libchaber’s does also not hold when the transport properties depend on temperature (i.e. non-Oberbeck-Boussinesq convection).

### 3.3 Conservation of the average plume density in spherical shells

As demonstrated in the previous section, none of the hypotheses classically employed accurately account for the temperature drops and the boundary layer asymmetry observed in spherical shells. We must therefore find a dynamical quantity that could be possibly identified between the two boundary layers.

Figure 6 shows visualisations of the thermal boundary layers for three selected numerical models with different radius ratios and gravity profiles. The isocontours displayed in panels (a-c) reveal the intricate plume structure. Long and thin sheet-like structures form the main network of plumes. During their migration along the spherical surfaces, these sheet-like plumes can collide and convolute with each other to give rise to mushroom-type plumes (see Zhou & Xia, 2010b; Chillà & Schumacher, 2012). During this morphological evolution, mushroom-type plumes acquire a strong radial vorticity component. These mushroom-type plumes are particularly visible at the connection points of the sheet plumes network at the inner thermal boundary layer (red isosurface in figure 6a-c). Figure 6(d-f) shows the corresponding equatorial and radial cuts of the temperature fluctuation . These panels further highlight the plume asymmetry between the inner and the outer thermal boundary layers. For instance, the case with and (top panels) features an outer boundary layer approximately 4.5 times thicker than the inner one. Accordingly, the mushroom-like plumes that depart from the outer boundary layer are significantly thicker than the ones emitted from the inner boundary. This discrepancy tends to vanish in the thin shell case (, bottom panels) in which curvature and gravity variations play a less significant role ( in that case).

Puthenveettil & Arakeri (2005) and Zhou & Xia (2010b) performed statistical analysis of the geometrical properties of thermal plumes in experimental RB convection. By tracking a large number of plumes, their analysis revealed that both the plume separation and the width of the sheet-like plumes follow a log-normal probability density function (PDF).

To further assess how the average plume properties of the inner and outer thermal boundary layers compare with each other in spherical geometry, we adopt a simpler strategy by only focussing on the statistics of the plume density. The plume density per surface unit at a given radius is expressed by

 ρp∼N4πr2, (20)

where is the number of plumes, approximated here by the ratio of the spherical surface area to the mean inter-plume area :

 N∼4πr2¯S. (21)

This inter-plume area can be further related to the average plume separation via .

An accurate evaluation of the inter-plume area for each thermal boundary layer however requires to separate the plumes from the background fluid. Most of the criteria employed to determine the location of the plume boundaries are based on thresholds of certain physical quantities (see Shishkina & Wagner, 2008, for a review of the different plume extraction techniques). This encompasses threshold values on the temperature fluctuations (Zhou & Xia, 2002), on the vertical velocity (Ching et al., 2004) or on the thermal dissipation rate (Shishkina & Wagner, 2005). The choice of the threshold value however remains an open problem. Alternatively, Vipin & Puthenveettil (2013) show that the sign of the horizontal divergence of the velocity might provide a simple and threshold-free criterion to separate the plumes from the background fluid

 \boldmath∇H⋅\boldmathu=1rsinθ∂∂θ(sinθuθ)+1rsinθ∂uϕ∂ϕ=−1r2∂∂r(r2ur).

Fluid regions with indeed correspond to local regions of positive vertical acceleration, expected inside the plumes, while the fluid regions with characterise the inter-plume area.

To analyse the statistics of , we thus consider here several criteria based either on a threshold value of the temperature fluctuations or on the sign of the horizontal divergence. This means that a given inter-plume area at the inner (outer) thermal boundary layer is either defined as an enclosed region surrounded by hot (cold) sheet-like plumes carrying a temperature perturbation ; or by an enclosed region with . To further estimate the possible impact of the chosen threshold value on , we vary between and . This yields

 S(r)≡r2∮Tsinθdθdϕ, (22)

where the physical criterion () to extract the plume boundaries at the inner (outer) boundary layer is given by

 Ti={T′(riλ,θ,ϕ)≤t,t∈[σ(riλ),σ(riλ/2),σ(riλ/4)],\boldmath∇H⋅\boldmathu≥0, (23) To={T′(roλ,θ,ϕ)≥t,t∈[σ(riλ),σ(riλ/2),σ(riλ/4)],\boldmath∇H⋅\boldmathu≥0,

where () for the inner (outer) thermal boundary layer.

Figure 7 shows an example of such a characterisation procedure for the inner thermal boundary layer of a numerical model with , , . Panel (b) illustrates a plume extraction process when using to determine the location of the plumes: the black area correspond to the inter-plume spacing while the white area correspond to the complementary plume network location. The fainter emerging sheet-like plumes are filtered out and only the remaining “skeleton” of the plume network is selected by this extraction process. The choice of is however arbitrary and can influence the evaluation of the number of plumes. The insets displayed in panels (c-e) illustrate the sensitivity of the plume extraction process on the criterion employed to detect the plumes. In particular, using the threshold based on the largest temperature fluctuations can lead to the fragmentation of the detected plume lanes into several isolated smaller regions. As a consequence, several neighbouring inter-plume areas can possibly be artificially connected when using this criterion. In contrast, using the sign of the horizontal divergence to estimate the plumes location yields much broader sheet-like plumes. As visible on panel (e), the plume boundaries frequently correspond to local maxima of the thermal dissipation rate (Shishkina & Wagner, 2008).

For each criterion given in (23), we then calculate the area of each bounded black surface visible in figure 7(b) to construct the statistical distribution of the inter-plume area for both thermal boundary layers. Figure 8 compares the resulting PDFs obtained by combining several snapshots for a numerical model with , and . Besides the criterion which yields PDFs that are slightly shifted towards smaller inter-plume spacing areas, the statistical distributions are found to be relatively insensitive to the detection criterion (23). We therefore restrict the following comparison to the criterion only.

Figure 9 shows the PDFs for the three numerical models of figure 6. For the two cases with and (panels b-c), the statistical distributions for both thermal boundary layers nearly overlap. This means that the inter-plume area is similar at both spherical shell surfaces. In contrast, for the case with (panel a), the two PDFs are offset relative to each other. However, the peaks of the distributions remain relatively close, meaning that once again the inner and the outer thermal boundary layers share a similar average inter-plume area. Puthenveettil & Arakeri (2005) and Zhou & Xia (2010b) demonstrated that the thermal plume statistics in turbulent RB convection follow a log-normal distribution (see also Shishkina & Wagner, 2008; Puthenveettil et al., 2011). The large number of plumes in the cases with and (figure 6b-c) would allow a characterisation of the nature of the statistical distributions. However, this would be much more difficult in the case (figure 6a) in which the plume density is significantly weaker. As a consequence, no further attempt has been made to characterise the exact nature of the PDFs visible in figure 9, although the universality of the log-normal statistics reported by Puthenveettil & Arakeri (2005) and Zhou & Xia (2010b) likely indicates that the same statistical distribution should hold here too.

The inter-plume area statistics therefore reveals that the inner and the outer thermal boundary layers exhibit a similar average plume density, independently of the spherical shell geometry and the gravity profile. Assuming would allow us to close the system of equations (13-14) and thus finally estimate , and . This however requires us to determine an analytical expression of the average inter-plume area or equivalently of the mean plume separation that depends on the boundary layer thickness and the temperature drop.

Using the boundary layer equations for natural convection (Rotem & Claassen, 1969), Puthenveettil et al. (2011) demonstrated that the thermal boundary layer thickness follows

 λi,oT(x)∼x(Rai,ox)1/5, (24)

where is the distance along the horizontal direction and is a Rayleigh number based on the lengthscale and on the boundary layer temperature jumps . As shown on figure 10, using (Puthenveettil & Arakeri, 2005; Puthenveettil et al., 2011) then allows to establish the following relation for the average plume spacing

 λT¯ℓ∼1Ra1/5ℓ. (25)

which yields

 ¯ℓi∼√αgiΔTiλiT5νκ,¯ℓo∼√αgoΔToλoT5νκ, (26)

for both thermal boundary layers. We note that an equivalent expression for the average plume spacing can be derived from a simple mechanical description of the equilibrium between production and coalescence of plumes in each boundary layer (see Parmentier & Sotin, 2000; King et al., 2013).

Equation (26) is however expected to be only valid at the scaling level. The vertical lines in figure 9 therefore correspond to the estimated average inter-plume area for both thermal boundary layers using (26) and . The predicted average inter-plume area is in good agreement with the peaks of the statistical distributions for the three cases discussed here. The expression (26) therefore provides a reasonable estimate of the average plume separation (Puthenveettil & Arakeri, 2005; Puthenveettil et al., 2011; Gunasegarane & Puthenveettil, 2014). The comparable observed plume density at both thermal boundary layers thus yields

 ρip=ρop→αgiΔTiλiT5νκ=αgoΔToλoT5νκ. (27)

Using Eqs. (13-14) then allows us to finally estimate the temperature jumps and the ratio of the thermal boundary layer thicknesses in our dimensionless units:

 ΔTi=11+η5/3χ1/6g,ΔTo=Tm=η5/3χ1/6g1+η5/3χ1/6g,λoTλiT=χ1/6gη1/3. (28)

Figure 11 shows the ratios , and the temperature jumps and . In contrast to the previous criteria, either coming from the marginal stability of the boundary layer (16, figure 4) or from the identity of the temperature fluctuations at mid-shell (28, figure 5), the ratio of the average plume separation now falls much closer to the unity line. Some deviations are nevertheless still visible for spherical shells with and (orange circles). The comparable average plume density between both boundary layers allows us to accurately predict the asymmetry of the thermal boundary layers and the corresponding temperature drops for the vast majority of the numerical cases explored here (solid lines in panels b-d).

As we consider a fluid with , the viscous boundary layers should show a comparable degree of asymmetry as the thermal boundary layers. (28) thus implies

 λoUλiU=λoTλiT=χ1/6gη1/3. (29)

Figure 12 shows the ratio of the viscous boundary layer thicknesses for the different setups explored in this study. The observed asymmetry between the two spherical shell surfaces is in a good agreement with (29) (solid black lines).

### 3.4 Thermal boundary layer scalings

Using (28) and the definition of the Nusselt number (6), we can derive the following scaling relations for the thermal boundary layer thicknesses:

 λiT=η1+η5/3χ1/6g1Nu,λoT=η2/3χ1/6g1+η5/3χ1/6g1Nu. (30)

Figure 13(a) demonstrates that the boundary layer thicknesses for the numerical simulations of Table 1 ( and ) are indeed in close agreement with the theoretical predictions. To further check this scaling for other spherical shell configurations, we introduce the following normalisation of the thermal boundary layer thicknesses