Project topics and materials
Get Complete Project Material File(s) Now! »
Introduction
Thermal convection is the main heat transfer mechanism in the interiors of planets and many stars (Phillips, 2013) 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, 1989; Busse and Carrigan, 1976). It is also one of the most studied example of pattern–forming system (Cross and Hohenberg, 1993).
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 (Chandrasekhar, 1961; Rayleigh, 1916) 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 Bénard (1900). 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 (Busse and Riahi, 1980). 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 Dirichlet and Neumann boundary conditions, respectively. In-deed, in some physical situations, the boundary conditions is more accurately described by a thermal flux proportional to the temperature itself, leading to β T + ∂T = 0 (2.1) ∂r with β a proportionality factor homogeneous to an inverse length-scale. Using a reference length-scale L, one can construct the non-dimensional Biot number Bi = βL. (2.2)
The boundary condition (2.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 (Busse and Riahi, 1980) or conductive-radiative empirical Newton’s law (O’Sullivan, 1990). Note that linearised black body radiative equilibrium are often interpreted as Newton’s law (for example in meteorological models, Savijärvi (1994)). Thermal contact between two layers of distinct length and thermal conductivity can also be modelled by a Robin boundary condition as done by Guillou and Jaupart (1995).
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 (Labrosse et al., 2015). Following the current scientific consensus, Earth may have undergone giant im-pacts during its formation (Solomatov, 2007). Such impacts were able to melt a large part of the nowadays solid mantle, producing what is called a magma ocean. This phys-ical object is not known with extreme precision but was characterised, at its formation, by large depth (∼ 1 × 106 m), low viscosity (∼ 0.1 Pa s) and rapid flows (∼ 10 m s−1) (Solomatov, 2007). 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 interac-tion 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 Sparrow et al. (Sparrow et al., 1964), 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, Chapman and Proctor (1980) 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 Busse and Riahi (1980) who showed that square-pattern convection planforms were preferred to two-dimensional rolls produced in the case of FT boundary conditions. Later, Ishiwatari et al. (1994) explored, by numerical integration in a Cartesian two-dimensional 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 (Chandrasekhar, 1961).
Adding the influence of global rotation, Takehiro et al. (2002) 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, Falsaperla et al. (2010) conducted studies on a Cartesian system with Robin boundary conditions influenced by weak rotation and also observed a monotonu-ous transition from FF to FT when varying the Biot number.
Influence of Robin boundary conditions has been mainly studied in Cartesian geom-etry. 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 η = ri/re = 2/3 (all the parameters and fields used are listed in the appendix A.2.1), with re and ri the external and internal radii, respectively, subjected to rotation of vector Ωez and filled with a Newtonian fluid of constant kinematic viscosity ν, thermal conductivity k, reference density ρ0, specific heat capacity cp and thermal expansivity α.
The system is subjected to a radial gravity field g = −ger, g being uniform. Note that this uniform value is relevant for thin shells (if the shell is thin enough, any variation of gravity inside can be neglected) but also generally speaking for Earth mantle studies. Indeed, due to the density distribution within the Earth, gravity at the bottom of the mantle is about 10.7 m.s−2 with variations below 10% up to the surface (Dziewonski and Anderson, 1981). Since our system of interest is the liquid equivalent of the current Earth mantle we have chosen to use a constant gravity. Nevertheless other mass and density distributions in other differentiated rocky bodies (Moon, Mars, Mercury …) may lead to different gravity profile in their respective mantle. Consequently, other radial dependencies of the gravity field can be used in other contexts (e.g. g ∝ r for uniform density models, or g ∝ r12 for centrally condensed mass).
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 ρ = 1 − α(T0 − T0), (2.3) ρ0
ρ being the density, while ρ0 and T0 are reference values of density and temperature and T 0 the temperature anomaly around T0.
The reference state is motionless (uniform zero velocity) and follows a radial con-ductive temperature profile Tref between the two boundaries with given temperatures (FT) ri re + T0,
Tref (r) = ΔT − 1 (2.4) re − ri r where ΔT is the total temperature difference of the reference profile across the shell. In the Boussinesq approximation, the value of T0 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 re as length scale, ν as velocity scale, re Δ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 = νρ0cp/k, the Rayleigh number Ra = gαΔT r3 ρ c p and the Ekman number Ek = ν .
The dimensionless conservation equations are then ∂Θ Tref 1 + u • r Θ + = r2Θ (2.5) ∂t ΔT Pr for the energy, Ek ∂u • ru + 2ez × u + rΠ = Ekr2u + (RaEk/Pr)Θer, + u (2.6) ∂t for the momentum and r • u = 0, (2.7) 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 (ur = 0), stress-free (i.e. tangential stress is null, that is with [σ] the tensor of viscous stresses of the fluid, t the tangential and n normal vector to the interface) ([σ] • n) • t = 0 at the top (r = re) and no-slip (u = 0) at the bottom (r = ri). For the thermal boundary conditions, we impose a fixed temperature at the bottom (FT, i.e. Θ(ri) = 0) and a Robin boundary condition at the top (i.e. −βreΘ(re) = ∂∂rΘ (re)).
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 (re) and surface heat flux −k∂rT (re), and an isothermal atmosphere at temperature Ta leads to an energy balance, σ(Ta4 − T 4(re)) = k∂rT (re), (2.8) 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 Tref such that its value at the surface is set by σ(Ta4 − Tref4 (re)) = k ∂Tref (re), (2.9) ∂r and the depth-dependent temperature anomaly θ = T − Tref . Assuming θ T0, equation (2.8) can be expanded to give − 4σTref3 (re)θ = k∂rθ(re). (2.10)
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 r0 can be expressed as Bi Y + ∂lY = 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 ex-changes can be modelled by a thermal Newton’s law with h a coefficient of the order of 10 W m−2 K−1 (experimentally determined). 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. Kreith (2000)).
In the case described in eq. (2.10), h = 4σTa3 so that 4σT3L Bi = a . (2.11)
Table 2.1 – Dimensionless numbers characterizing a terrestrial global magma ocean, maximal and minimal values are given, computed from Solomatov (2007), Maas and Hansen (2015) and Snyder et al. (1994) (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.
Here we choose to use the thickness of the shell L = re − ri as the typical length scale to compute and express the Biot number.
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 dimen-sionless 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. 2.1.
These parameters constitute extreme conditions for numerical simulations (for ex-ample in terms of Ek) which cannot be considered with current computational resources. 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. 2.1).
In addition to radiative equilibrium, Robin boundary condition can also model ther-mal contact between two layers of distinct conductivity (λ1 and λ2) and thickness (L1 and L2). In that case, following Guillou and Jaupart (1995), a Robin boundary condi-tion can be written at the frontier. It links linearly dimensionless temperature and flux by a factor Bi = λ2L1 . Such a model could be used for internal oceans of icy satellites of Jupiter (Sotin and Tobie, 2004) or Saturn (Collins and Goodman, 2007) with a layer of liquid water lies below an ice crust. In that case, the ratio of the conductivities would be typically around 3 (Kreith (2000)). 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; Monville et al., 2019; Vidal and Schaeffer, 2015) in order to determine the critical Rayleigh number of the system, RaC , and, in terms of spherical harmonics, the order of azimuthal 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 in-compressible and stratified fluids in spherical cavities. SINGE uses the SLEPc solver (Vidal and Schaeffer, 2015), finite difference (2nd order scheme with Nr points) in the radial direction and spherical harmonic decomposition, with degree l and azimuthal wave number m truncated at (lmax, mmax). 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 Kaplan et al. (2017) or Guervilly et al. (2019)). 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 (Schaeffer, 2013). 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 obtain a 3D-distribution of the anomalies of the main fields (temperature, velocity etc., see fig. 2.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 (with respect to the equatorial plane) are always the preferred unstable modes at the linear onset of convection (i.e.their RaC is the smallest).
In fig. 2.2 we plot RaC 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. Gastine et al. (2016)): in the first one (Ek → +∞) rotation does not affect RaC while this number increases when Ek decreases in the second regime. For small enough Ek, the expected power law RaC ∝ Ek is observed (Chandrasekhar, 1953). Note that, although RaC 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 RaC as a function of Ek is not affected by the value of Pr but its influence on the exact value of RaC 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 non-rotating Rayleigh-Bénard convection for which the onset of convection is known to be independent of the value of the Prandtl number (e.g. Chandrasekhar, 1961). 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. 2.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 Falsaperla et al. (2010) 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.
Table of contents :
1 Introduction
1.1 Une brève histoire du concept d’océan de magma
1.1.1 Premiers modèles physiques de la naissance du système solaire
1.1.2 Le refroidissement de la Terre selon Buffon
1.1.3 Une Terre primitive liquide ?
1.1.4 L’apport des roches lunaires
1.1.5 Impact géant et océan de magma lunaire
1.1.6 Quelques modélisations géophysiques de l’océan de magma terrestre
1.1.7 Des témoignages géochimiques plus ambigus
1.1.8 Bilan
1.2 Enjeux, motivations et méthodes
1.2.1 Des océans de magma bien mélangés : le modèle standard
1.2.2 Impact de la cristallisation fractionnée et motivation de notre contribution
1.2.3 Le choix d’une approche hydrodynamique
1.2.4 L’adimensionnement du problème, une démarche préalable indispensable
1.2.5 Bilan
2 Influence of Robin Boundary Condition on Thermal Rotating Convection
2.1 Résumé
2.2 Abstract
2.3 Erratum
2.4 Introduction
2.5 Methods
2.5.1 Physical System
2.5.2 Numerical Tools
2.6 Linear Stability Analysis
2.6.1 Differences between FF and FT cases with respect to Ekman number
2.6.2 Transition from FF to FT owing to Robin boundary condition
2.6.3 Possible dependence on the Prandtl number
2.7 Non Linear Simulations
2.7.1 Behaviour of purely FT or FF setups
2.7.2 FF to FT transition
2.8 Discussion
3 Thermo-Solutal Convection in Magma Oceans
3.1 Résumé
3.2 Abstract
3.3 Introduction
3.3.1 The Concept of Magma Ocean
3.3.2 Terrestrial Case
3.3.3 Thermal Convection & Partial Crystallization
3.3.4 Thermo-solutal Convection
3.4 Methods
3.4.1 Physical System
3.4.2 Reference Fields
3.4.3 Equations of Conservation
3.4.4 Boundary Conditions
3.4.5 Parameters
3.4.6 Numerical Tools & Diagnoses
3.5 Heat and Mass Transport Regimes Observed
3.5.1 Fully Convective Regime
3.5.2 Erosion Regime
3.5.3 Enduring Stratification Regime
3.6 Regimes Diagram
3.6.1 Dependence with other parameters
3.6.2 Dependence on initial conditions
3.7 Discussion
4 Conclusion
4.1 Résumé
4.2 Perspectives
Références
