Issue 
A&A
Volume 578, June 2015



Article Number  A106  
Number of page(s)  15  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201425125  
Published online  12 June 2015 
Overshooting by differential heating^{⋆}
Max Planck Institute for Astrophysics, KarlSchwarzschildstr. 1, 85748 Garching, Germany
email: robo@mpagarching.mpg.de
Received: 8 October 2014
Accepted: 8 January 2015
On the long nuclear time scale of stellar mainsequence evolution, even weak mixing processes can become relevant for redistributing chemical species in a star. We investigate a process of “differential heating”, which occurs when a temperature fluctuation propagates by radiative diffusion from the boundary of a convection zone into the adjacent radiative zone. The resulting perturbation of the hydrostatic equilibrium causes a flow that extends some distance from the convection zone. We study a simplified differentialheating problem with a static temperature fluctuation imposed on a solid boundary. The astrophysically relevant limit of a high Reynolds number and a low Péclet number (high thermal diffusivity) turns out to be interestingly nonintuitive. We derive a set of scaling relations for the stationary differential heating flow. A numerical method adapted to a high dynamic range in flow amplitude needed to detect weak flows is presented. Our twodimensional simulations show that the flow reaches a stationary state and confirm the analytic scaling relations. These imply that the flow speed drops abruptly to a negligible value at a finite height above the source of heating. We approximate the mixing rate due to the differential heating flow in a star by a heightdependent diffusion coefficient and show that this mixing extends about 4% of the pressure scale height above the convective core of a 10 M_{⊙} zeroage main sequence star.
Key words: convection / stars: evolution
Appendix A is available in electronic form at http://www.aanda.org
© ESO, 2015
1. Introduction
Our lack of understanding of (magneto)hydrodynamic transport processes in stars has hampered progress in developing the stellar evolution theory since its earliest beginnings. One particular aspect of the problem is the mixing in the boundary layers between convection and radiative zones in stellar interiors, which is also known as the problem of convective overshooting. Despite the indisputable advance in numerical simulations, the problem remains extremely challenging owing to the extreme range of the length and time scales involved in it.
The set of physical mechanisms that provide mixing at a convective/stable interface very likely depends on the type of convection zone involved. Because it is exposed to outer space at the top, a convective envelope is driven by the cold plumes originating in the photosphere. It is quite possible that the plumes span the whole convection zone and even provide mixing at its bottom boundary (cf. Andrássy & Spruit 2013, and references therein). A convective core or shell is, on the other hand, fully embedded in the star, its stratification is much weaker, and the temperature fluctuations within it are much smaller. Consequently, a different set of physical mechanisms may dominate mixing at its boundary.
It has long been known that the kinetic energy of the lowMachnumber flow in a convective core (or shell) is so low that the convective motions are stopped within about one per cent of the pressure scale height once they enter the steep entropy gradient of the radiative zone (Roxburgh 1965; Saslaw & Schwarzschild 1965). The motions can reach much farther, though, if they are vigorous enough to flatten the radiative entropy gradient above the convective core. In this case, known as the process convective penetration, the motions gradually “erode” the radiative stratification on the thermal time scale until radiative diffusion stops any further advance of the erosion front (Shaviv & Salpeter 1973; van Ballegooijen 1982; Zahn 1991). Finally, the fluid parcels hitting the stable stratification always generates a spectrum of internal gravity waves, which may also provide a certain amount of mixing (Press 1981; Garcia Lopez & Spruit 1991; Schatzman 1996).
Several of the processes mentioned may operate at the convection zone’s boundary at the same time. Their effects on long time scales and at long distances from the boundary are very different. In full numerical hydrodynamic simulations, the restrictions on time scales that can be covered makes it difficult to disentangle these effects. Physical insight developed by different means is needed to extrapolate them to longer time scales and distances.
We take a closer look at one specific process operating at a convective/stable interface in the interior of a star. Thermal diffusion causes temperature fluctuations from the convection zone’s boundary to spread into the stable stratification. Temperature differences on surfaces of constant pressure set up a flow even in the absence of momentum transport by hydrodynamic stress. We call this process “differential heating”, explore the physics of it in an idealised setup, and estimate what amount of mixing it could cause in the stellar interior.
2. The differential heating problem
2.1. Problem formulation and simplification
Consider a horizontal, solid surface with a stablystratified fluid overlying it^{1}. A temperature fluctuation imposed at the surface propagates into the fluid by a diffusive process and upsets the hydrostatic equilibrium. We investigate what the properties of the resulting flow are.
By replacing the convective/stable interface by a solid wall, we eliminate all the phenomena related to the inertia of the convective flows and the shear induced by them. This allows us to study the physics of differential heating in isolation. The upper boundary is taken far enough not to influence the flow. Next we introduce further assumptions to facilitate the mathematical description and the subsequent analysis of the problem:
 (1)
The flow is confined to a layer that is significantly thinner than the pressure scale height.
 (2)
The fluid is a chemically homogeneous, ideal gas.
 (3)
The BruntVäisälä (buoyancy) frequency of the stratification is constant.
 (4)
Thermal diffusivity is constant.
 (5)
The gravitational field is homogeneous.
 (6)
The differentially heated surface is flat and horizontal.
 (7)
The flow is constrained to two spatial dimensions.
The Boussinesq equations are (Spiegel & Veronis 1960)
where u is the fluid velocity, D/Dt = ∂/∂t + u·∇ the Lagrangian time derivative, ρ_{m} and T_{m} are the mean density and temperature, respectively, p′ and T′ the pressure and temperature perturbations, respectively, g is the gravitational acceleration, k a unit vector pointed in the vertical direction, ν the kinematic viscosity, N the BruntVäisälä frequency, w the vertical velocity component, and ϰ the thermal diffusivity.
Equations (1)–(3) still contain several dimensional parameters. It is crucial to realise that there is a natural system of units for the differential heating problem that makes the equations dimensionless. The flow in this problem is set off by thermal diffusion in a stably stratified medium, therefore the inverse of the BruntVäisälä frequency, 1 /N (or a multiple of it), is a natural unit of time. Having made this choice, we can define a natural unit of distance as , which is a typical thermaldiffusion length scale on the time scale 1 /N. The dimensionless Boussinesq equations are then
where we omit any symbol to indicate the new units. We have also introduced a new pressurelike variable p = p′/ρ_{m} and the buoyancy acceleration ϑ = gT′/T_{m}, which we continue to call the “temperature fluctuation” in the rest of the paper, because that is the central concept in the differential heating process. The Prandtl number Pr = ν/ϰ now becomes a measure of kinematic viscosity, because the new unit of diffusivity is ϰ. Equations (4)−(6) are particularly well suited to theoretical studies since their solution is fully determined by the Prandtl number, the initial, and the boundary conditions.
The distance unit is rather short in stellar interiors, and it only weakly depends on the stratification. To see this, we express the BruntVäisälä frequency in terms of the more common stellarstructure parameters, (7)where H_{p} is the pressure scale height, ∇_{ad} the adiabatic temperature gradient, and ∇ the actual temperature gradient. Close to a convection zone’s boundary, we can write (8)where α ≈ 10^{1} is a coefficient of proportionality and z the distance from the boundary (z> 0 in the stable stratification). When using Eqs. (7) and (8), the unit of distance can be expressed as (9)which is about 10^{7} cm for values typical of a point close to the convective/stable interface (z ≈ 10^{2}H_{p}) in the core of a massive (10 M_{⊙}), mainsequence star (ϰ≈ 10^{10} cm^{2} s^{1}, α ≈ 10^{1}, g ≈ 10^{5} cm s^{2}, H_{p} ≈ 10^{10} cm).
Two distinct regimes of differential heating can be expected, depending on the amplitude and the spatial scale of the temperature fluctuation imposed on the differentially heated surface. If the heating is strong enough, the heat transport is advectiondominated (i.e. the flow’s Péclet number is high), and the flow is generally unsteady. A similar phenomenon takes place right at the point where the convective flow leaves the unstable stratification, still retaining some positive temperature fluctuation. It quickly cools down as it rises in the stable medium, its temperature fluctuation turns negative, and the flow is brought to a halt. This is the place where we can impose a lower boundary condition for a much weaker kind of differentialheatinginduced flow, in which diffusive heat transport plays a major role (i.e. the flow’s Péclet number is low). The latter case is the main focus of this paper. We show in Sect. 3 that such a flow is generally smooth and reaches a stationary state (to be specified in Sect. 3.1) even at rather high values of the Reynolds number, up to Re = 4 × 10^{3}. This allows us to gain some insight into the problem by exploring the scaling properties of the stationary differentialheating equations, which we do in the next section.
2.2. Analytical considerations
The stationary differentialheating problem is described in two dimensions by the set of equations (cf. Eqs. (4)–(6))
where x and z are the horizontal and vertical coordinates, respectively, with the z axis pointed against the gravitational acceleration vector, u is the horizontal velocity component, and w the vertical one. In what follows, we show how the characteristic properties of the stationary flow depend on the typical amplitude Θ and the typical horizontal length scale L of the heating applied.
Assume that there is a welldefined vertical length scale H in the differential heating flow pattern. Let us denote the typical horizontal and vertical velocities by U and W, respectively, and the typical pressure fluctuation by P. We then introduce a new set of variables , ẑ, û, ŵ, , and , which all reach values of the order of unity close to the differentially heated surface, and Upon making these substitutions in Eq. (10), we obtain (20)which implies the approximate relation (21)The horizontal momentum equation (Eq. (11)) attains the form (22)where Eq. (21) has been used, so the equality is only approximate. Nonetheless, we can see that the viscous terms are of the order of 1/Re_{x} ≡ Pr/(UL) and 1/Re_{z} ≡ Pr/(WH), where Re_{x} and Re_{z} are Reynoldslike numbers associated with horizontal and vertical motions, respectively. We introduce this unusual notation to characterise the relative contributions of the two viscous terms in the case of L ≫ H. We focus on this limit because it turns out to be the relevant one in stellar interiors (see Sect. 5). From now on, we assume Re_{x} ≫ 1 and Re_{z} ≫ 1. Equation (22) shows that pressure fluctuations are of the order of U^{2} in this highReynoldsnumber limit, so that we can estimate (23)The vertical momentum equation (Eq. (12)), with the substitutions defined above and Eqs. (21) and (23), becomes (24)and implies a close balance between the vertical component of the pressure gradient and the buoyancyacceleration term provided that L ≫ H in addition to Re_{x} ≫ 1 and Re_{z} ≫ 1. This allows us to estimate (25)which is a plain, orderofmagnitude equality of the characteristic kinetic and potential energies. Finally, the energy equation (Eq. (13)) becomes (26)The diffusion terms in Eq. (26) are of the order of 1/Pe_{x} ≡ 1/(UL) and 1/Pe_{z} ≡ 1/(WH), where Pe_{x} and Pe_{z} are Pécletlike numbers associated with horizontal and vertical motions, respectively. We introduce them for the very same reason as we did in the case of Re_{x} and Re_{z}. Making use of Eqs. (21) and (25), we can put Eq. (26) into the form (27)which can be greatly simplified in the double limit of Pe_{x} ≫ Pe_{z} and Pe_{z} ≪ 1. In that case, the two terms in the parentheses on the righthand side have to closely balance one another, so that we can estimate (28)and Eq. (27) becomes linear, (29)Equation (29) is a special case of the energy equation in the lowPécletnumber approximation of Lignières (1999).
Using Eq. (28), we eliminate H from Eq. (25) to get an estimate of U(Θ,L) and, with Eq. (21), also an estimate of W(Θ,L). The resulting relations also enable us to express Re_{x}, Re_{z}, Pe_{x}, and Pe_{z} as functions of Θ, L, and Pr. This way we obtain
One might be tempted to estimate the time scale τ of flow acceleration towards the stationary state directly from the buoyancy acceleration Θ provided by the temperature fluctuation imposed on the bottom boundary. It is crucial to realise that, as Eq. (24) shows, the buoyancy acceleration is almost completely compensated for by the vertical component of the pressure gradient in the case L ≫ H. It is only their difference that contributes to the vertical acceleration. We can, however, consider the horizontal acceleration provided by the horizontal component of the pressure gradient and write U/τ ≈ P/L ≈ U^{2}/L (see Eq. (23)). Using Eq. (30) we obtain (36)Finally, we would like to point out that the characteristic thermaldiffusion length scale corresponding to the time scale τ is τ^{1/2} ≈ Θ^{− 2/7}L^{3/7}, which scales with Θ and L in quite a different way than H does (see Eq. (28)). This comes about because our estimates take the back reaction of the flow on the temperature distribution into account.
2.3. Numerical solutions
The orderofmagnitude estimates derived in the preceding section assume that the flow is stationary and that there is a welldefined vertical length scale in the flow pattern. We performed a series of timedependent, numerical simulations of the differential heating problem to confirm these assumptions and to determine how the solutions depend on the Reynolds number and how they decrease with height.
We have developed a specialised code dedicated to the study of the differentialheating problem, because the problem places rather high demands on the numerical scheme. For instance, it has to tackle the highly diffusive nature of the flow and its high aspect ratio and resolve a wide dynamic range within a single simulation box. The code is of the finitedifference type, and it solves the differentialheating equations on a collocated grid using a variant of the MacCormack integration scheme. The Poisson equation for pressure, which can be derived from Eqs. (4) and (5) (or (37), see below), is solved by a spectral method. Heatdiffusion terms are treated implicitly, again by a spectral method. In what follows, we discuss a few selected issues related to the numerical solution of the differentialheating equations that need to be borne in mind when interpreting our results. The reader interested in the details of the numerical scheme is referred to Appendix A.
We use periodic boundaries in the horizontal direction and force the shear stress and the vertical velocity component to vanish at the lower and upper boundaries of the computational domain. One could also use nonslip boundaries, but these are hardly more akin to the physical reality that motivated this study in the first place, so we omit this case. We impose a temperature fluctuation in the form ϑ(x,0) = Θsin(πx/L) at the bottom boundary and force the temperature fluctuation to vanish at the upper boundary. The parameters Θ and L can be identified with the same symbols as introduced in Sect. 2.2.
The high thermal diffusivity in the differentialheating problem forces us to use long implicit time steps for the heatdiffusion terms, which might have an adverse effect on the accuracy of the results. To show that this is not the case, we recomputed the simulations sr03, sr30, and Re1024 (see Tables 1 and 2 and Sect. 3), decreasing the time step by a factor of ten. This brings about a change in the velocity field, which is of the order of 0.1% in the cases sr03 and sr30 and of the order of 1% in the case of Re1024 (measured well away from the field’s zeroes). The reason for this insensitivity to the time step is the low Péclet number of the flow. Lignières (1999) shows that in the lowPe regime, the energy equation can be approximated by a Poisson equation for the temperature fluctuation with w as a source term (see also our Eq. (29)). We do not use this approximation to make our code more versatile; instead, we naturally obtain a close equilibrium between the terms ∇^{2}ϑ and w in Eq. (6) when the Péclet number is low. This equilibrium is reached so quickly that details of the evolution of ϑ towards the equilibrium become irrelevant.
Parameters of the series of simulations sampling a patch of the parameter space { Θ,L } at the constant value of Re = 2.6 × 10^{2}.
Parameters of the series of simulations with Re increasing at the fixed values of Θ = 10^{3} and L = 10^{0}.
It is a wellknown fact that any numerical advection scheme either requires adding a socalled artificialviscosity term to guarantee stability or it involves some viscous behaviour implicitly. In either case, the effective Reynolds number does not even come close to the astrophysical regime with current computing facilities. The artificial viscosity (be it explicit or implicit) thus exceeds the physical one by a wide margin, so it demands special attention.
Suppose we include an explicit viscous term as in Eq. (5) to model the artificial viscosity. Equations (32) and (33) show that for L ≈ 10^{3} (equivalent to ~H_{p} in the astrophysical case mentioned in Sect. 2.1) we have Re_{z} ≈ 10^{4} Re_{x} as a consequence of H ≪ L. Using equidistant grids with up to 10^{3} grid points in each direction, we can achieve Re_{x} ≈ 10^{3}. It follows that Re_{z} ≲ 10^{1} and the vertical momentum transport is dominated by the artificialviscosity term. A value Re_{z} ≫ 1 is, however, expected in stellar interiors. We use a simple workaround, replacing the viscous term Pr∇^{2}u by the anisotropic form Pr_{x}∂^{2}u/∂x^{2} + Pr_{z}∂^{2}u/∂z^{2}. The coefficients Pr_{x} and Pr_{z} are recomputed at each time step from the relations and , where h_{x} and h_{z} are the horizontal and vertical grid spacings, respectively, and Re_{grid} is the Reynolds number on the grid scale. We performed a few numerical tests of the code on a convection problem to determine that the value Re_{grid} = 4 is a conservative tradeoff between the amount of viscous dissipation and the code’s stability, so we use this value in all the simulations presented here.
The anisotropic form of artificial viscosity enables us to reach Re_{x} ≫ 1 and Re_{z} ≫ 1 at the same time on a grid of reasonable size. We show in Sect. 3 that the solutions with constant values of Pr_{x} and Pr_{z} decay exponentially with height. The effective, local Reynolds number decreases in proportion to the flow speed, and the solutions quickly become dominated by the artificial viscosity. This would also happen in an (otherwise idealised) stellar interior at some point but that point would be much farther from the convection zone’s boundary. Therefore, we generalise the artificialviscosity terms, and the momentum equation (Eq. (5)) in 2D becomes (37)where are Pr_{x}(z) and Pr_{z}(z) are proportional to e^{− ηz} with η being an adjustable parameter. We set η = 0 when we are not interested in the precise vertical profiles and use η> 0 to suppress the viscous terms when examining how the solutions decrease with height. The latter case, η> 0, is a rather touchy problem because one has to increase η in a few steps, always using the (almost) stationary flow from a previous run as an initial condition for the next run. Overestimating the value of η can lead to a lack of viscous dissipation in some parts of the computational domain and a numerical instability ensues. Finally, the very goal that we want to achieve by this treatment, i.e. the flow dynamics’ being dominated by inertial terms at great heights, becomes an issue since such a flow evolves on the extremely long time scale corresponding to its low speed.
3. Results
3.1. The stationarity and structure of the flow
Our numerical investigation of the diffusiondominated differential heating problem reveals that the flow reaches a stationary state at all values of the Reynolds number that we have been able to achieve (up to Re = 4 × 10^{3}). We use the rate of change of the quantity u_{max} = max  u  (taken over the whole simulation box) as a convergence monitor and a basis of our criterion for deciding the flow’s stationarity. We show in Sect. 2.2 that the relevant dynamical time scale near the differentially heated boundary should be close to τ given by Eq. (36) (confirmed a posteriori, see below). We pronounce the flow stationary and stop the simulation once the condition (38)has been met at least for one time scale τ. A direct implementation of this condition would involve extrapolation from the time scale of one time step, Δt, to a much longer time scale τ, which would amplify the roundoff noise by a factor of τ/ Δt ≫ 1. Instead, we approximate Eq. (38) by (39)where is the Eulertimestepped solution of the equation . Thus, is a smoothed version of u_{max}, lagging behind it approximately by τ in time.
The flow in all of the simulations presented in this paper is composed of several layers of overturning cells with flow speed rapidly decreasing from one layer to the next (see Figs. 4 or 6). We characterise the flow properties close to the differentially heated surface by a vertical length scale H, defined as the height above the hottest spot at which the flow first turns over (i.e. w(L/ 2,H) = 0), and the typical horizontal and vertical velocity components and , respectively, where the maxima are taken over the whole simulation box. The symbols H, U, and W can be identified with the same symbols as used in Sect. 2.2. The flow always reaches its maximal horizontal speed at the bottom boundary and the maximal vertical speed in the first overturning cell above the hot spot. The flow pattern is asymmetric, with the maximum downward flow speed (reached above the cold spot), max( − w), always lower than the maximal upward flow speed, max(w). We define the characteristic numbers Re_{x}, Re_{z}, Pe_{x}, and Pe_{z} in an analogous way to what is used in Sect. 2.2 with the difference that now we have two Prandtllike numbers Pr_{x} and Pr_{z} instead of one Prandtl number Pr.
Fig. 1 Dependence of the global flow characteristics on the heating amplitude Θ and length scale L. Circles show the values derived from numerical simulations (Table 1). Solid lines show the scaling relations (Eqs. (28), (30), (31), and (36)), normalised to fit all but the four simulations at L = 10^{1}, which are expected to deviate from the scaling relations. 

Open with DEXTER 
3.2. Scaling relations
We computed a grid of 16 simulations to verify the analytical relations derived in Sect. 2.2. All of these simulations, summarised in Table 1, have a resolution of 256 × 512, and the vertical grid spacing was adjusted so as to obtain Re_{x} = Re_{z} ≡ Re = 2.6 × 10^{2}. The decision to fix the value of Re is motivated by the fact that the flow pattern turns out to be scalable over a large part of the parameter space provided that Re = const. In other words, while changing the heating parameters Θ and L at Re = const. does change the amplitude and the vertical scale of the flow, the structure of the flow, as seen in a system of normalised coordinates x/L and z/H, remains unchanged (see Fig. 4). We show in Fig. 1 our numerical results as compared with the scaling relations fitted to all but the four data points at L = 10^{1}. The excluded data points do not comply well with the premise L ≫ H and are thus expected not to follow the scaling relations. Allowing only the constants of proportionality to change in the fitting process, we obtain The unexpectedly good fit is a result of the flow’s scalability.
Fig. 2 Dependence of the maximal horizontal velocity on the Reynolds number. Simulation data (circles) are connected by the solid line to guide the eye. The scaling law u_{max} ∝ Re^{0.054} is shown by the dashed line for comparison. 

Open with DEXTER 
The constants of proportionality in Eqs. (40)–(43), as well as the structure of the flow, depend on Re. We illustrate this dependence by computing a series of simulations with resolution ranging from 32^{2} to 4096^{2}. This way, we cover about two orders of magnitude in the Reynolds number from 3 × 10^{1} to 4 × 10^{3} (again with Re_{x} = Re_{z} ≡ Re). We perform this experiment at L = const. and Θ = const., so any change in Re reflects a change in Pr_{x} and Pr_{z}. Nevertheless, we present the dependence on Re, because the scalability of the flow shows that the absolute values of Pr_{x} and Pr_{z} do not matter. Ideally, we should choose the heating parameters so as to have Θ ≪ 1 and L ≫ 1 as the scaling relations hold true in this limit (see Sect. 2.2).
Equation (43) shows, however, that the flow’s dynamical time scale becomes extremely long in the same limit, thus making any highresolution computation unfeasible. Therefore we use Θ = 10^{3} and L = 10^{0}, which still keeps the energy equation approximately linear (since Pe_{x} ≈ Pe_{x} ≈ 10^{2}), but we forgo having L ≫ H. Nevertheless, we expect the changes in the flow with increasing Re in this case to be similar to those that would be seen in a simulation with L ≫ H because of the energy equation’s being linear in both cases. All of this series of simulations, summarised in Table 2, reach the stationary state as defined by Eq. (39). Figure 2 shows that the maximum horizontal velocity in the computational domain slowly increases in proportion to Re^{0.054} in the highRe regime. The flow also becomes increasingly asymmetric, as shown by the ratio of the maximum upward and downward flow speeds plotted as a function of Re in Fig. 3. The seemingly asymptotic trend changes at the highest Reynolds number considered, but we do not know the reason for this change.
3.3. Flow at great heights
The flow speed in all our simulations quickly decreases with height. Figure 5 compares the vertical profiles of the rootmeansquare (rms; computed in the x direction) vertical velocity component, w_{rms}(z), in four simulations with widely disparate heating parameters (sr00, sr03, sr30, and sr33). We find that w_{rms} decreases approximately as e^{−βwz/H} in a global sense with β_{w} ≐ 1.5 almost independently of Θ and L. We use the values H(Θ,L) given by Eq. (40) instead of those measured in the simulations to normalise the z coordinate, because this brings the slopes much closer to one another. We have to keep in mind, though, that these flows are reasonably close to a stationary state only up to z/H ≐ 2.5, because our convergence criterion (Eq. (39)) is ignorant of the weak flow in the upper part of the simulation box, and consequently, that part of the flow is still slowly evolving when the computation is stopped.
Fig. 3 Dependence on the Reynolds number of the flow’s asymmetry, characterised by the ratio of the maximum upward and downward flow speeds. 

Open with DEXTER 
Fig. 4 Comparison of the flow structure in two simulations with very disparate heating parameters. The left panel shows simulation sr30 (Θ = 10^{3}, L = 10^{1}). The right panel shows simulation sr03 (Θ = 10^{0}, L = 10^{4}). In both cases, the vertical velocity component w, normalised to its maximal absolute value, is plotted on a split logarithmic colour scale. The length of the velocity vectors (arrows) is scaled in a nonlinear way to aid visualisation. The spatial coordinates are normalised by the characteristic length scales defined in Sects. 2.3 and 3. 

Open with DEXTER 
The simulations discussed so far use constant artificialviscosity parameters Pr_{x} and Pr_{z}, which leads to a rapid decrease in the local Reynolds number with height (in proportion to the decreasing flow speed). We computed another two simulations, this time with Θ = 10^{4} and L = 10^{1}. In the first one, we set Pr_{x} = const. and Pr_{z} = const. (the constantPr case hereinafter), just as we have done so far. In the other one, we set Pr_{x} ∝ e^{− ηz} and Pr_{z} ∝ e^{− ηz} as described in Sect. 2.3 to keep a local version of the Reynolds number approximately constant (the constantRe case hereinafter). We increased the slope η from 0 in a few steps in order to make the ratio of the rms advection terms to the rms viscous terms, i.e. the local Reynolds number, as independent of height as possible; η = 2 turns out to be a good compromise in this case. There is a largescale, residual variation by about a factor of four in the local Reynolds number, because the simple exponential profile of the artificial viscosity is not flexible enough to compensate for it. Using Fig. 2 we estimate that this variation can change the velocities by ~0.1 dex at most. Since our usual stopping condition, Eq. (39), cannot “sense” the weak flow at great heights, we judge the stationarity of the flow by comparing the rms values of the ∂/∂t terms to the rms values of all other terms that appear in Eq. (37) and require the former to be significantly smaller than the latter. This way, we obtain the results summarised in Figs. 6−8. The constantPr flow can be considered stationary over the whole range shown, whereas the constantRe flow is only stationary for z/H ≲ 3.8, because the topmost part of that flow evolves so slowly that a global oscillation develops before it has reached equilibrium (see Sect. 3.4 for details).
Figure 7 illustrates that the flow is somewhat faster at z/H> 1 in the constantRe case, as could be expected from the massive increase in the local Reynolds number by as much as two orders of magnitude at z ≐ 3. Much more interesting is, however, that the overturning cells in the constantRe case become apparently thinner with increasing height, hence with decreasing local temperature fluctuation. This observation suggests that the scaling relations derived in Sect. 2.2 could be used locally (see the dependence of H on Θ in Eq. (28)). Another piece of evidence for this hypothesis is shown in Fig. 8, in which we compare the relative rates of decrease in ϑ_{rms}(z), u_{rms}(z), and w_{rms}(z). The envelope of ϑ_{rms}(z) can be approximated well by the function ϑ_{e}(z) ∝ e^{− βϑz/H} with β_{ϑ} = 1.7 for z/H ≲ 3. We then regard ϑ_{e}(z) as an estimate of the local temperature fluctuation and rewrite the scaling relations for the velocity components, Eqs. (30) and (31), to obtain their local versions, where u_{e}(z) and w_{e}(z) are expected to be good envelope models of u_{rms}(z) and w_{rms}(z).
Fig. 5 Decline with height of the relative rms vertical velocity in four simulations with very different heating parameters. The flows are reasonably close to a stationary state only up to z/H ≐ 2.5. 

Open with DEXTER 
In other words, we expect u_{e}(z) ∝ e^{− βuz/H} and w_{e}(z) ∝ e^{− βwz/H} with and . Indeed, these scalings turn out to be correct, as shown in Fig. 8. Similarly, we can produce a local version of Eq. (28), (46)where h(z) is a local, heightdependent estimate of a vertical length scale analogous to H. As a result, we expect h(z) ∝ e^{− βhz/H} with , i.e. a slow thinning of the overturning cells with increasing height, similar to what we observe in Figs. 6−8. This seemingly innocuous phenomenon has very grave consequences for the flow at great heights. Instead of fading out exponentially, it decreases even faster (see Fig. 8). We expand on this in Sect. 4.1 and derive a better model for the flow’s decline with height to show that the flow speed drops dramatically above a certain point.
Fig. 6 Effect of two different artificialviscosity prescriptions on the flow structure. The constantPr flow (Pr_{x}, Pr_{z} = const.; left panel) is compared with the constantRe flow (Pr_{x}, Pr_{z} ∝ e^{− 2z}; right panel). The vertical velocity component w is in both cases plotted on a splitlogarithmic colour scale. The length of the velocity vectors (arrows) is scaled in a nonlinear way to aid visualisation. 

Open with DEXTER 
3.4. Latetime evolution of the flow
Having continued some of our simulations for as much as 10^{4}τ, we discover an intriguing phenomenon. At first, a horizontal mean shear flow develops on top of the differentialheating flow. Its amplitude grows, and the shear flow begins to oscillate at some point. Finally, the oscillation saturates at an amplitude ranging from ~10^{3} to ~10^{0} of the differential heating flow’s amplitude, depending on the parameters of the simulation. The oscillation’s period and development time strongly decrease with increasing Reynolds number. They do not seem to have an upper limit but approach 10τ at Re ≈ 10^{3}. This phenomenon most likely has a physical origin because decreasing the time step by a factor of ten does not affect the shear flow or its behaviour significantly. Any detailed study of this phenomenon is certainly beyond the scope of this paper, but our preliminary research suggests that it is unlikely to be a cumulative effect induced by internal gravity waves since it (1) also occurs in very small computational boxes, in which all internalwave modes are overdamped by radiative diffusion; and (2) the temporal spectra of the average horizontal velocity are featureless at periods significantly shorter than that of the shear flow oscillation.
Fig. 7 Effect of two different artificialviscosity prescriptions on the rms vertical velocities. Note that the constantPr flow can be considered stationary over the whole range shown whereas the constantRe flow is stationary only up to z/H ≐ 3.8. 

Open with DEXTER 
Fig. 8 Comparison of the rms velocities and temperature fluctuations in the constantRe case with two models approximating their global behaviour. Solid lines show u_{rms} (top), w_{rms} (middle), and ϑ_{rms} (bottom). Dashed lines show the model, in which u_{e}(z) ∝ e^{− βuz/H}, w_{e}(z) ∝ e^{− βwz/H}, and ϑ_{e}(z) ∝ e^{− βϑz/H} with , , and β_{ϑ} = 1.7. Dotted lines show the improved model given by Eqs. (52)−(54) with γ = 1.3. The coefficients of proportionality have been adjusted for each variable independently. 

Open with DEXTER 
4. Interpretation of the results
4.1. Improving the model at great heights
After picking up the threads of Sect. 3.3, we presently find that the diffusiondominated, highRe differentialheating flow actually decreases faster than exponentially with height. To see this, we make use of two results from Sect. 3.3. First, that the scaling relations derived in Sect. 2.2 have their local analogues, which hold within the flow (compare Eqs. (28), (30), and (31) with Eqs. (46), (44), and (45), respectively). Second, that the envelope of ϑ_{rms}(z) can be approximated well by the function ϑ_{e}(z) ∝ e^{− βϑz/H} at low heights, where β_{ϑ} is independent of Θ and L. This allows us to write (47)The characteristic vertical scale H is linked to the heating amplitude Θ by Eq. (28) and is thus relevant close to the differentially heated surface, where the typical temperature fluctuations are of the order of Θ. A straightforward generalisation of Eq. (47) is obtained by replacing H by the local, heightdependent estimate h(z) given by Eq. (46). Upon doing so, we have (48)where we have introduced the new variables , z′ = z/H, and h′(z) = h(z) /H. By Eqs. (28) and (46) we have (49)and Eq. (48) becomes (50)where β may differ slightly from β_{ϑ}, because we have used an orderofmagnitude relation in the last step. Equation (50) shows that decreases with a fairly constant slope over a few orders of magnitude, but the slope starts to change as soon as a wider dynamic range is considered. Since the slope is proportional to , Eq. (50) describes a runaway process. Indeed, the solution is (51)and vanishes at a finite height of . The constant must be very close to unity as and we can simplify Eq. (51) to obtain (52)where we have returned to the nonprimed variables and introduced a new constant , which is a parameter to be adjusted to fit the numerical data. Using the local scaling relations, Eqs. (44)−(46), we derive the functional dependencies
The functions ϑ_{e}(z), u_{e}(z), and w_{e}(z) are shown in Fig. 8. The constants of proportionality in Eqs. (52)–(54) have been adjusted independently, but all three functions share the value γ = 1.3. The good fit indicates that our line of reasoning is probably correct.
Can we conclude that the flow stops at the finite height we have just derived? No, since the scaling relations only work in the highRe regime. Provided that Re is high close to z = 0, the flow speed quickly decreases according to Eqs. (53) and (54) until Re ≈ 1 is achieved at some height z_{1} < 7H/γ. The weak flow above this point is supported by viscosity and gradually vanishes as z → ∞.
4.2. Allowing for a buoyancyfrequency gradient
So far, we have assumed that the flow occurs in a particularly simple type of thermal stratification – one characterised by a typical buoyancy frequency N_{typ} = const. Nevertheless, we aim to apply our results to the immediate vicinity of a convection zone, i.e. to a medium, in that the buoyancy frequency rises continuously from zero to a finite value. In this section, we first show how to estimate the value of N_{typ} in such a setting and then reapply the techniques developed in Sect. 4.1 to demonstrate how the varying buoyancy frequency affects the global flow field.
To do this, we have to recover the dependence of all the relevant flow properties on N_{typ} by returning to a system of physical units. We recall that we use 1 /N_{typ} as a unit of time and (ϰ/N_{typ})^{1/2} as a unit of distance, which implies that the unit of velocity is (ϰN_{typ})^{1/2} and the unit of acceleration (hence of ϑ) is . We use these conversion factors throughout this section without mentioning them further. The height of the bottommost overturning cell is by Eq. (40) (56)where we have introduced the index “ph” to indicate the use of physical units for quantities that are dimensionless in the rest of our analysis. The buoyancy frequency N is now an increasing function of z_{ph} and can be approximated by Eqs. (7) and (8), (57)The overall scale of the flow pattern is given by the bottommost overturning cell, which is thus the most important. Therefore we estimate N_{typ} = N(H_{ph}/ 2), i.e. (58)and combine Eqs. (56) and (58) to obtain (59)where we have also expanded Θ_{ph} = gΔT/T_{m} to emphasise the dependence on the imposed temperature fluctuation ΔT/T_{m}. We use the sign ≐ in Eq. (59) and also in Eqs. (60)−(62) below to indicate that we do not expect these estimates to be off by more than a few tens of percent. The dependence of H_{ph} on the heating amplitude and length scale is somewhat weaker in Eq. (59) compared with Eq. (56), because Eq. (59) takes into account that any gain in the flow’s vertical extent brings about an increase in the typical buoyancy frequency, which in turn makes further penetration harder. This effect can also be seen when we express the characteristic velocity components and the flow’s dynamical time scale in physical units, where the exponents have slightly changed compared with Eqs. (41)–(43).
The spatial variation of N brings on a firstorder effect, too; that is, the stratification offers less resistance to overturning in the bottom part of the flow field compared with the rest of it. We mimic this effect by using the flow’s excellent scaling properties under the assumption that the flow behaves locally as if N was constant. Our goal is to improve upon the envelope models of ϑ_{rms}(z), u_{rms}(z), and w_{rms}(z) derived in Sect. 4.1 by taking the dependence of N on height into account.
Our starting point is Eq. (48) with the difference that now we define , z′ = z_{ph}/H_{ph} and h′ = h_{ph}/H_{ph}. We caution the reader that ϑ_{e,ph} refers to a model with N = N(z) and not to a direct translation of ϑ_{e} that appears in Eq. (48) to physical units. The local vertical length scale of the flow, h(z) given by Eq. (46), can be translated to physical units directly, (63)This equation, together with Eq. (56), implies (64)where N′ = N/N_{typ} = (2z′)^{1/2} (see Eqs. (57) and (58)). It is evident that h′ diverges for N′ → 0^{+}, i.e. z′ → 0^{+}. This effect is purely artificial because the divergence occurs within the bottommost overturning cell of the flow, and the largescale model we are developing here cannot capture such local phenomena. We ignore the divergence for now because only h^{′ − 1} appears in Eq. (48) and use the same procedure as in Sect. 4.1 to derive a generalised version of Eq. (50), (65)where the parameter β has absorbed all coefficients of the order of unity. Its value should still be of the order of unity, but it may be different in this model compared with the model developed in Sect. 4.1. By analogy to the derivation in Sect. 4.1, we can write the solution to Eq. (65) in the form (66)where we have also returned to the nonprimed quantities, and is a parameter of the order of unity. The typical velocity components and the typical vertical vertical length scale can be estimated using the local scaling relations, Eqs. (44)−(46), and (66). We obtain where an explicit dependence on N appears after the transition to physical units. These expressions diverge for z → 0^{+} where N → 0^{+} (see Eq. (57)), which is just another illustration of the envelope models’ inability to capture local phenomena (see also the discussion above). The bottommost part of the flow should in reality behave approximately as if it was in a medium with N = N_{typ} = const., so we can cut off the problematic part of the N(z) profile and use, for example, the function (70)instead of N(z) in practical calculations. Doing so makes the righthand sides of Eqs. (67)−(69) converge to unity as z_{ph} → 0^{+}.
Just as the results of Sect. 4.1 do not mean that the flow vanishes at a finite height, neither the results of this section mean that. Again, the sudden drop in the typical velocities predicted by Eqs. (67) and (68) only signifies that the flow undergoes a transition to the lowRe regime at a relatively low height. Equations (67) and (68) cease to be usable from that point on and the weak flow supported by viscosity gradually vanishes as z_{ph} → ∞.
5. Application to stellar conditions
The flow in a layer of thickness h_{ph}(z_{ph}) and vertical velocity w_{e,ph}(z_{ph}) at distance z_{ph} from the boundary overturns a passive tracer in it on a time scale τ_{m} = h_{ph}/w_{e,ph}. This suggests an effective diffusion coefficient D_{eff} ≈ h_{ph}w_{e,ph}. For the first layer above the boundary, this is (71)At distance z_{ph}, Eqs. (68) and (69) give (72)where N_{typ} is given by Eq. (58) and we have replaced N(z_{ph}) in Eqs. (68) and (69) by given by Eq. (70) as discussed in Sect. 4.2. The constant γ is of the order of unity but cannot be constrained further by our present analysis. It determines the maximum height z_{max,ph} that the mixing process can reach, z_{max,ph} = (9 /γ)^{7/9}H_{ph}.
For a specific example, consider the boundary of the core convection zone in a 10 M_{⊙} zero age main sequence star. This environment is characterised by α = d(∇_{ad} − ∇)/d(z_{ph}/H_{p}) = 0.14, a thermal diffusivity ϰ= 5.9 × 10^{10} cm^{2} s^{1}, a gravitational acceleration g = 1.1 × 10^{5} cm s^{2} and a pressure scale height H_{p} = 2.9 × 10^{10} cm. A mixinglength estimate for convection in the core produces temperature fluctuations ΔT/T_{m} ≈ 10^{6} on a horizontal length scale L_{ph} ≈ H_{p}. Equation (59) then predicts that the typical vertical length scale is H_{ph} ≈ 2 × 10^{8} cm = 7 × 10^{3}H_{p}. The typical vertical velocity (Eq. (61)) is W_{ph} ≈ 5 × 10^{1} cm s^{1}. These numbers imply Pe_{z} = (W_{ph}H_{ph})/ϰ ≈ 2; i.e., the bottom part of the flow is located right at the transition between the regions of advectiondominated and diffusiondominated heat transport. This is not a coincidence, because we are modelling the region where heat leaks from the convective eddies, allowing them to turn over and sink back to the convection zone. Such a flow has to have Pe_{z} ≈ 1. Therefore, the effective diffusivity close to the convection zone, D_{eff}(0) in Eq. (72), is of the same order as the diffusivity of heat ϰ. Diffusivities that are several orders of magnitude smaller than ϰ can be important on the long nuclear time scale. The maximum height reached by the differential heating process on this time scale can thus be approximated by z_{max,ph}. Assuming γ = 1 we obtain z_{max,ph} ≈ 4 × 10^{2}H_{p}.
Equation (72) is likely to be somewhat of an overestimate of the actual mixing rate of the differentialheating process. The layers mix on the hydrodynamic time scale in their interiors, but as long as they are stationary, transport of the tracer between layers takes place by diffusion. As in the case of semiconvective layering (cf. Spruit 2013), this reduces the effective mixing rate to the geometric mean of the microscopic diffusion coefficient κ_{t} of the tracer and the estimate (72).
More significantly, the picture is complicated by the time dependence of the convective heat source. For the 10 M_{⊙} example, only the bottommost part of the flow can approach the stationary flow speed before the heating pattern changes because the dynamical time scale τ_{ph} ≈ 5 × 10^{6} s (Eq. (62) with the parameter values stated above) is of the same order as the convective overturning time scale in the core. This is likely to lead to some form of averaging, reducing the effective amplitude of the source. The level of this effect can probably be investigated with a timedependent simulation.
6. Summary
Various observations show that there is a need for some additional mixing at the interfaces between the convective and radiative layers of stars. Even processes that are too weak to be detectable in numerical hydrodynamic simulations need to be considered as candidate sources of this mixing, because the nuclear time scale on the main sequence is so much longer than the dynamical time scale of convection, and cumulative effects are likely to play an important role.
In this work, we have investigated one such weak process, which we call “differential heating”. The differential heating process occurs when radiative diffusion transports a temperature fluctuation from the boundary of a convection zone into the neighbouring stable stratification. The resulting perturbation of hydrostatic equilibrium triggers a weak flow, which may provide mixing up to some distance from the convection zone. We investigated the flow that is driven by a static temperature fluctuation varying sinusoidally along the solid horizontal boundary of a stably stratified, thin layer of gas. This lowPéclet number problem (i.e. a slow flow dominated by thermal diffusion) turns out to be intrinsically nonlinear, in the sense that the horizontal structure of the flow is asymmetric. Even for symmetric boundary conditions, the upflow is narrower than the downflowing part for the flow, and the shape of the flow pattern is nearly independent of the amplitude of the driving temperature perturbation.
A few additional assumptions (Sect. 2) allow us to describe the problem by a set of dimensionless equations, the solution to which depends (apart from the boundary and initial conditions) only on the Prandtl number. We analysed these differentialheating equations for their scaling properties under the assumption that the flow is stationary (Sect. 2.2). An astrophysically interesting corner of the parameter space is characterised by Re_{x} ≫ 1, Re_{z} ≫ 1, Pe_{x} ≫ Pe_{z}, and Pe_{z} ≪ 1. (The x and z directions have to be distinguished because such flow has a high aspect ratio.) In this limit we derive a set of simple relations (Eqs. (28) and (30)–(36)) to describe how the global flow properties depend on the heating amplitude Θ and length scale L. We find, in particular, that the characteristic vertical length scale H depends only weakly on the heating parameters (Eq. (28)).
We developed a dedicated numerical code to solve the equations. The main difficulties are related to the highly diffusive nature of the flow, its high aspect ratio, and the need to resolve a wide dynamic range in the flow amplitude within the computational box (as much as five orders of magnitude). The flow in our twodimensional, timedependent simulations reaches a stationary state at all values of the Reynolds number that we have been able to achieve (up to Re ≡ Re_{x} = Re_{z} = 4 × 10^{3}). The flow is always composed of several layers of overturning cells, the shape of which depends only on the Reynolds number and not on the heating length scale L and amplitude Θ. This property makes the flow scaleable in the sense that the flow field corresponding to some heating parameters L_{1}, Θ_{1} can be stretched in space and scaled in amplitude to get a good approximation of the flow field corresponding to a different set of heating parameters L_{2}, Θ_{2} provided that Re is in both cases the same. This is also the reason the scaling relations derived in Eq. (2.2) fit the simulation data remarkably well at Re = const. (see Fig. 1). Increasing the Reynolds number has little influence on the flow speed, but it makes the flow pattern increasingly asymmetric.
We decrease the artificialviscosity coefficients in the code with height in order to keep the Reynolds number approximately the same in every layer of flow cells. The numerical data show that the global scaling relations derived in Sect. 2.2 have their local analogues, which can be used within the flow. The flow speed’s decrease with height, being locally exponential, steepens with the decreasing flow amplitude according to the local scaling relations. Based on this we derive a model of the flow’s dependence on height, which closely fits the numerical data over the whole dynamic range that we have been able to cover (as much as five orders of magnitude, see Fig. 8). The model shows that the flow speed drops abruptly to a negligible value at a finite height. The local scaling relations also allow us to generalise our results to the more realistic case, in which the buoyancy frequency N increases with height (see Sect. 4.1).
We illustrated the typical scales associated with the stellar differentialheating process with the example of the convective core of a 10 M_{⊙} zeroage main sequence star (see Sect. 5). We approximate the mixing due to the differentialheating flow by an “effective” diffusion coefficient D_{eff}, which is of the order of the diffusivity of heat near the convection zone and decreases with height according to Eq. (72). The mixing relevant for stellar evolution extends about 4% of the pressure scale height above the convection zone.
The main findings of the paper are:

(1)
The flow has a cellular structure and reachesa stationary state at all values of the Reynoldsnumber that we have been able to achieve (up toRe = 4 × 10^{3}).

(2)
Both global and local properties of the flow can be described by a set of simple analytical relations.

(3)
The flow speed drops abruptly to a negligible value at a finite height above the source of heating.

(4)
The mixing relevant for stellar evolution extends about 4% of the pressure scale height above the convection zone of a 10 M_{⊙} zeroage main sequence star.
Acknowledgments
We would like to thank Achim Weiss for the 10 M_{⊙} model, Ewald Müller and Maxime Viallet for enlightening discussions on numerical hydrodynamics, and the anonymous referee for critical comments that improved the overall presentation of the text.
References
 Andrássy, R., & Spruit, H. C. 2013, A&A, 559, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Frigo, M., & Johnson, S. G. 2005, Proc. IEEE, 93, 216 [CrossRef] [Google Scholar]
 Garcia Lopez, R. J., & Spruit, H. C. 1991, ApJ, 377, 268 [NASA ADS] [CrossRef] [Google Scholar]
 Lignières, F. 1999, A&A, 348, 933 [NASA ADS] [Google Scholar]
 Press, W. H. 1981, ApJ, 245, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Roxburgh, I. W. 1965, MNRAS, 130, 223 [NASA ADS] [Google Scholar]
 Saslaw, W. C., & Schwarzschild, M. 1965, ApJ, 142, 1468 [NASA ADS] [CrossRef] [Google Scholar]
 Schatzman, E. 1996, J. Fluid Mech., 322, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Shaviv, G., & Salpeter, E. E. 1973, ApJ, 184, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Spiegel, E. A., & Veronis, G. 1960, ApJ, 131, 442 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Spruit, H. C. 2013, A&A, 552, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Ballegooijen, A. A. 1982, A&A, 113, 99 [NASA ADS] [Google Scholar]
 Zahn, J.P. 1991, A&A, 252, 179 [NASA ADS] [Google Scholar]
Online material
Appendix A: Numerical methods
Appendix A.1: Integration scheme
We have adapted the standard MacCormack method to suit our specific problem. In the simplest case of a onedimensional vector q of conserved quantities being advected on and equidistant grid with a spacing of Δx, MacCormack’s method can be written as where is the value of q at the kth grid point and the nth time step, Δt the time step, f(q) the flux function, and we use the convention that any parenthesised upper index refers to a substep of the method instead of a timestep index. The method is linearly stable provided that the CFL condition Δt ≤ Δx/ρ(A) is met, where A is the Jacobian matrix of the flux vector and ρ(A) is the largest characteristic value of A. Nonlinear stability typically requires the addition of some form of artificial viscosity. MacCormack’s method is secondorder accurate both in space and time.
We discretise Eqs. (4), (6), and (37) on a collocated, twodimensional grid of M × N cells with constant cell spacing (Δx, Δz). The two spatial dimensions and the presence of source terms in the equations forces us to significantly extend the basic MacCormack scheme. We begin by advecting the vector of variables q = (u,w,ϑ) in both spatial directions using Strang splitting, where we have written out the explicit form of the flux terms. The indices k and l refer to the position along the x and z axes, respectively. We proceed by adding the source terms to the momentum equations, where we use secondorderaccurate central differences to keep up with the order of accuracy of the advection scheme, ν_{l} = Pr_{x}(z_{l}) and μ_{l} = Pr_{z}(z_{l}) are the coefficients of our anisotropic artificialviscosity prescription (see Sect. 2.3), and μ_{l + 1/2} = (μ_{l} + μ_{l + 1})/2. The new velocity field is, in general, slightly divergent. We correct for this divergence by subtracting the gradient of a pressurecorrection field, u^{(1)} = u^{(1d)} − Δt∇(Δp)^{(1)}. The condition ∇·u^{(1)} = 0 leads to a Poisson equation for the pressure correction, (A.9)Since we use central differences to compute partial derivatives, the discrete form of the Laplace operator in Eq. (A.9) should be derived by applying the central differences twice. That would, however, lead to a sparse operator and cause oddevendecoupling problems on our collocated grid. Therefore we use the standard compact Laplacian and solve the approximate pressurecorrection equation (A.10)Equation (A.10) is solved by a spectral solver, see Sect. A.3. Having computed the pressure correction, we apply it to the velocity field, The approximate nature of the pressurecorrection equation (Eq. (A.10)) causes ∇·u^{(1)} to be small, but nonzero. Practical experience has shown that the residual divergence is negligibly small in the flows analysed in this paper provided that the boundary conditions are treated properly, see Sect. A.2. We should also write at this point, but our numerical tests have shown that the residual divergence in the velocity field becomes much smaller if we set , so we use the latter form. The next step is to integrate the remaining two terms in the energy equation. We begin by adding the −w term, (A.13)where its latest available value, −w^{(1)}, has been used. The diffusion substep is given by the implicit equation (A.14)which is also solved by a spectral solver, see Sect. A.3. We have thus completed the first step of the MacCormack scheme, analogous to Eq. (A.1), and obtained the new variables u^{(1)}, w^{(1)}, p^{(1)}, and ϑ^{(1)}. The second step, which we do not do not go into detail on, differs from the first one at two points. First, advection is done using backwardspace flux differencing, as in Eq. (A.2) (compare with Eq. (A.1)). Second, the pressure field is updated in this step, i.e. . The final step of the MacCormack’s scheme, Eq. (A.3), is used in the same form, with q = (u,w,ϑ). We also update the pressure field in the same way, , so that we obtain an estimate of the pressure field for the next time step.
Fig. A.1 Effect of three different methods of treating the pressure at the solid top and bottom boundaries. In all three panels, the vertical velocity component, w, is plotted on a split logarithmic colour scale. We use Θ = 10^{3}, L = 10^{1}, constant kinematic viscosity and set the resolution to only 16 × 64 to make the spurious oscillations visible. We obtain the result plotted in the left panel using the simple symmetry conditions for pressure (Eqs. (A.19) and (A.20)). Preceding the pressuregradient computation by thirdorder pressure extrapolation to the ghost cells reduces the oscillations’ amplitude by a factor of ~100 (middle panel). Increasing the extrapolation order to six brings about another decrease by a factor of ~30 in the oscillations’ amplitude (right panel). The pressure gradient is in all three cases computed by the secondorder central differences in the whole computational domain. 

Open with DEXTER 
Finally, there is a simple way of increasing the accuracy of the scheme at a given grid resolution, which we use. The MacCormack method contains a builtin asymmetry: Eqs. (A.1) and (A.2) show that it always starts with forwardspace flux differencing and continues with backwardspace flux differencing. The two fluxdifferencing methods can be reversed, obtaining a “reverse” MacCormack method, without decreasing the order of accuracy of the overall scheme. We compute every time step using both the “direct” and the “reverse” methods and use the arithmetic average of the estimates given by the two methods.
Appendix A.2: Boundary conditions
The treatment of boundaries is restricted by our decision to use spectral solvers, which do not allow changing the differentiation operators anywhere in the computation domain. We use the ghostcell technique for this reason. The boundary conditions we impose on the differentialheating flow are summarised in Sect. 2.3. The periodic boundaries in the horizontal direction are trivial to implement. The solid boundaries on the top and bottom of the computational domain, however, require much more care. We implement them using reflective boundary conditions for the velocity vector, so that the imaginary walls are located at l = −1/2 and at l = N − 1/2. The conditions imposed on u also eliminate any shear on the boundary. The pressure field is required to be symmetric with respect to the solid boundaries, The conditions imposed by Eqs. (A.15)–(A.20) can easily be shown to be consistent with the pressurecorrection equation (Eq. (A.10); sum both sides over k = 0, 1,...,M and l = 0, 1,...,N). They typically do, however, bring about a cusp in the pressure field along the normal to the walls. The resulting discontinuity in the vertical pressure gradient then propagates to the rest of the domain and can be seen as a lowamplitude oscillatory field superimposed on the true pressure field (see the left panel of Fig. A.1). We tried to cure this problem by changing the discretisation of the verticalgradient operator at the walls, so that the ghost cells would not be used when computing the pressure gradient. This solution has met with very little success, most likely because the abrupt change in the operator brings about an abrupt change in the discretisation error so the problem remains. Quite surprisingly, preceding the pressuregradient computation by highorder pressure extrapolation to the ghost cells has turned out to be an effective solution, able to eliminate nearly all of the spurious oscillations (see the middle and right panels of Fig. A.1). We therefore use sixthorder extrapolation in the simulations with constant artificial viscosity and increase the extrapolation order to ten when we let the artificial viscosity decrease with height. This technique cannot be viewed, however, as an allpurpose solution, because it is likely to be too unstable to be useful when computing highly turbulent flows.
We require the temperature fluctuation ϑ to have a fixed sinusoidal profile at the bottom boundary and to vanish at the upper boundary, which translates into
Appendix A.3: Spectral solvers
We use spectral methods to solve the two equations involving the Laplace operator, the Poisson equation for the pressurecorrection equation (Eq. (A.10)) and the implicit heatdiffusion equation (Eq. (A.14)). We express both the knowns and unknowns as linear combinations of the Laplacian’s eigenfunctions that comply with the desired boundary conditions. The solution procedure is then much simplified and effective, provided that the transform to the eigenfunction basis can be computed efficiently.
In case of the pressurecorrection equation (Eq. (A.10)), we use the linear transform (A.23)and its inverse (A.24)to transform any field f_{k,l} to an array of complex amplitudes and back. We can see that the basis functions in Eq. (A.24) are periodic in k and even around l = −1/2 and l = N − 1/2; i.e., they comply with our boundary conditions on the pressure field (see Sect. A.2). Upon using the spectral decomposition defined by Eq. (A.24) on both sides of the pressurecorrection equation (Eq. (A.10)), we readily obtain its solution in the wavenumber space, (A.25)where we have omitted the upper indices because the expression applies to both steps of the MacCormack scheme, Ŝ_{k,l} is the transformed righthand side of Eq. (A.10). The eigenvalues λ_{m,n} of the Laplacian are (A.26)and can be precomputed. We set λ_{0,0} to a large number to prevent division by zero and make the undetermined component vanish.
In case of the heatdiffusion equation (Eq. (A.14)), we use the linear transform (A.27)\newpage
and its inverse (A.28)to transform any field g_{k,l} to an array of complex amplitudes ĝ_{m,n} and back. We can see that the basis functions in Eq. (A.28) are periodic in k and odd around l = −1/2 and l = N − 1/2; i.e., they comply with our boundary conditions on the temperature field in case of a vanishing heating amplitude (see Sect. A.2). To allow for an arbitrary heating profile at the bottom boundary, we take out the known boundary term from the Laplacian on the righthand side of Eq. (A.14) and treat it as a source term. One can show that it is the same as replacing the diffusion equation ∂ϑ/∂t = ∇^{2}ϑ by the equivalent equation ∂(ϑ − ζ) /∂t = ∇^{2}(ϑ − ζ), where ζ is the static solution to the diffusion equation ∂ζ/∂t = ∇^{2}ζ with the desired boundary conditions (ζ can be precomputed for a fixed heating profile). The boundary conditions on the difference ϑ − ζ are then identically zero, and the spectral decomposition defined by Eq. (A.28) can be used. This way we obtain an explicit expression for the solution of the implicit Eq. (A.14) in the wavenumber space, (A.29)where the eigenvalues Λ_{m,n} of the Laplacian are (A.30)and can be precomputed. An equation analogous to Eq. (A.29) relates to .
In the practical implementation, we use the FFTW library (Frigo & Johnson 2005) to compute the transforms in Eqs. (A.23), (A.24), (A.27), and (A.28). We combine standard, onedimensional transforms of different kinds to obtain the nonstandard, twodimensional transforms that we need. Namely, Eq. (A.23) is implemented as a series of DCTII transforms over the rows of the input array, after which the columns of the resulting array are transformed by a series of DTF transforms. The backward transform (Eq. (A.24)) is then computed by a series of DFTs followed by a series of DCTIIIs. The transforms for the diffusion equation (Eqs. (A.27) and (A.28)) are implemented in the same way, but simply replacing the DCTIIs by DSTIIs and DCTIIIs by DSTIIIs. The transforms from the FFTW library do not include the normalisation factor (2MN)^{1}.
All Tables
Parameters of the series of simulations sampling a patch of the parameter space { Θ,L } at the constant value of Re = 2.6 × 10^{2}.
Parameters of the series of simulations with Re increasing at the fixed values of Θ = 10^{3} and L = 10^{0}.
All Figures
Fig. 1 Dependence of the global flow characteristics on the heating amplitude Θ and length scale L. Circles show the values derived from numerical simulations (Table 1). Solid lines show the scaling relations (Eqs. (28), (30), (31), and (36)), normalised to fit all but the four simulations at L = 10^{1}, which are expected to deviate from the scaling relations. 

Open with DEXTER  
In the text 
Fig. 2 Dependence of the maximal horizontal velocity on the Reynolds number. Simulation data (circles) are connected by the solid line to guide the eye. The scaling law u_{max} ∝ Re^{0.054} is shown by the dashed line for comparison. 

Open with DEXTER  
In the text 
Fig. 3 Dependence on the Reynolds number of the flow’s asymmetry, characterised by the ratio of the maximum upward and downward flow speeds. 

Open with DEXTER  
In the text 
Fig. 4 Comparison of the flow structure in two simulations with very disparate heating parameters. The left panel shows simulation sr30 (Θ = 10^{3}, L = 10^{1}). The right panel shows simulation sr03 (Θ = 10^{0}, L = 10^{4}). In both cases, the vertical velocity component w, normalised to its maximal absolute value, is plotted on a split logarithmic colour scale. The length of the velocity vectors (arrows) is scaled in a nonlinear way to aid visualisation. The spatial coordinates are normalised by the characteristic length scales defined in Sects. 2.3 and 3. 

Open with DEXTER  
In the text 
Fig. 5 Decline with height of the relative rms vertical velocity in four simulations with very different heating parameters. The flows are reasonably close to a stationary state only up to z/H ≐ 2.5. 

Open with DEXTER  
In the text 
Fig. 6 Effect of two different artificialviscosity prescriptions on the flow structure. The constantPr flow (Pr_{x}, Pr_{z} = const.; left panel) is compared with the constantRe flow (Pr_{x}, Pr_{z} ∝ e^{− 2z}; right panel). The vertical velocity component w is in both cases plotted on a splitlogarithmic colour scale. The length of the velocity vectors (arrows) is scaled in a nonlinear way to aid visualisation. 

Open with DEXTER  
In the text 
Fig. 7 Effect of two different artificialviscosity prescriptions on the rms vertical velocities. Note that the constantPr flow can be considered stationary over the whole range shown whereas the constantRe flow is stationary only up to z/H ≐ 3.8. 

Open with DEXTER  
In the text 
Fig. 8 Comparison of the rms velocities and temperature fluctuations in the constantRe case with two models approximating their global behaviour. Solid lines show u_{rms} (top), w_{rms} (middle), and ϑ_{rms} (bottom). Dashed lines show the model, in which u_{e}(z) ∝ e^{− βuz/H}, w_{e}(z) ∝ e^{− βwz/H}, and ϑ_{e}(z) ∝ e^{− βϑz/H} with , , and β_{ϑ} = 1.7. Dotted lines show the improved model given by Eqs. (52)−(54) with γ = 1.3. The coefficients of proportionality have been adjusted for each variable independently. 

Open with DEXTER  
In the text 
Fig. A.1 Effect of three different methods of treating the pressure at the solid top and bottom boundaries. In all three panels, the vertical velocity component, w, is plotted on a split logarithmic colour scale. We use Θ = 10^{3}, L = 10^{1}, constant kinematic viscosity and set the resolution to only 16 × 64 to make the spurious oscillations visible. We obtain the result plotted in the left panel using the simple symmetry conditions for pressure (Eqs. (A.19) and (A.20)). Preceding the pressuregradient computation by thirdorder pressure extrapolation to the ghost cells reduces the oscillations’ amplitude by a factor of ~100 (middle panel). Increasing the extrapolation order to six brings about another decrease by a factor of ~30 in the oscillations’ amplitude (right panel). The pressure gradient is in all three cases computed by the secondorder central differences in the whole computational domain. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.