Extended subadiabatic layer in simulations of overshooting convection
Petri J. Käpylä, Matthias Rheinhardt, Axel Brandenburg, Rainer Arlt, Maarit J. Käpylä, Andreas Lagg, Nigul Olspert, Jörn Warnecke
AA STROPHYS . J. 845, L23 (2017)
Preprint typeset using L A TEX style emulateapj v. 08/22/09
EXTENDED SUBADIABATIC LAYER IN SIMULATIONS OF OVERSHOOTING CONVECTION P ETRI
J. K ¨
APYL ¨ A , , , M ATTHIAS R HEINHARDT , A XEL B RANDENBURG , , , , R AINER A RLT , M AARIT
J. K ¨
APYL ¨ A , ,A NDREAS L AGG , N IGUL O LSPERT , AND
J ¨
ORN W ARNECKE Leibniz-Institut f¨ur Astrophysik, An der Sternwarte 16, D-14482 Potsdam, Germany ReSoLVE Centre of Excellence, Department of Computer Science, P.O. Box 15400, FI-00076 Aalto, Finland Max-Planck-Institut f¨ur Sonnensystemforschung, Justus-von-Liebig-Weg 3, D-37077 G¨ottingen, Germany NORDITA, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-10691 Stockholm, Sweden Department of Astronomy, AlbaNova University Center, Stockholm University, SE-10691 Stockholm, Sweden JILA and Department of Astrophysical and Planetary Sciences, Box 440, University of Colorado, Boulder, CO 80303, USA Laboratory for Atmospheric and Space Physics, 3665 Discovery Drive, Boulder, CO 80303, USA (Dated: Received 2017 March 20; revised 2017 July 27; accepted 2017 July 29; published 2017 August 22)
Astrophys. J. 845, L23 (2017)
ABSTRACTWe present numerical simulations of hydrodynamic overshooting convection in local Cartesian domains. Wefind that a substantial fraction of the lower part of the convection zone (CZ) is stably stratified according tothe Schwarzschild criterion while the enthalpy flux is outward directed. This occurs when the heat conductionprofile at the bottom of the CZ is smoothly varying, based either on a Kramers-like opacity prescription as afunction of temperature and density or a static profile of a similar shape. We show that the subadiabatic layerarises due to nonlocal energy transport by buoyantly driven downflows in the upper parts of the CZ. Analysis ofthe force balance of the upflows and downflows confirms that convection is driven by cooling at the surface. Wefind that the commonly used prescription for the convective enthalpy flux being proportional to the negativeentropy gradient does not hold in the stably stratified layers where the flux is positive. We demonstrate theexistence of a non-gradient contribution to the enthalpy flux, which is estimated to be important throughout theconvective layer. A quantitative analysis of downflows indicates a transition from a tree-like structure wheresmaller 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 substantialsubadiabatic layer is present in the lower part of the CZ.
Subject headings:
Hydrodynamics — convection — turbulence INTRODUCTION
Convection plays a vital role in stellar activity by generatingturbulence that, together with the star’s overall rotation, leadsto differential rotation (e.g. R¨udiger 1989) and dynamo action(e.g. Krause & R¨adler 1980). Energy transport due to con-vection is important for almost all stars during some stages oftheir evolution (e.g. Kippenhahn et al. 2012). Hence, a properparameterization of convection is crucial for stellar structureand evolution in one-dimensional models.Mixing length theory (MLT) continues to be a popular de-scription of stellar convection. The formulation of MLT, asit is used today, goes back to the seminal work of Vitense(1953), where the properties of convection are related to thelocal value of the superadiabatic gradient ∇ − ∇ ad , with ∇ = d ln T /d ln p and ∇ ad being the actual and adiabaticlogarithmic temperature gradients, respectively. Here, T and p are temperature and pressure, respectively, while overbarsdenote horizontal averaging. Convection is supposed to oc-cur only if the horizontally averaged temperature stratifica-tion is superadiabatic, ∇ > ∇ ad , which is equivalent to theSchwarzschild criterion for convection, d s/dz < , where s is the specific entropy.The MLT has deeply influenced the way in which three-dimensional ab initio convection models are constructed: of-ten a fixed profile of radiative heat conductivity K is chosen,producing a superadiabatic layer of fixed depth (e.g. Hurlburtet al. 1986). Alternatively, the static thermodynamic back-ground state is taken from an MLT-based stellar evolutionmodel (e.g. Brun et al. 2011; Kitiashvili et al. 2016). In such setups, convection is driven at the largest scale available. Thisis a possible cause for the discrepancy in the convective veloc-ities at large horizontal scales between simulations and time–distance helioseismology (Hanasoge et al. 2012). Smallerlength scales could instead be imprinted by convection drivensolely by the surface layers (Cossette & Rast 2016), whichleads to a topology change of the downdrafts from a tree-likestructure near the surface to strong plumes penetrating intodeeper layers as cool entropy rain ; see Spruit (1997). The lat-ter can take part in the convective flux in these layers through anon-gradient contribution known as Deardorff flux (Deardorff1966; Brandenburg 2016).We present simulations in which we use either a physicallymotivated heat conduction formulation based on a Kramers-like opacity (Brandenburg et al. 2000) or two types of fixedheat conductivity profiles to study their effect on the struc-ture of the convection zone (CZ). Furthermore, we demon-strate the existence of a non-gradient contribution to the en-thalpy flux and use a quantitative analysis to study the topol-ogy change of the downflow and upflow structures in the sim-ulations. THE MODEL
Basic equations
We solve the equations of compressible hydrodynamics: D ln ρDt = − ∇ · u , (1) a r X i v : . [ a s t r o - ph . S R ] A ug D u Dt = g − ρ ( ∇ p − ∇ · νρ S ) , (2) T DsDt = − ρ [ ∇ · ( F rad + F SGS )] + 2 ν S , (3)where D/Dt = ∂/∂t + u · ∇ is the advective derivative, ρ is the density, u is the velocity, g = − g ˆ e z , is the gravita-tional acceleration with g > , and ν is the constant kine-matic viscosity. F rad and F SGS are the radiative and subgridscale (SGS) fluxes, respectively, and S is the traceless rate-of-strain tensor with S ij = ( u i,j + u j,i ) − δ ij ∇ · u . Weconsider an optically thick, fully ionized gas. Thus, radia-tion is taken into account through the diffusion approxima-tion, and the ideal gas equation of state p = R ρT applies,where R = c P − c V is the gas constant and c P , V are the spe-cific heats at constant pressure and volume, respectively. Theradiative flux is given by F rad = − K ∇ T , where K is theradiative heat conductivity. It has either a fixed profile K ( z ) or it is a function of density and temperature, K ( ρ, T ) , givenby K = 16 σ SB T / κρ , where σ SB is the Stefan–Boltzmannconstant and κ = κ ( ρ/ρ ) a ( T /T ) b is the opacity with coef-ficient κ , exponents a and b , and reference values of densityand temperature, ρ , T . These relations combine into K ( ρ, T ) = K ( ρ/ρ ) − ( a +1) ( T /T ) − b . (4)Here we use a = 1 and b = − / , corresponding to Kramersopacity law for free-free transitions (Weiss et al. 2004).The radiative diffusivity χ = K/ρc P can vary by severalorders of magnitude as a function of depth in the Kramersopacity case. In order to keep the simulations numericallystable, we apply additional turbulent SGS diffusivities in theentropy equation, F SGS = − ρT (cid:16) χ (0)SGS ∇ s + χ (1)SGS ∇ s (cid:48) (cid:17) , (5)where s (cid:48) = s − s is the fluctuation of specific entropy, and χ (0)SGS acts on the mean entropy, is non-zero only near the sur-face, while χ (1)SGS is constant and acts on the entropy fluctua-tions. We use the P ENCIL C ODE . Geometry, initial and boundary conditions
The computational domain is rectangular with − ≤ ( x, y ) /d ≤ and − . ≤ z/d ≤ , where d is the depthof the initially isentropic layer. The initial stratification con-sists of two polytropic layers with indices n = 3 . in − . ≤ z/d < and n = 1 . in ≤ z/d ≤ . The for-mer is the same as in the special case where the temperaturegradient in the corresponding hydrostatic state is constant; seeBarekat & Brandenburg (2014).The horizontal boundaries are periodic whereas the verticalboundaries are impenetrable and stress free for the flow. Weset the temperature gradient at the bottom according to ∂ z T = − F tot /K bot , where F tot is a fixed input flux and K bot is thevalue of the heat conductivity at the bottom of the domain. Onthe upper boundary we assume, for simplicity, a fixed gradientof specific entropy such that ( d/c P ) ∂ z s ≡ (cid:103) ∂ z s = − . Thiscondition allows the density and temperature to vary locally. https://github.com/pencil-code Table 1
Definition of the Zones. F enth ∇ − ∇ ad Zone Label Lower Limit Thickness > > Buoyancy BZ z BZ d BZ > < Deardorff DZ z DZ d DZ < < Overshoot OZ z OZ d OZ ≈ < Radiative RZ · · · · · ·
Control parameters and diagnostics
Our models are fully defined by choosing the values of ν , g , a , b , (cid:103) ∂ z s , F tot , K , ρ , T , the SGS Prandtl numbers Pr (0)SGS = ν/χ (0)SGS ( z/d = 1) and Pr (1)SGS = ν/χ (1)SGS , the z -dependent profile of χ (0)SGS , and the initial normalized pres-sure scale height at the surface, ξ ≡ H (top)p /d = R T top /gd .The value of K is fixed by assuming F rad = F tot at thebottom of the domain. The normalized input flux is givenby F n = F tot /ρ bot c , bot , where ρ bot and c s , bot are densityand sound speed, respectively, at z/d = − . in the initialnon-convecting state. We also quote the Reynolds number, Re = u rms /νk , where u rms is the volume-averaged rms ve-locity and k = 2 π/d .Dominant contributions to the mean vertical energy flux are F rad = − K ∂ z T , F kin = ρ u u z , (6) F enth = c P ( ρu z ) (cid:48) T (cid:48) , F (0)SGS = − χ (0)SGS ρT ∂ z s. (7)The viscous energy flux is negligible. No mean flows are gen-erated, hence primes on u are dropped. RESULTS
Here, we describe the results of three simulations where theheat conductivity is either based on Kramers’ law (Run K),or it has a fixed profile that either coincides with the Kramersconductivity (P) in the initial state of Run K or a piecewiseconstant profile (S); see, e.g., Hurlburt et al. (1994).
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 thesuperadiabaticity, ∇ − ∇ ad . Depending on the signs of en-thalpy flux and superadiabaticity, four different regimes andcorresponding zones can be identified; see Table 1. The topthree layers, BZ, DZ, and OZ, are efficiently mixed while inthe lowermost (radiative) layer (RZ), mixing is inefficient. In z > z BZ , we have ∇ ≥ ∇ ad , whereas in z DZ < z < z BZ we have ∇ ≤ ∇ ad and yet F enth > . In z < z DZ , wehave F enth < , while in z < z OZ we also have | F enth | ≤ . F tot . The positions and thicknesses of the respective lay-ers, d BZ , d DZ , and d 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 CZand OZ, respectively, while the DZ has no counterpart in theusual paradigm of convection. Here, we consider the layerswhere F enth > , i.e. the combination of BZ and DZ, as therevised CZ.We identify a DZ in Runs K and P. In Run K, F enth remainspositive down to z/d ≈ , although the superadiabaticity al-ready turns negative at z/d = 0 . . Runs P and K are similar,but the BZ is somewhat shallower in P. This is a consequenceof the static profile of the heat conductivity as opposed to thedynamic formulation of Run K, where the depth of the MZ is Figure 1.
Solid lines: radiative (purple), enthalpy (blue), kinetic energy(magenta), and SGS (green) fluxes for Runs K (a) and S (b). Red: ∇ − ∇ ad .Dashed lines in (a): corresponding data from Run P. The thick horizontallines on the abscissae mark the extent of the MZ. Table 2
Summary of the runs, All with meshpoints.Run K Re z BZ z DZ z OZ d BZ d DZ d OZ K Kramers
27 0 .
26 0 . − .
27 0 .
74 0 .
26 0 . P profile
25 0 .
34 0 . − .
19 0 .
66 0 .
24 0 . S step
26 0 .
02 0 . − .
28 0 .
98 0 .
02 0 . Note . — Column ‘ K ’: heat conduction scheme. Remaining columns: depthsand thicknesses of the zones; see Table 1. Pr (0)SGS = 0 . , Pr (1)SGS = 1 , F n ≈ . · − , and ξ = 0 . . not known a priori . In Run S with a fixed step profile for K ,the difference is more striking: even though the depths of theMZ and CZ are the same as in Run K, the DZ is negligiblythin; see Figure 1(b). This is due to the fact that the constantheat conductivity above z/d = 0 forces radiative diffusionto transport a certain fraction of the flux (Brandenburg et al.2005) and the abrupt change of K around z = 0 prevents asmooth transition to a stable stratification beneath.It is remarkable that in Runs K and P, the lower ∼ of the MZ is stably stratified according to the Schwarzschildcriterion. In these runs, the mixed, but stably stratified layeris roughly equally divided into DZ and OZ. This is similar tothe results of Brandenburg et al. (2000), who were the firstto use a Kramers-based heat conductivity. Roxburgh & Sim-mons (1993) used a temperature-dependent heat conductivityin earlier two-dimensional simulations and found a subadia-batic convective layer at the base of the CZ. However, simu-lations of a M (cid:12) red giant star, employing a heat conductionprofile based on OPAL opacities, did not indicate a DZ (Vial-let et al. 2013). An extended subadiabatic convective layerwas also reported by Chan & Gigas (1992) from a large-eddy Figure 2. (a) Temperature fluctuation (solid, left axis) and density fluctua-tion (dashed, right axis), averaged separately over upflows (red) and down-flows (blue). (b) Horizontally averaged force F z = ρDu z /Dt (solid lines,left axis), and the accelerating power P z = u z F z of those forces (dashed,right axis). P ↑ z is scaled up by a factor of five. F visc and P visc are thecorresponding viscous force and its power. (c) Averaged enthalpy flux (solidlines) along with parameterizations according to Equation (8) (dashed). Theinset shows τ to and τ bu as functions of depth in units of (cid:112) d/g . Data forRun K. simulation, although they applied a prescribed step functionfor the heat conductivity. However, their models had low res-olution and their subsequent works did not mention similarfindings. More recent studies reported subadiabatic convec-tive layers in different contexts (e.g. Tremblay et al. 2015;Korre et al. 2017; Hotta 2017). Why does an extended DZ exist?
We show in Figure 2(a) for Run K that the fluid in the up-flows (downflows) is lighter (heavier) and warmer (cooler)than average in almost all of the DZ, indicating that both con-tribute to positive F enth . Furthermore, Figure 2(b) shows that,in the BZ, the total force F z = ρDu z /Dt is negative for thedownflows and changes sign at z BZ , while for the upflows, Figure 3. (a) Vertical velocity and (b) volume rendering of downflows at the periphery and from depths z/d = (0 . , , , . , − . corresponding to anear-surface layer, and the middles of BZ, DZ, and OZ, respectively, in Run K. Tildes refer to normalization with ˜ u z = u z / √ gd . (c) Filling factor of downflows(dash–dotted lines), a proxy filling factor for IDPs (solid), IGLs (dashed, only occurring at z/d (cid:38) . ) and their sum (dotted), (d) average widths of IDPs andIGLs with the standard deviations from the temporal evolution (shaded areas), and (e) their numbers as functions of z from Runs K (black), P (blue), and S (red). F z is positive everywhere except very near the surface. Theassociated power, P z = u z F z (blue dashed), shows that thedownflows gain energy mostly near the surface while the up-flows (red dashed) are accelerated throughout the MZ exceptnear the surface. The viscous force is non-negligible only nearthe surface (dotted and dash–dotted). We speculate that therethe upflows are decelerated by viscous momentum exchangewith the downflows.Based on these data, we interpret the DZ as an overshoot-ing phenomenon , but with an upward enthalpy flux, and withthe resulting stable stratification being nearly adiabatic. Weexplain the appearance of the DZ such that cool fluid ele-ments that originate near the surface are accelerated down-ward by gravity and gain enough momentum to penetrate theconvectively stable layer beneath the BZ. There they are pro-gressively decelerated and heated up. If this proceeds fastenough, fluid elements, having kept a sufficient part of theirmomentum, can continue moving downward, but now with an excess of entropy thus forming the OZ. The upflows in the sta-bly stratified OZ cannot be due to convective instability, butare driven by the pressure excess exerted by the matter in thedownflows. In the Schwarzschild stable DZ, the upflows arelighter than their surroundings, which is an important propertyof 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. In-stead, we argue that they are pressure driven, just like in theOZ. We conclude that the Deardorff layer is associated with nonlocal transport of heat as the downdrafts of cool matteroriginating from the strongly cooled surface propagate notonly through the BZ, but further on to the bottom of the DZ.This process is called entropy rain which characterizes stellarconvection as being driven by radiative cooling at the surface(Stein & Nordlund 1989).In an attempt to quantify the different contributions to theenthalpy flux, we compare the numerical results for F enth with a mean-field parameterization that takes into accountthe non-gradient contribution introduced by Deardorff (1961,1966). Here, we use the expression derived by Brandenburg(2016): F MFenth = τ red ρ T (cid:0) g s (cid:48) /c P − u z ∂ z s (cid:1) ≡ F D + F G , (8)where τ red is a reduced relaxation time, taking into accountradiative cooling and turbulent energy transport. The firstterm, F D , describes the non-gradient Deardorff flux, whichis positive irrespective of the sign of the entropy gradient,whereas the latter term, F G , is the traditional mean-fielddescription of the (gradient) enthalpy flux. In Figure 2(c),we show F D , F G , and their sum for Run K, assuming that τ red ( z ) = c τ τ ( z ) , where τ ( z ) = min( τ to , τ bu ) is the min-imum of the convective turnover time τ to = H p /u rms andthe buoyancy timescale τ bu = ( c P /g )( u z /s (cid:48) ) / , and c τ is afree parameter, for the best fit set to 0.73. Figure 2(c) showsthat Equation (8) provides a good description in the BZ andeven suggests that the Deardorff term is the dominant contri-bution to the heat transport. However, the expression in Equa-tion (8) breaks down in the DZ and OZ. Further, we separate F enth in Figure 2(c) into contributions from upflows ( F ↑ enth )and downflows ( F ↓ enth ) for Run K. The downflows dominatethe enthalpy flux with F ↓ enth ≈ F ↑ enth . We find that gradientand Deardorff contributions match F ↑ enth and F ↓ enth , respec-tively, in the BZ. However, F D and F G contain contributionsfrom both upflows and downflows, and the correspondence to F ↑ enth and F ↓ enth is likely coincidental. This conjecture is sup-ported by Figure 2(b), which suggests that the downflows feelthe local entropy gradient. The generality of these results willbe investigated elsewhere using wider parameter studies.Two recent studies (Korre et al. 2017; Hotta 2017) have re-ported subadiabatic layers from convection simulations. Theformer authors studied “weakly non-Boussinesq” convectionin spherical coordinates where the radial dependence of thesuperadiabatic temperature gradient of the background statewas varied. They found that a subadiabatic layer appeared inregions where convective transport was efficient (or radiativediffusion inefficient). They also studied the contributions ofupflows and downflows separately and found that the upflowsin these cases contributed to downward flux of heat. Thisis qualitatively different to our cases where F ↑ enth is alwaysnearly zero (OZ, lower part of DZ) or positive (upper part ofDZ, BZ). In their case, the upflows are clearly pressure driveneven in the bulk of the CZ – in contrast to our simulations.The study of Hotta (2017) is a close parallel to ours in thata density-stratified, initially piecewise polytropic setup wasused to study overshooting in fully compressible simulations,although with magnetic fields. The main difference to ourruns is that our density contrast is larger ( in the CZ ofRun 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 suba-diabatic layer at the base of the CZ as in our Runs K andP. Moreover, Hotta (2017) analyzed the vertical force balance(his Figure 13) and came to the conclusion that the total buoy-ancy force switches sign roughly at the same location as theentropy gradient. We conclude that the mechanism forminga subadiabatic layer in the runs of Hotta (2017) is very likelythe same as in our models. Structure of convection
Given the appearance of an extended DZ, we analyze theflow in detail to find out whether its topology in the DZis altered in comparison to the BZ. We adopt the approachof Brandenburg (2016), where the structure of convection ischaracterized 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 downflowplumes (IDPs) and intergranular lanes (IGLs), compute theirsizes and numbers, and thereby their filling factors for such ananalysis. In Case I of Brandenburg (2016), corresponding toforest-like downflow structures, number, size, and filling fac-tor of the downflows are independent of depth. His MLT de-scription 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 ofthe vertical velocity of Run K from BZ, DZ, and OZ. They arequalitatively similar to those found in numerous other studies of stratified convection (e.g. Stein & Nordlund 1989; Steinet al. 2009; Hotta et al. 2014; K¨apyl¨a et al. 2016; Kitiashviliet al. 2016): upwelling granules with downflows along a net-work of connected IGLs near the surface. Deeper down, thecellular structure disintegrates and IDPs appear. In the DZand OZ only a few IDPs survive in the midst of much largerscale upflows.The filling factors f ↓ of all downflows, and of IGLs andIDPs separately (for the latter two a proxy using average sizesand numbers to ease comparison with panels (d), (e)), areshown in Figure 3(c). We find that f ↓ of all downflows isalmost independent of depth in the redefined CZ (= BZ + DZ)for all runs. The filling factors of IDPs and IGLs reveal thatfor Run S in the bulk of the CZ ( . (cid:46) z/d (cid:46) . ), the dom-inant IGLs have roughly a constant filling factor while theirsize increases and their number decreases, consistent with thetree-like structure of Case III. For Runs K and P, the structureof convection is clearly distinct from S, with the IGL networkstarting to disappear at much smaller depths, and the IDPsalready taking over at z/d ≈ . .After the IDPs take over as the dominating structure of con-vection, we observe another difference between Run S andRuns K and P. In the latter cases, after a smooth transition,the IDP filling factor attains a constant value, while theirsize is constant and the number is mildly increasing. Thesedata correspond most closely to Case I. This holds roughly atdepths (cid:46) z/d (cid:46) . , encompassing the DZ and the lowerparts of the BZ. Thus, the tree-like picture roughly holds un-til z/d ≈ . , below which a depth-independent number ofIDPs persist. In Run S, the IDP filling factor also tends to aconstant in between − . (cid:46) z/d (cid:46) . , but this is accompa-nied by an increase of the number and decrease of the size ofIDPs, incompatible with Cases I and III. These results suggestthat the structure of the downflows in the DZ and the bottompart of the BZ is qualitatively different from that in the upperparts of the BZ (forest-like instead of tree-like). CONCLUSIONS
We have shown that, when a smoothly varying heat con-duction profile is used, a substantial part of the lower CZ isweakly subadiabatic although F enth > . A smooth tran-sition can also be expected to occur in deep stellar interiorswhere a Kramers-based conductivity is valid. Furthermore,with such a heat conduction law, the depth of the CZ is anoutcome of the simulation and cannot be determined a priori .We have shown that the subadiabatic layer can still lead toan upward enthalpy flux due to downflows bringing low en-tropy material from near the surface to the stably stratifiedlayers below. We also found that the upflows in the over-shoot zone are driven by the pressure excess due to the matterbrought down by the downflows. Except for the lowermostparts, the ascending matter in the Deardorff layer is lighterthan the surroundings. Yet, we argue that also in the Dear-dorff layer the upflows are pressure driven. There is no buoy-ant acceleration, and the downflows are instead decelerated inaccordance with the Schwarzschild criterion. Our results con-firm that convection is highly nonlocal and driven by coolingat the surface resulting in cool entropy rain. The traditionalmean-field expression of the enthalpy flux fails in the suba-diabatic 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 byBrandenburg (2016), for the first time from numerical simula-tions. A geometric analysis shows a transition from a tree-liketo a forest-like structure in the deep parts of the CZ. Energeticconsiderations reveal that the downflows provide the domi-nant contribution to the enthalpy flux everywhere in the CZ.Furthermore, as they are largely responsible for driving theupflows, their importance for the overall convective structureis indeed crucial.The current simulations may have too low resolution tofully capture the driving of strong downflows near the sur-face so that their effect in more realistic conditions can beeven more pronounced. This can lead to a further reduction ofthe depth of the BZ, which, in conjunction with the topologychange, could alleviate the discrepancy between helioseismicand simulation-based estimates of convective velocities (Mi-esch et al. 2012). Furthermore, our results are obviously atodds with MLT, which is widely used in stellar structure mod-els, calling for more advanced one-dimensional models (see,e.g. Kupka 1999; Snellman et al. 2015). Another implicationof a subadiabatic layer above the base of the CZ comes frompotentially breaking the Taylor–Proudman balance of the so-lar rotation profile (Rempel 2005). These questions will beaddressed elsewhere.The authors thank an anonymous referee for his/her con-structive comments and criticism on the manuscript. The sim-ulations were performed using the supercomputers hosted byCSC – IT Center for Science Ltd. in Espoo, Finland, whoare administered by the Finnish Ministry of Education. Spe-cial 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 fundingfrom the People Programme (Marie Curie Actions) of the Eu-ropean Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant agreement No. 623609.are driven by the pressure excess due to the matterbrought down by the downflows. Except for the lowermostparts, the ascending matter in the Deardorff layer is lighterthan the surroundings. Yet, we argue that also in the Dear-dorff layer the upflows are pressure driven. There is no buoy-ant acceleration, and the downflows are instead decelerated inaccordance with the Schwarzschild criterion. Our results con-firm that convection is highly nonlocal and driven by coolingat the surface resulting in cool entropy rain. The traditionalmean-field expression of the enthalpy flux fails in the suba-diabatic 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 byBrandenburg (2016), for the first time from numerical simula-tions. A geometric analysis shows a transition from a tree-liketo a forest-like structure in the deep parts of the CZ. Energeticconsiderations reveal that the downflows provide the domi-nant contribution to the enthalpy flux everywhere in the CZ.Furthermore, as they are largely responsible for driving theupflows, their importance for the overall convective structureis indeed crucial.The current simulations may have too low resolution tofully capture the driving of strong downflows near the sur-face so that their effect in more realistic conditions can beeven more pronounced. This can lead to a further reduction ofthe depth of the BZ, which, in conjunction with the topologychange, could alleviate the discrepancy between helioseismicand simulation-based estimates of convective velocities (Mi-esch et al. 2012). Furthermore, our results are obviously atodds with MLT, which is widely used in stellar structure mod-els, calling for more advanced one-dimensional models (see,e.g. Kupka 1999; Snellman et al. 2015). Another implicationof a subadiabatic layer above the base of the CZ comes frompotentially breaking the Taylor–Proudman balance of the so-lar rotation profile (Rempel 2005). These questions will beaddressed elsewhere.The authors thank an anonymous referee for his/her con-structive comments and criticism on the manuscript. The sim-ulations were performed using the supercomputers hosted byCSC – IT Center for Science Ltd. in Espoo, Finland, whoare administered by the Finnish Ministry of Education. Spe-cial 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 fundingfrom the People Programme (Marie Curie Actions) of the Eu-ropean Union’s Seventh Framework Programme (FP7/2007-2013) under REA grant agreement No. 623609.