Compressible turbulent convection in highly stratified adiabatic background

Abstract

Buoyancy-driven turbulent convection leads to a fully compressible flow with a prominent top-down asymmetry of first- and second-order statistics when the adiabatic equilibrium profiles of temperature, density and pressure change very strongly across the convection layer. The growth of this asymmetry and the formation of an increasingly thicker stabilized sublayer with a slightly negative mean convective heat flux $J_c(z)$ at the top of the convection zone is reported here by a series of highly resolved three-dimensional direct numerical simulations beyond the Oberbeck–Boussinesq and anelastic limits for dimensionless dissipation numbers, $0.1 \le D\le 0.8$, at fixed Rayleigh number $Ra=10^6$ and superadiabaticity $\epsilon =0.1$. The highly stratified compressible convection regime appears for $D > D_{crit}\approx 0.65$, when density fluctuations collapse to those of pressure; it is characterized by an up to nearly 50 % reduced global turbulent heat transfer and a sparse network of focused thin and sheet-like thermal plumes falling through the top sublayer deep into the bulk.

1. Introduction

Solar convection close to the surface proceeds in the presence of an extremely stratified adiabatic background equilibrium which is characterized by scale heights, the distances across which temperature $T$, density $\rho$ or pressure $p$ drop by an order of magnitude, as small as $\ell \sim 10^3$ km amounting to only $5\ {\rm per}\ {\rm mile}$ of the total depth of the solar convection layer, $H\sim 2\times 10^5$ km (Christensen-Dalsgaard 2002; Nordlund, Stein & Asplund 2009; Schumacher & Sreenivasan 2020). It is one prominent example for a fully compressible convection flow which is highly asymmetric when up- and down-welling thermal plumes are compared. Driven by the strong radiative cooling at the top surface, cold plasma sinks into the interior ($u_z<0$, where $u_{z}$ is the vertical velocity) with nearly the speed of sound at the narrow edges of the granules, the convection cells in the solar case. Plasma rises moderately to the surface across the granule interior with a typical diameter of approximately $10^3$ km (Magic et al. 2013). It is yet open, how deep the narrow plumes reach into the highly stratified interior (Brandenburg 2016; Cossette & Rast 2016; Anders, Lecoanet & Brown 2019) and how much they get focused by compressibility effects to overcome turbulent dispersion. Answers to these questions would also provide insights on the vertical extent of supergranules, the larger-scale convection cells in the background with an extension of approximately 30 granule diameters (Rincon & Rieutord 2018; Hanson et al. 2020) and on the related solar conundrum of anomalously low velocity-fluctuation amplitudes from helioseismology in comparison with models (Hanasoge, Duvall Jr & Sreenivasan 2012). All this sets the physical motivation for the present study.

In this work, we want to study fully compressible convection in the presence of highly stratified adiabatic equilibrium profiles for $T$, $\rho$ and $p$. We thus analyse a series of direct numerical simulations (DNS) towards the high stratification limit of $D\to 1$ for the dissipation number $D$ (defined below). Our goal is to analyse the genuine compressibility effects isolated from complex multi-physics, as in the solar example above, and from temperature dependencies of dynamic viscosity and thermal conductivity. Our three-dimensional DNS beyond the Oberbeck–Boussinesq (Chillà & Schumacher 2012) and anelastic regimes (Ogura & Phillips 1962; Verhoeven, Wiesehöfer & Stellmach 2015; Alboussière & Ricard 2017; Jones, Mizerski & Kessar 2022) lead to strongly differing boundary-layer dynamics at the top and bottom, such as a sublayer formation at the top with a slightly negative mean convective heat flux, which causes a reduction of the global turbulent heat transfer by almost 50 % for $D=0.8$. In this layer, narrow focused thermal plumes form at the top and fall deeply into the highly stratified bulk of the layer. Our study thus demonstrates that some typical phenomena in natural convection flows can be understood in significantly simpler compressible flow configurations even though being far away from the natural flow case in terms of Rayleigh and Prandtl numbers.

2. Simulation model and parameters of compressible convection

The equations of motion for compressible convection are given by

(2.1)\begin{gather} \partial_t \rho + \partial_i (\rho u_{i}) = 0, \end{gather}

(2.2)\begin{gather}\partial_t (\rho u_{i}) + \partial_j (\rho u_{i} u_{j}) ={-}\partial_i p + \partial_j \sigma_{ij} - \rho g \delta_{i,3} , \end{gather}

(2.3)\begin{gather}\partial_t (\rho e) + \partial_j (\rho e u_{j} ) ={-}p \partial_i u_{i} + \partial_i(k \partial_i T) + \sigma_{ij} S_{ij}, \end{gather}

(2.4)\begin{gather}p = \rho R T \quad\text{where }\ R= C_{p} -C_{v} . \end{gather}

These equations correspond to mass, momentum and energy conservation laws along with the ideal gas equation of state. Here, $\rho$, $\rho u_{i}$, $p$, $\rho e$, $T$ are the mass density, momentum density components, pressure, internal energy density and temperature, respectively. The acceleration due to gravity is g, viscous stress tensor is ${\boldsymbol {\sigma }}=2\mu {\boldsymbol S}+ 2\mu {\boldsymbol I}({\boldsymbol {\nabla }}\boldsymbol {\cdot } {\boldsymbol u})/3$ with the Kronecker tensor ${\boldsymbol I}$ and the rate of strain tensor ${\boldsymbol S} = ({\boldsymbol {\nabla }} {\boldsymbol u}+{\boldsymbol {\nabla }} {\boldsymbol u}^{\rm T})/2$. Here, $\boldsymbol{u}$ stands for the velocity field in vector notation. The dynamic viscosity $\mu$ is assumed to be constant in these simulations. The thermal conductivity $k$ is related to the viscosity through the Prandtl number, $k= \mu C_{p}/Pr$. In our DNS, $Pr= 0.7$; $C_{p}$ and $C_{v}$ correspond to specific heat at constant pressure and volume, respectively. Their ratio, $\gamma = C_{p}/C_{v}= 1.4$ for a diatomic gas. The internal energy is defined as $e=C_{v}T$.

A uniform grid is used in $x$- and $y$-directions along with periodic boundary conditions. In wall-normal $z$-direction, a non-uniform grid with a point clustering near the walls is taken, which follows a hyperbolic tangent stretching function. Spatial derivatives are calculated by a 6th-order compact scheme for all points except near the walls (Lele 1992; Baranwal, Donzis & Bowersox 2022); 4th- and 3rd-order compact schemes are used at the last two grid points near the wall. No-slip, isothermal boundary conditions are applied at the top and bottom. The boundary condition for $p$ is evaluated using the $z$-component of the momentum equation at $z=0, H$, $\partial p/\partial z = \left( \partial \sigma_{iz}/\partial x_{i} \right) - \rho g$. The fields are advanced in time by a low-storage 3rd-order Runge–Kutta with a Courant number of ${\rm CFL}=0.5$.

Incompressible Rayleigh–Bénard convection is determined by the Rayleigh and Prandtl numbers, $Ra$ and $Pr$. In compressible convection, additional dimensionless parameters need to be introduced. The first one is the dissipation number $D$ (Verhoeven et al. 2015; Jones et al. 2022) which is given by

(2.5)\begin{equation} D= \frac{gH}{C_{p}T^{B}}. \end{equation}

The parameter $D$ can be considered as a rescaled dry adiabatic lapse rate as it determines differences in the thermodynamic properties across the layer depth in a purely isentropic process. The adiabatic temperature ratio between top (T) and bottom (B) is $T^{T}_{a} /T^{B}_{a}=( 1 - D)$, the corresponding density ratio follows to $\rho ^{T}_{a} /\rho ^{B}_{a} = [ T^{T}_{a} /T^{B}_{a} ]^{1/ ( \gamma -1 ) }$. The subscript ‘$a$’ corresponds to the adiabatic equilibrium without convection. The equilibrium profiles for all state variables and all values of $D$ are shown in figure 1. In order to initiate convective motion, the actual temperature gradient across the layer must be greater than that due to a purely adiabatic process (Jeffreys 1930), i.e. $T^{B} - T^{T} > T_{a}^{B} - T_{a}^{T}$. In our study, the bottom plate is taken as the reference, thus $T^{B}= T_{a}^{B}$ is a constant for all cases; the instability criterion reduces to $T^{T} < T_{a}^{T}$.

Figure 1. Adiabatic equilibrium profiles of the thermodynamic state variables. (a) Temperature $T$, (b) density $\rho$ and (c) pressure $p$. The difference between top and bottom values for all three variables increases with $D$, see table 1 for the colour code description.

The second additional compressible convection parameter is the superadiabaticity $\epsilon$ (Verhoeven et al. 2015; Jones et al. 2022) , which is given by

(2.6)\begin{equation} \epsilon= \frac{T_{a}^{T} - T^{T}}{T^{B}}. \end{equation}

Superadiabaticity is the excess relative temperature gradient with respect to the adiabatic equilibrium state, or, in other words, a measure of the driving of convection. If one takes the characteristic free-fall velocity, $U_{f}= \sqrt {\epsilon g H}$, then a free-fall Mach number,

(2.7)\begin{equation} M_{f}(z) = \frac{U_{f}}{\sqrt{\gamma R \langle T(z) \rangle_{A,t}}} = \sqrt{\frac{\epsilon D}{\gamma - 1}} \sqrt{\frac{T^{B}}{\langle T(z)\rangle_{A,t}}} , \end{equation}

follows. The notation $\langle \boldsymbol {\cdot }\rangle _{A,t}$ corresponds to a combined average over the horizontal directions and time which is given for a field $X$ as

(2.8)\begin{equation} \langle X (z) \rangle_{A,t} = \frac{ \displaystyle\sum^{N_{t}}_{m=1} \displaystyle\sum^{N_{x}}_{i=1} \sum^{N_{y}}_{j=1} X(x_i, y_j,z,t_m)}{N_{x}N_{y}N_{t}}, \end{equation}

where $N_{x}$, $N_{y}$ and $N_{t}$ are the grid point numbers with respect to $x$-, $y$-directions and number of snapshots, respectively. The strength of compressibility is thus a function of both $\epsilon$ and $D$. From (2.7), it is expected that compressibility at the top is higher than at the bottom. For $\epsilon \rightarrow 0$, we get the anelastic approximation which includes first-order terms in $\epsilon$ only (Ogura & Phillips 1962; Gough 1969; Lantz & Fan 1999). This results in a filtering of high-frequency acoustic waves. Rayleigh–Bénard convection in the Oberbeck–Boussinesq (OB) limit follows for a subsequent $D \rightarrow 0$ with $D$ decreasing faster than $\epsilon$. On the basis of the new parameters, we will use $Ra= \epsilon g H^{3}/(\nu ^{B} \kappa ^{B})$ and Prandtl number, $Pr= \nu ^B/\kappa ^B$, where $\kappa ^B= k/(C_{p} \rho _a^B )$ is the temperature diffusivity and $\nu ^B=\mu /\rho _a^B$ the kinematic viscosity. The Rayleigh number is here defined with respect to the bottom variables and thus constant for all cases. Unlike $Ra$, which varies across depth due to density dependence, Prandtl number $Pr$ is a constant in the entire domain, since the density drops out for $\nu$ and $\kappa$. In this work, we focus on a constant and low but finite superadiabaticity, $\epsilon =0.1$ for a range of $D$ from $0.1$ to $0.8$. All other governing parameters such as $Ra \approx 10^6$, $Pr= 0.72$ and aspect ratio, $\varGamma = 4$ are kept constant, see table 1. In all figures that follow, unless specified otherwise, we use the plotting style given in the last column. From table 1, we observe that in our study, the strength of density stratification $(\rho _{a}^{B}/\rho _{a}^{T})$ ranges from 1.3 at $D= 0.1$ to 55.9 at $D=0.8$. The degree of stratification in previous compressible convection simulations is thus significantly smaller than what we will study here. The strongest stratification ratios in Verhoeven et al. (2015), Jones et al. (2022), and Tilgner (2011) are 20.1, 10 and 4.6, respectively.

Table 1. List of the DNS. All cases have $Pr= 0.7$, $Ra \approx 10^6$, $\epsilon = 0.1$, $\gamma = C_{p}/C_{v}= 1.4$, and an aspect ratio, $\varGamma =L/H=4$ resolved by $N_x:N_y:N_z=512: 512: 256$ grid points. Here, L is the horizontal domain size. We list dissipation number $D$, bottom-to-top ratios of adiabatic temperatures and densities, the global maximum turbulent Mach number $M_t^{max}$, the Nusselt and Reynolds numbers, $Nu$ and $Re$. The error bars of $Nu$ and $Re$ are less than $2\,\%$ for all cases. Furthermore, we display boundary-layer scales for density, convective heat flux and density–temperature correlations. In all figures that follow, coloured lines are given as in this table.

3. Growing asymmetry with increased stratification

The adiabatic temperature profile, $T_{a}(z)=T^B_a(1-Dz)$, is the equilibrium profile. Thus the part of temperature corresponding to convection is the superadiabatic temperature defined as

(3.1)\begin{equation} T_{sa} ({\boldsymbol x},t)= T({\boldsymbol x},t) - T_{a}(z). \end{equation}

In figure 2, we compare the contour visualization of $\ln (| \boldsymbol {\nabla } T_{sa}| )$ for the two extreme cases in our series with $D= 0.1$ (low) and $D= 0.8$ (high). Panel (a) corresponds to case 1, closest to the OB limit. Not surprisingly, the plumes from the top and bottom are more or less symmetric in shape and frequency. However, this top-bottom symmetry disappears completely for the highly stratified case 8 in panel (b). One observes a network of slender thin plumes (Rast 1998) originating from the top boundary layer that fall deep into the bulk. Recently, John & Schumacher (2023) showed that thinner plumes from the top boundary are a common characteristic of compressible convection regardless of the relative boundary-layer thickness, but here the asymmetry is very pronounced.

Figure 2. Visualization of the plume structure of the superadiabatic temperature field $T_{sa}$. Contours of $\ln |{\boldsymbol \nabla } T_{sa} |$ are shown at two instants, (a) for case 1 with $D=0.1$ and (b) for case 8 with $D=0.8$. Minimum/maximum contour levels correspond to $-7.0/-2.4$ in (a) and $-9.5/-4.0$ in (b). The top contour surface is slightly below $z=H$.

Similar asymmetry is observed for the mean profile of the superadiabatic temperature, $\langle T_{sa} (z) \rangle _{A,t}$. In figure 3(a), we plot the normalized superadiabatic temperature, $\langle T_{sa}(z)\rangle _{A,t}/ \varDelta \langle T_{sa} \rangle _{A,t}$ as a function of depth $z$ where $\varDelta \langle T_{sa} \rangle _{A,t}=\langle T_{sa}(z=0) \rangle _{A,t} - \langle T_{sa}(z=H) \rangle _{A,t}$. Clearly one observes that the asymmetry, i.e. the offset from 0.5 and thickness of the top boundary layer in comparison with the bottom one, increases considerably with growing dissipation number. More detailed, as $D$ increases, the bulk temperature is increasingly closer to the bottom temperature, consistent with Tilgner (2011), Verhoeven et al. (2015), Jones et al. (2022) or Curbelo et al. (2019) for $Pr\to \infty$. Differently, the experiments of Wu & Libchaber (1991) showed that the bulk temperature is closer to the prescribed top plate value with $(\nu ^{T}, \kappa ^{T} ) < ( \nu ^{B}, \kappa ^{B} )$. This results in a thinner top boundary layer and a bulk temperature closer to that at the top. Due to strong density stratification, we get $(\nu ^{T}, \kappa ^{T} ) > ( \nu ^{B}, \kappa ^{B} )$ and thus thicker layers at the top, which is detailed later. A trend similar to Wu & Libchaber (1991) is observed for $\epsilon / D \gg 1$ in John & Schumacher (2023).

Figure 3. Asymmetry of the (a) normalized superadiabatic temperature and (b) the turbulent Mach number $M_t$. The inset in (a) corresponds to mean thermal boundary-layer thickness $\lambda ^{B,T}(D)$. Symbols $+$ and ${\rm o}$ are for the top and bottom boundary-layer thickness, respectively. The inset in (b) shows $M_f(D)$ from (2.7). The red, blue and green lines are for $z/H = 0, 1$ and 0.5, respectively. The black squares are the maximum Mach numbers as given table 1, which are observed in the simulations over the whole space and time. The legend in (b) corresponds to the inset figure.

We can quantify this thickness as the mean of a local thermal boundary-layer thickness (Scheel & Schumacher 2014) which is given by

(3.2) \begin{equation} \tilde{\lambda}^{B,T}(x,y)= \frac{1}{H}\left | \langle T_{sa} \left( 0,H \right) - T_{sa} \left( H/2 \right) \rangle_{A,t} \right| \left | \frac{\partial {{T_{sa}} _{A,t}} }{\partial z} \right|^{{-}1} _{z=0,H}. \end{equation}

We plot $\lambda ^{B,T}=\langle \tilde {\lambda }^{B,T}\rangle _{A,t}$ as a function of $D$ in the inset of figure 3(a). The difference between the boundary-layer thickness at the bottom and top plate increases with $D$ in agreement with the differences observed in the distributions of the instantaneous, local thermal boundary-layer thicknesses (John & Schumacher 2023). The top boundary-layer thickness increases drastically with $D$ whereas the bottom boundary-layer thickness is almost invariant. This, along with increased asymmetry in the bulk towards the bottom plate, as seen in the normalized $T_{sa}$-profiles, suggests that at high $D$ the mixing is not as efficient as in OB and compressibility effects are mainly localized near the top boundary.

An asymmetry similar to that of $\langle T_{sa}(z)\rangle _{A,t}$ is observed for the turbulent Mach number too with significant localized compressibility seen near the top boundary. The turbulent Mach number is defined as

(3.3)\begin{equation} M_{t}(z) = \frac{u_{ rms}(z)}{\sqrt{\gamma R \langle T(z) \rangle_{A,t}}}\quad \mbox{with} \quad u_{ rms}(z)= \sqrt{ \sum_{i=1}^3 \left [\langle u_{i}^{2}(z)\rangle_{A,t} - \langle u_{i}(z) \rangle^{2}_{A,t}\right ]}. \end{equation}

The turbulent Mach number at the top increases monotonically with $D$, whereas the behaviour at the bottom boundary is non-monotonic with $D$. At the bottom, $M_{t}$ increases up to $D= 0.5$ and then decreases with $D$ as a result of the increased asymmetry between the boundary layers, see figure 3(b).

In the inset of figure 3(b), we compare the resulting maximum turbulent Mach number in the flow with the estimate from (2.7). The black squares are the maximum local turbulent Mach number for each $D$ encompassing all time steps and spatial locations. The red and blue lines correspond to $M_{f}$ estimates using the bottom and top temperatures, $T^B$ and $T^{\rm T}$, respectively. One finds that the resulting maxima of $M_t$ are between these estimates. This is due to no-slip boundary conditions at the wall; thus, the velocity fluctuations peak only at the outer region of the boundary layer and in the bulk. Since the (mean) bulk temperature is closer to that at the bottom plate, the values of $M_t$ are closer to $M_{f} (T^{B} )$. A really close agreement for the simulations is seen when the Mach number is estimated using the bulk temperature, green line in the inset figure. This suggests that a simple estimate such as in (2.7) can give us a good idea about the level of compressibility in the flow.

We would like to point out that in compressible turbulence $M_{t}$ is not sufficient to characterize the system; an additional parameter $\delta$, the ratio of dilatational to solenoidal root-mean-square velocities, needs to be included as discussed in Donzis & John (2020), John, Donzis & Sreenivasan (2019, 2020, 2021). This is not surprising as in compressible convection two extremely different conditions, $\epsilon /D \ll 1$ and $\epsilon /D \gg 1$, can lead to the same $M_t$ despite different dynamics. John & Schumacher (2023) showed that an increasing $\epsilon$ leads to enhanced heat transfer, whereas it is decreased for increasing $D$ at low $\epsilon$. Even though $M_t$ would be comparable for both cases, $\delta$ will differ in these limits. A similar conclusion, that the Mach number along with a measure of dilatation in the flow field is required for the statistical analysis, was reported for compressible turbulent channel flows (Baranwal, Donzis & Bowersox 2023a,b).

4. Fluctuations of state variables and regime transition

A major consequence of compressibility are the fluctuations of $p$, $\rho$ and $T$. In figure 4, we compare these fluctuations as a function of depth for the two extreme cases of our DNS series, $D= 0.1$ and $0.8$. John & Schumacher (2023) showed that $D\le 1-\epsilon$, which would result in a maximally possible theoretical value of $D=0.9$ in our DNS series. The relative fluctuations are defined as $X_{ rms}(z)/\langle X(z) \rangle _{A,t}$ with $X_{ rms} = (\langle X^2 \rangle _{A,t} - \langle X \rangle _{A,t}^2)^{1/2}$. Panel (a) at $D=0.1$ shows that pressure fluctuations are negligible. The magnitude of the density and temperature fluctuations are of the same order of magnitude and the profiles of both quantities are practically identical. These findings are consistent with the OB limit at $\epsilon, D \rightarrow 0$.

Figure 4. Relative fluctuations of thermodynamic quantities for (a) $D = 0.1$ and (b) $D= 0.8$. (c) Maximum of the relative thermodynamic fluctuations vs dissipation number, $D$. For all figures, the red, blue and green correspond to temperature, density and pressure, respectively. Inset of (c) shows four different distances $\lambda _{X}$ as a function of $D$. For the definition of the specific $\lambda _{X}$, we refer to the text in § 4 and the legend of the inset.

In figure 4(b), the relative fluctuations of $p$, $\rho$ and $T$ at $D= 0.8$ are shown. As a result of high compressibility at $D=0.8$ at the top, see figure 3(b) for $M_t(z)$, the thermodynamic fluctuations are considerably higher compared with the bottom one. Close to the top boundary, the pressure and density fluctuations dominate the temperature fluctuations. This is due to the constant temperature imposed at the boundaries which results in an isothermal condition. For an assumed ideal gas, $p\propto \rho$. Indeed this is evident for $D= 0.8$ at the top boundary as the relative fluctuations of density and pressure perfectly match at the top boundary. This would be also true at the bottom boundary and even for $D=0.1$, but it is not evident due to low $M_{t}$. Compressibility effects imply in the following discussion that we face a large magnitude of pressure fluctuations.

As we move away from the wall, fluctuations of both $p$ and $\rho$ decrease, but are not perfectly correlated anymore due to the growing temperature fluctuations. After a critical distance from the wall, pressure fluctuations continue to decrease; the density fluctuations reach a minimum at $z=z_{\ast }$. This local minimum of $\rho _{ rms}(z)/\langle \rho (z)\rangle _{A,t}$ exists when pressure and temperature fluctuations become comparable in magnitude. For the highest dissipation number $D= 0.8$, density fluctuations at $z< z_{\ast }$ are driven by temperature fluctuations in contrast to $z>z_{\ast }$. One also notes that the location of the secondary maximum of $\rho _{ rms}(z)/\langle \rho (z)\rangle _{A,t}$ does not match the local maxima of $T_{ rms}(z)/\langle T(z)\rangle _{A,t}$. This implies that pressure fluctuations, although not dominant, remain significant and thus compressibility effects too (Jones et al. 2022).

These features are absent for $D= 0.1$ in figure 4(a). Thus, there seems to exist a critical $D$ where a transition to the compressibility dominated convection regime exists that is qualitatively different to the OB or anelastic limits. In order to identify the critical $D$, we plot the maximum relative thermodynamic fluctuations with respect to $D$ in figure 4(c). We find a $D_{ crit}\approx 0.65$ for the present $(Ra,\epsilon )$. For $D < D_{ crit}$, maxima of density and temperature fluctuations dominate over those of the pressure fluctuations. As $D$ increases, the maximum pressure fluctuations start to increase as well which is connected with increasing differences in the maxima of density and temperature fluctuations. For $D \ge D_{ crit}$, the maxima of pressure and density fluctuations collapse and are larger than those for the temperature fluctuations.

Clearly, there are multiple physical processes and thus scales at play near the top boundary at high $D$. This seems to be similar to the turbulent boundary layer with its inner and outer scales (Smits, McKeon & Marusic 2011). The behaviour of density can be analysed to demarcate various regimes inside the boundary layers in compressible convection. We thus introduce the following boundary layer scales: (1) $\lambda _{T_{ max}}$ as the distance between the top wall and the location of the local maximum of the relative temperature fluctuations; (2) $\lambda _{\rho _{ max}}$ as the distance between the top wall and the maximum of the relative density fluctuations. Note that in OB and even anelastic convection, $\lambda _{T_{ max}} = \lambda _{\rho _{ max}}$ (Jones et al. 2022). Furthermore, (3) $\lambda _{\rho _{ min}}$ as the distance between the top wall and the location of the local minimum of the relative density fluctuations. This distance corresponds to the location where we start to observe significant compressibility effects. Finally, (4) $\lambda _{\rho _{{ max},T}}$ is the distance between the top wall and the local maximum of the relative density fluctuations, corresponding to the maximum of local temperature fluctuations.

Next, we display all four scales $\lambda _X$ vs $D$ in the inset of figure 4(c). Here $\lambda _{T_{ max}}$ increases with $D$ which is consistent with the growth of the thermal boundary-layer thickness at the top, $\lambda ^{T}$, see figure 3(a). We see a sudden transition of $\lambda _{\rho _{ max}}$ at $D_{crit}= 0.65$ indicating that the maximum for $D>D_{ crit}$ is at $z=H$. Note also that for $D < 0.65$, $\lambda _{\rho _{ max}} \approx \lambda _{\rho _{{ max}, T}}$. Both scales increase with $D$, but start to deviate from $\lambda _{T_{ max}}$ at $D < D_{ crit}= 0.65$. This deviation starts for $D \ge 0.5$; it suggests that pressure fluctuations are important already, but not as dominant as discussed before. It might thus indicate that the anelastic approximation is strictly speaking not applicable even for $D \ge 0.5$. Finally, $\lambda _{\rho _{ min}}$ starts from nearly zero at $D = 0.1$ with negligible compressibility. With increasing dissipation number, $\lambda _{\rho _{ min}}$ increases. This is a new top boundary-layer regime induced by the increasing stratification. It is in line with an increase of $\lambda _{T_{ max}}$ and $\lambda ^{T}$; see again the inset of figure 3(a). We discuss the implications of this sublayer for the turbulent heat transfer in the next section.

5. Stratification impact on turbulent heat transfer

A systematic increase of the level of stratification, controlled by $D$, at relatively small superadibaticity of $\epsilon = 0.1$ has a significant impact on the normalized stress profiles as seen in figure 5. We show the variation of the normalized convective heat flux profiles $J_c(z)$ in panel (a) and the normalized density–temperature correlation profiles $J_{\rho }(z)$ in panel (b) for all eight runs. They are qualitatively similar, but sign-flipped behaviour is observed for both second-order correlations in panels (a,b). An increasing sublayer with an average ‘anti-convection’-behaviour, i.e. negative convective heat flux and positive density–temperature correlation can be detected. Note that $J_{\rho }(z)$ is only present in fully compressible convection and neglected in the anelastic case. Exactly, this term becomes significant for the highest $D$ and is expected to further increase when $D\to 0.9$ (which is computationally very hard). For $z\lesssim 0.9$ ($D= 0.8$) both mean profiles correspond to unstable convection.

Figure 5. Profiles of plane-time averaged stresses vs $D$. (a) Normalized convective heat flux $J_{c}(z)= \langle u_z T^{\prime }(z)\rangle _{A,t}/[u_{z,{rms}}(z)\langle T(z)\rangle _{A,t}]$. (b) Normalized density–temperature correlation $J_{\rho }(z)=\langle \rho ^{\prime } T^{\prime }(z)\rangle _{A,t}/[\langle \rho (z)\rangle _{A,t}\langle T(z)\rangle _{A,t}]$. (c) Normalized density averaged convective heat flux where $J^{\rho }_{c}(z)=\langle \rho u_z T^{\prime }(z)\rangle _{A,t}$ and (d) Nusselt number $Nu$ vs dissipation number $D$, see (5.1). Inset of (d) displays $1/Nu$ vs $\lambda _{\rho _{min}}$. The dashed line corresponds to a fit $A + B \lambda _{\rho _{ min}}$ with coefficients $A=0.12$ and $B=1.05$. Line colouring corresponds to table 1. Note that $T^{\prime }$ follows by a standard Reynolds decomposition, $T({\boldsymbol x},t)=\langle T(z)\rangle _{A,t}+T'({\boldsymbol x},t)$.

Next, we relate the depth of these stabilized regions to $\lambda _{\rho _{ min}}$ and thus to the regime transition at $D=D_{ crit}$. We therefore define three further distances or scales: (1) $\delta _{u_{z}T}$ as the width of negative normalized heat flux with respect to the top boundary; (2) $\lambda ^{max}_{\rho 'T'}$ as the distance from the top boundary to the location of the positive maximum of the normalized density–temperature correlation; (3) $\delta _{\rho 'T'}$ as the width of the region of positive density–temperature correlations.

These length scales along with $\lambda _{\rho _{ min}}$ are listed in table 1 for all DNS. We observe that for all runs, except for $D= 0.1$, the region of positive correlation exactly matches with the compressibility dominated region, i.e. $\delta _{\rho 'T'}=\lambda _{\rho _{ min}}$. Thus, there is also a direct proportionality between regions dominated by compressibility and anti-convection with $\delta _{{u_z}T} < \lambda _{\rho _{ min}}$. For $D \le 0.70$, we also find that $\delta _{{u_z}T} \approx \lambda ^{max}_{\rho 'T'}$, the location of the maximum positive correlation. For $D > 0.7$, we get $\delta _{{u_z}T} > \lambda ^{max}_{\rho 'T'}$, suggesting an increasing influence of compressibility. Furthermore, $\delta _{{u_z}T} < \lambda _{\rho _{ min}}$ implies that on average weak convective motions can occur at heights, $\lambda _{\rho _{ min}}\gtrsim z\gtrsim \delta _{{u_z}T}$, even though there are already significant compressibility effects. We recall for this discussion that even for the highest $D$, we still observe focused cold thermal plumes which detach from the top thus indicating locally positive convective flux. All these processes take places in the region where density fluctuations are tightly correlated to pressure fluctuations as we discussed in figure 4(c). We also report that there is no trend of the location of maximum turbulent Mach number, $M^{\ast }_{t}=\max _z M_t(z)$, see figure 3(b), with $D$ even though $M_{t}^{\ast }$ increases with $D$. For all the cases considered, the distance of this maximum to the top wall is $H-z_{\ast } \in [0.16, 0.174]$, significantly larger than $\lambda _{\rho _{ min}}$. This again suggests the inadequacy of $M_{t}$ alone to characterize general compressibility conditions (Donzis & John 2020). Particularly for the strongly asymmetric cases at highest $D$, the behaviour of compressible convection is completely different for same $M_{t}$ magnitude at either side of $z_{\ast }$.

In figure 5(c) we also plot the plane-averaged heat flux, $J^{\rho }_{c}(z)=\langle \rho ' u_z T'\rangle _{A,t}$, normalized by the maximum value. We observe a region of practically zero convective heat flux near the top that increases with $D$. For $z>H-\delta _{u_{z}T}$, the flux $J^{\rho }_{c}$ is less than 1 % of its maximum. There are other ways to define convective flux in terms of entropy or enthalpy (Currie & Browning 2017). The entropy changes can written as $\rho T dS = \rho C_{p} \,{\rm d}T - {\rm d}P$. The pressure term is neglected when boundary layers are incompressible. At high $D$, the boundary layer near the top is, however, strongly compressible. We consider the contributions due to pressure changes as a genuine compressibility effect; they become important and we analysed them in this work separated from the convective flux.

Finally, we plot the Nusselt number defined as

(5.1)\begin{equation} Nu = \frac{1}{2} \left[\left . - \frac{H}{\epsilon T_{B}} \frac{{\rm d} \langle T_{sa} \rangle _{A,t}}{{\rm d}z} \right|_{z=0} \left . - \frac{H}{\epsilon T_{B}} \frac{{\rm d} \langle T_{sa} \rangle _{A,t}}{{\rm d}z} \right|_{z=H}\right]. \end{equation}

It is shown as a function of $D$ in figure 5(d). The Nusselt number is calculated as the average of the heat flux at the top and bottom. Since the flows are statistically steady, the Nusselt number at the top and bottom are the same and the maximum mean Nusselt number difference between the plate is found to be less than $3\,\%$ for all cases. Clearly visible is that the Nusselt number $Nu$ monotonically decreases with $D$. At critical $D= 0.65$, there seems to be a slight change in behaviour but it is not significant. A decrease of $Nu$ with $D$ is again consistent with previous studies (Verhoeven et al. 2015; Jones et al. 2022). We plot the inverse of the Nusselt number vs $\lambda _{\rho _{ min}}$ in the inset of panel (d) and we observe that the trend can be reasonably approximated by a linear fit. Thus, we conclusively demonstrate that an increasingly inefficient heat transfer is connected to structural changes in the top boundary layer.

6. Final discussion

The objective of the present study was to take the least complex, but fully compressible turbulent convection system to study the genuine effects of compressibility isolated from other major contributions to a non-Boussinesq behaviour, such as temperature-dependent material properties, and to go beyond the OB and anelastic limits. We took a configuration in the limit of strongly decreasing mean profiles of the thermodynamic state variables – the high stratification regime. It corresponds to the dissipation number of $D\to 1$ for small superadiabaticity $\epsilon$. One prominent example for such a configuration might be the convection dynamics close to the solar surface. Even though the aim of our work is not to provide a realistic modelling of particularly this complex multi-physics flow, that would involve the coupling to magnetic fields and radiation transfer, the present simulations provide an insight into generic physical mechanisms which lead to highly asymmetric thermal plume motions, the fundamental local process in any convection flow that triggers fluid turbulence and transports heat. The limit cases of our simplified configuration lead already to a sparse network of thin plume structures which resembles similarities to the granule network at the Sun; see again figure 2(b). These coherent structures ‘leak’ through the stratified top region and thus maintain a significantly reduced global turbulent heat transport. Future studies will aim at an increase of the aspect ratio and Rayleigh number together with a decrease of the Prandtl number and $T$-dependent material parameters.