Extended subadiabatic layer in simulations of overshooting convection

TL;DR

This study uses numerical simulations to reveal that a significant subadiabatic layer can form in overshooting convection zones due to nonlocal energy transport, challenging traditional flux-gradient assumptions.

Contribution

It demonstrates the formation of an extended subadiabatic layer in overshooting convection driven by nonlocal buoyant downflows, with implications for stellar convection modeling.

Findings

Subadiabatic layers form when heat conduction varies smoothly.

Enthalpy flux includes a non-gradient component.

Flow topology shifts with subadiabatic layer presence.

Abstract

We present numerical simulations of hydrodynamic overshooting convection in local Cartesian domains. We find that a substantial fraction of the lower part of the convection zone (CZ) is stably stratified according to the Schwarzschild criterion while the enthalpy flux is outward directed. This occurs when the heat conduction profile at the bottom of the CZ is smoothly varying, based either on a Kramers-like opacity prescription as a function of temperature and density or a static profile of a similar shape. We show that the subadiabatic layer arises due to nonlocal energy transport by buoyantly driven downflows in the upper parts of the CZ. Analysis of the force balance of the upflows and downflows confirms that convection is driven by cooling at the surface. We find that the commonly used prescription for the convective enthalpy flux being proportional to the negative entropy gradient does not hold in the stably stratified layers where the flux is positive. We demonstrate the existence of a non-gradient contribution to the enthalpy flux, which is estimated to be important throughout the convective layer. A quantitative analysis of downflows indicates a transition from a tree-like structure where smaller downdrafts merge into larger ones in the upper parts to a structure in the deeper parts where a height-independent number of strong downdrafts persist. This change of flow topology occurs when a substantial subadiabatic layer is present in the lower part of the CZ.

Tables

Table 1. Table 1 Definition of the Zones.

${\overline{F}}_{\mathrm{enth}}$	$\nabla -{\nabla }_{\mathrm{ad}}$	Zone	Label	Lower Limit	Thickness
$>0$	$>0$	Buoyancy	BZ	$z_{\mathrm{BZ}}$	$d_{\mathrm{BZ}}$
$>0$		Deardorff	DZ	$z_{\mathrm{DZ}}$	$d_{\mathrm{DZ}}$
		Overshoot	OZ	$z_{\mathrm{OZ}}$	$d_{\mathrm{OZ}}$
$\approx 0$		Radiative	RZ	$\cdots$	$\cdots$

Table 2. Table 2 Summary of the runs, All with 288 3 superscript 288 3 288^{3} meshpoints.

Run	$K$	$\mathrm{Re}$	$z_{\mathrm{BZ}}$	$z_{\mathrm{DZ}}$	$z_{\mathrm{OZ}}$	$d_{\mathrm{BZ}}$	$d_{\mathrm{DZ}}$	$d_{\mathrm{OZ}}$
K	Kramers	$27$	$0.26$	$0.00$	$-0.27$	$0.74$	$0.26$	$0.27$
P	profile	$25$	$0.34$	$0.10$	$-0.19$	$0.66$	$0.24$	$0.29$
S	step	$26$	$0.02$	$0.00$	$-0.28$	$0.98$	$0.02$	$0.28$

Equations

$$\frac{D\ln \rho }{Dt}$$

$$\frac{D\mathbf{u}}{Dt}$$

$$T\frac{Ds}{Dt}$$

$$K(\rho ,T)=K_0(\rho /{\rho }_0{)}^{-(a+1)}(T/T_0{)}^{3-b}.$$

$${\mathbf{F}}_{\mathrm{SGS}}=-\rho T\left({\chi }_{\mathrm{SGS}}^{(0)}\nabla \overline{s}+{\chi }_{\mathrm{SGS}}^{(1)}\nabla s^{\prime }\right),$$

$${\overline{F}}_{\mathrm{rad}}=-\overline{K}{\partial }_z\overline{T},$$

$${\overline{F}}_{\mathrm{enth}}=c_{\mathrm{P}}\overline{(\rho u_z{)}^{\prime }{T}^{\prime }},$$

$$\overline{F}_{\rm enth}^{\rm MF}=\tau_{\rm red}\,\overline{\rho}\,\overline{T}\big{(}g\,\overline{{s^{\prime}}^{2}}/c_{\rm P}-\overline{u_{z}^{2}}\,\partial_{z}\overline{s}\big{)}\equiv\overline{F}_{\rm D}+\overline{F}_{\rm G},$$

Full text

1. Introduction

Convection plays a vital role in stellar activity by generating turbulence that, together with the star’s overall rotation, leads to differential rotation (e.g. Rüdiger 1989) and dynamo action (e.g. Krause & Rädler 1980). Energy transport due to convection is important for almost all stars during some stages of their evolution (e.g. Kippenhahn et al. 2012). Hence, a proper parameterization of convection is crucial for stellar structure and evolution in one-dimensional models.

Mixing length theory (MLT) continues to be a popular description of stellar convection. The formulation of MLT, as it is used today, goes back to the seminal work of Vitense (1953), where the properties of convection are related to the local value of the superadiabatic gradient $\nabla-\nabla_{\rm ad}$, with $\nabla=d\ln\overline{T}/d\ln\overline{p}$ and $\nabla_{\rm ad}$ being the actual and adiabatic logarithmic temperature gradients, respectively. Here, $T$ and $p$ are temperature and pressure, respectively, while overbars denote horizontal averaging. Convection is supposed to occur only if the horizontally averaged temperature stratification is superadiabatic, $\nabla>\nabla_{\rm ad}$, which is equivalent to the Schwarzschild criterion for convection, $d\,\overline{s}/dz<0$, where $s$ is the specific entropy.

The MLT has deeply influenced the way in which three-dimensional ab initio convection models are constructed: often a fixed profile of radiative heat conductivity $K$ is chosen, producing a superadiabatic layer of fixed depth (e.g. Hurlburt et al. 1986). Alternatively, the static thermodynamic background state is taken from an MLT-based stellar evolution model (e.g. Brun et al. 2011; Kitiashvili et al. 2016). In such setups, convection is driven at the largest scale available. This is a possible cause for the discrepancy in the convective velocities at large horizontal scales between simulations and time–distance helioseismology (Hanasoge et al. 2012). Smaller length scales could instead be imprinted by convection driven solely by the surface layers (Cossette & Rast 2016), which leads to a topology change of the downdrafts from a tree-like structure near the surface to strong plumes penetrating into deeper layers as cool entropy rain; see Spruit (1997). The latter can take part in the convective flux in these layers through a non-gradient contribution known as Deardorff flux (Deardorff 1966; Brandenburg 2016).

We present simulations in which we use either a physically motivated heat conduction formulation based on a Kramers-like opacity (Brandenburg et al. 2000) or two types of fixed heat conductivity profiles to study their effect on the structure of the convection zone (CZ). Furthermore, we demonstrate the existence of a non-gradient contribution to the enthalpy flux and use a quantitative analysis to study the topology change of the downflow and upflow structures in the simulations.

2. The Model

2.1. Basic equations

We solve the equations of compressible hydrodynamics:

[TABLE]

[TABLE]

where $D/Dt=\partial/\partial t+{\bm{u}}\cdot\bm{\nabla}$ is the advective derivative, $\rho$ is the density, ${\bm{u}}$ is the velocity, $\bm{g}=-g\hat{\bm{e}}_{z}$, is the gravitational acceleration with $g>0$, and $\nu$ is the constant kinematic viscosity. ${\bm{F}}_{\rm rad}$ and ${\bm{F}}_{\rm SGS}$ are the radiative and subgrid scale (SGS) fluxes, respectively, and $\bm{\mathsf{S}}$ is the traceless rate-of-strain tensor with $\mathsf{S}_{ij}={\textstyle{1\over 2}}(u_{i,j}+u_{j,i})-{\textstyle{1\over 3}}\delta_{ij}\bm{\nabla}\cdot{\bm{u}}$. We consider an optically thick, fully ionized gas. Thus, radiation is taken into account through the diffusion approximation, and the ideal gas equation of state $p={\cal R}\rho T$ applies, where ${\cal R}=c_{\rm P}-c_{\rm V}$ is the gas constant and $c_{\rm P,V}$ are the specific heats at constant pressure and volume, respectively. The radiative flux is given by ${\bm{F}}_{\rm rad}=-K\bm{\nabla}T$, where $K$ is the radiative heat conductivity. It has either a fixed profile $K(z)$ or it is a function of density and temperature, $K(\rho,T)$, given by $K=16\sigma_{\rm SB}T^{3}/3\kappa\rho$, where $\sigma_{\rm SB}$ is the Stefan–Boltzmann constant and $\kappa=\kappa_{0}(\rho/\rho_{0})^{a}(T/T_{0})^{b}$ is the opacity with coefficient $\kappa_{0}$, exponents $a$ and $b$, and reference values of density and temperature, $\rho_{0}$, $T_{0}$. These relations combine into

[TABLE]

Here we use $a=1$ and $b=-7/2$, corresponding to Kramers opacity law for free-free transitions (Weiss et al. 2004).

The radiative diffusivity $\chi=K/\rho c_{\rm P}$ can vary by several orders of magnitude as a function of depth in the Kramers opacity case. In order to keep the simulations numerically stable, we apply additional turbulent SGS diffusivities in the entropy equation,

[TABLE]

where $s^{\prime}=s-\overline{s}$ is the fluctuation of specific entropy, and $\chi_{\rm SGS}^{(0)}$ acts on the mean entropy, is non-zero only near the surface, while $\chi_{\rm SGS}^{(1)}$ is constant and acts on the entropy fluctuations. We use the Pencil Code111https://github.com/pencil-code.

2.2. Geometry, initial and boundary conditions

The computational domain is rectangular with $-2\leq(x,y)/d\leq 2$ and $-0.5\leq z/d\leq 1$, where $d$ is the depth of the initially isentropic layer. The initial stratification consists of two polytropic layers with indices $n_{1}=3.25$ in $-0.5\leq z/d<0$ and $n_{2}=1.5$ in $0\leq z/d\leq 1$. The former is the same as in the special case where the temperature gradient in the corresponding hydrostatic state is constant; see Barekat & Brandenburg (2014).

The horizontal boundaries are periodic whereas the vertical boundaries are impenetrable and stress free for the flow. We set the temperature gradient at the bottom according to ${\partial_{z}T}=-{F_{\rm tot}}/{K_{\rm bot}}$, where $F_{\rm tot}$ is a fixed input flux and $K_{\rm bot}$ is the value of the heat conductivity at the bottom of the domain. On the upper boundary we assume, for simplicity, a fixed gradient of specific entropy such that $({d}/{c_{\rm P}})\,{\partial_{z}s}\equiv\widetilde{\partial_{z}s}=-10$. This condition allows the density and temperature to vary locally.

2.3. Control parameters and diagnostics

Our models are fully defined by choosing the values of $\nu$, $g$, $a$, $b$, $\widetilde{\partial_{z}s}$, $F_{\rm tot}$, $K_{0}$, $\rho_{0}$, $T_{0}$, the SGS Prandtl numbers ${\rm Pr}_{\rm SGS}^{(0)}={\nu}/{\chi_{\rm SGS}^{(0)}}(z/d=1)$ and ${\rm Pr}_{\rm SGS}^{(1)}={\nu}/{\chi_{\rm SGS}^{(1)}}$, the $z$-dependent profile of $\chi_{\rm SGS}^{(0)}$, and the initial normalized pressure scale height at the surface, $\xi_{0}\equiv H_{\rm p}^{\rm(top)}/d={{{\cal R}}\overline{T}_{\rm top}}/{gd}$. The value of $K_{0}$ is fixed by assuming $F_{\rm rad}=F_{\rm tot}$ at the bottom of the domain. The normalized input flux is given by $F_{n}=F_{\rm tot}/\rho_{\rm bot}c_{\rm s,bot}^{3}$, where $\rho_{\rm bot}$ and $c_{\rm s,bot}$ are density and sound speed, respectively, at $z/d=-0.5$ in the initial non-convecting state. We also quote the Reynolds number, ${\rm Re}={u_{\rm rms}}/{\nu k_{1}}$, where $u_{\rm rms}$ is the volume-averaged rms velocity and $k_{1}=2\pi/d$.

Dominant contributions to the mean vertical energy flux are

[TABLE]

The viscous energy flux is negligible. No mean flows are generated, hence primes on $\bm{u}$ are dropped.

3. Results

Here, we describe the results of three simulations where the heat conductivity is either based on Kramers’ law (Run K), or it has a fixed profile that either coincides with the Kramers conductivity (P) in the initial state of Run K or a piecewise constant profile (S); see, e.g., Hurlburt et al. (1994).

3.1. Revising the CZ structure

As a basis for our analysis, we show in Figures 1(a) and (b) the energy fluxes, defined by Equations (6) and (7), and the superadiabaticity, $\nabla-\nabla_{\rm ad}$. Depending on the signs of enthalpy flux and superadiabaticity, four different regimes and corresponding zones can be identified; see Table 1. The top three layers, BZ, DZ, and OZ, are efficiently mixed while in the lowermost (radiative) layer (RZ), mixing is inefficient. In $z>z_{\rm BZ}$, we have $\nabla\geq\nabla_{\rm ad}$, whereas in $z_{\rm DZ}<z<z_{\rm BZ}$ we have $\nabla\leq\nabla_{\rm ad}$ and yet $\overline{F}_{\rm enth}>0$. In $z<z_{\rm DZ}$, we have $\overline{F}_{\rm enth}<0$, while in $z<z_{\rm OZ}$ we also have $|\overline{F}_{\rm enth}|\leq 0.03F_{\rm tot}$. The positions and thicknesses of the respective layers, $d_{\rm BZ}$, $d_{\rm DZ}$, and $d_{\rm OZ}$ (see Table 1), are listed in Table 2. We refer to the union of BZ, DZ, and OZ as the ‘mixed zone’ (MZ). Our BZ and OZ are in the traditional parlance the CZ and OZ, respectively, while the DZ has no counterpart in the usual paradigm of convection. Here, we consider the layers where $\overline{F}_{\rm enth}>0$, i.e. the combination of BZ and DZ, as the revised CZ.

We identify a DZ in Runs K and P. In Run K, $\overline{F}_{\rm enth}$ remains positive down to $z/d\approx 0$, although the superadiabaticity already turns negative at $z/d=0.26$. Runs P and K are similar, but the BZ is somewhat shallower in P. This is a consequence of the static profile of the heat conductivity as opposed to the dynamic formulation of Run K, where the depth of the MZ is not known a priori. In Run S with a fixed step profile for $K$, the difference is more striking: even though the depths of the MZ and CZ are the same as in Run K, the DZ is negligibly thin; see Figure 1(b). This is due to the fact that the constant heat conductivity above $z/d=0$ forces radiative diffusion to transport a certain fraction of the flux (Brandenburg et al. 2005) and the abrupt change of $K$ around $z=0$ prevents a smooth transition to a stable stratification beneath.

It is remarkable that in Runs K and P, the lower $\sim 40\%$ of the MZ is stably stratified according to the Schwarzschild criterion. In these runs, the mixed, but stably stratified layer is roughly equally divided into DZ and OZ. This is similar to the results of Brandenburg et al. (2000), who were the first to use a Kramers-based heat conductivity. Roxburgh & Simmons (1993) used a temperature-dependent heat conductivity in earlier two-dimensional simulations and found a subadiabatic convective layer at the base of the CZ. However, simulations of a $5M_{\odot}$ red giant star, employing a heat conduction profile based on OPAL opacities, did not indicate a DZ (Viallet et al. 2013). An extended subadiabatic convective layer was also reported by Chan & Gigas (1992) from a large-eddy simulation, although they applied a prescribed step function for the heat conductivity. However, their models had low resolution and their subsequent works did not mention similar findings. More recent studies reported subadiabatic convective layers in different contexts (e.g. Tremblay et al. 2015; Korre et al. 2017; Hotta 2017).

3.2. Why does an extended DZ exist?

We show in Figure 2(a) for Run K that the fluid in the upflows (downflows) is lighter (heavier) and warmer (cooler) than average in almost all of the DZ, indicating that both contribute to positive $\overline{F}_{\rm enth}$. Furthermore, Figure 2(b) shows that, in the BZ, the total force $\overline{\mathscr{F}}_{z}=\overline{\rho Du_{z}/Dt}$ is negative for the downflows and changes sign at $z_{\rm BZ}$, while for the upflows, $\overline{\mathscr{F}}_{z}$ is positive everywhere except very near the surface. The associated power, $\overline{P}_{z}=\overline{u_{z}\mathscr{F}_{z}}$ (blue dashed), shows that the downflows gain energy mostly near the surface while the upflows (red dashed) are accelerated throughout the MZ except near the surface. The viscous force is non-negligible only near the surface (dotted and dash–dotted). We speculate that there the upflows are decelerated by viscous momentum exchange with the downflows.

Based on these data, we interpret the DZ as an overshooting phenomenon, but with an upward enthalpy flux, and with the resulting stable stratification being nearly adiabatic. We explain the appearance of the DZ such that cool fluid elements that originate near the surface are accelerated downward by gravity and gain enough momentum to penetrate the convectively stable layer beneath the BZ. There they are progressively decelerated and heated up. If this proceeds fast enough, fluid elements, having kept a sufficient part of their momentum, can continue moving downward, but now with an excess of entropy thus forming the OZ. The upflows in the stably stratified OZ cannot be due to convective instability, but are driven by the pressure excess exerted by the matter in the downflows. In the Schwarzschild stable DZ, the upflows are lighter than their surroundings, which is an important property of the DZ; see Figure 1(d) of Brandenburg (2016). However, the force on the downflows is decelerating; see Figure 2(b). Therefore, the upflows in the DZ are not buoyancy-driven. Instead, we argue that they are pressure driven, just like in the OZ.

We conclude that the Deardorff layer is associated with nonlocal transport of heat as the downdrafts of cool matter originating from the strongly cooled surface propagate not only through the BZ, but further on to the bottom of the DZ. This process is called entropy rain which characterizes stellar convection as being driven by radiative cooling at the surface (Stein & Nordlund 1989).

In an attempt to quantify the different contributions to the enthalpy flux, we compare the numerical results for $\overline{F}_{\rm enth}$ with a mean-field parameterization that takes into account the non-gradient contribution introduced by Deardorff (1961, 1966). Here, we use the expression derived by Brandenburg (2016):

[TABLE]

where $\tau_{\rm red}$ is a reduced relaxation time, taking into account radiative cooling and turbulent energy transport. The first term, $\overline{F}_{\rm D}$, describes the non-gradient Deardorff flux, which is positive irrespective of the sign of the entropy gradient, whereas the latter term, $\overline{F}_{\rm G}$, is the traditional mean-field description of the (gradient) enthalpy flux. In Figure 2(c), we show $F_{\rm D}$, $F_{\rm G}$, and their sum for Run K, assuming that $\tau_{\rm red}(z)=c_{\tau}\tau(z)$, where $\tau(z)=\min(\tau_{\rm to},\tau_{\rm bu})$ is the minimum of the convective turnover time $\tau_{\rm to}=H_{\rm p}/u_{\rm rms}$ and the buoyancy timescale $\tau_{\rm bu}=(c_{\rm P}/g)(\overline{u_{z}^{2}}/\overline{{s^{\prime}}^{2}})^{1/2}$, and $c_{\tau}$ is a free parameter, for the best fit set to 0.73. Figure 2(c) shows that Equation (8) provides a good description in the BZ and even suggests that the Deardorff term is the dominant contribution to the heat transport. However, the expression in Equation (8) breaks down in the DZ and OZ. Further, we separate $\overline{F}_{\rm enth}$ in Figure 2(c) into contributions from upflows ($\overline{F}_{\rm enth}^{\uparrow}$) and downflows ($\overline{F}_{\rm enth}^{\downarrow}$) for Run K. The downflows dominate the enthalpy flux with $\overline{F}_{\rm enth}^{\downarrow}\approx 3\overline{F}_{\rm enth}^{\uparrow}$. We find that gradient and Deardorff contributions match $\overline{F}_{\rm enth}^{\uparrow}$ and $\overline{F}_{\rm enth}^{\downarrow}$, respectively, in the BZ. However, $F_{\rm D}$ and $F_{\rm G}$ contain contributions from both upflows and downflows, and the correspondence to $\overline{F}_{\rm enth}^{\uparrow}$ and $\overline{F}_{\rm enth}^{\downarrow}$ is likely coincidental. This conjecture is supported by Figure 2(b), which suggests that the downflows feel the local entropy gradient. The generality of these results will be investigated elsewhere using wider parameter studies.

Two recent studies (Korre et al. 2017; Hotta 2017) have reported subadiabatic layers from convection simulations. The former authors studied “weakly non-Boussinesq” convection in spherical coordinates where the radial dependence of the superadiabatic temperature gradient of the background state was varied. They found that a subadiabatic layer appeared in regions where convective transport was efficient (or radiative diffusion inefficient). They also studied the contributions of upflows and downflows separately and found that the upflows in these cases contributed to downward flux of heat. This is qualitatively different to our cases where $\overline{F}_{\rm enth}^{\uparrow}$ is always nearly zero (OZ, lower part of DZ) or positive (upper part of DZ, BZ). In their case, the upflows are clearly pressure driven even in the bulk of the CZ – in contrast to our simulations.

The study of Hotta (2017) is a close parallel to ours in that a density-stratified, initially piecewise polytropic setup was used to study overshooting in fully compressible simulations, although with magnetic fields. The main difference to our runs is that our density contrast is larger ($37$ in the CZ of Run K) compared to about 6–7 in Hotta (2017). Furthermore, he used a fixed profile of $K$ that is smoother than in our Run S, but steeper than in our Run P. This results in a similar subadiabatic layer at the base of the CZ as in our Runs K and P. Moreover, Hotta (2017) analyzed the vertical force balance (his Figure 13) and came to the conclusion that the total buoyancy force switches sign roughly at the same location as the entropy gradient. We conclude that the mechanism forming a subadiabatic layer in the runs of Hotta (2017) is very likely the same as in our models.

3.3. Structure of convection

Given the appearance of an extended DZ, we analyze the flow in detail to find out whether its topology in the DZ is altered in comparison to the BZ. We adopt the approach of Brandenburg (2016), where the structure of convection is characterized by the number and filling factor of downflows (cf. his Figure 2). We have developed a dedicated algorithm (to be reported on elsewhere) to detect isolated downflow plumes (IDPs) and intergranular lanes (IGLs), compute their sizes and numbers, and thereby their filling factors for such an analysis. In Case I of Brandenburg (2016), corresponding to forest-like downflow structures, number, size, and filling factor of the downflows are independent of depth. His MLT description including the Deardorff flux is closest to his Case III, with a tree-like structure, where the filling factor is constant, but the number (size) of the downflows decreases (increases).

In Figures 3(a) and (b), we show representative patterns of the vertical velocity of Run K from BZ, DZ, and OZ. They are qualitatively similar to those found in numerous other studies of stratified convection (e.g. Stein & Nordlund 1989; Stein et al. 2009; Hotta et al. 2014; Käpylä et al. 2016; Kitiashvili et al. 2016): upwelling granules with downflows along a network of connected IGLs near the surface. Deeper down, the cellular structure disintegrates and IDPs appear. In the DZ and OZ only a few IDPs survive in the midst of much larger scale upflows.

The filling factors $f_{\downarrow}$ of all downflows, and of IGLs and IDPs separately (for the latter two a proxy using average sizes and numbers to ease comparison with panels (d), (e)), are shown in Figure 3(c). We find that $f_{\downarrow}$ of all downflows is almost independent of depth in the redefined CZ (= BZ + DZ) for all runs. The filling factors of IDPs and IGLs reveal that for Run S in the bulk of the CZ ($0.3\lesssim z/d\lesssim 0.8$), the dominant IGLs have roughly a constant filling factor while their size increases and their number decreases, consistent with the tree-like structure of Case III. For Runs K and P, the structure of convection is clearly distinct from S, with the IGL network starting to disappear at much smaller depths, and the IDPs already taking over at $z/d\approx 0.6$.

After the IDPs take over as the dominating structure of convection, we observe another difference between Run S and Runs K and P. In the latter cases, after a smooth transition, the IDP filling factor attains a constant value, while their size is constant and the number is mildly increasing. These data correspond most closely to Case I. This holds roughly at depths $0\lesssim z/d\lesssim 0.35$, encompassing the DZ and the lower parts of the BZ. Thus, the tree-like picture roughly holds until $z/d\approx 0.35$, below which a depth-independent number of IDPs persist. In Run S, the IDP filling factor also tends to a constant in between $-0.1\lesssim z/d\lesssim 0.2$, but this is accompanied by an increase of the number and decrease of the size of IDPs, incompatible with Cases I and III. These results suggest that the structure of the downflows in the DZ and the bottom part of the BZ is qualitatively different from that in the upper parts of the BZ (forest-like instead of tree-like).

4. Conclusions

We have shown that, when a smoothly varying heat conduction profile is used, a substantial part of the lower CZ is weakly subadiabatic although $\overline{F}_{\rm enth}>0$. A smooth transition can also be expected to occur in deep stellar interiors where a Kramers-based conductivity is valid. Furthermore, with such a heat conduction law, the depth of the CZ is an outcome of the simulation and cannot be determined a priori.

We have shown that the subadiabatic layer can still lead to an upward enthalpy flux due to downflows bringing low entropy material from near the surface to the stably stratified layers below. We also found that the upflows in the overshoot zone are driven by the pressure excess due to the matter brought down by the downflows. Except for the lowermost parts, the ascending matter in the Deardorff layer is lighter than the surroundings. Yet, we argue that also in the Deardorff layer the upflows are pressure driven. There is no buoyant acceleration, and the downflows are instead decelerated in accordance with the Schwarzschild criterion. Our results confirm that convection is highly nonlocal and driven by cooling at the surface resulting in cool entropy rain. The traditional mean-field expression of the enthalpy flux fails in the subadiabatic part of the CZ and a non-gradient term is required. Our work demonstrates the existence of such a contribution, introduced by Deardorff (1961) and applied to stellar MLT by Brandenburg (2016), for the first time from numerical simulations. A geometric analysis shows a transition from a tree-like to a forest-like structure in the deep parts of the CZ. Energetic considerations reveal that the downflows provide the dominant contribution to the enthalpy flux everywhere in the CZ. Furthermore, as they are largely responsible for driving the upflows, their importance for the overall convective structure is indeed crucial.

The current simulations may have too low resolution to fully capture the driving of strong downflows near the surface so that their effect in more realistic conditions can be even more pronounced. This can lead to a further reduction of the depth of the BZ, which, in conjunction with the topology change, could alleviate the discrepancy between helioseismic and simulation-based estimates of convective velocities (Miesch et al. 2012). Furthermore, our results are obviously at odds with MLT, which is widely used in stellar structure models, calling for more advanced one-dimensional models (see, e.g. Kupka 1999; Snellman et al. 2015). Another implication of a subadiabatic layer above the base of the CZ comes from potentially breaking the Taylor–Proudman balance of the solar rotation profile (Rempel 2005). These questions will be addressed elsewhere.

The authors thank an anonymous referee for his/her constructive comments and criticism on the manuscript. The simulations were performed using the supercomputers hosted by CSC – IT Center for Science Ltd. in Espoo, Finland, who are administered by the Finnish Ministry of Education. Special Grand Challenge allocation NEOCON is acknowledged. Financial support from the Academy of Finland grant No. 272157 to the ReSoLVE Centre of Excellence (P.J.K., M.R., M.J.K., N.O.) is acknowledged. J.W. acknowledges funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant agreement No. 623609.