Weakly non-Boussinesq convection in a gaseous spherical shell

TL;DR

This paper investigates weakly compressible convection in a spherical shell with asymmetric boundary conditions, revealing pressure-dominated flows, anomalous heat transport, and the development of subadiabatic layers despite vigorous convection.

Contribution

It demonstrates how asymmetry in boundary conditions causes a transition from buoyancy- to pressure-dominated convection and reveals non-local effects in spherical shell convection.

Findings

Flow becomes pressure-dominated under asymmetric conditions

Anomalous heat transport occurs via upflows

Subadiabatic layers can develop with vigorous convection

Abstract

We examine the dynamics associated with weakly compressible convection in a spherical shell by running 3D direct numerical simulations using the Boussinesq formalism [1]. Motivated by problems in astrophysics, we assume the existence of a finite adiabatic temperature gradient $\nabla T_{\rm{ad}}$ and use mixed boundary conditions for the temperature with fixed flux at the inner boundary and fixed temperature at the outer boundary. This setup is intrinsically more asymmetric than the more standard case of Rayleigh-B\'{e}nard convection in liquids between parallel plates with fixed temperature boundary conditions. Conditions where there is substantial asymmetry can cause a dramatic change in the nature of convection and we demonstrate that this is the case here. The flows can become pressure- rather than buoyancy- dominated leading to anomalous heat transport by upflows. Counter-intuitively, the background temperature gradient $\nabla\bar{T}$ can develop a subadiabatic layer (where $\boldsymbol{g}\cdot\nabla\bar{T}<\boldsymbol{g}\cdot\nabla T_{\rm{ad}}$, where $\boldsymbol{g}$ is gravity) although convection remains vigorous at every point across the shell. This indicates a high degree of non-locality.

Figures

Tables

Table 1. Table 1: Table with all the different model configurations and the input parameters used in each case. The resolution is provided both in number of equivalent meshpoints N r subscript 𝑁 𝑟 N_{r} , N θ subscript 𝑁 𝜃 N_{\theta} , N ϕ subscript 𝑁 italic-ϕ N_{\phi} , as well as in the number of spherical harmonics used in the horizontal directions, L m ​ a ​ x subscript 𝐿 𝑚 𝑎 𝑥 L_{max} and M m ​ a ​ x subscript 𝑀 𝑚 𝑎 𝑥 M_{max} .

	Model	$\chi$	Rao	$N_r$	$N_{\theta }$	$N_{\varphi }$	$L_{max}$	$M_{max}$
	A	0.1	${10}^7$	250	402	480	268	134
	A	0.5	${10}^7$	220	346	384	230	120
	A	1	${10}^7$	220	346	384	230	120
	B	0.001	${10}^7$	200	288	320	192	96
	B	0.01	${10}^6$	200	192	192	128	64
	B	0.01	${10}^7$	200	288	320	192	96
	B	0.01	${10}^8$	300	516	640	344	172
	B	0.1	${10}^7$	200	288	320	192	96
	B	0.1	${10}^8$	300	516	640	344	172
	B	0.5	${10}^7$	200	288	320	192	96
	C	0.01	${10}^7$	200	288	320	192	96
	C	0.1	${10}^7$	200	288	320	192	96

Equations

$$\kappa {\nabla }^2{T}_{\mathrm{rad}}(r)=0\Rightarrow \kappa r^2\frac{d{T}_{\mathrm{rad}}}{dr}=\mathrm{const},$$

$$-\kappa\frac{dT_{\rm rad}}{dr}\bigg{|}_{r=r_{i}}=F_{\rm{rad}},$$

$$\frac{d{T}_{\mathrm{rad}}}{dr}=-\frac{F_{\mathrm{rad}}}{\kappa }{\left(\frac{r_i}{r}\right)}^2,$$

$$T_{\mathrm{rad}}(r)=\frac{F_{\mathrm{rad}}{r}_i^2}{\kappa }\left(\frac{1}{r}-\frac{1}{r_o}\right)+T_o.$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}=-\frac{1}{{\rho }_m}\nabla p+\alpha \Theta g{\mathbf{e}}_{\mathbf{r}}+\nu {\nabla }^2\mathbf{u},$$

$$\displaystyle\frac{\partial\Theta}{\partial t}+\bm{{u}}\cdot\nabla\Theta+u_{r}\left(\frac{dT_{\rm{rad}}}{dr}-\frac{dT_{\rm{ad}}}{dr}\right)=\kappa\nabla^{2}\Theta{\color[rgb]{0,0,0},}$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}=-\nabla p+\frac{{\text{Ra}}_o}{\text{Pr}}\Theta {\mathbf{e}}_{\mathbf{r}}+{\nabla }^2\mathbf{u},$$

$$\frac{\partial \Theta }{\partial t}+\mathbf{u}\cdot \nabla \Theta +\beta (r)u_r=\frac{1}{\mathrm{Pr}}{\nabla }^2\Theta .$$

$$\text{Pr}=\frac{\nu }{\kappa }\text{and}{\text{Ra}}_o=\frac{\alpha g\left|\frac{d{T}_o}{dr}-\frac{d{T}_{\mathrm{ad}}}{dr}\right|r_o^4}{\kappa \nu },$$

$$\beta (r)=\frac{\frac{d{T}_{\mathrm{rad}}}{dr}-\frac{d{T}_{\mathrm{ad}}}{dr}}{\left|\frac{d{T}_o}{dr}-\frac{d{T}_{\mathrm{ad}}}{dr}\right|}=\frac{{\left(\frac{1}{r}\right)}^2\frac{d{T}_o}{dr}-\frac{d{T}_{\mathrm{ad}}}{dr}}{\left|\frac{d{T}_o}{dr}-\frac{d{T}_{\mathrm{ad}}}{dr}\right|}.$$

$$\beta (r)=-\frac{\mathrm{Ra}(\mathrm{r})}{{\mathrm{Ra}}_{\mathrm{o}}},$$

$$\displaystyle\text{Ra}(r)=\frac{\alpha g\left|\displaystyle\frac{dT_{\rm{rad}}}{dr}-\frac{dT_{\rm{ad}}}{dr}\right|r_{o}^{4}}{\kappa\nu}{\color[rgb]{0,0,0}=\frac{\alpha g\left|\displaystyle-\frac{F_{\rm{rad}}}{\kappa}\left(\frac{r_{i}}{r}\right)^{2}-\frac{dT_{\rm{ad}}}{dr}\right|r_{o}^{4}}{\kappa\nu}}.$$

$$\chi=\displaystyle\left|\frac{\displaystyle\frac{dT_{o}}{dr}-\displaystyle\frac{dT_{\rm{ad}}}{dr}}{\displaystyle\frac{dT_{o}}{dr}}\right|{\color[rgb]{0,0,0}=\displaystyle\left|1+\dfrac{\kappa\dfrac{dT_{\rm{ad}}}{dr}r_{o}^{2}}{F_{\rm{rad}}r_{i}^{2}}\right|},$$

$$\beta (r)=\frac{1-\chi -(1/r{)}^2}{\chi }.$$

$$\bar{q}(r)=\displaystyle\frac{1}{{\color[rgb]{0,0,0}4\pi(t_{2}-t_{1})}}\int_{t_{1}}^{t_{2}}\int_{0}^{2\pi}\int_{0}^{\pi}q(r,\theta,\phi,t)\sin\theta d\theta d\phi.$$

$$\displaystyle{\color[rgb]{0,0,0}{\bar{E}}_{k}(r)=\displaystyle\frac{1}{2}(\overline{{u_{r}^{2}+u_{\theta}^{2}+u_{\phi}^{2})}}}.$$

$${\color[rgb]{0,0,0}\bar{N}^{2}(r)}=\alpha g\left(\frac{d\bar{T}}{dr}-\frac{dT_{\rm{ad}}}{dr}\right)\frac{r_{o}^{4}}{\nu^{2}}=\left(\beta(r)+\frac{d\bar{\Theta}}{dr}\right)\frac{\text{Ra}_{o}}{\text{Pr}},$$

$$\text{Ra}(r)=-\frac{\alpha (r)g\left(\frac{d{T}_{\mathrm{rad}}}{dr}-\frac{d{T}_{\mathrm{ad}}}{dr}\right)r_o^4}{\kappa \nu }=-{\text{Ra}}_o\frac{\alpha (r)}{{\alpha }_o}\beta (r)\equiv {\text{Ra}}_o$$

$$\displaystyle\frac{\partial{\bm{{u}}}}{\partial{t}}+{\bm{{u}}}\cdot{\nabla}{\bm{{u}}}=-{\nabla}{p}+\frac{\alpha(r)}{\alpha_{o}}\frac{\text{Ra}_{o}}{\text{Pr}}{\Theta}\bm{e_{r}}+{\nabla}^{2}{\bm{{u}}}{\color[rgb]{0,0,0},}$$

$$\displaystyle{\color[rgb]{0,0,0}{\rm{{Ra}}}_{b}}=\frac{\int_{r_{i}}^{r_{o}}{\rm{Ra}}(r)r^{2}dr}{\int_{r_{i}}^{r_{o}}r^{2}dr}.$$

$${\bar{N}}^2(r)=\frac{\alpha (r)}{{\alpha }_o}\left(\beta (r)+\frac{d\bar{\Theta }}{dr}\right)\frac{{\text{Ra}}_o}{\text{Pr}}.$$

$${\bar{F}}_h-\frac{1}{\mathrm{Pr}}\frac{d\bar{\Theta }}{dr}=0.$$

$${\bar{N}}^2(r)=\frac{\alpha (r)}{{\alpha }_o}(\beta (r)+\mathrm{Pr}{\bar{F}}_h)\frac{{\mathrm{Ra}}_{\mathrm{o}}}{\mathrm{Pr}},$$

Full text

I INTRODUCTION

Convection is a ubiquitous physical process in geophysical fluid dynamics, which has been extensively studied analytically, experimentally and numerically because of the vital role it plays in the global dynamics of the Earth’s mantle (e.g. Bercovici et al. (1989a); Glatzmaier and Schubert (1993); Breuer et al. (2013) and for a review see Schubert (1992) and references therein), oceans Marshall and Schott (1999) and atmosphere Emanuel (1994). Convection is also important in astrophysical settings such as the gaseous interiors of stars and planets where the convective zones are usually global, either spanning the entire object or at least a deep spherical shell. By contrast with geophysical convection, relatively little is known about convection in astrophysical objects. Observationally speaking, a limited amount of information can be obtained either through direct imaging of the surface (e.g. see Rieutord and Rincon (2010) for a review), or indirectly using asteroseismology to infer, for instance, the mean temperature profile within the convection zone Gough et al. (1996). Meaningful physical experiments are almost impossible to design because the governing parameters appropriate to the interiors of stars and planets are vastly different from those achievable in a laboratory. In particular, the Prandtl number, which is the ratio of the kinematic viscosity to the thermal diffusivity, is much smaller than unity in astrophysical plasmas (e.g. \sim\color[rgb]{0,0,0}10^{-2} in giant planets, and $\sim 10^{-6}$ in stars) whereas it is usually of order unity or much larger in geophysical applications. Among other things, this implies that the ordering of the relevant dynamical timescales is different in the two regimes, and that the effects of the inertial terms in astrophysical convection are much larger than in geophysical convection.

In this paper, we attempt to shed new light on the subject by investigating the dynamics of convection in weakly compressible gaseous spherical shells in the low Prandtl number parameter regime using direct numerical simulations (DNSs) with mixed temperature boundary conditions (here meant to imply fixed flux at the inner boundary and fixed temperature at the outer boundary). This setup is designed to capture some of the most salient features of convection in stars and giant planets, and differs in significant ways from most studies of convection to date.

Arguably, the most commonly studied form of convection is thermal Rayleigh-Bénard convection (RBC thereafter) between two parallel plates where a Boussinesq liquid (in the original Boussinesq sense Boussinesq (1903); Oberbeck (1879)) is heated from below and cooled from above, and the two rigid boundaries are held at constant temperatures. For sufficiently strong driving, as measured by the Rayleigh number, buoyancy forces overcome thermal and viscous damping and turbulent heat transport by convection dominates conduction. This highly symmetric idealized model setup has been studied extensively in both 2D and 3D Veronis (1966); Moore and Weiss (1973); Shraiman and Siggia (1990); Julien et al. (1996); Guervilly et al. (2014) (also for a general review of RBC see Ahlers et al. (2009) and references therein).

When studying geophysical problems, several extensions of this basic model are usually considered depending on the specific application. Studies of mantle convection usually adopt a spherical shell geometry and consider the limit of infinite Prandtl number Chandrasekhar (1961); Busse (1975); Zebib et al. (1980); Schubert and Zebib (1980); Zebib et al. (1985); Machetel et al. (1986); Bercovici et al. (1989b); Vangelov and Jarvis (1994); Jarvis et al. (1995); Choblet (2012). More generally, geophysically motivated studies of convection in spherical shells sometimes include the effect of rotation or allow for a finite Prandtl number Dormy et al. (2004); Tilgner (1996); Feudel et al. (2011); Gastine et al. (2015, 2016) but have so far nearly always used fixed temperature boundary conditions. The majority of these investigations have focused on the derivation of scaling laws for global quantities such as the heat flux or the total kinetic energy as functions of input parameters, as well as developing models for the boundary layers (see Grossmann and Lohse (2000) for a review).

In astrophysical applications, on the other hand, the fluid is generally compressible. Solving the compressible Navier-Stokes equations requires the resolution of timescales associated with fast sound waves, as well as the much slower timescales associated with global thermal or viscous adjustment. This stiffness is a severe impediment to simulation and filtering out the fast sonic dynamics is often desirable. One way of accounting for weak compressibility in astrophysical convection is through the anelastic approximation Batchelor (1953); Ogura and Phillips (1962); Gough (1969); Latour et al. (1976); Gilman and Glatzmaier (1981), which filters out sound waves while allowing for strong variations in the background density. This is the more commonly adopted formalism for the study of solar convection (Miesch (2005) and references therein) and stellar convection (e.g. Browning et al. (2004); Brown et al. (2008); Augustson et al. (2012); Brown et al. (2012)) but it has significant drawbacks. First of all, there are numerous formulations of the approximation and there is some debate about their relative validity Vasil et al. (2013); Verhoeven et al. (2015). Secondly, the anelastic approximation is usually based on the assumption of small departures from adiabaticity which is not guaranteed in all reasonable problems.

Another commonly used approximation under which sound waves are filtered out is the Boussinesq approximation for gases Spiegel and Veronis (1960). It is important to note that the standard Boussinesq approximation Boussinesq (1903); Oberbeck (1879) should not be used in astrophysical applications because of the compressibility of the gas (although it is still sometimes used for simplicity Gilman (1976, 1977, 1978a); Christensen (2002); Heimpel et al. (2005); Heimpel and Aurnou (2007); Aurnou et al. (2007)). However, Spiegel & Veronis (1960) Spiegel and Veronis (1960) (SV thereafter) showed that it is possible to generalize the approximation to take into account some effects of compressibility, thereby allowing its use in modeling convection in gaseous systems, such as the Earth’s atmosphere or the interiors of stars and planets. Assuming that the size of the convective region is much smaller than any scale height of the system (including the local radius, if the convection zone is a spherical shell), and that the fluid motions are much slower than the local speed of sound, they showed that the only effect of compressibility is to heat or cool a parcel of fluid as it shrinks or expands to adjust to the ambient hydrostatic pressure. As a result, their formulation contains an additional term in the temperature equation which is proportional to the local adiabatic temperature gradient (which is non-zero for gases) but is otherwise identical to the traditional Boussinesq approximation. Studies of astrophysical convection in Cartesian geometry or in a very thin spherical shell using the SV Boussinesq (SVB thereafter) approximation were presented for instance in Cattaneo et al. (2003) and Miesch (2001).

In this work we propose to study 3D DNSs of low Prandtl number convection in a spherical geometry using the SVB equations with particular applications to solar and stellar convection in mind. We therefore deviate from the usual assumption of fixed temperature boundary conditions and instead, consider a more astrophysically realistic setup with fixed flux at the inner boundary and fixed temperature at the outer boundary. Indeed, in stars like the Sun for instance, the flux through the base of an outer convection zone is fixed by the luminosity of the star, which in turn is set by the nuclear generation rate within the core.

These four elements (spherical geometry, weak compressibility, mixed temperature boundary conditions and low Prandtl number) have never, to the authors’ knowledge, been used in conjunction and yet are all crucial elements of astrophysical convection. As we demonstrate in this work, their combined effect is to create substantial asymmetry between the upflows and downflows, which in turn transforms the nature of convection near the lower boundary, with surprising repercussions throughout the entire convection zone.

The paper is organized as follows: Sections II and III present two different model setups that both possess all four properties listed above, and report on the results of a suite of numerical simulations for varying Rayleigh number. Surprising new dynamics are observed, which are then analyzed and explained in detail in Section IV. In Section V, we explore a third model setup which more closely resembles the Sun (although still simplistically), in order to test the robustness of our results and to assess the applicability of what we have found to the circumstances that most interest us. Finally, in Section VI, we summarize our results and discuss the possible limitations of the SVB approximation.

II Boussinesq convection in a weakly compressible spherical shell

II.1 Mathematical formulation

We begin our systematic investigation of the effects of mixed temperature boundary conditions and weak compressibility on the dynamics of Rayleigh-Bénard convection in a spherical shell by constructing the simplest possible model with these properties. In this model, and in all of the ones that follow, we consider a spherical shell located between an inner sphere of radius $r_{i}$ and an outer sphere of radius $r_{o}$. For simplicity, we assume constant thermal expansion coefficient $\alpha$, viscosity $\nu$, thermal diffusivity $\kappa$, adiabatic temperature gradient $dT_{\rm{ad}}/dr$ and gravity $g$. In the absence of fluid motion and when the system is in a steady state, the background radiative temperature gradient is obtained by solving

[TABLE]

where $r$ is the local radius. The inner fixed flux boundary condition implies that

[TABLE]

where $F_{\rm{rad}}$ is the temperature flux per unit area through the inner boundary, whereas the outer fixed temperature boundary condition is $T(r_{o})=T_{o}$. Then, solving Eq. (1) using the first boundary condition implies that

[TABLE]

which, along with the second boundary condition, gives

[TABLE]

We clearly see that, in contrast to the Cartesian case, the radiative temperature gradient in a spherical geometry is not constant but depends on the radius. This implies in turn that $dT_{\rm{rad}}/dr-dT_{\rm{ad}}/dr$ also varies with depth. Note that for the SVB approximation to be valid, $\Delta T=T_{\rm{rad}}(r_{i})-T_{\rm{rad}}(r_{o})$ must be much smaller than, say, $T_{o}$. This is true either for small enough $r_{o}-r_{i}$ (thin layer) given $F_{\rm{rad}}$, or for small enough $F_{\rm{rad}}$ given $r_{o}-r_{i}$.

We now let $T(r,\theta,\phi,t)=T_{\rm{rad}}(r)+\Theta(r,\theta,\phi,t)$ where $\Theta$ is the temperature perturbation to the radiative background. We also assume a linear relationship between the temperature and density perturbations consistent with the SVB approximation, $\rho/\rho_{m}=-\alpha\Theta$, where $\rho_{m}$ is the mean density of the background fluid. With these assumptions, the governing SVB equations are:

[TABLE]

[TABLE]

and

[TABLE]

where $\bm{{u}}=(u_{r},u_{\theta},u_{\phi})$ is the velocity field, and $p$ is the pressure. We non-dimensionalize the problem by using 111We numerically solve the non-dimensional Boussinesq equations in which we have used the outer radius as the lengthscale. If we wanted to compare spherical numerical simulations with Cartesian ones, we would have to non-dimensionalize the problem using the thickness of the shell $[l]=r_{o}-r_{i}=L$ such that both the problems could have the same effective Rayleigh number for accurate comparison. $[l]=r_{o}$, $[t]=r_{o}^{2}/\nu$, $[u]=\nu/r_{o}$ and $[T]=|dT_{o}/dr-dT_{\rm{ad}}/dr|r_{o}$ as the unit length, time, velocity and temperature, where ${dT_{o}}/{dr}\equiv{dT_{\rm{rad}}}/{dr}|_{r=r_{o}}$. Then, we can write the non-dimensional equations as:

[TABLE]

[TABLE]

and

[TABLE]

All the variables and parameters are now implicitly non-dimensional, which introduces the Prandtl number Pr and the Rayleigh Rao defined as

[TABLE]

and the non-dimensional superadiabaticity

[TABLE]

Another way to interpret $\beta$ is to note that it is minus the ratio of the local Rayleigh number Ra$(r)$ to Rao i.e.

[TABLE]

where

[TABLE]

Finally, note that while $\beta$ seems to depend on two dimensional quantities $dT_{\rm{rad}}/dr$ and $dT_{\rm{ad}}/dr$ (see Equation (12)), it can be rewritten in this simple model just in terms of a single non-dimensional parameter $\chi$, defined as

[TABLE]

so that

[TABLE]

Note that $\beta(1)=-1$ and $\beta(r_{i}/r_{o})=(1-\chi-(r_{o}/r_{i})^{2})/\chi$.

Figure 1 illustrates how $\beta$, and therefore the local Rayleigh ratio Ra$(r)$/Rao, depends on $\chi$. Note that for $\chi=0.5$ the local Rayleigh number at the inner boundary is about 3 times larger than Rao, whereas for $\chi=0.1$, it is 11 times larger, illustrating that a small $\chi$ implies a stronger variation of the local Rayleigh number across the shell. In the limit of a very thin shell $(r_{i}/r_{o}\rightarrow 1)$ on the other hand (which recovers the case of convection between infinite parallel plates), $\beta(r)$ tends to the constant $-1$ regardless of $\chi$. The functional form of $\beta$ therefore depends on our choice of boundary conditions and on the fact that we are operating in an appreciably deep spherical shell (see Equation (16)).

II.2 Numerical simulations

In order to study the influence of weak compressibility and sphericity (which manifest themselves in a variable $\beta(r)$), and of mixed temperature boundary conditions on the model dynamics, we have run 3D DNSs solving Equations (8)-(10) in a spherical shell, exactly as outlined above, using the PARODY code Aubert et al. (2008). The boundary conditions for the temperature perturbations $\Theta$ are such that we have fixed flux at the inner boundary, $\partial\Theta/\partial r|_{r_{i}}=0$ and fixed temperature at the outer boundary, $\Theta(r_{o})=0$. The velocity boundary conditions are stress-free at both the inner and outer boundaries. The simulations discussed in this section are referred to as “Model A” simulations. Table 1 summarizes our various runs in this setup, as well as those later discussed in Sections III and V. Note that $r_{i}/r_{o}=0.7$ and Pr $=0.1$ for all the simulations presented in this paper.

We now examine the qualitative and quantitative properties of our simulations, focusing on three typical cases with varying $\chi$ ($\chi=0.1$, $\chi=0.5$ and $\chi=1$) for fixed Ra${}_{o}=10^{7}$. A simple way of visualizing the turbulent motions due to convection is to look at snapshots of the velocity components $u_{r}$, $u_{\theta}$ or $u_{\phi}$ at a typical time after saturation of the linear instability. Figure 2 shows snapshots of $u_{r}$. In each panel, the left hemisphere shows the velocity field on a spherical shell close to the upper boundary, illustrating the convective motions near the surface. The right hemisphere is a meridional slice showing the radial velocity as a function of depth and latitude, for a selected longitude. Figure 2(a) is for $\chi=0.1$, while Figure 2(b) is for $\chi=0.5$. We notice that the $\chi=0.1$ case appears somewhat more turbulent than the $\chi=0.5$ case, as visualized by stronger eddies with a wider range of scales.

To see more clearly any difference among the runs, we turn to more quantitative measures. Figure 3(a) shows the kinetic energy per unit volume $E_{k}$ within the shell as a function of time for the Model A simulations (solid lines). We clearly observe the initial development of the convective instability, visible as a large spike in the interval $t\in[0,0.01]$, followed by its nonlinear saturation. Note that $E_{k}(t)$ reaches a stationary state very fast but reaching thermal equilibrium is a much slower process. We estimate that a simulation has reached thermal equilibrium when $\partial\Theta/\partial r|_{r=r_{o}}$ is statistically stationary and equal to zero. This happens around $t\approx 0.02$ for the $\chi=1$, $\chi=0.5$ and $\chi=0.1$ simulations. In all that follows, we only present the results of simulations once they have achieved thermal equilibrium.

Figure 3(a) shows that $E_{k}(t)$ is much larger for the $\chi=0.1$ run than for cases with larger $\chi$, confirming our rapid visual inspection of Figure 2. To understand why this may be the case, recall that for smaller values of $\chi$ the local Rayleigh number Ra$(r)$ increases more with depth than for larger $\chi$ (Fig. 1). A higher Rayleigh number near the lower boundary drives convection more vigorously, which increases the overall kinetic energy.

Throughout the paper, we define the time- and spherical- average of a quantity as

[TABLE]

Figure 3(b) shows the non-dimensional kinetic energy profiles $\bar{E}_{k}(r)$ given by

[TABLE]

The forms of these profiles look similar for $\chi=0.1$, $\chi=0.5$ and $\chi=1$, taking their highest value at the top of the convection zone and then decreasing inward to a plateau from approximately $r=0.95$ down to $r=0.75$. Below $r=0.75$, there is a small increase in the kinetic energy associated with the inner boundary layer. As we already saw in Figure 3(a), the $\chi=0.1$ case has significantly higher kinetic energy than the other runs. Somewhat surprisingly, however, the kinetic energy is larger everywhere even though Ra$(r)$/Rao is only larger near the inner boundary. This could be explained by the fact that the convection in this model is a highly non-local process i.e. that stronger driving deeper down implies strong upflows and downflows throughout the domain.

In Figure 4, we plot the square of the non-dimensional buoyancy frequency

[TABLE]

(solid line) for $\chi=0.1$, $\chi=0.5$, $\chi=1$ and Ra${}_{o}=10^{7}$. We also show the square of the background buoyancy frequency {\color[rgb]{0,0,0}N^{2}_{\rm{rad}}(r)}=\beta(r)Ra${}_{o}/$Pr as a dashed line for reference. As expected, we find that convective motions outside the boundary layers generally mix potential temperature and drive the mean radial temperature gradient towards an adiabatic state where $\bar{N}^{2}\approx 0$. However, subtle differences arise when $\chi$ decreases, which manifest themselves in two different ways. Firstly, note that for lower $\chi$, $|N^{2}_{\rm{rad}}|$ is much larger, consistent with stronger convective driving. Nonetheless, even though the kinetic energy per unit volume is larger, we see that the interior is not mixed as well for $\chi=0.1$ as for $\chi=0.5$ and $\chi=1$. Secondly, for $\chi=0.1$, we observe the surprising emergence of a slightly subadiabatic region ($\bar{N}^{2}>0$) just below the upper boundary layer.

This remarkable behavior, namely the emergence of a layer in the flow that is subadiabatic and therefore ostensibly convectively stable, only occurs for the lowest value of $\chi$ we were able to simulate. Proceeding to lower $\chi$ to test the robustness of this observation would be an obvious path, but one that is numerically difficult. For example, using $\chi=0.01$ would require the Rayleigh function Ra$(r)$ to reach values of approximately $100$Rao at the inner boundary. Such a range is hard to resolve. For this reason, and furthermore to elicit which of the physics elements are responsible for the unexpected emergence of a subadiabatic layer, we now switch to a different model setup (Model B).

III Spherical shell with a constant Rayleigh function

In the Model A simulations discussed in the previous section, both $\beta(r)$ and the local Rayleigh number Ra$(r)$ vary with depth proportionally to one another. As a result, it is difficult to determine what may be the direct cause of some of the interesting features we observe. Thus, we construct a second model (called “Model B”) where $\beta(r)$ is the same as in Section II, but where Ra$(r)$ is constant across the convection zone. We can achieve this by varying the thermal expansion coefficient $\alpha$ with radius in order to compensate for the radial variation of $dT_{\rm{rad}}/dr-dT_{\rm{ad}}/dr$. Continuing to assume that $\kappa$, $\nu$ and $g$ are constant, we now choose $\alpha(r)$ such that Ra$(r)$=Rao. That is,

[TABLE]

as long as {\alpha(r)}/{\alpha_{o}}\equiv-{1}/{\beta(r)}{\color[rgb]{0,0,0},} where $\alpha_{o}=\alpha(r_{o})$. In this new setup, the non-dimensional momentum equation is

[TABLE]

while the thermal energy equation remains unchanged, and is given by Eq. (10).

Figure 5 shows the variation of $\alpha$ needed to keep the Rayleigh function constant for our fiducial values of $\chi$. In all cases, $\alpha(r)$ decreases with depth, and $\alpha(r_{i})/\alpha_{o}$ is smaller for smaller $\chi$. Physically speaking, this implies that the effective buoyancy of fluid elements with fixed temperature perturbation $\Theta$ decreases with depth. Note that, in contrast with Model A, Model B now explicitly violates the conditions of use of the SVB approximation when $\chi$ is small. This will be discussed in Section VI, but lends credence to our use of the paper title “weakly non-Boussinesq convection”.

We now compare the convective dynamics of the Model A and B setups in order to try and understand the respective roles of Ra$(r)$ and $\beta(r)$ in driving convection and mixing. To do so, we have run numerical simulations using Model B for three different values of $\chi$, for a fixed Ra$(r)=$Rao=$10^{7}$. In these constant Rayleigh function runs, we were able to achieve values of $\chi$ down to $0.001$.

Figure 3 compares the energetics of Model A (solid lines) and Model B (dashed lines) runs. In Figure 3(a), we observe that the saturation level of the kinetic energy per unit volume $E_{k}$ varies much less with $\chi$ in Model B than in Model A. This might be expected since both Ra and Pr are now constant at all radii and in all configurations of Model B presented. In Figure 3(b), we see that the kinetic energy profiles $\bar{E}_{k}(r)$ of the various Model B runs almost coincide in the bulk of the convection zone, showing that $\beta$ alone does not influence this quantity much.

Figure 6 shows the mean kinetic energy ${E}$ (i.e. the time average of $E_{k}$) against the bulk Rayleigh number defined as

[TABLE]

We see that the mean kinetic energy depends solely on the bulk Rayleigh number (for fixed Pr and $r_{o}-r_{i}$) and not on the setup used or on the value of $\chi$. This is a very interesting finding, since it illustrates that the mean kinetic energy is model-independent and can be predicted as long as the bulk Rayleigh number {\color[rgb]{0,0,0}{\rm{{Ra}}}_{b}} of the problem is known. Fitting the available data, we find that E=(3.7\pm 2.6){\color[rgb]{0,0,0}{\rm{{Ra}}}_{b}}^{0.72\pm 0.04}. Note that we expect the exponent to be universal, but, the prefactor likely depends on Pr or on $(r_{o}-r_{i})$.

In order to compare the efficiency of mixing in this new system, we again look at the square of the non-dimensional buoyancy frequency, defined for Model B as

[TABLE]

In Figure 7, we plot $\bar{N}^{2}(r)$Pr/Rao compared with the background $N^{2}_{\rm{rad}}$Pr/Ra${}_{o}=(\alpha(r)/\alpha_{o})\beta(r)$ for Model B. Note that by construction in this setup $N^{2}_{\rm{rad}}$Pr$/$Ra${}_{o}=-1$ regardless of $\chi$. We see that as $\chi$ decreases, $\bar{N}^{2}$ increases and for $\chi\leq 0.1$ a subadiabatic region does indeed emerge as in Model A. This unusual effect is much more pronounced at $\chi=0.01$. Overall, this conclusively shows that the appearance of the subadiabatic region is not model-dependent, but instead, a fairly generic property of simulations that combine mixed temperature boundary conditions with varying superadiabaticity.

In order to determine more precisely how the emergence of a subadiabatic region depends on the model parameters, we ran additional simulations at Ra${}_{o}=10^{6}$ and Ra${}_{o}=10^{8}$ for $\chi=0.01$, as well as a simulation with Ra${}_{o}=10^{8}$ for $\chi=0.1$. Figure 8 shows the square of the buoyancy frequency profiles for these comparative runs. We observe that, for a given value of $\chi$, there is a threshold value of Rao above which the subadiabatic region appears, and that the size and subadiabaticity of that region increases with Rao beyond that threshold. For fixed Rao we see a similar behavior with decreasing $\chi$. These considerations suggest that the subadiabatic layer appears only for sufficiently vigorous convection (high Rayleigh number) and/or in systems with sufficiently large radial variations in the background superadiabaticity (here generated by low $\chi$).

Interestingly, convection appears unaffected by the emergence of the subadiabatic layer, and proceeds as if it did not exist. This can be seen both in snapshots of the velocity field (Figure 9) and in the kinetic energy profiles as a function of radius (Figure 10). Figure 9 shows snapshots of $u_{\phi}$ as a function of radius and latitude, for a fixed longitude, for $\chi=0.01$ and the three different Rayleigh numbers used in that case. As the Rayleigh number increases, the convective eddies are more pronounced and the turbulent motions are apparently stronger. However, none of the simulations show any obvious indication of a non-convective or “dead” zone due to the subadiabatic layer (which is present for the Ra${}_{o}=10^{7}$ and Ra${}_{o}=10^{8}$ cases). The same can be seen more quantitatively in Figure 10, which shows the kinetic energy profiles for the same three cases (Ra${}_{o}=10^{6}$, $10^{7}$ and $10^{8}$, and $\chi=0.01$). As in Model A, we find that they have roughly the same shape, but that the kinetic energy increases with Rao. Crucially, there is no sign of any dip in the kinetic energy profiles at the locations of the subadiabatic layers in the Ra${}_{o}=10^{7}$ and Ra${}_{o}=10^{8}$ runs, which proves that convection is efficient everywhere across the shell. All the above provide strong indications that convection in these models is a very non-local process.

IV Interpretation of the results

Having established that the emergence of a subadiabatic layer is a robust phenomenon in these models, we now proceed to explain the observed dynamics more quantitatively. As we shall demonstrate, the phenomenon is directly related to the use of mixed temperature boundary conditions, and is facilitated by the presence of a strongly-varying background superadiabaticity which act together to create strong asymmetries between upflows and downflows.

We note that in traditional Boussinesq Rayleigh-Bénard convection between parallel plates with fixed temperature boundary conditions (BRBC thereafter), an initially superadiabatic mean temperature profile becomes relaxed via convective motions to a state where upper and lower superadiabatic boundary layers are joined by an adiabatic interior. Upflows and downflows are driven by buoyancy forces in the boundary layers, and interact non-linearly in the bulk of the fluid, mixing it towards an adiabat. The intrinsic up-down symmetry of BRBC implies that upflows and downflows contribute equally to the upward heat transport. As we now demonstrate, our findings here are very different.

In Figure 11(a), we show the temperature perturbation profiles $\bar{\Theta}(r)$ in Model B runs with $\chi=0.01$, and Ra${}_{o}=10^{6}$, Ra${}_{o}=10^{7}$ and Ra${}_{o}=10^{8}$. This quantity is proportional to the term (Ra$/$Pr)$\bar{\Theta}(r)$ which is the time- and spherically-averaged non-dimensional buoyancy force in the statistically stationary state (see Equation (21)). We notice that $\bar{\Theta}$ is negative almost everywhere so the average buoyancy force is downwards. That is quite different from what happens in the BRBC case (with or without adiabatic temperature gradient), where $\bar{\Theta}$ would be found to be positive in the bottom half of the domain, and negative in the top half. This raises two questions: (1) why is $\bar{\Theta}$ almost entirely negative in our setup, and (2) how are the upflows driven if the average buoyancy force is downwards?

To answer the first question, we look at the behavior of an individual fluid parcel. Ignoring thermal diffusion, the evolution of the temperature within the parcel is given by the Lagrangian derivative $D\Theta/Dt=-\beta u_{r}$, which can be re-written as {\color[rgb]{0,0,0}D\Theta/Dr}=-\beta since the radial velocity of the parcel is $u_{r}=Dr/Dt$. The temperature of the parcel at the boundary is given by $\Theta(r_{o})=0$, therefore by integration inwards, the temperature within downward flowing parcels is $\Theta_{\rm{ad}}(r)=-\int_{r_{o}}^{r}\beta(r^{\prime})dr^{\prime}$. This quantity represents the temperature perturbation within a fluid parcel at radius $r$ that is moving downward adiabatically and without any mixing with its surroundings. Comparing $\Theta_{\rm{ad}}(r)$ to the mean temperature perturbation profiles $\bar{\Theta}(r)$ of Model B simulations with varying Rao in Figure 11(a), we see that they coincide throughout most of the domain except when approaching the inner boundary. This suggests that indeed the fluid parcels travel downwards more or less adiabatically, and that is what controls $\bar{\Theta}$ over much of the spherical shell.

In answer to the second question, there are two possibilities: either the mean hides information about the upflows and perhaps the return flows from the lower boundary are very rare but strong, arising from the tail of the distribution of the temperature perturbations $\Theta$ around the mean $\bar{\Theta}$, or pressure dominates over buoyancy and pushes the parcels back upwards. Figure 11(b) shows $\Theta$ profiles at specific example locations (fixed longitude but various latitudes) and a specific time, and indicates that there really are no fluid parcels with positive $\Theta$ anywhere in the lower part of the domain. We thus conclude that the upflows are not buoyantly driven, and must be pressure dominated.

This idea can be confirmed by looking at the turbulent temperature fluxes and the respective contributions from the upflows and downflows directly. In quintessential convection such as BRBC, the direction of the force acting on a fluid parcel is given by the sign of $\Theta$. As a result, aside from short transients, there is a very strong correlation between the sign of the temperature perturbation and the sign of the vertical velocity of the fluid parcel (with $\Theta>0$ corresponding to $u_{r}>0$, and $\Theta<0$ corresponding to $u_{r}<0$). Therefore, in both cases, one might expect the turbulent (or convective) temperature flux $F_{h}=u_{r}\Theta$ to be positive, if motion is due to buoyancy. With this in mind, we therefore examine the contributions from upflows and downflows to the total turbulent temperature flux separately. We investigate by looking firstly at the mean flux (Fig. 12) and then secondly at the pointwise flux (Fig. 13), for a characteristic Model B simulation with $\chi=0.01$, and Ra${}_{o}=10^{7}$ at a typical time when the system is in a statistically stationary and thermally relaxed state. In Figure 12, we plot $\bar{F}_{h}$ together with the corresponding mean turbulent temperature flux carried by the upflows only, $\bar{F}_{up}$ (given by $\bar{F}_{h}$ using only those points where $u_{r}>0$) and by the downflows only, $\bar{F}_{down}$ (given by $\bar{F}_{h}$ using only those points where $u_{r}<0$). We notice that the average flux carried by the downflows is positive, which means that they transport relatively cold material downward as expected. However, the average flux carried by the upflows is negative, indicating that they are carrying cold material up, contrary to expectations. In standard BRBC for instance, both downflows and upflows would on average transport heat upwards, i.e. $\bar{F}_{h}$, $\bar{F}_{up}$ and $\bar{F}_{down}$ would all be positive, with downflows carrying relatively cold fluid and upflows carrying relatively warm fluid.

To check whether this odd behavior of the upward fluxes is true only in the average, or applies to all fluid parcels, we now look at the pointwise turbulent temperature flux, $F_{h}=u_{r}\Theta$. Figure 13 shows the pointwise flux at every point on two spherical shells, located close to the inner boundary (at $r\approx 0.75$, in Figure 13(a)) and in the bulk of the domain (at $r\approx 0.85$, in Figure 13(b)) respectively at a representative time. As mentioned above, one might more normally (in BRBC for example) expect that nearly all points would lie in the two upper quadrants. Figures 13(a) and 13(b) confirm our findings from Figure 12, but provide more detail. While the downflows all appear to be working productively at transporting heat upwards (i.e. cold downwards; upper left quadrant), it is clear that a considerable number of upflows are working counter-productively at both depths. At $r\approx 0.85$, in the middle of the domain, some of the upflows work productively to transport heat upwards (upper right quadrant) but many more work counter-productively to transport heat downwards (lower right quadrant), leading to the negative mean flux in the upflows seen in Figure 12. Closer to the inner boundary all of the upflows are counterproductive to heat transport.

All of these results point to the same conclusion, namely that the upflows are not buoyantly driven in the lower part of the domain. Since the only other force in the system is the pressure gradient, we conclude that the upflows must be pressure-driven. The dynamic pressure is really only a manifestation of the divergence condition in the Boussinesq approximation, thus, another equivalent interpretation is that the upflows are merely an inertial continuation of the downflows which are forced to turn around at the lower boundary. Indeed, if a downflow driven by its negative temperature perturbation simply “rebounds” off the lower boundary without changing its heat content, it becomes a counter-productive upflow, i.e. fluid moving upwards but with the same negative temperature perturbation. It is perhaps surprising that this result persists so high into the shell. This appears to be a feature of low $\chi$, low Pr and high Rao convection.

The above results are implicitly related to the choice of mixed boundary conditions for $\Theta$. The fixed flux at the inner boundary is a source of strong asymmetry in the dynamics of the problem since it allows the temperature perturbations $\Theta$ to be negative there, roughly following the adiabatic profile $\Theta_{\rm{ad}}(r)$. In a system with fixed temperature conditions on the other hand, $\Theta$ and therefore $\bar{\Theta}$ would be forced to be zero at both bottom and top boundaries, and the system would be much more symmetric (though not perfectly because of the sphericity and the non-zero constant adiabatic temperature gradient which are additional sources of asymmetry). This would then guarantee the existence of temperature perturbations of both signs near the inner boundary and therefore some buoyantly-driven upflows there.

The use of mixed boundary conditions has a second very important impact on the convective dynamics, namely that the total perturbed temperature flux through the system (turbulent + diffusive) must be equal to that at the inner boundary, and thus zero everywhere (Fig. 14). Non-dimensionally, this is expressed as

[TABLE]

In thermal equilibrium, the diffusive and non-diffusive contributions to the perturbed temperature flux must therefore cancel out exactly. The magnitude of the temperature perturbations depends on $\chi$ (through the increasingly negative values of $\Theta_{\rm{ad}}$ as $\beta(r)$ decreases rapidly with $r$ for low $\chi$) as in Figure 11(a). Furthermore, the rms velocity of the convective eddies increases substantially with Rao (see Figure 6). Thus for low $\chi$ and high Rao, the turbulent flux increases and the diffusive flux of the temperature perturbations must follow accordingly. This is crucial, and causes the emergence of the subadiabatic layer in our simulations as follows.

Using Equation (24) we can re-write $\bar{N}^{2}$ as

[TABLE]

where we carefully note that $\beta(r)<0$ while $\bar{F}_{h}>0$. Since $\bar{F}_{h}$ increases monotonically with increasing Rao (because of the increase in $u_{rms}$) or decreasing $\chi$ (because $\Theta_{rms}$ is larger), there exists a region of parameter space where $(\beta(r)+{\rm{Pr}}\bar{F}_{h})$ becomes positive, at least somewhere within the shell, leading to a positive $\bar{N}^{2}(r)$.

Using the results we have obtained so far, we can in fact provide an order-of-magnitude estimate for $\bar{F}_{h}$ as a function of $\chi$ and of the bulk Rayleigh number {\color[rgb]{0,0,0}{\rm{{Ra}}}_{b}} given in (22). The typical amplitude of the temperature perturbations $\Theta_{rms}$ can be estimated from $\Theta_{\rm{ad}}$, which is proportional to $1/\chi$ for low enough $\chi$ . The rms velocity of the flow $u_{rms}$ can be estimated from {\color[rgb]{0,0,0}{\rm{{Ra}}}_{b}} using $u_{rms}=\sqrt{2{E}}$. In Section III, we found that E=(3.7\pm 2.6){\color[rgb]{0,0,0}{\rm{{Ra}}}_{b}}^{0.72\pm 0.04}, so $u_{rms}$ approximately scales as {\color[rgb]{0,0,0}{\rm{{Ra}}}_{b}}^{0.36}. Combining these two estimates suggests that the turbulent temperature flux should scale as {\color[rgb]{0,0,0}{\rm{{Ra}}}_{b}}^{0.36}/\chi for low enough $\chi$. In Figure 15, we plot $\bar{F}_{h}\chi/$Ra${}_{b}^{0.36}$ versus $r$ for Model B runs at Ra${}_{o}=10^{7}$ and for four different values of $\chi$. The predicted scaling seems to work well for $\chi\leq 0.1$. We conclude that the emergence of a subadiabatic layer is a generic result which occurs for large bulk Rayleigh numbers and/or low values of $\chi$, as we have found in our simulations. Note that the scaling \bar{F}_{h}\sim{\color[rgb]{0,0,0}{\rm{{Ra}}}_{b}}^{0.36}/\chi suggests that the subadiabatic layer could in fact appear even when $\chi=1$ provided {\color[rgb]{0,0,0}{\rm{{Ra}}}_{b}} is large enough. In that sense, it might even be realized in the limit of a Cartesian RBC system as long as mixed thermal boundary conditions are used, though the Rayleigh number may need to be extremely large in that limit to exhibit the desired dynamics.

V A more solar $\beta(r)$ profile: Setup and numerical results

Until now we have used a profile for $\beta(r)$ dictated by the geometry and the boundary conditions of our model setup. To see whether our findings have any bearing on the dynamics of convection in stars, we now compute the equivalent $\beta(r)$ profile from a standard solar model (Model S, Christensen-Dalsgaard et al. (1996)). To do so, we evaluate the difference between $dT_{\rm{rad}}/dr=-3\kappa\rho L/(64\pi r^{2}\sigma T^{3})$ (where the Model S is used to extract the density $\rho$, the luminosity $L$, the temperature $T$, and the opacity $\kappa$, and where $\sigma$ is the Stefan-Boltzmann constant), and the adiabatic temperature gradient $dT_{\rm{ad}}/dr=-g/c_{p}$. The results are shown in Figure 16. We see that, by contrast with the models we have been using so far, $|\beta|$ decreases inwards instead of increasing inwards.

In the light of this information, we conduct a final set of numerical experiments where we construct a more solar-like profile for $\beta(r)$ choosing $\beta(r)=\chi/(1-\chi-(1/r^{2}))$ in order to ensure that $|\beta(r)|$ decreases inward, and letting $\alpha(r)/\alpha_{o}=-1/\beta(r)$ as before to have Ra$(r)=$Rao. Note that in this model, $\chi$ does not have the same physical meaning as in Equation (15), but it is merely used as a parameter to describe a family of functions $\beta(r)$ with a “solar-like” profile. Figure 17 illustrates the $\beta(r)$ functions thus created for two different values of $\chi$. Note that $\beta$ lies in the range (0,1] but crucially the ratio of the inner to outer values is large and approximately equal to $11$ for $\chi=0.1$ and $105$ for $\chi=0.01$.

We have run two simulations, for two different values of $\chi$ ($\chi=0.01$ and $\chi=0.1$) at Ra${}_{o}=10^{7}$. In the previous models, these cases led to the emergence of a subadiabatic region close to the outer boundary of the convection zone. Looking at the square of the non-dimensional buoyancy frequency profile for this model (Fig. 18), we observe that a slightly subadiabatic region does indeed appear, this time close to the inner boundary of the convection zone.

Note that the general mechanism for the appearance of this layer is the same as before, although the quantitative details differ. In this model setup, the mean kinetic energy is again controlled only by the bulk Rayleigh number (see Fig. 6), hence the velocity fluctuations remain large. However because $\beta(r)$ varies between [math] and $1$ the typical amplitude of the temperature perturbations $\Theta_{rms}$ is much smaller (i.e. this time $\Theta_{rms}\propto\chi$). This results in a much smaller total turbulent temperature flux $\bar{F}_{h}$ compared with Models A and B. As shown in Equation (25), however, whether a subadiabatic layer appears or not depends on the relative amplitude of $\bar{F}_{h}$ compared to $\beta(r)$. Since $\beta(r)$ is close to [math] near the inner boundary, a small turbulent temperature flux is indeed sufficient to create a subadiabatic layer there according to Eq. (24).

Figure 19 shows a snapshot of $u_{r}$ for the $\chi=0.01$ case, in which the subadiabatic region is the deepest observed so far. Notice that convection is still vigorous throughout, again supporting our conclusions from previous models that this type of convection is highly non-local.

VI Discussion

We have studied convection in a weakly compressible gaseous spherical shell, assuming a constant adiabatic temperature gradient as well as mixed temperature boundary conditions (fixed flux at the inner boundary and fixed temperature at the outer boundary). In Sections II, III and V, we presented results from three different model setups, that all have the same remarkable properties for sufficiently large Rayleigh number Ra, and sufficiently large variations in the superadiabaticity across the shell (measured by $\chi$). All these simulations showed substantial asymmetry between upflows and downflows, as well as the emergence of a subadiabatic layer which is still fully mixed by non-local convection. In Section IV, we explained these findings as follows:

Asymmetry between upflows and downflows: As in standard convection downward-traveling parcels are heated up by adiabatic compression, but remain cooler than the background temperature and therefore proceed to sink. In BRBC, these parcels would eventually have to warm up to match the temperature at the lower boundary, but in the case of fixed flux inner boundary condition, this is not the case and the parcels remain cooler than the surroundings as they reach the bottom of the convection zone. There they rebound from the boundary essentially being diverted (i.e. pressure-driven) into cool upflows. We have found that for large enough Rao, all upflows are pressure-driven. This asymmetric driving mechanism of upflows and downflows persists in much of the domain and proves that convection is highly non-local. These results are similar, though not identical, to those recently reported in Käpylä et al. (2017).

Emergence of a subadiabatic layer: The fixed flux boundary condition at the inner boundary has a second consequence, namely that of tying the turbulent temperature flux to the perturbed diffusive temperature flux. Hence, for sufficiently large turbulent temperature flux, the diffusive temperature flux must also become large and can cause the background temperature gradient to exceed the adiabatic one, which results in a subadiabatic stratification.

Two natural questions hence arise: What are the minimal necessary conditions for these dynamics to manifest themselves and why have these never been reported before in other Boussinesq studies? As an answer to the first question, we argue that the necessary conditions are (1) mixed temperature boundary conditions with fixed flux at one boundary and fixed temperature at the other, and (2) sufficiently turbulent flows. The importance of (1) should be clear from the description above. Condition (2) is necessary for the turbulent fluxes to be large enough and for inertia to be strong enough. These conditions are necessary, but do not have to be met necessarily in exactly the same way they are created as in Models A, B or C. For instance, we believe that with a sufficiently deep shell, it may be possible to achieve this dynamical regime even if $\chi$ were closer to one. Although not a strictly necessary condition, we have also found that having a radially-varying superadiabaticity $\beta(r)$ more easily creates a large enough contrast across the domain and therefore lowers the threshold in Ra for the emergence of the subadiabatic layer. As such, it might even be possible to be in this unusual regime in a strictly Boussinesq Cartesian RBC, though the Rayleigh number may need to be extremely high, or one may need to invoke additional non-Boussinesq effects to generate a rapidly-varying $\beta$ and the therefore significant $\beta$ contrast (such as varying gravity, for example, which causes a varying $dT_{\rm{ad}}/dz$, or by varying the diffusivities, or by adding internal heat sources within the fluid). By this reasoning, we conjecture that these dynamics may be found in high Pr number convection for large enough Ra.

Given these necessary conditions, we can now easily answer the second question. This kind of convection has not been previously observed in other Boussinesq studies because the vast majority of investigations to date have used fixed temperature boundary conditions or fixed flux at both boundaries (e.g. Otero et al. (2002); Johnston and Doering (2009)). There are certain studies in which mixed temperature boundary conditions have been implemented, notably in Gilman (1978b). There, low and intermediate Ra were investigated, with a Prandtl number equal to one and the flows were likely insufficiently turbulent (see condition 2) for the subadiabatic layer to appear. In Verzicco and Sreenivasan (2008) (see also Stevens et al. (2011)), turbulent convection in the high Ra regime using fixed flux at the bottom and fixed temperature at the top was also studied but there was no report of any subadiabatic layer. Interestingly however, they indeed found larger rms temperatures near the lower boundary, similar to our results. They also noticed that the plumes were less buoyant and cooler and as a result carried less heat compared with cases where the temperature was fixed at both boundaries. In both cases however, the fluid was incompressible with $dT_{\rm{ad}}/dz=0$ and $dT_{\rm{rad}}/dz$ constant, so that $\beta(z)=-1$ and there was much less imposed asymmetry in their system. Hence, although some prior studies have considered the effects of mixed temperature boundary conditions, ours appears to be the first to report the emergence of a subadiabatic layer. That implies that a combination of both conditions described above has to be satisfied.

Finally, it is important to recall that we have used the SVB approximation for weakly compressible gases even though our model setups do not necessarily satisfy all the requirements of this approximation. Indeed, the two fundamental assumptions entering the SVB approximation are:

1. that any perturbations of the thermodynamic quantities $\rho$ and $T$ about their domain mean $\rho_{m}$ and $T_{m}$ should be small, and 2) that the flow velocities be much smaller than the sound speed. These approximations then justify the use of equations (5) to (7). Note that SV also used a Cartesian domain, and further assumed, for simplicity, that $\alpha$, $g$, $\kappa$ and $\nu$, as well as the radiative and adiabatic temperature gradients were constant, but these are not necessary conditions for the applicability of their equations. However, in a spherical geometry, or if these quantities vary with depth (as in our own models), one should verify that both the background state and the perturbations continue to satisfy conditions (1) and (2) a posteriori. In Model A, as discussed in Section II-A, the sphericity of the domain implies that $T_{\rm{rad}}$ must vary with depth, even if everything else is held constant. As such, the SVB approximation can only be used if $\Delta T\ll T_{m}$, or equivalently, if $F_{\rm{rad}}(r_{i}^{2}/\kappa)(\frac{1}{r_{i}}-\frac{1}{r_{o}})\ll T_{m}$. To satisfy this condition and maintain a large Rayleigh number at the same time can be achieved by letting $\nu\rightarrow 0$ for instance. In Model B on the other hand, the validity of the SVB approximation is definitely violated, as $\alpha(r)$ varies quite significantly across the domain when $\chi$ is small. This invalidates the linearization of the equation of state used. A similar statement applies to Model C.

One may therefore rightfully question whether the new dynamics discovered here are artifacts associated with breaking the bounds of validity of the SVB equations, or whether they would indeed occur in a more realistic, fully compressible model setup as well. Based on published work and our own unpublished recent findings, the latter statement is likely true. Indeed, Chan and Gigas (1992) reported the emergence of a subadiabatic layer close to the outer boundary of 3D, large-eddy simulations of compressible convection in an effectively plane parallel domain. However the resolution used was very low, shedding doubt at the time on the robustness of their results. Other studies have also noted the appearance of a subadiabatic mixed region in compressible hydrodynamic simulations (e.g. Mocák et al. (2009); Herwig et al. (2006)) but the setup in these cases was not based purely on a convection zone, but rather, on a convection zone embedded in between two stable regions. Recently, Käpylä et al. (2017) also reported a subadiabatic layer in their fully compressible 3D simulations of overshooting convection. Nevertheless, a complete and definitive answer to this question requires the solution of the fully compressible equations, which is the subject of a future publication. Preliminary results obtained by the authors with 3D fully compressible DNSs in a Cartesian box with mixed temperature boundary conditions do indeed show the appearance of a subadiabatic region, suggesting that the most salient feature of this kind of convection is robust (Brummell et al. (2017), in prep.). We, however, expect the details to differ substantially, since the SVB approximation inherently suppresses some essential compressible dynamics, in particular the dynamic role of pressure in compressional heating.

ACKNOWLEDGMENTS

The authors would like to thank Céline Guervilly for her help with the PARODY code and Hugues Faller for his useful comment. This work was financially supported by NASA NNX14AG08G. The numerical simulations were performed on the Hyades supercomputer at UCSC purchased using an NSF MRI grant.