Open Access
Volume 667, November 2022
Article Number A96
Number of page(s) 17
Section Stellar structure and evolution
Published online 11 November 2022

© F. Kupka et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1. Introduction

From the early work of Biermann (1932) onwards, research on convection has remained a major challenge in stellar astrophysics, in particular as convection turned out to be one of the most important mechanisms of energy transport and mixing in stars. As is described, for instance, in Canuto et al. (2009), when we compare the spatial scales of viscous processes derived from the results of Chapman (1954) on fully ionised gases with the spatial scales of convective flow observed at stellar surfaces (Kupka & Muthsam 2017), stellar convection is characterised by very high Reynolds numbers. Stellar convective flows are thus highly turbulent, even though the direct detection of turbulence is difficult due to the nature and resolution of the observational methods available to us (cf. Kupka & Muthsam 2017).

Modelling this class of flows poses serious challenges for stellar structure and evolution models (for introductions and reviews see, e.g. Weiss et al. 2004; Canuto et al. 2009; Kupka & Muthsam 2017; Kupka et al. 2020). Due to the extreme range of scales in space and time, numerical hydrodynamical simulations cannot be used directly in stellar evolution calculations (cf. the estimates given in Kupka & Muthsam 2017). Consequently, turbulent convection has to be modelled in a framework affordable for direct coupling into one-dimensional (1D) stellar models of stellar evolution. The turbulent convection models (TCMs) used in this approach differ widely in computational costs, physical completeness, and general principles considered in their derivation, from completely phenomenological to more systematic approaches based on turbulence theory (see Kupka & Muthsam 2017 for an overview).

One methodological way to derive TCM equations that are suitable for stellar evolution calculations is the Reynolds stress approach. The splitting of variables in turbulent flow into a mean and a fluctuating component was first introduced by Reynolds (1894), followed by the suggestion of Keller & Friedmann (1925) to consider this Reynolds splitting for a moment expansion approach that was first completed by Chou (1945). Dynamical variables such as velocity v, density ρ, or entropy s, for example, can be subject to such splitting:

Strictly speaking these are ‘ensemble averages’ over different initial conditions. In practice, the variables are also subject to spatial averaging, in 1D stellar models typically over the θ and ϕ directions, to which the overbar in the above notation refers to whereas the component with a prime refers to the fluctuating part of each quantity.

Due to their immediate physical meaning, the higher order combinations of the fluctuating parts which appear in such Reynolds stress models of turbulent convection are also model predictions of direct astrophysical interest. The second order moment of velocity fluctuations, characterising the turbulent kinetic energy (TKE) of the convective flow, is directly related to the highly efficient chemical mixing induced by convection. In stars with nuclear burning in convective cores, this has a direct impact on the luminosity and the lifetime of the nuclear burning phase. Similarly, the second order moment of velocity and entropy fluctuations, related to the convective flux, determines the energy transported by convection. Computing the convective flux allows one to predict the temperature gradient in convective regions. Recently, the temperature gradient in core boundary layers of an intermediate-mass main-sequence star was probed using asteroseismology (Michielsen et al. 2021), an observation that can directly be compared to results from a TCM.

Presently, the most commonly used theory to describe convection in stellar structure and evolution models is still the mixing length theory (Böhm-Vitense 1958, MLT). However, MLT is not able to describe the convective boundary in a physically accurate way. Observations have shown that chemical mixing beyond the boundary of convectively unstable regions, commonly known as overshooting, is required (see, for example, Maeder & Mermilliod 1981; Bressan et al. 1981; Pietrinferni et al. 2004). In stellar models using MLT, parametrised ad hoc mixing beyond the boundary is introduced to achieve this. Likewise, the temperature structure of an overshooting region cannot be predicted by MLT. These examples highlight the need for more physical theories of convection, such as TCM, being included in stellar structure and evolution models.

A large number of TCM have been developed (Xiong et al. 1997; Canuto 1992, 1993; Canuto & Dubovikov 1998; Li & Yang 2001, 2007; Kuhfuß 1986, 1987) which differ in the set of variables used and the set of approximations and assumptions made (see Canuto 1993 and Kupka & Muthsam 2017 for comparisons and a review). Among other physical effects, the dissipation of TKE requires a careful discussion in the context of TCM. Acting as a sink term for TKE in overshooting layers, the dissipation rate has a direct impact on the extent of convectively mixed regions. Assuming a Kolmogorov spectrum of turbulence, the dissipation rate of TKE can conveniently be computed by a local expression involving a dissipation length scale with a single constant parameter. This expression is, however, inapplicable in non-local situations, encountered in layers adjacent to convectively unstable zones. To treat the dissipation of TKE in non-local convection models, a physically more complete description of the dissipation rate is required (Zeman & Tennekes 1977; Canuto & Dubovikov 1998).

We begin this paper by discussing local and non-local descriptions of the dissipation of TKE in Sect. 2. From the dissipation rate in non-local convection theories we derive a model to account for the dissipation of TKE by buoyancy waves in overshooting layers in Sect. 3. In Sect. 4 we then discuss implications of the improved dissipation model when applied to stellar models. For the computation of the stellar models we use the TCM derived by Kuhfuß (1987) implemented into the GARching STellar Evolution Code (GARSTEC, Weiss & Schlattl 2008). The key assumptions and approximations of the Kuhfuß (1987) model are reviewed in Appendix A. Using the local expression for the dissipation rate of the TKE, we find an excessive overshooting extent beyond convective cores. When including the dissipation by buoyancy waves, this overshooting is limited to a physically more reasonable range. This allows us to predict the convective core sizes and temperature structures of stars with different masses. We present our conclusions in Sect. 5. A detailed discussion of the results obtained from the improved TCM can be found in Ahlborn et al. (2022; Paper II in the following).

2. On the dissipation rate ϵ of turbulent kinetic energy

The necessity to account for the dissipation rate of turbulent kinetic energy, ϵ, in models of convection stems from the fact that it is not a negligibly small quantity. Indeed, the expression from which ϵ is computed is proportional to the kinematic viscosity ν. The latter is small in stars compared to the radiative diffusivity χ which results in the small values of the Prandtl number Pr = ν/χ typical for stars. Energy conservation requires ϵ to remain finite and non-negligible even in the limit of small viscosity (see Canuto 1997b). Neglecting compressibility (for its modelling cf. Canuto 1997a) we can compute ϵ from


in case of a locally isotropic, homogeneous flow. Here, ui is the i-th velocity component, xi is the i-th component of location, and E(k) is the spectrum of turbulent kinetic energy as a function of wavenumber k1. Although convection is neither isotropic nor homogeneous on those large scales on which its contribution to energy transport is maximal, Eq. (1) is a sufficient approximation to explain some basic properties of turbulent convection2. In a quasi-stationary state where the amount of kinetic energy injected into the system per unit of time equals ϵ, the enstrophy Ω of the flow increases, if ν decreases. The latter follows from the vorticity ω through . Thus, ϵ is constrained by energy conservation and quantifies the amount of kinetic energy converted into thermal one.

If for a flow both the first and second Kolmogorov hypotheses hold (Pope 2000), then there exists a range of length scales  = π/k for which ϵ is independent of ν and independent of the details of the large scale input of kinetic energy into the flow. This region is known as the Kolmogorov inertial range. In that region ϵ is solely described by the exchange of energy between larger and smaller scales. If this exchange peaks between neighbouring scales (see Lesieur 2008), which is assumed to hold for turbulent flows except for corrections due to intermittency (see also Pope 2000), it can be modelled as a flux in k-space. This is one of the basic inputs for the turbulence model of Canuto & Dubovikov (1996) used in Canuto & Dubovikov (1998) to justify the mathematical form and the constants involved in closure relations derived for their one-point closure Reynolds stress model of convection (see their Eq. (9c)).

One important consequence for Eq. (1) is the following one: if an inertial range exists, it can be shown to require


to hold, i.e. a Kolomogorov spectrum to exist. Here, Ko is the Kolomogorov constant which also turns out to equal 5/3 in the model of Canuto & Dubovikov (1996, 1998), just as the power law index for k in the spectrum, Eq. (2). Recalling Eq. (1) this demands that the contributions of small scales k to ϵ increase with k1/3 and a region where the Kolmogorov inertial range no longer holds, just around the dissipation scale kd = π/d, would have to be characterised more accurately than through Eq. (2) for a direct computation of ϵ from the spectral energy distribution E(k) and Eq. (1). Within a one-point closure model and thus in any of the prescriptions used in astrophysics to compute the convective flux inside a stellar structure code, this is not feasible and a different approach is required to compute ϵ.

2.1. Computation in local models: spectra and local limits

One way around computing the spectrum ϵ(k) is to just compute its integral value ϵ from a model of E(k) as follows. Assume that ν is negligibly small. In the limit of vanishing ν the latter can ensure that its product with ∫k2E(k) dk remains finite even though Ω might increase indefinitely for arbitrarily large k. Hence, as in the derivation of Eq. (5b) of Canuto & Dubovikov (1998) and as also in their Sect. 6.4, assume that E(k) is given by Eq. (2) from a certain value k0 onwards, i.e. the entire energy spectrum is given by a Kolmogorov spectrum with an energy cutoff for k < k0. Thus, E(k) = 0 for k < k0 and E(k)∼k−5/3 for arbitrarily large k with k ≥ k0. In this case, it is easy to first obtain K, the turbulent kinetic energy (TKE), from integration of E(k) over all wavenumbers:


and with Eq. (2) we obtain from Eq. (3) that


For 0 = π/k0 this can be quickly rearranged to yield


as shown in Canuto & Dubovikov (1998)3. This is also the standard ‘local’ or ‘mixing length’ prescription for the computation of ϵ. It assumes maximum separation of the energy carrying scales around 0 and the Kolmogorov dissipation scale d (assumed to be negligibly small, and not to be confused with Λ). Moreover, it assumes validity of the inertial range as if all the energy input were at one length scale only, i.e. at 0, here set to be equal to Λ. All other scales for which  ≳ 0 behave as if they were unaffected by the very small scales (scale separation) and also by the details of the energy input. Thus, a perfect energy cascade is assumed. Mixing length theory (MLT) in addition replaces E(k) with a δ-function peaked at l0 such that its integral yields Eq. (5). It is thus a ‘one-eddy approximation’ where all the energy transport due to convection occurs on the critical (mixing) length scale Λ which has to be computed for each layer. Either way, the challenge of computing ϵ turns into the challenge of prescribing Λ.

Evidently, this cannot be an accurate model, since at least a range of scales spanning easily an order of magnitude (consider different granule sizes as an example) is expected to contribute to energy transport at convective stellar surfaces such as those of our Sun. Thus, Eq. (5) can at most be an estimate of order O(1). For some flows such as a shear flow in a pipe (Poiseuille flow), for which the mixing length formalism to compute the turbulent viscosity had originally been proposed by L. Prandtl (see Pope 2000, for example), this length can be fairly well constrained from geometrical arguments. Not surprisingly this is the application for which this prescription is most reliable. For compressible convection on the other hand this length is much more difficult to constrain and the standard choice is to assume that


where α is the MLT-parameter or mixing length parameter and Hp is the local pressure scale height in the convective zone. This situation has motivated Canuto & Mazzitelli (1991, 1992) to suggest a new convection model in which ϵ is computed directly from Eq. (3). That removes the uncertainties introduced by the one-eddy approximation, but a scale length Λ is still introduced in this model. It compares the geometric size of flow features which transport most of the energy with the length scales dominated by dissipation. This has permitted easy implemention into existing stellar evolution codes based on MLT. The same approach was used by Canuto et al. (1996).

But that concept collapses if an overshooting zone has to be modelled. In such a region, located just underneath or above a convectively unstable zone, the convective flow is fundamentally non-local: the only way to sustain a non-vanishing solution is transport of kinetic and potential energy from the adjacent convective zone (cf. Sect. 10 in Canuto & Dubovikov 1998). For such a region there is no reason to assume that the prescription of Eq. (6) with an α independent of vertical location can still hold.

Thus, even if other equations in a convection model are treated non-locally, the continued use of Eq. (6) with Eq. (5) along with a constant α even just within a single object may lead to inconsistent or unphysical results, a fact long acknowledged in the atmospheric sciences by much more advanced modelling (see, for instance, Zeman & Tennekes 1977). As we show below, this is exactly the problem one encounters when using the 3-equation Kuhfuß theory (Kuhfuß 1987), and it motivated the present work on how to proceed and improve the computation of ϵ in such a case.

2.2. Computation in non-local models: the dissipation rate equation

A common starting point for non-local models of convection is the dynamical equation for turbulent kinetic energy:


as given in Canuto (1993), for example. In the Boussinesq case though, the corrections due to compressibility given by the term Cii are zero. For the case of a low Prandtl number and if there are no contributions by a mean shear or rotation, we obtain (Canuto 1992)


which within the Boussinesq approximation is an exact equation, though yet unclosed. Here, ∂t and ∂z are partial derivatives with respect to time t and vertical (radial) coordinate z (in case of radial coordinates ∂z → r−2rr2). This is a prognostic equation for the second order moment with derived directly from the Boussinesq approximation of the Navier-Stokes equations through ensemble averaging. The non-local transport includes the flux of kinetic energy (in the Boussinesq approximation given by with w as the fluctuating component of vertical velocity) and of pressure fluctuations, . This is to be balanced by local production, , and the local sink given by −ϵ. Through the cross-correlation the production is readily linked to the convective (enthalpy) flux . The latter is exact in the Boussinesq approximation and can be generalised to a compressible flow. The quantities g, αv, cp, and ρ are the local (vertical) gravitational acceleration, the volume expansion coefficient, the specific heat at constant pressure, and mass density.

To solve Eq. (8) we need to know ϵ. The exact evolution equation for ϵ was first derived by Davidov (1961). In their Sect. 3, Hanjalić & Launder (1972) emphasised4 why it is difficult to close this equation. But in the same paper they also point out how to proceed to derive a new equation which models the transport of ϵ. One term (diffusional transport due to pressure fluctuations) is argued to be small on general grounds compared to other contributions while others are modelled such that the ensuing closure constants can be determined in the case of simple flows directly from experiments: decaying turbulence behind a grid and a constant-stress layer adjacent to a wall. Their model equation for ϵ eventually reads


where P means production of dissipation (due to shear or buoyancy or both). The term ∂z(νzϵ) is only relevant at moderate or low Reynolds numbers and can always be neglected for small Prandtl numbers as is the case for stars. The term Df(ϵ) was suggested to be parametrised as


where νt requires a model for turbulent viscosity such as5νt = CμK2/ϵ with a closure constant Cμ. Although this term is mainly relevant for moderate to low Reynolds numbers, it must be kept and modelled, since this is just what we also encounter in the case of overshooting zones. This is in contrast with terms only relevant for moderate to large Prandtl numbers (i.e. only in a non-stellar case) or which are small independently of the parameter space considered: those we can safely neglect for our applications. We emphasise that contrary to Eq. (8) all contributions to Eq. (9) contain closure approximations. Hence, Eq. (9) is essentially a model for ∂tϵ and not an exact evolution equation.

Equation (9) was reconsidered by Canuto et al. (1994) and Canuto & Dubovikov (1998), who also suggested the additional contribution to Eq. (9) introduced in Zeman & Tennekes (1977):


Here, β  =   − ((∂T/∂z)−(∂T/∂z)ad) is the superadiabatic gradient. In addition to c1  =  1.44 and c2  =  1.92, which is close to the middle of the typical range of values in earlier work (Tennekes & Lumley 1972; Hanjalić & Launder 1976), Canuto & Dubovikov (1998) suggested Cμ = 0.08 from their turbulence model (Canuto & Dubovikov 1996), which they obtained using Eq. (2).

Before quantifying the new term more closely, the physical origin of the contributions to Eq. (11) requires some explanation. The first term on the right hand side provides a closure for the production of dissipation by buoyancy (Hanjalić & Launder 1972). The second term was discussed already in detail by Hanjalić & Launder (1972) and represents a closure for the combined effects of the exact terms describing the generation of vorticity fluctuations through self-stretching in turbulent flows and the decay of turbulence due to viscosity. For the exact term of diffusion of ϵ by velocity fluctuations, Df(ϵ), both a down-gradient closure (Hanjalić & Launder 1972) and a direct closure based on the flux of turbulent kinetic energy (Canuto 1992) have been proposed. The viscous diffusion term ∂z(νzϵ) is also part of the exact expression for diffusional transport and is suggested to be kept when modelling flows in the regime of low to moderately high Reynolds numbers, especially in the case of moderate to high Prandtl numbers (see Hanjalić & Launder 1976).

For buoyancy driven flows Eq. (9) requires several changes in comparison with Hanjalić & Launder (1972, 1976). We refer the reader to the work by Zeman & Lumley (1976) and Zeman & Tennekes (1977) which eventually allowed the derivation of Eq. (11). What follows from their and similar considerations is that, irrespectively of the detailed physical nature of increased local dissipation in the overshooting zone, a separately parametrised loss term that involves the superadiabatic temperature gradient β, or actually, the Brunt-Väisälä frequency, , is needed. With hindsight gravity waves are expected to play the most important role as a source of ϵ. As argued by Zeman & Tennekes (1977), this involves a characteristic length scale which can be computed from the ratio of flow velocity w2 and . It can also be viewed as the distance which eddies of a certain size that penetrate into the stable layer with a certain lapse rate can travel until their potential energy is fully converted into kinetic energy. It turns out that this yields the same expression as the parametrisation of dissipation by internal gravity waves: their contributions may differ in magnitude, but their functional form remains the same.

Hence, Canuto et al. (1994) suggested that this term should indeed be added to the standard form of Eq. (9). As they pointed out, this contribution also allows to maintain stationarity in homogeneous, stratified turbulence as confirmed by data from direct numerical simulations of shear turbulence by Holt et al. (1992). Thus, Canuto et al. (1994) suggested c3 = 0.3 for stably stratified layers and c3 = 0 elsewhere to complete Eq. (11). Canuto & Dubovikov (1998) followed that proposal.

Clearly though, among all the parametrisations which appear in Eq. (11), remains the most uncertain one, but yet it is also crucial. Its choice requires to be tested carefully. Otherwise, the width of convective overshooting may turn out to be sensitive to the detailed calibration of its parameters. In Appendix B we discuss more recent suggestions to further improve the physical content of Eq. (11).

3. A new model for the dissipation rate in non-local convection models in GARSTEC

3.1. The problem: Overshooting zones of convective cores growing unlimitedly during main-sequence stellar evolution

The Garching Stellar Evolution Code (GARSTEC; see Weiss & Schlattl 2008) offers several models to compute the contributions of convection to energy transport and mixing in stellar evolution calculations (including those of Böhm-Vitense 1958; Canuto & Mazzitelli 1991; Kuhfuß 1987). In particular, the model of Kuhfuß (1987) has been implemented (Flaskamp et al. 2002; Flaskamp 2003) in GARSTEC both in its 1-equation version, i.e. with an additional differential equation for turbulent kinetic energy, K, and in its full, 3-equation version (Kuhfuß 1987; for a brief discussion of this model see Appendix A). The latter features differential equations for the TKE, the squared fluctuations of entropy, Φ, and for the turbulent flux of entropy fluctuations, Π. Those three equations are essentially equivalent to the dynamical equations for the TKE, the squared fluctuations of temperature , and for the cross correlation between velocity and temperature fluctuations, denoted here by . The latter can be derived from the phyically more complete model of Canuto & Dubovikov (1998) by assuming (i) an isotropic velocity distribution, (ii) a local prescription to compute the distribution of the dissipation rate ϵ, (iii) the diffusion approximation for the non-local fluxes, and (iv) some minor simplifications in the closures used in the dynamical equations6. As a variant, the 3-equation model may be used with local limit expressions for the non-local transport terms for as well as J. As a theoretical analysis shows (see Kupka et al. 2020 and references therein) only a full 3-equation model can feature a countergradient or ‘Deardorff’ layer where J is positive, while the superadiabatic gradient β is negative. Only in such a model both quantities can change their sign independently (the key to a positive convective flux in a countergradient stratification is the non-local transport of as originally shown by Deardorff 1961, 1966). However, in both the fully non-local and the local limit of the 3-equation model variant as described above, overshooting gradually mixes the entire star in a stellar evolution calculation for a 5 solar mass (B-type) main-sequence star. In Fig. 1 we show the profile of the TKE as a function of fractional mass in this calculation. It can be seen that the energy extends substantially beyond the Schwarzschild boundary, reaching very close to the surface of the star. Due to the high efficiency of convective mixing the whole star would become essentially homogeneous which is unrealistic, because the star would evolve from the hydrogen to the helium main-sequence, i.e. to the left in the colour-magnitude diagram, contrary to all observations (see Kippenhahn et al. 2012, Chap. 23.1). This problem was originally identified in the PhD thesis of Flaskamp (2003).

thumbnail Fig. 1.

TKE as a function of the fractional mass for the original Kuhfuß model. The formal Schwarzschild boundary, defined by ∇rad = ∇ad, is indicated by a dashed black line.

To solve this problem Flaskamp (2003) suggested to give up the assumption of isotropy of TKE of the model of Kuhfuß (1987) in the overshooting (OV) zone and let the ratio of vertical to horizontal kinetic energy tend to zero. This limits the mixing efficiency in the outer layers of the OV zone, located above the stellar convective core, and avoids its unphysical growth throughout main-sequence evolution. If this simulation were plausible, also a more realistic model for the anisotropy of the convective velocity field, derived, for instance, from the stationary limit of Eq. (19d) of Canuto & Dubovikov (1998), should solve this problem. Both variants of this approach are discussed below in Sect. 4.1.

3.2. A comparison with a fully non-local Reynolds stress model

A progressive growth of the overshooting zone with time is not observed in 3D radiation hydrodynamical simulations of overshooting in DA white dwarfs (Kupka et al. 2018) either. Since the extension of the different zones in that case (Schwarzschild unstable convective zone with J > 0 and β > 0, countergradient region with J > 0 and β < 0, plume dominated region with J < 0 and β < 0, and wave dominated region with J ≈ 0 and β < 0) compare quite well with results from the non-local Reynolds stress model of Canuto & Dubovikov (1998) solved in Montgomery & Kupka (2004) for the same type of stars, the latter can provide a guideline for the behaviour of variables such as ϵ as a function of depth. The overall structure of the OV zones and the behaviour of the convection related variables described in Montgomery & Kupka (2004) is very similar to that one which had already been found for A-type main-sequence stars in Kupka & Montgomery (2002) which in turn had been found in qualitative agreement with earlier 2D radiation hydrodynamical simulations of Freytag et al. (1996).

We hence use the Reynolds stress convection model calculations of Kupka & Montgomery (2002) in Fig. 2 to illustrate the convective flux, the root mean square vertical velocity, and the dissipation rate as a function of depth. The left panels show results for the full third order moment model while the right panels show results computed using the downgradient approximation (for the latter cf. also Appendix A.3). For Teff = 8000 K and a log g value slightly below the main sequence (see Kupka & Montgomery 2002 for further details) we find two convective zones, an upper one due to ionisation of neutral hydrogen and a lower one caused by double-ionisation of helium. They are connected by an overshooting region at a radius of ∼931 Mm and there is another overshooting region underneath the lower convective zone at ∼926.5 Mm. For this setting we compare the computation of dissipation rates from Eq. (11), the full equation of Canuto & Dubovikov (1998), with the standard mixing length prescription for a range of bulk convective and overshooting layers. Clearly, the dissipation rate ϵ becomes much larger than the value computed from the MLT prescription as soon as the plume region of the OV zones (with J < 0 and β < 0) is reached, and which can be determined from the behaviour of the convective flux. At the bottom of the lower overshooting zone, ϵ becomes even order(s) of magnitudes larger than the oversimplified MLT prescription would predict. Note that if the downgradient (diffusion) approximation is used to compute third order moments such as in the model of Canuto & Dubovikov (1998; the non-local fluxes of K, J, , and ), a smaller overshooting is obtained in comparison with the complete third order moment model used in Kupka & Montgomery (2002). Hence, the two convection zones become separated at Teff = 8000 K which allows observing this behaviour of ϵ even between the two convective zones. At lower Teff, for example at 7500 K, convection and overshooting are stronger also for the downgradient approximation of third order moments and the same behaviour is recovered as for the physically more complete third order moment model at Teff = 8000 K. For that latter model the two convective zones become more tightly coupled and the increase of ϵ compared to the MLT prescription is eventually restricted to the lower overshooting zone only, for instance, for models with Teff = 7200 K.

thumbnail Fig. 2.

Left panels: convective flux in units of total flux, root mean square vertical velocity in units of km s−1, and dissipation rate ϵ from Eq. (11) relative to a value computed from Eqs. (5) and (6) with α as given in the figure legend. Right panels: same quantities as left panels, however, the downgradient approximation (DGA) is used to compute third order moments instead of the full model used in Kupka & Montgomery (2002). The results are for one of the A-star envelope models discussed in Kupka & Montgomery (2002).

We hence can draw the following conclusions from solutions of the Reynolds stress model of Canuto & Dubovikov (1998) for convective envelopes of A-type stars: irrespective of the various situations described above, deep inside the plume-dominated region characterised by J < 0 and β < 0 the MLT prescription to compute ϵ begins to fail by entirely missing out the drastic increase in dissipation in that region. However, the proper computation of ϵ is essential to determine the extent of the mixed region, since it drains kinetic energy from the overshooting flow. From Eq. (5) one can immediately conclude that underestimating ϵ in the MLT framework can be easily caused by overestimating the mixing length Λ or 0.

3.3. Reducing the mixing length in the OV zone

There is also a physical argument why the mixing length must be limited and even gradually shrink in the OV zone on top of a stellar convective core. Taking Λ to be about a pressure scale height at the convective core boundary results in a very large length scale. This is essentially the size of the convective core itself. The claim that such a large structure penetrates into the radiative zone makes no sense, both from the viewpoint of available potential energy and from the viewpoint of the typical size of a convective structure. We note here that existing numerical simulations of convective cores are actually for extremely different physical parameter regimes, featuring mostly Pr ≳ 1 or even Pr ≫ 1 (see, for instance, Rogers et al. 2013; Rogers 2015; Edelmann et al. 2019a). They are unable to reproduce the very small levels of superadiabaticity (β > 0, but |β/(∂T/∂r)ad|≪1) at realistic stellar luminosities. This inevitably leads to excessive numerical heat diffusion and unrealistically small effective Peclet numbers (see Kupka & Muthsam 2017 for a discussion). Numerical simulations of convective cores may hence also be affected by the convective conundrum problem reported for the Sun (cf. Gizon & Birch 2012; Hanasoge et al. 2016). Probably, they are not as reliable for guiding us as numerical simulations are in the case of convective overshooting near stellar surfaces (cf. Freytag et al. 1996; Tremblay et al. 2015; Kupka et al. 2018, and many others). We return to the problem of comparing results on convective cores from stellar evolution models with 3D hydrodynamical simulations of convective cores in Sect. 4.3. In the following, we thus use a different chain of arguments to derive an improved estimate of Λ.

As a very first step, one could let Λ decay to zero within the OV zone, either linearly or exponentially, from the value it has at the boundary of the convective zone. This ad hoc ‘fix’ has been implemented into GARSTEC. The exponential decay model was chosen and indeed this easily stops the growth of the overshooting zone as a function of stellar evolution time. The enhanced dissipation rate introduced this way can be seen in Fig. 3 at the outer edge of the convective region. The model including the exponential decay has a central hydrogen abundance of 0.6. The stellar model computed with the original Kuhfuß model was chosen to have the same maximum TKE in the convection zone to make the dissipation rates comparable.

thumbnail Fig. 3.

Dissipation rate as a function of fractional mass for the original Kuhfuß model and the Kuhfuß model including an ad hoc exponential decay of the dissipation length, shown with a grey dotted and a blue continuous line, respectively. The ad hoc exponential decay of the dissipation length leads to an increased dissipation rate at the beginning of the overshooting zone, indicated by the local maximum beyond the Schwarzschild boundary, followed by a sharp drop due to the rapid decay of TKE. The models have been chosen to have the same maximum TKE.

Physically plausible extensions of the OV zone can be obtained from a “reduction factor”, which forces an e-folding extent of the “decay” of the mixing length of 2% to 6% of the mass of the Schwarzschild-unstable region. In a 5 M main-sequence star this limits the OV zone to contain about 12% to 29% in terms of the Schwarzschild core mass. The relative extent of the overshooting region in terms of the Schwarzschild core mass remains mostly constant along the main-sequence. For an e-folding extent of 4% the overshooting region contains about 5% of the stellar mass at the beginning of the main-sequence while it is shrinking to about 2% of the total mass at the end of the main-sequence. The procedure introduces a free parameter, but it is sufficient as a proof of concept: a physically more complete model of ϵ constrains the OV contrary to earlier, alternative explanations that require unphysical parameter values to do so (such as which is at variance with Kupka et al. 2018, see Sect. 4.1 below).

3.4. Boundary conditions and regularity constraints

As a prerequisite to derive an improved estimate for Λ, we first discuss its asymptotic behaviour in the centre of a convective core. Regularity properties of non-local models of convection at the centre of stellar cores are a rather delicate issue which has been analysed in Roxburgh et al. (2007b). Under the assumption that non-zero convective motions can also occur at the centre of a convective core, for the second order moments they demonstrated that , K, and are all positive and have an even order expansion in r just like the gas pressure P. Moreover, from their Eq. (11), the horizontal component of TKE has to balance the vertical component in the sense that for the velocity components in spherical polar coordinates (r, θ, ϕ). Hence, and the flow is isotropic. Clearly, also ϵ has to be positive in this case.

Thus, if the relation ϵ = cϵK3/2/Λ is used, a positive Λ guarantees positivity of ϵ. An appropriate prescription which ensures this property is to use the curvature of the pressure profile to define a local scale height, since its gradient vanishes at the centre. This has been worked out in Roxburgh et al. (2007a) where the scale height at the centre is defined from and the subscript c denotes the value of the local scale height, Hc, of pressure, Pc, and density, ρc, at the centre (and G is, of course, the gravitational constant). For the centre, Λ = αHc and in general Λ = αmin(Hp, Hc). Roxburgh et al. (2007a) suggest a smooth interpolation between the limit at the core centre and the expression for Hp ≫ Hc.

For reasons of regularity and energy conservation, Fconv → 0 at the centre in that case, which is fulfilled by the above prescription of the mixing length. An expansion in odd powers of r is found for the Reynolds stress equation for J and thus for Fconv (Roxburgh et al. 2007b). This implies non-trivial constraints on closures for the third order moments. Roxburgh et al. (2007b) demonstrate that the downgradient closure forces the core centre to be convectively neutral (J ∝ r3 instead of J ∝ r) while other closures have to be modified to ensure regularity of the solution.

In GARSTEC, the Wuchterl (1995) prescription for Λ is used by default. This requires a different approach at the core centre, as it assumes Λ → 0 for r → 0. Thus, in GARSTEC, it is ensured by power series expansions that the convective variables are not forced to zero while the temperature gradient at the centre is the adiabatic one. As a result, the convective quantities become small in the central region (see Paper II). It is hence important to clarify whether the differences between these two rather different prescriptions for the convective variables in the centre of a convective core are relevant for applications. Fortunately, it turns out that they remain constrained to less than the innermost 10% of the stellar core. In either case, the stellar core is predicted to be fully mixed and has a temperature gradient close to the adiabatic one. For this study, we hence prefer to stay within the standard setup used for GARSTEC, i.e. we start from the prescription for Λ proposed by Wuchterl (1995).

3.5. Some input from the dissipation rate equation

We now show how it is possible to carry over some of the physics contained in Eqs. (9) or (11) into a local model for ϵ and still avoid the solution of an additional differential equation. If we model the non-local transport of TKE in Eq. (8) by a downgradient approximation, the closure relates to in Eq. (11). The same behaviour is found for a direct downgradient closure for (i.e., computing it from ∂zϵ) as for example Eq. (10). Let us hence assume a local approximation for Df(ϵ), the non-local flux of ϵ, which replaces the derivatives of the outer divergence operator and the gradient operator in Eq. (11) by a product of reciprocal length scales, 1/2. Inspecting Eq. (11), for the sake of simplicity, it appears desirable to model as many contributions as possible by expressions of type ϵ2/K ∝ ϵ/τ. Instead of a diffusion length scale () we hence use the characteristic transport time scale τ = 2K/ϵ to approximate Df(ϵ)∝ − αϵϵ/τ. The same can be done also in the case of Eq. (9). If we furthermore assume the local limit of Eq. (8), P = Pb = ϵ, i.e. production of TKE by buoyancy equals its dissipation, and if we also assume c3 = 0, we obtain the following approximation for both Eqs. (9) and (11):


To remain consistent with gαJ = ϵ we have to require that αϵ = 2c2 − 2c1 if ϵ itself is computed from Eqs. (5)–(6). In this case we obtain a completely local model for the computation of ϵ.

We can use Eq. (12) to understand some implications from the different physical contributions which its physically more complete counterpart, Eq. (9), would instead account for. To this end let us relax the requirement Pb = ϵ in Eq. (8) somewhat. In this case, whether the 1-equation or the 3-equation version of the Kuhfuß (1987) model is used (cf. Appendix A), due to the non-locality of the flux of kinetic energy in Eq. (8), , there is always a point where J = 0 (cf. Chap. 5 in Kupka et al. 2020). At such a point, αϵ = 2c2 is required from Eq. (12) for a non-vanishing dissipation rate ϵ. Right next to such a point, where ϵ > 0 with J < 0, a value of αϵ > 2c2 would be required whereas αϵ < 2c2 where J > 0. So αϵ would have to be a function that has to be fine-tuned to obtain consistent results from Eq. (12) in the vicinity of J = 0. Moreover, because of the downgradient closure for also constraints on would be imposed.

Such constraints appear unphysical: Eq. (12) does not provide a good starting point for a local model capable to capture at least the main gist of either Eq. (9) or Eq. (11). To proceed we need a physically more complete model for ϵ, i.e. we either have to abandon the mixing length prescription altogether or we need a more complete model equation than Eq. (9) to start from. Let us hence first have a look at Eq. (11), i.e. we no longer impose c3 = 0 everywhere. The sibling of Eq. (12) which accounts for the production of dissipation by gravity waves in stably stratified fluid then reads:


If we were to combine this equation with the 1-equation model of Kuhfuß (1986), β and J change sign at the same point so the perfect balancing constraint between Df(ϵ) and −2 c2ϵ/τ reappears. In the region where J < 0, more freedom of how Df(ϵ) behaves is permitted. This changes once we switch to the 3-equation model of Kuhfuß (1987): since β and J then change sign at different locations, αϵ is no longer forced by c2 at any point. In the end, the contribution decouples both Df(ϵ) and from peculiar constraints required to be fulfilled at where β = 0 or where J = 0.

On the other hand, now there is an efficient local source for ϵ also where β < 0. This is particularly important for the 3-equation model which through its countergradient layer permits much larger enthalpy (and hence also TKE) fluxes in this region: considering that property it is understandable that the 3-equation model can be prone to large overshooting, unless the latter is limited by efficient dissipation. And this is just what gravity waves can provide.

3.6. Deriving a local model for ϵ with enhanced dissipation

For the sake of physical completeness it would be preferable to switch to Eq. (11) and give up the local model Eqs. (5)–(6) altogether. However, as a first step into that direction we can aim at modifying the computation of Λ for the stably stratified layers by guiding the necessary physical input through Eq. (11) and in particular through its local approximation, Eq. (13). In a local framework we cannot accurately account for Df(ϵ). Hence, we first express τ in terms of Λ in the local limit,


from which we obtain that


To proceed we can now rewrite as follows:


Following the analysis in the previous subsection we now compare Eq. (16) with


In the stationary, local limit and assuming that we can absorb the contribution from αϵϵ/τ + 2 c1gαJ/τ into −2 c2ϵ/τ for sufficiently small J and we obtain from Eqs. (13), (16), and (17) that


where the numerical value is obtained from setting c2 = 1.92 and c3 = 0.3. Contributions absorbed into the −2 c2ϵ/τ term could be accounted for by small change of c2. As inspection of the full Reynolds stress models solved in Kupka & Montgomery (2002) demonstrates this is well justified since the two terms compared in Eq. (18) completely dominate where J < 0.

This motivates the idea to also scale Λ, which according to Eq. (15) is proportional to τ, by a contribution . In GARSTEC the mixing length required for the turbulent convection model of Kuhfuß (1987) is computed following the prescription of Wuchterl (1995),


where βs is a factor chosen to be 1 in convectively unstable layers, where β = −(dT/dr − (dT/dr)ad) > 0 and thus ∇ − ∇ad > 0, and βs is possibly less than 1 elsewhere. We now account for the effect of enhanced dissipation by gravity waves through reducing βs to values less than 1. To this end we can interpolate between the two asymptotic cases and through


where Mschw is the mass of the convectively unstable core and thus identifies the mass shell for which ∇ = ∇ad and λs is a model parameter. Comparisons with solutions of the non-local Reynolds stress model of Canuto & Dubovikov (1998) for A-type stars (Kupka & Montgomery 2002) show that τb ≈ 0.1 τ where Fconv reaches its negative minimum. This range of values for τb is what we also expect from Eq. (18) for a moderate variation of c2.

The results of Kupka & Montgomery (2002) can hence provide a rough guideline for the choice of λs and imply that Λ is rapidly reduced by an order of magnitude already within the countergradient region from the value it has at the Schwarzschild stability boundary (see Fig. 4). This value is then maintained throughout the remainder of the countergradient region and the entire region where Fconv < 0, in agreement with the suggested by Canuto (2011c) in his Eq. (5h). The preceding arguments and the analysis in the previous subsection show how this relation is connected with the full Eq. (11) and how this result can be implemented into a physically motivated reduction factor for the mixing length through Eqs. (19) and (20). Since the rough constancy of τ/τb (or the “dominance” of the term in Eq. (11)) also causes the linear decay of the root mean square velocity as a function of distance in the results of Kupka & Montgomery (2002) and Montgomery & Kupka (2004), and because the latter has also been recovered from 3D radiation hydrodynamical simulations (Kupka et al. 2018) for just those layers, the entire procedure is at least indirectly supported by this physically much more complete modelling. Similar results are not yet available for convective cores, however.

thumbnail Fig. 4.

Ratio of τ/τb as a function of convective stability from a solution of the non-local Reynolds stress model as presented in Kupka & Montgomery (2002) assuming the downgradient approximation for third order moments. The time scale τb is computed from where the absolute value of β is taken. Sign changes are hence indicated by spikes. Both the overshooting zones below and above the lower and the upper convectively unstable zone show the same increase of τ/τb from 0 to more than 10 (the finite grid resolution prevents τ/τb from becoming actually zero).

In spite of its simplicity the disadvantage of Eq. (20) is the fact that λs is a dimensional parameter. It hence has to be determined separately for each stellar evolution model by numerical experiments which yield the value it has to have for a sufficient reduction of Λ by an order of magnitude. For stars of different mass this may have to be changed and for modelling later stages of stellar evolution it is even less convenient. What we need instead here is an estimate for τ. Without solving Eq. (11) this is akin to a hen and egg problem, since in the end this would require just the quantity Λ we are up to compute: λs = (25/320) τ with τ computed from Eq. (15). We could simplify this by setting τ = (2/cϵ)(αHpK−1/2) or τ = (2/cϵ)(rK−1/2), as this formula is to be used only for r > 0 and Hp < ∞ anyway. However, this has the disadvantage that near the outer edge of the overshooting zone where K → 0 one obtains τ → ∞. From standard calculus applied to Eq. (19) we then obtain that Λ ≈ αHp right there which is exactly not what we want. But we can rewrite Eq. (20) into




for which we have used cϵ = π(2/(3 Ko))3/2 ≈ 0.7948 ≈ 0.8 with Ko = 5/3 from Canuto & Dubovikov (1998)7. This is achieved by realising that which is just Eq. (18) and where we have used Eq. (15) for the last step. Equation (21) is equivalent to Eq. (20) and also interpolates between the two asymptotic cases, the transition between locally stable to unstable stratification () as well as the overshooting region far away from the convective zone, where flow motions are dominated by waves (). Equation (19) combined with Eqs. (21)–(22) can be rewritten into a quadratic equation for Λ for which the positive branch can be taken or which can be solved implicitly, for instance, by an iterative scheme (the former will be done in Paper II). In principle, the parameter c4 could be adjusted to achieve the goal of τb ≈ 0.1 τ or rather Λ(min(Fconv)) ≈ 0.1Λ(Mr = Mschw) which mimics the result discussed in Fig. 4 and in the previous paragraphs. However, we prefer to assume sufficient generality of Eq. (11) and its parameters and therefore use them without further adjustments. Some numerical experiments on the effects of varying c4 can be found in Appendix B of Paper II. In the next section we show that this procedure also leads to a finite overshooting layer which does not (notably) grow during stellar evolution.

4. Discussion: Kuhfuß 3-equation model with enhanced dissipation

4.1. Flow anisotropy instead of enhanced dissipation

A very important difference between the Kuhfuß (1987) and the Canuto & Dubovikov (1998) model is the set of convective variables considered. In addition to the TKE Canuto & Dubovikov (1998) also solve for the vertical TKE. This means that the ratio of is not fixed a priori but is an outcome of the theory. Kuhfuß (1987) on the other hand assumes full isotropy in the whole convection zone which translates to a fixed ratio of . Furthermore, the Kuhfuß model uses an isotropic estimate of the radial velocity in the non-local terms. Hence, these terms are potentially overestimated by overestimating the ratio of vertical to total kinetic energy. This could result in an unreasonably large overshooting zone. The treatment of the flow anisotropy is especially problematic at convective boundaries where the flow turns over. In the convective boundary layers the motions change from being predominantly radial to becoming predominantly horizontal. This means that the ratio of vertical to total kinetic energy should drop from the isotropic value to smaller values.

To study the impact of anisotropy we mimic the change of the flow pattern by introducing an artificial anisotropy factor . This anisotropy factor is set to a value of in the bulk of the convection zone and then linearly decreases to a value of zero from the Schwarzschild boundary outwards. This is most probably not a very physical behaviour but just meant for illustrative purposes. The profile of this artificial anisotropy factor is shown in the upper panel of Fig. 5. The profile of the TKE computed with this anisotropy factor is shown in the lower panel of the same figure. The black dashed line indicates the Schwarzschild boundary. It can be seen that an overshooting zone beyond the Schwarzschild boundary emerges, which has, however, a clearly limited extent. As intended, a limitation of the anisotropy could solve the problems observed with the original version of the 3-equation model. The description requires another free parameter which is the slope of the linear function. The slope parameter directly controls the overshooting distance which is very similar to other ad hoc descriptions of convective overshooting. Also, the functional form of ξ has not been determined by physical arguments but has been chosen arbitrarily.

thumbnail Fig. 5.

Artificial anisotropy factor ξ and TKE as a function of fractional mass in the upper and lower panel, respectively. The black dashed line indicates the Schwarzschild boundary.

This unfavourable situation should be avoided by a physically motivated estimate for the anisotropy factor. This requires to compute the vertical kinetic energy. To obtain an estimate of the distribution of the turbulent kinetic energy in the Kuhfuß (1987) model we start from the fourth equation of the Canuto & Dubovikov (1998) model:


which solves for the vertical turbulent kinetic energy . Not solving for implies that also is unknown. A reasonable way to compute this quantity from the Kuhfuß (1987) model is again to assume an isotropic distribution of the fluxes: . By rearranging and neglecting the time-dependence in Eq. (23) we can define an anisotropy factor:


All quantities in Eq. (24) can be computed within the Kuhfuß 3-equation model.

We have computed the anisotropy factor according to Eq. (24) for a stellar model which used the original version of the Kuhfuß 3-equation model. The result is shown in Fig. 6. In the bulk of the convection zone within the Schwarzschild boundary the estimated anisotropy points towards a radially dominated flow. Directly beyond the Schwarzschild boundary the estimated anisotropy factor drops below the isotropic value of 2/3. This can be attributed to the negative convective flux in the overshooting zone which according to Eq. (24) reduces the ratio of vertical to total kinetic energy. Further out in mass coordinate the estimated anisotropy increases again slightly above a value of 2/3 and remains to a good approximation constant over the region in which positive kinetic energy is observed (see Fig. 1).

thumbnail Fig. 6.

Estimate of the anisotropy factor according to Eq. (24) for a 3-equation model without limited dissipation length-scale Λ. The profile of the turbulent kinetic energy of this model is shown in Fig. 1.

Introducing this anisotropy factor into the Kuhfuß 3-equation model would not substantially reduce the estimate of the radial velocity. On the contrary, over large parts of the model the value of the radial velocity would be even larger than the current estimate as we find an anisotropy factor above the isotropic value of 2/3. To finally settle the question of the flow anisotropy in Reynolds stress models one also has to solve the respective equation for the vertical kinetic energy (Eq. (23) shown here, as taken from the Canuto & Dubovikov 1998 model) self-consistently coupled to the non-local convection model. However, since such a more realistic anisotropy factor cannot resolve the problem of excessive mixing found in the original Kuhfuß 3-equation model and because its implementation as an additional differential equation increases the complexity of the model, we first perform a thorough analysis of the improved 3-equation model in Paper II and postpone the extension of this new model to future work.

4.2. Dissipation in the Kuhfuß 1- and 3-equation model

We have implemented the enhanced dissipation mechanism, developed in Sect. 3.6, into GARSTEC. For the details of the implementation we here refer to Paper II. With this implementation we solve the stellar structure equations and the convective Eqs. (A.4)–(A.6) self-consistently. We note that for consistency and to simplify the comparison between the 1-equation and the 3-equation model, we set cϵ = CD (see Appendix A), whence it follows that c4 ≈ 0.072 in those calculations. As an example we show here the TKE in a 5 M main-sequence star in Fig. 7. The Schwarzschild boundary is indicated with a black dashed line. In this model the convective energy extends slightly beyond the Schwarzschild boundary which means that an overshooting zone emerges consistently from the solution of the model equations. However, in contrast to Fig. 1 the energy does no longer extend throughout the whole star but has a clearly limited extent as one would expect for this kind of star in this evolutionary phase.

thumbnail Fig. 7.

Convective energy as a function of the fractional mass for the Kuhfuß model including the improved dissipation mechanism. The Schwarzschild boundary is indicated by a dashed black line.

This shows already that the enhanced dissipation mechanism proposed above is able to solve the problems observed in the original version of the 3-equation Kuhfuß convection model. The detailed structure and the behaviour of stellar models with different initial masses is discussed in Paper II.

The results obtained from the different versions of the Kuhfuß model can be interpreted by studying the individual terms of the TKE equation (Eq. (A.4)) in more detail. In Fig. 8 we show the three terms of the TKE equation–buoyant driving, dissipation and non-local flux–with a corresponding red, black, and blue line respectively for the 1-equation model (panel a), the original 3-equation model (panel b) and the improved 3-equation model (panel c).

thumbnail Fig. 8.

Comparison of the different terms in the TKE equation (Eq. (A.4)) in the Kuhfuß 1-equation (panel a), original 3-equation (panel b) and improved 3-equation (panel c) model. The buoyant driving term, the dissipation term and the non-local flux term are shown with a red, black, and blue line here.

Stellar models applying the non-local 1-equation theory posses a clearly bounded convective region with a reasonable extent. However, this is achieved by suppressing the countergradient layer and artificially coupling the sign of the convective flux to that one of the superadiabatic gradient.

When using the 3-equation model in its original version this welcome property vanishes and the stellar models become fully convective. As discussed in Appendix A the 3-equation model does not approximate the convective flux by a local model but rather solves an additional differential equation for it. This reduces the coupling of the different convective variables. Intuitively one would expect this model to be physically more complete than the 1-equation model and to yield physically improved models (see the discussion in Sect. 5 of Kupka et al. 2020). However, the stellar models computed with the 3-equation model look physically unreasonable, as the existence of fully convective B-stars with 5 M is excluded from the lack of stars hotter than the hydrogen main-sequence.

This rises the question why a seemingly physically more complete model leads to worse results. It can be illustrated by comparing the TKE terms in the 1- and original 3-equation models shown in panels a) and b) in Fig. 8. In the 1-equation model the buoyant driving term which is proportional to the convective flux shows negative values in the overshooting zone, which is expected due to the buoyant braking in the stable layers. The buoyant term even exceeds the actual dissipation term in magnitude. This means that in the 1-equation model it is not the dissipation term but rather the buoyant driving term which acts as the main sink term in the overshooting zone. When applying the 3-equation model the buoyant term is still negative in the overshooting zone. The values are, however, much smaller in magnitude compared to the 1-equation model. The dissipation and non-local flux term have about the same magnitude in the overshooting zone as obtained with the 1-equation model, because their functional form did not change. Considering that it was the buoyant driving term which was acting as the main sink term, the 3-equation model in its original form is lacking a sink term in the overshooting zone. This naturally explains the excessive overshooting distance found for this model.

To understand how the dissipation by buoyancy waves can mitigate this problem it is worth to recall the approximation for the convective flux in the 1-equation model. Kuhfuß (1987) has approximated this to be Π ∝ (∇−∇ad). As the convective flux is the major sink term in the overshooting zone in the 1-equation model one possibility is to introduce a dissipation term which has the same dependence, ϵ ∝ (∇−∇ad). A process with this dependence would be, for example, the dissipation by buoyancy waves as proposed above. We have demonstrated that the enhanced dissipation by buoyancy waves reduces the overshooting distance again to a more reasonable extent for the TKE (see Fig. 7). The related terms of the TKE equation are shown in Fig. 8 in panel c). In the overshooting zone the magnitude of the dissipation term is now substantially larger than the negative buoyancy term such that it acts as the dominant sink term. Also the shape of the dissipation profile has changed compared to the original 3-equation case. The transition from finite to zero values looks smoother for the improved 3-equation model because the temperature gradient which has readjusted differs in comparison with the 1-equation model.

This comparison shows why the original version of the 3-equation model results in fully convective stars. The fact that a sink term is missing points again at the importance of a dissipation term which is proportional to (∇−∇ad). On a first glance, a negative convective flux with larger magnitude in the overshooting zone could also increase the sink term in the TKE equation. But the following line of arguments shows that this hypothesis leads to unplausibly large non-local fluxes.

Here, we consider Eqs. (A.4)–(A.6)8. Let us assume that J becomes larger, or, equivalently, Π in Eqs. (A.4)–(A.6) becomes larger in magnitude in the region where it is negative. Then, the buoyant driving term shown in panel b of Fig. 8 changes towards more negative values. This permits the dissipation rate of kinetic energy to become smaller. However, in that case the buoyant driving term (containing Π) also becomes larger in Eq. (A.6), which predicts the magnitude of entropy fluctuations.

Since the vertical velocities have to become smaller, when the dissipation rate of kinetic energy becomes smaller (and we assume a constant anisotropy in this thought experiment), the squared fluctuations of entropy, Φ, or of temperature, , have to become larger instead. But for Π < 0 in the region we consider here, both −Π/τrad and (2∇adT/Hp)Φ act as sources which are boosted in Eq. (A.5). Unless we would consider a large rate of change in the non-local transport of convective flux and entropy fluctuations, the only way to obtain an equilibrium solution in this model is to increase velocities and thus also the flux of kinetic energy. This is exactly the solution observed in panel b of Fig. 8 with its excessively extended overshooting. The closure used in Canuto (1993) and Canuto & Dubovikov (1998), which also accounts for buoyancy contributions to the correlation between fluctuations of temperature and the pressure gradient (the −Π/τrad term in Eq. (A.5)) does not change this argument. But a scenario that builds up large fluctuations of entropy in the overshooting region, where radiative cooling should efficiently smooth them while it has to suppress high velocities, appears unphysical. Thus, this alternative can be excluded.

Since extensive overshooting, which eventually mixes the entire B-star, is ruled out by observations, we are left with flow anisotropy or enhanced dissipation due to the generation of waves as physical mechanisms to limit overshooting in the 3-equation framework. Because extreme levels of flow anisotropy are neither found in solar observations nor in numerical simulations of overshooting in white dwarfs (Kupka et al. 2018), nor in solutions of the model of Canuto & Dubovikov (1998) for A-stars (Kupka & Montgomery 2002) or white dwarfs (Montgomery & Kupka 2004), there is hardly evidence for this idea. On the contrary, the enhanced energy dissipation rate is contained in the full model of Canuto & Dubovikov (1998) which yields at least some qualitative agreement with numerical simulations of several scenarios of stellar overshooting (see Kupka & Montgomery 2002 and Montgomery & Kupka 2004 and compare with Kupka et al. 2018 for the latter). This makes the improved computation of the dissipation rate of kinetic energy the most plausible improvement of the 3-equation model to remove the deficiency the model has had in its original version proposed by Kuhfuß (1987).

4.3. Comparing the Kuhfuß 3-equation model with overshooting models and numerical simulations

Viallet et al. (2015) have reviewed several models suitable for parametrisation of overshooting above stellar convective cores. One of them is the model proposed by Freytag et al. (1996) based on 2D hydrodynamical simulations of thin convective zones which appear in the atmosphere and upper envelope of stars. The simulations had to be restricted to low Peclet (Pe) numbers where highly efficient radiative diffusion competes with convective energy transport. The simple exponential decay law for velocity as a function of distance from the convection zone has been particularly attractive for stellar evolution modelling and the model is available in most actively used stellar evolution codes including GARSTEC. As the velocity scales with the pressure scale height, this model requires an additional cut-off to prevent diverging overshoot from very small convective cores as found in stars with less than two solar masses. We discuss this issue in detail in Paper II. Additionally, Kupka et al. (2018) have pointed out that within the countergradient and plume dominated regions of convective overshooting zones exponential decay rates for velocity work only within a limited spatial range (see also Montgomery & Kupka 2004). Cunningham et al. (2019) argued for different decay rates for the plume dominated and the wave dominated regime. Such distinctions are, however, not made in applications of that model. We refer to Paper II to a detailed comparsion of convective core sizes between the Kuhfuß 3-equation model and the exponential overshooting model, and here just emphasise that the energy loss of turbulent flows due to waves is readily built into the improved Kuhfuß 3-equation model.

Another model, suitable for a higher Pe regime, where penetrative convection due to plumes occurs, is the one originally suggested by Zahn (1991). When applied to convective cores his model had to rely on invoking Roxburgh’s integral constraint (Roxburgh 1989) for self-consistent predictions which effectively turns it into a model similar in complexity to the 1-equation model by Kuhfuß (1987). We recall here that the 3-equation model with enhanced dissipation has a built-in dependence on Pe by accounting for radiative losses in its dynamical equations. A detailed discussion on the role of Pe in the 1-equation and 3-equation models can also be found in Paper II. The latter model is also not subject to the simplifications made in Roxburgh (1989) concerning the treatment of the dissipation rate ϵ.

Finally, for the very high Pe regime of convective entrainment Viallet et al. (2015) considered a model based on estimates relying on the variation of the inverse buoyancy time scale in the stably stratified layer next to a convective zone and the kinetic energy available at the boundary of the convective zone. The Kuhfuß 3-equation model can also deal with this case since it is the regime in which heat conduction is negligibly small. Hence, instead of relying on physically different models which have not been designed to be compatible among each other, the new 3-equation model can deal with the different regimes discussed in Viallet et al. (2015) within a single formalism and without the necessity of fine tuning for these different cases.

Comparing the predictions of the new model with those derived from 3D hydrodynamical simulations of convection is more difficult: as already mentioned in Sect. 3.3 they currently have to be restricted to a different parameter range. In Kupka et al. (2020) it is explained why the effective (numerical) heat conductivity in the simulations has to be higher than the physical one which leads to values of Pe several orders of magnitudes smaller than those found in stars. As the numerical diffusion of momentum and heat in high Pe simulations have to remain comparable to each other, we have to expect differences in flow structures and overshooting distances when compared to the actual, stellar parameter range (see, for instance, Scheel & Schumacher 2017 and Käpylä 2019). Nevertheless, it is a very important finding for the veracity of the enhanced Kuhfuß 3-equation model that the 3D simulation results concerning convective cores by Browning et al. (2004), Gilet et al. (2013), Rogers et al. (2013), Augustson et al. (2016), and Edelmann et al. (2019b), among others, and the related simulations of convective shells by Meakin & Arnett (2007), all show that convective zones excited by nuclear burning are subject to convective entrainment and penetration, respectively, depending on the specific setup, and in each case gravity waves are excited which extend throughout the radiative stellar envelope. This supports the theoretical analysis of Linden (1975) and Zeman & Tennekes (1977) for the equivalent scenario in meteorology which led to the non-local dissipation rate equation proposed in Canuto et al. (1994) and generalised to applications in stellar convection by Canuto & Dubovikov (1998), see Eq. (11), and which has been the starting point for our investigations we detail in this paper.

5. Conclusions

The original model by Kuhfuß (1987) was shown by Flaskamp (2003) to lead to convective overshooting zones on top of convective cores that fully mix the entire object on a fraction of its main sequence life time. We verified that the ad hoc cure to reduce the ratio of vertical to total TKE to zero no longer works once realistic models for that quantity are used. From a physical point of view the ad hoc cure is hence ruled out as an explanation for this deficiency of the model by Kuhfuß (1987). In this paper a physically motivated modification of the mixing length has hence been suggested which takes into account that the dissipation rate of TKE has been underestimated by the original 3-equation model of Kuhfuß (1987). In Paper II we present more detailed tests of the improved 3-equation model proposed in this paper based on stellar evolution tracks for A- and B-type main sequence stars of different masses.

One conclusion from these analyses appears to be that the minimum physics to obtain realistic models of overshooting layers require to account for non-locality of the fluxes of kinetic energy and potential temperature (as intended by Kuhfuß 1987) and in addition to account for the variation of the anisotropy of turbulent kinetic energy as a function of local stability and non-local transport. If the latter is done in a realistic way, it becomes also clear that a physically more complete model of the dissipation rate of TKE is needed. All these features are already provided by the model of Canuto & Dubovikov (1998) which in its most simple form accounts for non-locality with the downgradient approximation (as in the model of Kuhfuß 1987). The present simplification is an attempt to carry over the most important features of the more complete model by Canuto & Dubovikov (1998) into the Kuhfuß (1987) model which is already coded within GARSTEC.

Switching to more complex non-local convection models in a stellar evolution code is not an easy task. This requires that the model and its implementation also account for the following:

  1. Realistic, mathematically self-consistent boundary conditions. This is taken care of in the current implementation of the Kuhfuß (1987) model in GARSTEC.

  2. A fully implicit, relaxation based numerical solver for the resulting set of equations. This is fulfilled by GARSTEC as well. Adding further differential equations always means some non-trivial work on this side.

  3. A stable, monotonic interpolation scheme for the equation of state. Again this is fulfilled in GARSTEC (Weiss & Schlattl 2008). If this is not fulfilled, β cannot be computed correctly and any closure depending on its sign becomes uncertain, since oscillations may be fed into its computation.

  4. A robust formulation of the dynamical equations which avoids cancellation errors introduced through a nearly perfectly adiabatic stratification. This is realised in the implementation of the Kuhfuß (1987) model in GARSTEC indicated by the smoothness of the equation terms in Fig. 8. This can be attributed to the fact that the implementation uses Eq. (A.7) to compute the temperature gradient instead of numerical derivatives.

Naturally, as discussed in Ireland & Browning (2018), in Augustson & Mathis (2019), and in Korre & Featherstone (2021), among others, rotation and magnetic fields influence convection and convective overshooting. A path towards including rotation in non-local convection models has been investigated by Canuto (1998) and by Canuto (2011a), for example, but such extensions have to be left for future work: the present model is only a first step beyond MLT-like models.

If the modified mixing length Eqs. (19) and (20) and even more so Eq. (19) with Eqs. (21)–(22) turns out to produce stable, physically meaningfully evolving overshooting zones with GARSTEC, further tests of this approach are highly warranting. These may also motivate the implementation of fully non-local Reynolds stress models at the complexity level of Canuto & Dubovikov (1998) which completely avoid the introduction of a mixing length with all its shortcomings.


Concerning notation the convention of summation over equally named indices is assumed.


In real world systems the spectra of turbulent kinetic energy, E(k), usually depend on location r and in the most general sense an averaging over directions in k-space would have to be performed, i.e. E(k) becomes a two-point correlation function E(k, r) and would also have to account for density fluctuations.


It is important to note that there is a typo in their Eq. (5c) which should have the constant Ko in the denominator.


In the literature the model discussed here is known as K − ϵ model or “Imperial College model” since there the model had been developed by Hanjalić & Launder (1972).


Note that this definition is different from Canuto & Dubovikov (1998), Eq. (24c), which appears to have a typo.


Entropy gradients in turn are numerically easier to compute than the small differences between temperature and adiabatic temperature gradients during stellar evolution calculations.


If we used the value of cϵ ≈ 2.18 suggested in Kuhfuß (1987) we would instead obtain that c4 ≈ 0.07. However, in the product cϵK3/2/Λ the constant cϵ to some extent cancels out, hence, the overshooting distance is only weakly depending on this parameter. We discuss this further in Appendix B of Paper II.


We point out that exactly the same sequence of arguments applies to the equivalent three equations for the turbulent kinetic energy , the squared fluctuation of the difference between local temperature and its Reynolds average, , and the cross correlation between velocity and temperature fluctuations, , as they appear in the model of Canuto & Dubovikov (1998) and discussed in Kupka et al. (2020).


F. Kupka is thankful to the Austrian Science Fund FWF for support through projects P29172-N and P33140-N and support from European Research Council (ERC) Synergy Grant WHOLESUN #810218. F. Ahlborn thanks Martin Flaskamp for his pioneering work on the 3-equation non-local model by R. Kuhfuß.


  1. Ahlborn, F., Kupka, F., Weiss, A., & Flaskamp, M. 2022, A&A, 667, A97 (Paper II) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Augustson, K. C., & Mathis, S. 2019, ApJ, 874, 83 [Google Scholar]
  3. Augustson, K. C., Brun, A. S., & Toomre, J. 2016, ApJ, 829, 92 [Google Scholar]
  4. Biermann, L. 1932, Z. Astrophys., 5, 117 [NASA ADS] [Google Scholar]
  5. Böhm-Vitense, E. 1958, Z. Astrophys., 46, 108 [Google Scholar]
  6. Bressan, A. G., Chiosi, C., & Bertelli, G. 1981, A&A, 102, 25 [Google Scholar]
  7. Browning, M. K., Brun, A. S., & Toomre, J. 2004, ApJ, 601, 512 [Google Scholar]
  8. Canuto, V. M. 1992, ApJ, 392, 218 [Google Scholar]
  9. Canuto, V. M. 1993, ApJ, 416, 331 [NASA ADS] [CrossRef] [Google Scholar]
  10. Canuto, V. M. 1997a, ApJ, 482, 827 [Google Scholar]
  11. Canuto, V. M. 1997b, ApJ, 489, L71 [NASA ADS] [CrossRef] [Google Scholar]
  12. Canuto, V. M. 1998, ApJ, 508, 767 [NASA ADS] [CrossRef] [Google Scholar]
  13. Canuto, V. M. 2009, in Interdisciplinary Aspects of Turbulence, eds. W. Hillebrandt, & F. Kupka, Lecture Notes in Physics, 756, 107 (Berlin: Springer) [NASA ADS] [CrossRef] [Google Scholar]
  14. Canuto, V. M. 2011a, A&A, 528, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Canuto, V. M. 2011b, A&A, 528, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Canuto, V. M. 2011c, A&A, 528, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Canuto, V. M., & Mazzitelli, I. 1991, ApJ, 370, 295 [NASA ADS] [CrossRef] [Google Scholar]
  18. Canuto, V. M., & Mazzitelli, I. 1992, ApJ, 389, 724 [NASA ADS] [CrossRef] [Google Scholar]
  19. Canuto, V. M., & Dubovikov, M. 1996, Phys. Fluids, 8, 571 [NASA ADS] [CrossRef] [Google Scholar]
  20. Canuto, V. M., & Dubovikov, M. 1998, ApJ, 493, 834 [NASA ADS] [CrossRef] [Google Scholar]
  21. Canuto, V. M., Minotti, F., Ronchi, C., Ypma, R. M., & Zeman, O. 1994, J. Atmos. Sci., 51, 1605 [NASA ADS] [CrossRef] [Google Scholar]
  22. Canuto, V. M., Goldman, I., & Mazzitelli, I. 1996, ApJ, 473, 550 [NASA ADS] [CrossRef] [Google Scholar]
  23. Canuto, V. M., Cheng, Y., & Howard, A. 2001, J. Atmos. Sci., 58, 1169 [NASA ADS] [CrossRef] [Google Scholar]
  24. Canuto, V. M., Cheng, Y., & Howard, A. M. 2010, J. Atmos. Sci., 67, 1678 [NASA ADS] [CrossRef] [Google Scholar]
  25. Chapman, S. 1954, ApJ, 120, 151 [NASA ADS] [CrossRef] [Google Scholar]
  26. Chou, P. Y. 1945, Quart. Appl. Math., 3, 38 [CrossRef] [Google Scholar]
  27. Cunningham, T., Tremblay, P.-E., Freytag, B., Ludwig, H.-G., & Koester, D. 2019, MNRAS, 488, 2503 [NASA ADS] [CrossRef] [Google Scholar]
  28. Daly, B. J., & Harlow, F. H. 1970, Phys. Fluids, 13, 2634 [NASA ADS] [CrossRef] [Google Scholar]
  29. Davidov, B. I. 1961, Dokl. Akad. Nauk. SSSR, 136, 47 [Google Scholar]
  30. Deardorff, J. W. 1961, J. Meteor. (later on: J. Atmos. Sci.), 18, 540 [NASA ADS] [Google Scholar]
  31. Deardorff, J. W. 1966, J. Atmos. Sci., 23, 503 [NASA ADS] [CrossRef] [Google Scholar]
  32. Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019a, ApJ, 876, 4 [NASA ADS] [CrossRef] [Google Scholar]
  33. Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019b, ApJ, 876, 4 [NASA ADS] [CrossRef] [Google Scholar]
  34. Flaskamp, M. 2002, in IAU Colloq. 185: Radial and Nonradial Pulsationsn as Probes of Stellar Physics, eds. C. Aerts, T. R. Bedding, & J. Christensen-Dalsgaard, ASP Conf. Ser., 259, 456 [NASA ADS] [Google Scholar]
  35. Flaskamp, M. 2003, Ph.D. Thesis, Max-Planck-Institut für Astrophysik, Technische Universität München [Google Scholar]
  36. Freytag, B., Ludwig, H.-G., & Steffen, M. 1996, A&A, 313, 497 [Google Scholar]
  37. Gilet, C., Almgren, A. S., Bell, J. B., et al. 2013, ApJ, 773, 137 [NASA ADS] [CrossRef] [Google Scholar]
  38. Gizon, L., & Birch, A. C. 2012, Proc. Natl. Acad. Sci., 109, 11896 [Google Scholar]
  39. Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Annu. Rev. Fluid Mech., 48, 191 [Google Scholar]
  40. Hanjalić, K., & Launder, B. E. 1972, J. Fluid Mech., 52, 609 [CrossRef] [Google Scholar]
  41. Hanjalić, K., & Launder, B. E. 1976, J. Fluid Mech., 74, 593 [CrossRef] [Google Scholar]
  42. Holt, S. E., Koseff, J. R., & Ferziger, J. H. 1992, J. Fluid Mech., 237, 499 [NASA ADS] [CrossRef] [Google Scholar]
  43. Ireland, L. G., & Browning, M. K. 2018, ApJ, 856, 132 [NASA ADS] [CrossRef] [Google Scholar]
  44. Käpylä, P. 2019, A&A, 631, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Keller, L. V., & Friedmann, A. A. 1925, in Proceedings of the First International Congress of Applied Mechanics, eds. C. Biezeno, & J. Burgers, (Delft: J. Waltman Jr.) 395 [Google Scholar]
  46. Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution (Springer Verlag) [Google Scholar]
  47. Kolmogorov, A. N. 1962, J. Fluid Mech., 13, 82 [Google Scholar]
  48. Kolmogorov, A. N. 1968, Sov. Phys. Usp., 10, 734 [NASA ADS] [CrossRef] [Google Scholar]
  49. Korre, L., & Featherstone, N. A. 2021, ApJ, 923, 52 [NASA ADS] [CrossRef] [Google Scholar]
  50. Kuhfuß, R. 1986, A&A, 160, 116 [NASA ADS] [Google Scholar]
  51. Kuhfuß, R. 1987, Ph.D. Thesis, Max-Planck-Institut für Astrophysik Technische Universität München [Google Scholar]
  52. Kumar, P., Talon, S., & Zahn, J.-P. 1999, ApJ, 520, 859 [Google Scholar]
  53. Kupka, F. 2020, in Multi-Dimensional Processes In Stellar Physics, eds. Y. L. M. Rieutord, & I. Baraffe (Les Ulis: EDP Sciences Proceedings), 69 [CrossRef] [Google Scholar]
  54. Kupka, F., & Montgomery, M. H. 2002, MNRAS, 330, L6 [NASA ADS] [CrossRef] [Google Scholar]
  55. Kupka, F., & Muthsam, H. J. 2007, in Convection in Astrophysics, eds. F. Kupka, I. W. Roxburgh, & K. L. Chan, Symposium, 239, 86 (Cambridge University Press) [NASA ADS] [Google Scholar]
  56. Kupka, F., & Muthsam, H. 2017, Liv. Rev. Comput. Astrophys., 3, 159 [Google Scholar]
  57. Kupka, F., Zaussinger, F., & Montgomery, M. H. 2018, MNRAS, 474, 4660 [NASA ADS] [CrossRef] [Google Scholar]
  58. Launder, B. E., Reece, G. J., & Rodi, W. 1975, J. Fluid Mech., 68, 537 [NASA ADS] [CrossRef] [Google Scholar]
  59. Lesieur, M. 2008, Turbulence in Fluids, 4th edn. (Dordrecht: Springer-Verlag) [CrossRef] [Google Scholar]
  60. Li, Y., & Yang, J.-Y. 2001, Chin. J. Astron. Astrophys., 1, 66 [NASA ADS] [CrossRef] [Google Scholar]
  61. Li, Y., & Yang, J. Y. 2007, MNRAS, 375, 388 [NASA ADS] [CrossRef] [Google Scholar]
  62. Linden, P. F. 1975, J. Fluid Mech., 71, 385 [NASA ADS] [CrossRef] [Google Scholar]
  63. Maeder, A., & Mermilliod, J. C. 1981, A&A, 93, 136 [NASA ADS] [Google Scholar]
  64. Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448 [NASA ADS] [CrossRef] [Google Scholar]
  65. Michielsen, M., Aerts, C., & Bowman, D. M. 2021, A&A, 650, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Montgomery, M. H., & Kupka, F. 2004, MNRAS, 350, 267 [NASA ADS] [CrossRef] [Google Scholar]
  67. Pietrinferni, A., Cassisi, S., Salaris, M., & Castelli, F. 2004, ApJ, 612, 168 [Google Scholar]
  68. Pope, S. B. 2000, Turbulent Flows (Cambridge: Cambridge University Press) [Google Scholar]
  69. Reynolds, O. 1894, Philos. Trans. R. Soc. London Ser. A, 186, 123 [NASA ADS] [Google Scholar]
  70. Rogers, T. M. 2015, ApJ, 815, L30 [Google Scholar]
  71. Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772, 21 [NASA ADS] [CrossRef] [Google Scholar]
  72. Roxburgh, I. W. 1989, A&A, 211, 361 [NASA ADS] [Google Scholar]
  73. Roxburgh, I. W., & Kupka, F. 2007a, in Convection in Astrophysics, eds. F. Kupka, & I. W. Roxburgh, Symposium, 239, 98 (Cambridge University Press) [NASA ADS] [Google Scholar]
  74. Roxburgh, I. W., & Kupka, F. 2007b, in Convection in Astrophysics, eds. F. Kupka, I. W. Roxburgh, & K. L. Chan, Symposium, 239, 77 (Cambridge University Press) [NASA ADS] [Google Scholar]
  75. Scheel, J. D., & Schumacher, J. 2017, Phys. Rev. Fluids, 2 [CrossRef] [Google Scholar]
  76. Tennekes, H., & Lumley, J. L. 1972, A First Course in Turbulence (M.I.T. Press) [Google Scholar]
  77. Tremblay, P.-E., Ludwig, H.-G., Freytag, B., et al. 2015, ApJ, 799, 142 [Google Scholar]
  78. Viallet, M., Meakin, C., Arnett, D., & Mocák, M. 2013, ApJ, 769, 1 [NASA ADS] [CrossRef] [Google Scholar]
  79. Viallet, M., Meakin, C., Prat, V., & Arnett, D. 2015, A&A, 580, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Weiss, A., & Schlattl, H. 2008, Ap&SS, 316, 99 [CrossRef] [Google Scholar]
  81. Weiss, A., Hillebrandt, W., Thomas, H.-C., & Ritter, H. 2004, Cox& Giuli’s Principles of Stellar Structure. Extended, 2nd edn. (Cambridge: Cambridge Scientific Publishers) [Google Scholar]
  82. Wuchterl, G. 1995, Comput. Phys. Commun., 89, 119 [NASA ADS] [CrossRef] [Google Scholar]
  83. Xiong, D. R. 1978, Chin. Astron., 2, 118 [NASA ADS] [CrossRef] [Google Scholar]
  84. Xiong, D. R., Cheng, Q. L., & Deng, L. 1997, ApJS, 108, 529 [NASA ADS] [CrossRef] [Google Scholar]
  85. Zahn, J.-P. 1991, A&A, 252, 179 [NASA ADS] [Google Scholar]
  86. Zeman, O., & Lumley, J. L. 1976, J. Atmos. Sci., 33, 1974 [NASA ADS] [CrossRef] [Google Scholar]
  87. Zeman, O., & Tennekes, H. 1977, J. Atmos. Sci., 34, 111 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: The Kuhfuß convection model

In this appendix we summarise the turbulent convection model developed by Kuhfuß (1987) who derived dynamical equations for three of the second order moments to model turbulent convection in the stellar interior: the turbulent kinetic energy, the turbulent convective flux, and the squared entropy fluctuations. The (specific) turbulent kinetic energy (TKE) is denoted by K in the main text. Here, we summarise those equations as used inside GARSTEC. They model entropy fluctuations instead of temperature fluctuations. To avoid confusion with other models and their implementation here we stick to the notation of Kuhfuß (1987): TKE is denoted by ω. The radial component of the turbulent convective flux is written as Π and is computed from entropy fluctuations, consistent with choosing the squared entropy fluctuations Φ as the third dynamical variable of the system. Hence, the Reynolds splitting is performed for

and the second order moments are computed from

As for any TCM a number of assumptions and approximations is required to obtain closed systems of equations that can actually be applied in stellar structure and evolution models. In the following we briefly review the key assumptions of the Kuhfuß model. By using only the total TKE ω the Kuhfuß (1987) model is not able to account for a variable distribution of the kinetic energy in radial and horizontal directions. Instead the distribution of kinetic energy in radial and horizontal directions is assumed to be isotropic at all radii, such that one third of the energy is attributed to each spatial direction. The Kuhfuß model further neglects turbulent pressure fluctuations. As pointed out by Viallet et al. (2013) pressure fluctuations play an important role for convection in envelopes, hence the Kuhfuß model is probably not suited to model envelope convection. Finally, Kuhfuß (1987) also made use of the Boussinesq approximation. In the current implementation suggested by Flaskamp (2003) we also neglect effects due to the chemical composition, e.g. composition gradients.

A.1. Viscous dissipation

In the Kuhfuß model, most terms containing the molecular viscosity are neglected because they are of minor importance compared to competing terms. Only the viscous dissipation term for the kinetic energy is considered to be non-negligible. Kuhfuß (1987) models the dissipation of the kinetic energy with a Kolmogorov-type term (Kolmogorov 1962, 1968):


where CD is a parameter. Kuhfuß (1987) suggests a value of to be compatible with MLT in the local limit of his model.

In the Kolmogorov picture kinetic energy is dissipated thanks to a cascade through which energy is transferred to smaller and smaller spatial scales. The rate at which this dissipation happens is dominated by the largest scales at which energy is fed into the cascade. In Eq. (A.1) the length-scale Λ refers to this largest scale of the turbulent cascade. As in the mixing length theory the length-scale is parametrised using the pressure scale height Hp and an adjustable parameter α: Λ = αHp. Problems with this parametrisation are discussed in the main text.

A.2. Radiative dissipation

Convective elements lose energy through radiation. This is considered in the energy conservation equation by including radiative fluxes as sink terms. In the Kuhfuß equations the radiative losses finally appear as dissipation terms:

where Kuhfuß (1987) models radiative dissipation by introduing the radiative time-scale τrad, which he defines as:

Here, γR is a parameter which Kuhfuß (1987) sets to , again to recover the MLT model in the local limit. Furthermore, cp refers to the specific heat capacity at constant pressure, κ to Rosseland opacity and σ to the Stefan-Boltzmann-constant. The variables T and ρ are temperature and density, as usual in stellar structure models.

A.3. Higher order moments

The Navier-Stokes equations contain non-linear advection terms. When constructing the equations for the second order moments these advection terms give rise to third order moments (TOMs). These higher order moments are the source of the non-local behaviour of the convection model. They can be cast into the form:

where a is a second order quantity. The closure of these TOMs is one of the main challenges of any TCM. Kuhfuß (1987) closes the system of equations at second order and describes each TOM using the so-called down-gradient approximation (e.g., Daly & Harlow 1970; Launder et al. 1975; Xiong 1978; Li & Yang 2007). In the down-gradient approximation the fluxes ja are modelled following Fick’s law:


This approximation is applied for the TOMs appearing in the equations for ω, Π, and Φ with a = v ′2/2, v ′s′, or s2/2. The parameters αa control the impact of the non-local terms. Kuhfuß (1987) suggests a default value of αω ≈ 0.25. The values for the parameters αΠ, Φ are calibrated to MLT in a local version of the Kuhfuß theory. However, no values for the non-local case are provided.

Alternatively, one could compute the TOMs by deriving equations for them in the same way as for the second order moments. This has been shown in Canuto (1992, 1993), Canuto & Dubovikov (1998), or Xiong et al. (1997), for example, and introduces fourth order moments which again have to be closed.

A.4. Final model equations

The above listed approximations are implemented in the derivation of the Kuhfuß model. The final set of partial differential equations reads:




where ∇ and ∇ad refer to the model and adiabatic temperature gradient, respectively. The substantial derivative is defined as . For more details about the derivation we refer to the original work by Kuhfuß (1987) and Flaskamp (2003).

Using the convective flux from the convection model one can compute the temperature gradient of the stellar model self-consistently from



where a and c denote the radiation constant and the speed of light, respectively. Here, we neglect the kinetic energy flux jω, which is assumed to be small compared to the convective flux. Equation (A.7) couples the convection model to the stellar structure equations. The self-consistent computation of the temperature gradient allows us to study its behaviour in the overshooting region. This is an advantage over ad hoc descriptions of overshooting in which the temperature gradient is set manually.

Appendix B: Alternatives to improve Eq. (11)

Eq. (11) is heavily parametrised. Canuto et al. (2009) hence discussed a number of simplified models used in geophysics for the computation of ϵ. They are based on modified mixing lengths which account for physical processes relevant to dissipation. However, those models are not directly applicable to stellar convection: some of them consider a solid wall as a boundary and none of them has been designed for the extreme density contrast of deep stellar envelopes or the peculiarities of convective cores in massive stars.

As the closures used to derive Eq. (9), which were modelled on the basis of turbulent channel flows and freely decaying turbulence, may not be universal, they should ideally be obtained from a more general framework. This approach has been taken in Canuto et al. (2010) who derived a dynamical equation for the TKE dissipation rate ϵ using the general turbulence model of Canuto & Dubovikov (1996). That requires the spectrum of the source driving turbulence to be known. For shear-driven flows power law spectra for the TKE and the Reynolds stress spectrum can readily be specified. Note that these concern scales k < k0, i.e. below the maximum of the TKE spectrum E(k). In addition, energy conservation is invoked which allows computing the non-local contribution to ϵ from the flux of turbulent kinetic energy. That closure was already used in Canuto (1992) (Eq. (37f)) and the non-local character it introduces into Eq. (11) was discussed in Sect. 11 of Canuto (1993). It was tested in Kupka et al. (2007) who found it to be one of the most robust ones among all the closures suggested for the Reynolds stress models of Canuto (1992), Canuto (1993), Canuto & Dubovikov (1998), Canuto et al. (2001), and in Canuto et al. (2009). It specifies that with . In practice, the accuracy of this closure is degraded, if can only be computed from a downgradient approximation, but even in this case it justifies that Df(ϵ) can be evaluated from Df(K) which is required anyway. Hence, Canuto et al. (2010) use the (exact) dynamical equation for the turbulent kinetic energy and a closure for Fkin to compute Df(ϵ). The equivalents of c1 and c2 of Eq. (9) are obtained from within the model, too. The resulting dissipation rate equation passes the same tests as the original Eq. (9) for turbulent channel flow and also two tests concerning the shear dominated planetary boundary of the Earth atmosphere. Unfortunately, this procedure is currently not feasible for the case of convection in stars, since this would require accurate knowledge of the turbulent kinetic energy spectrum over a large range of scales and as a function of depth throughout the star (see also the discussions in Gizon & Birch 2012 and Fig. 5 in Hanasoge et al. 2016 on difficulties in modelling the turbulent kinetic energy spectrum for the Sun).

The dissipation rate equation Eq. (9) has hence remained part of the Reynolds stress model of Canuto (2011a), whether for dealing with double-diffusive convection (Canuto 2011b) or overshooting (Canuto 2011c). The latter paper provides a detailed discussion of the computation of ϵ, which considers Eq. (9) and for non-local contributions. The role of gravity waves as a source of dissipation in the overshooting zone is emphasised, too. From earlier work of Kumar et al. (1999), it is concluded in Canuto (2011b) that ϵ ≈ 10−3 cm2 s−3. However, as also pointed out in Canuto (2011c), it is unclear how this result could be applied to overshooting zones other than the solar tachocline. Thus, in his Eq. (5h), Canuto (2011c) suggests to use to compute τ and hence via τ = 2 K/ϵ the dissipation rate ϵ in the overshooting region. This, however, is consistent with the claim that the term , neglected in the explicit form of the ϵ-equation in Canuto (2011a,b,c), actually dominates in the overshooting region.

Recalling Kupka & Montgomery (2002) and Montgomery & Kupka (2004), who had found the term to dominate the solution of Eq. (11) in their applications of the Reynolds stress model of Canuto & Dubovikov (1998) to overshooting in envelopes of A-stars and white dwarfs and taking into account the confirmation of their results for the case of a DA white dwarf by 3D radiation hydrodynamical simulations in Kupka et al. (2018), Eq. (11) is still the physically most complete model for the computation of ϵ available at the moment. It is thus used to guide the considerations in Sect. 3. Note that the dynamical equation for ϵ which is discussed here does not account for physical effects due to compressibility. Canuto (1997a) has presented several different models to extend Eq. (11) beyond its solenoidal (incompressible) component ϵs and account for a dilation (compressible) contribution ϵd (cf. Sect. 14 in that paper). For current modelling in stellar structure and evolution theory such extensions appear yet too advanced: the very first step is to give up the MLT approach to compute ϵ as specified by Eq. (5)–(6).

All Figures

thumbnail Fig. 1.

TKE as a function of the fractional mass for the original Kuhfuß model. The formal Schwarzschild boundary, defined by ∇rad = ∇ad, is indicated by a dashed black line.

In the text
thumbnail Fig. 2.

Left panels: convective flux in units of total flux, root mean square vertical velocity in units of km s−1, and dissipation rate ϵ from Eq. (11) relative to a value computed from Eqs. (5) and (6) with α as given in the figure legend. Right panels: same quantities as left panels, however, the downgradient approximation (DGA) is used to compute third order moments instead of the full model used in Kupka & Montgomery (2002). The results are for one of the A-star envelope models discussed in Kupka & Montgomery (2002).

In the text
thumbnail Fig. 3.

Dissipation rate as a function of fractional mass for the original Kuhfuß model and the Kuhfuß model including an ad hoc exponential decay of the dissipation length, shown with a grey dotted and a blue continuous line, respectively. The ad hoc exponential decay of the dissipation length leads to an increased dissipation rate at the beginning of the overshooting zone, indicated by the local maximum beyond the Schwarzschild boundary, followed by a sharp drop due to the rapid decay of TKE. The models have been chosen to have the same maximum TKE.

In the text
thumbnail Fig. 4.

Ratio of τ/τb as a function of convective stability from a solution of the non-local Reynolds stress model as presented in Kupka & Montgomery (2002) assuming the downgradient approximation for third order moments. The time scale τb is computed from where the absolute value of β is taken. Sign changes are hence indicated by spikes. Both the overshooting zones below and above the lower and the upper convectively unstable zone show the same increase of τ/τb from 0 to more than 10 (the finite grid resolution prevents τ/τb from becoming actually zero).

In the text
thumbnail Fig. 5.

Artificial anisotropy factor ξ and TKE as a function of fractional mass in the upper and lower panel, respectively. The black dashed line indicates the Schwarzschild boundary.

In the text
thumbnail Fig. 6.

Estimate of the anisotropy factor according to Eq. (24) for a 3-equation model without limited dissipation length-scale Λ. The profile of the turbulent kinetic energy of this model is shown in Fig. 1.

In the text
thumbnail Fig. 7.

Convective energy as a function of the fractional mass for the Kuhfuß model including the improved dissipation mechanism. The Schwarzschild boundary is indicated by a dashed black line.

In the text
thumbnail Fig. 8.

Comparison of the different terms in the TKE equation (Eq. (A.4)) in the Kuhfuß 1-equation (panel a), original 3-equation (panel b) and improved 3-equation (panel c) model. The buoyant driving term, the dissipation term and the non-local flux term are shown with a red, black, and blue line here.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.