Issue 
A&A
Volume 667, November 2022



Article Number  A96  
Number of page(s)  17  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202243125  
Published online  11 November 2022 
Stellar evolution models with overshooting based on 3equation nonlocal theories
I. Physical basis and the computation of the dissipation rate
^{1}
Dept. Applied Mathematics and Physics, Univ. of Applied Sciences, Technikum Wien, Höchstädtplatz 6, 1200 Wien, Austria
email: friedrich.kupka@technikumwien.at
^{2}
WolfgangPauliInstitute c/o Faculty of Mathematics, University of Vienna, OskarMorgensternPlatz 1, 1090 Wien, Austria
^{3}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
^{4}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStraße 1, 85748 Garching, Germany
Received:
15
January
2022
Accepted:
18
July
2022
Context. Mixing by convective overshooting has long been suggested to play an important role in the amount of hydrogen available for nuclear burning in convective cores of stars. The best way to model this effect is still debated.
Aims. We suggest an improved model for the computation of the dissipation rate of turbulent kinetic energy which can be used in nonlocal models of stellar convection and can readily be implemented and selfconsistently used in 1D stellar evolution calculations.
Methods. We review the physics underlying various models to compute the dissipation rate of turbulent kinetic energy, ϵ, in local models of convection in stellar astrophysics and particularly in nonlocal ones. The different contributions to the dissipation rate and their dependence on local stratification and on nonlocal transport are analysed and a new method to account for at least some of these physical mechanisms is suggested.
Results. We show how the new approach influences predictions of stellar models of intermediatemass mainsequence stars and how these changes differ from other modifications of the nonlocal convection model that focus on the ratio of horizontal to vertical (turbulent) kinetic energy.
Conclusions. The new model is shown to allow for a physically more complete description of convective overshooting and mixing in massive stars. Dissipation by buoyancy waves is found to be a key ingredient which has to be accounted for in nonlocal models of turbulent convection.
Key words: convection / turbulence / stars: evolution / stars: interiors
© F. Kupka et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), 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 SubscribetoOpen 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 onedimensional (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 intermediatemass mainsequence 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öhmVitense 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 nonlocal situations, encountered in layers adjacent to convectively unstable zones. To treat the dissipation of TKE in nonlocal 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 nonlocal descriptions of the dissipation of TKE in Sect. 2. From the dissipation rate in nonlocal 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 nonnegligible 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, u_{i} is the ith velocity component, x_{i} is the ith component of location, and E(k) is the spectrum of turbulent kinetic energy as a function of wavenumber k^{1}. 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 convection^{2}. In a quasistationary 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 kspace. 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 onepoint 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 k^{1/3} and a region where the Kolmogorov inertial range no longer holds, just around the dissipation scale k_{d} = π/ℓ_{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 onepoint 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 ∫k^{2} E(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 k_{0} onwards, i.e. the entire energy spectrum is given by a Kolmogorov spectrum with an energy cutoff for k < k_{0}. Thus, E(k) = 0 for k < k_{0} and E(k)∼k^{−5/3} for arbitrarily large k with k ≥ k_{0}. 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} = π/k_{0} 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 l_{0} such that its integral yields Eq. (5). It is thus a ‘oneeddy 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 MLTparameter or mixing length parameter and H_{p} 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 oneeddy 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 nonlocal: the only way to sustain a nonvanishing 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 nonlocally, 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 3equation 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 nonlocal models: the dissipation rate equation
A common starting point for nonlocal 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 C_{ii} 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^{−2}∂_{r} r^{2}). This is a prognostic equation for the second order moment with derived directly from the Boussinesq approximation of the NavierStokes equations through ensemble averaging. The nonlocal 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 crosscorrelation 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}, c_{p}, 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) emphasised^{4} 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 constantstress 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 D_{f}(ϵ) was suggested to be parametrised as
where ν_{t} requires a model for turbulent viscosity such as^{5}ν_{t} = C_{μ} K^{2}/ϵ 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 nonstellar 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 c_{1} = 1.44 and c_{2} = 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 selfstretching in turbulent flows and the decay of turbulence due to viscosity. For the exact term of diffusion of ϵ by velocity fluctuations, D_{f}(ϵ), both a downgradient 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 BruntVä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 w^{2} 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 c_{3} = 0.3 for stably stratified layers and c_{3} = 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 nonlocal convection models in GARSTEC
3.1. The problem: Overshooting zones of convective cores growing unlimitedly during mainsequence 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öhmVitense 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 1equation version, i.e. with an additional differential equation for turbulent kinetic energy, K, and in its full, 3equation 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 nonlocal fluxes, and (iv) some minor simplifications in the closures used in the dynamical equations^{6}. As a variant, the 3equation model may be used with local limit expressions for the nonlocal transport terms for as well as J. As a theoretical analysis shows (see Kupka et al. 2020 and references therein) only a full 3equation 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 nonlocal transport of as originally shown by Deardorff 1961, 1966). However, in both the fully nonlocal and the local limit of the 3equation model variant as described above, overshooting gradually mixes the entire star in a stellar evolution calculation for a 5 solar mass (Btype) mainsequence 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 mainsequence, i.e. to the left in the colourmagnitude 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).
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 mainsequence 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 nonlocal 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 nonlocal 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 Atype mainsequence 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 T_{eff} = 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 doubleionisation 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 nonlocal 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 T_{eff} = 8000 K which allows observing this behaviour of ϵ even between the two convective zones. At lower T_{eff}, 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 T_{eff} = 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 T_{eff} = 7200 K.
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 Astar 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 Atype stars: irrespective of the various situations described above, deep inside the plumedominated 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.
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 efolding extent of the “decay” of the mixing length of 2% to 6% of the mass of the Schwarzschildunstable region. In a 5 M_{⊙} mainsequence 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 mainsequence. For an efolding extent of 4% the overshooting region contains about 5% of the stellar mass at the beginning of the mainsequence while it is shrinking to about 2% of the total mass at the end of the mainsequence. 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 nonlocal 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 nonzero 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_{ϵ}K^{3/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, H_{c}, of pressure, P_{c}, and density, ρ_{c}, at the centre (and G is, of course, the gravitational constant). For the centre, Λ = αH_{c} and in general Λ = αmin(H_{p}, H_{c}). Roxburgh et al. (2007a) suggest a smooth interpolation between the limit at the core centre and the expression for H_{p} ≫ H_{c}.
For reasons of regularity and energy conservation, F_{conv} → 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 F_{conv} (Roxburgh et al. 2007b). This implies nontrivial 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 ∝ r^{3} 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 nonlocal 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 D_{f}(ϵ), the nonlocal 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 D_{f}(ϵ)∝ − α_{ϵ}ϵ/τ. The same can be done also in the case of Eq. (9). If we furthermore assume the local limit of Eq. (8), P = P_{b} = ϵ, i.e. production of TKE by buoyancy equals its dissipation, and if we also assume c_{3} = 0, we obtain the following approximation for both Eqs. (9) and (11):
To remain consistent with gαJ = ϵ we have to require that α_{ϵ} = 2c_{2} − 2c_{1} 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 P_{b} = ϵ in Eq. (8) somewhat. In this case, whether the 1equation or the 3equation version of the Kuhfuß (1987) model is used (cf. Appendix A), due to the nonlocality 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, α_{ϵ} = 2c_{2} is required from Eq. (12) for a nonvanishing dissipation rate ϵ. Right next to such a point, where ϵ > 0 with J < 0, a value of α_{ϵ} > 2c_{2} would be required whereas α_{ϵ} < 2c_{2} where J > 0. So α_{ϵ} would have to be a function that has to be finetuned 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 c_{3} = 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 1equation model of Kuhfuß (1986), β and J change sign at the same point so the perfect balancing constraint between D_{f}(ϵ) and −2 c_{2}ϵ/τ reappears. In the region where J < 0, more freedom of how D_{f}(ϵ) behaves is permitted. This changes once we switch to the 3equation model of Kuhfuß (1987): since β and J then change sign at different locations, α_{ϵ} is no longer forced by c_{2} at any point. In the end, the contribution decouples both D_{f}(ϵ) 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 3equation 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 3equation 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 D_{f}(ϵ). 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 c_{1}gαJ/τ into −2 c_{2}ϵ/τ for sufficiently small J and we obtain from Eqs. (13), (16), and (17) that
where the numerical value is obtained from setting c_{2} = 1.92 and c_{3} = 0.3. Contributions absorbed into the −2 c_{2}ϵ/τ term could be accounted for by small change of c_{2}. 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 M_{schw} 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 nonlocal Reynolds stress model of Canuto & Dubovikov (1998) for Atype stars (Kupka & Montgomery 2002) show that τ_{b} ≈ 0.1 τ where F_{conv} reaches its negative minimum. This range of values for τ_{b} is what we also expect from Eq. (18) for a moderate variation of c_{2}.
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 F_{conv} < 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.
Fig. 4. Ratio of τ/τ_{b} as a function of convective stability from a solution of the nonlocal 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_{ϵ})(αH_{p}K^{−1/2}) or τ = (2/c_{ϵ})(rK^{−1/2}), as this formula is to be used only for r > 0 and H_{p} < ∞ 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 Λ ≈ αH_{p} right there which is exactly not what we want. But we can rewrite Eq. (20) into
with
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 c_{4} could be adjusted to achieve the goal of τ_{b} ≈ 0.1 τ or rather Λ(min(F_{conv})) ≈ 0.1Λ(M_{r} = M_{schw}) 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 c_{4} 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ß 3equation 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 nonlocal 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 3equation 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.
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 timedependence in Eq. (23) we can define an anisotropy factor:
All quantities in Eq. (24) can be computed within the Kuhfuß 3equation model.
We have computed the anisotropy factor according to Eq. (24) for a stellar model which used the original version of the Kuhfuß 3equation 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).
Fig. 6. Estimate of the anisotropy factor according to Eq. (24) for a 3equation model without limited dissipation lengthscale Λ. The profile of the turbulent kinetic energy of this model is shown in Fig. 1. 
Introducing this anisotropy factor into the Kuhfuß 3equation 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) selfconsistently coupled to the nonlocal convection model. However, since such a more realistic anisotropy factor cannot resolve the problem of excessive mixing found in the original Kuhfuß 3equation 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 3equation model in Paper II and postpone the extension of this new model to future work.
4.2. Dissipation in the Kuhfuß 1 and 3equation 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) selfconsistently. We note that for consistency and to simplify the comparison between the 1equation and the 3equation model, we set c_{ϵ} = C_{D} (see Appendix A), whence it follows that c_{4} ≈ 0.072 in those calculations. As an example we show here the TKE in a 5 M_{⊙} mainsequence 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.
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 3equation 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 nonlocal flux–with a corresponding red, black, and blue line respectively for the 1equation model (panel a), the original 3equation model (panel b) and the improved 3equation model (panel c).
Fig. 8. Comparison of the different terms in the TKE equation (Eq. (A.4)) in the Kuhfuß 1equation (panel a), original 3equation (panel b) and improved 3equation (panel c) model. The buoyant driving term, the dissipation term and the nonlocal flux term are shown with a red, black, and blue line here. 
Stellar models applying the nonlocal 1equation 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 3equation model in its original version this welcome property vanishes and the stellar models become fully convective. As discussed in Appendix A the 3equation 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 1equation 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 3equation model look physically unreasonable, as the existence of fully convective Bstars with 5 M_{⊙} is excluded from the lack of stars hotter than the hydrogen mainsequence.
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 3equation models shown in panels a) and b) in Fig. 8. In the 1equation 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 1equation 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 3equation model the buoyant term is still negative in the overshooting zone. The values are, however, much smaller in magnitude compared to the 1equation model. The dissipation and nonlocal flux term have about the same magnitude in the overshooting zone as obtained with the 1equation 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 3equation 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 1equation model. Kuhfuß (1987) has approximated this to be Π ∝ (∇−∇_{ad}). As the convective flux is the major sink term in the overshooting zone in the 1equation 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 3equation case. The transition from finite to zero values looks smoother for the improved 3equation model because the temperature gradient which has readjusted differs in comparison with the 1equation model.
This comparison shows why the original version of the 3equation 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 nonlocal 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∇_{ad}T/H_{p})Φ act as sources which are boosted in Eq. (A.5). Unless we would consider a large rate of change in the nonlocal 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 Bstar, 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 3equation 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 Astars (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 3equation model to remove the deficiency the model has had in its original version proposed by Kuhfuß (1987).
4.3. Comparing the Kuhfuß 3equation 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 cutoff 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ß 3equation 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ß 3equation 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 selfconsistent predictions which effectively turns it into a model similar in complexity to the 1equation model by Kuhfuß (1987). We recall here that the 3equation model with enhanced dissipation has a builtin dependence on Pe by accounting for radiative losses in its dynamical equations. A detailed discussion on the role of Pe in the 1equation and 3equation 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ß 3equation 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 3equation 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ß 3equation 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 nonlocal 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 3equation model of Kuhfuß (1987). In Paper II we present more detailed tests of the improved 3equation model proposed in this paper based on stellar evolution tracks for A and Btype 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 nonlocality 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 nonlocal 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 nonlocality 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 nonlocal 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:

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

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 nontrivial work on this side.

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.

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 nonlocal 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 MLTlike 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 nonlocal Reynolds stress models at the complexity level of Canuto & Dubovikov (1998) which completely avoid the introduction of a mixing length with all its shortcomings.
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 kspace would have to be performed, i.e. E(k) becomes a twopoint correlation function E(k, r) and would also have to account for density fluctuations.
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.
If we used the value of c_{ϵ} ≈ 2.18 suggested in Kuhfuß (1987) we would instead obtain that c_{4} ≈ 0.07. However, in the product c_{ϵ}K^{3/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).
Acknowledgments
F. Kupka is thankful to the Austrian Science Fund FWF for support through projects P29172N and P33140N and support from European Research Council (ERC) Synergy Grant WHOLESUN #810218. F. Ahlborn thanks Martin Flaskamp for his pioneering work on the 3equation nonlocal model by R. Kuhfuß.
References
 Ahlborn, F., Kupka, F., Weiss, A., & Flaskamp, M. 2022, A&A, 667, A97 (Paper II) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Augustson, K. C., & Mathis, S. 2019, ApJ, 874, 83 [Google Scholar]
 Augustson, K. C., Brun, A. S., & Toomre, J. 2016, ApJ, 829, 92 [Google Scholar]
 Biermann, L. 1932, Z. Astrophys., 5, 117 [NASA ADS] [Google Scholar]
 BöhmVitense, E. 1958, Z. Astrophys., 46, 108 [Google Scholar]
 Bressan, A. G., Chiosi, C., & Bertelli, G. 1981, A&A, 102, 25 [Google Scholar]
 Browning, M. K., Brun, A. S., & Toomre, J. 2004, ApJ, 601, 512 [Google Scholar]
 Canuto, V. M. 1992, ApJ, 392, 218 [Google Scholar]
 Canuto, V. M. 1993, ApJ, 416, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V. M. 1997a, ApJ, 482, 827 [Google Scholar]
 Canuto, V. M. 1997b, ApJ, 489, L71 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V. M. 1998, ApJ, 508, 767 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Canuto, V. M. 2011a, A&A, 528, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Canuto, V. M. 2011b, A&A, 528, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Canuto, V. M. 2011c, A&A, 528, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Canuto, V. M., & Mazzitelli, I. 1991, ApJ, 370, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V. M., & Mazzitelli, I. 1992, ApJ, 389, 724 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V. M., & Dubovikov, M. 1996, Phys. Fluids, 8, 571 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V. M., & Dubovikov, M. 1998, ApJ, 493, 834 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V. M., Minotti, F., Ronchi, C., Ypma, R. M., & Zeman, O. 1994, J. Atmos. Sci., 51, 1605 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V. M., Goldman, I., & Mazzitelli, I. 1996, ApJ, 473, 550 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V. M., Cheng, Y., & Howard, A. 2001, J. Atmos. Sci., 58, 1169 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V. M., Cheng, Y., & Howard, A. M. 2010, J. Atmos. Sci., 67, 1678 [NASA ADS] [CrossRef] [Google Scholar]
 Chapman, S. 1954, ApJ, 120, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Chou, P. Y. 1945, Quart. Appl. Math., 3, 38 [CrossRef] [Google Scholar]
 Cunningham, T., Tremblay, P.E., Freytag, B., Ludwig, H.G., & Koester, D. 2019, MNRAS, 488, 2503 [NASA ADS] [CrossRef] [Google Scholar]
 Daly, B. J., & Harlow, F. H. 1970, Phys. Fluids, 13, 2634 [NASA ADS] [CrossRef] [Google Scholar]
 Davidov, B. I. 1961, Dokl. Akad. Nauk. SSSR, 136, 47 [Google Scholar]
 Deardorff, J. W. 1961, J. Meteor. (later on: J. Atmos. Sci.), 18, 540 [NASA ADS] [Google Scholar]
 Deardorff, J. W. 1966, J. Atmos. Sci., 23, 503 [NASA ADS] [CrossRef] [Google Scholar]
 Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019a, ApJ, 876, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019b, ApJ, 876, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Flaskamp, M. 2002, in IAU Colloq. 185: Radial and Nonradial Pulsationsn as Probes of Stellar Physics, eds. C. Aerts, T. R. Bedding, & J. ChristensenDalsgaard, ASP Conf. Ser., 259, 456 [NASA ADS] [Google Scholar]
 Flaskamp, M. 2003, Ph.D. Thesis, MaxPlanckInstitut für Astrophysik, Technische Universität München [Google Scholar]
 Freytag, B., Ludwig, H.G., & Steffen, M. 1996, A&A, 313, 497 [Google Scholar]
 Gilet, C., Almgren, A. S., Bell, J. B., et al. 2013, ApJ, 773, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., & Birch, A. C. 2012, Proc. Natl. Acad. Sci., 109, 11896 [Google Scholar]
 Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Annu. Rev. Fluid Mech., 48, 191 [Google Scholar]
 Hanjalić, K., & Launder, B. E. 1972, J. Fluid Mech., 52, 609 [CrossRef] [Google Scholar]
 Hanjalić, K., & Launder, B. E. 1976, J. Fluid Mech., 74, 593 [CrossRef] [Google Scholar]
 Holt, S. E., Koseff, J. R., & Ferziger, J. H. 1992, J. Fluid Mech., 237, 499 [NASA ADS] [CrossRef] [Google Scholar]
 Ireland, L. G., & Browning, M. K. 2018, ApJ, 856, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. 2019, A&A, 631, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 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]
 Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution (Springer Verlag) [Google Scholar]
 Kolmogorov, A. N. 1962, J. Fluid Mech., 13, 82 [Google Scholar]
 Kolmogorov, A. N. 1968, Sov. Phys. Usp., 10, 734 [NASA ADS] [CrossRef] [Google Scholar]
 Korre, L., & Featherstone, N. A. 2021, ApJ, 923, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Kuhfuß, R. 1986, A&A, 160, 116 [NASA ADS] [Google Scholar]
 Kuhfuß, R. 1987, Ph.D. Thesis, MaxPlanckInstitut für Astrophysik Technische Universität München [Google Scholar]
 Kumar, P., Talon, S., & Zahn, J.P. 1999, ApJ, 520, 859 [Google Scholar]
 Kupka, F. 2020, in MultiDimensional Processes In Stellar Physics, eds. Y. L. M. Rieutord, & I. Baraffe (Les Ulis: EDP Sciences Proceedings), 69 [CrossRef] [Google Scholar]
 Kupka, F., & Montgomery, M. H. 2002, MNRAS, 330, L6 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Kupka, F., & Muthsam, H. 2017, Liv. Rev. Comput. Astrophys., 3, 159 [Google Scholar]
 Kupka, F., Zaussinger, F., & Montgomery, M. H. 2018, MNRAS, 474, 4660 [NASA ADS] [CrossRef] [Google Scholar]
 Launder, B. E., Reece, G. J., & Rodi, W. 1975, J. Fluid Mech., 68, 537 [NASA ADS] [CrossRef] [Google Scholar]
 Lesieur, M. 2008, Turbulence in Fluids, 4th edn. (Dordrecht: SpringerVerlag) [CrossRef] [Google Scholar]
 Li, Y., & Yang, J.Y. 2001, Chin. J. Astron. Astrophys., 1, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Li, Y., & Yang, J. Y. 2007, MNRAS, 375, 388 [NASA ADS] [CrossRef] [Google Scholar]
 Linden, P. F. 1975, J. Fluid Mech., 71, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A., & Mermilliod, J. C. 1981, A&A, 93, 136 [NASA ADS] [Google Scholar]
 Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448 [NASA ADS] [CrossRef] [Google Scholar]
 Michielsen, M., Aerts, C., & Bowman, D. M. 2021, A&A, 650, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Montgomery, M. H., & Kupka, F. 2004, MNRAS, 350, 267 [NASA ADS] [CrossRef] [Google Scholar]
 Pietrinferni, A., Cassisi, S., Salaris, M., & Castelli, F. 2004, ApJ, 612, 168 [Google Scholar]
 Pope, S. B. 2000, Turbulent Flows (Cambridge: Cambridge University Press) [Google Scholar]
 Reynolds, O. 1894, Philos. Trans. R. Soc. London Ser. A, 186, 123 [NASA ADS] [Google Scholar]
 Rogers, T. M. 2015, ApJ, 815, L30 [Google Scholar]
 Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Roxburgh, I. W. 1989, A&A, 211, 361 [NASA ADS] [Google Scholar]
 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]
 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]
 Scheel, J. D., & Schumacher, J. 2017, Phys. Rev. Fluids, 2 [CrossRef] [Google Scholar]
 Tennekes, H., & Lumley, J. L. 1972, A First Course in Turbulence (M.I.T. Press) [Google Scholar]
 Tremblay, P.E., Ludwig, H.G., Freytag, B., et al. 2015, ApJ, 799, 142 [Google Scholar]
 Viallet, M., Meakin, C., Arnett, D., & Mocák, M. 2013, ApJ, 769, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Viallet, M., Meakin, C., Prat, V., & Arnett, D. 2015, A&A, 580, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weiss, A., & Schlattl, H. 2008, Ap&SS, 316, 99 [CrossRef] [Google Scholar]
 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]
 Wuchterl, G. 1995, Comput. Phys. Commun., 89, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Xiong, D. R. 1978, Chin. Astron., 2, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Xiong, D. R., Cheng, Q. L., & Deng, L. 1997, ApJS, 108, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J.P. 1991, A&A, 252, 179 [NASA ADS] [Google Scholar]
 Zeman, O., & Lumley, J. L. 1976, J. Atmos. Sci., 33, 1974 [NASA ADS] [CrossRef] [Google Scholar]
 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 nonnegligible. Kuhfuß (1987) models the dissipation of the kinetic energy with a Kolmogorovtype term (Kolmogorov 1962, 1968):
where C_{D} 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 lengthscale Λ refers to this largest scale of the turbulent cascade. As in the mixing length theory the lengthscale is parametrised using the pressure scale height H_{p} and an adjustable parameter α: Λ = αH_{p}. 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 timescale τ_{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, c_{p} refers to the specific heat capacity at constant pressure, κ to Rosseland opacity and σ to the StefanBoltzmannconstant. The variables T and ρ are temperature and density, as usual in stellar structure models.
A.3. Higher order moments
The NavierStokes equations contain nonlinear 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 nonlocal 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 socalled downgradient approximation (e.g., Daly & Harlow 1970; Launder et al. 1975; Xiong 1978; Li & Yang 2007). In the downgradient approximation the fluxes j_{a} 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 s′^{2}/2. The parameters α_{a} control the impact of the nonlocal 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 nonlocal 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 selfconsistently from
with
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 selfconsistent 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 sheardriven flows power law spectra for the TKE and the Reynolds stress spectrum can readily be specified. Note that these concern scales k < k_{0}, i.e. below the maximum of the TKE spectrum E(k). In addition, energy conservation is invoked which allows computing the nonlocal contribution to ϵ from the flux of turbulent kinetic energy. That closure was already used in Canuto (1992) (Eq. (37f)) and the nonlocal 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 D_{f}(ϵ) can be evaluated from D_{f}(K) which is required anyway. Hence, Canuto et al. (2010) use the (exact) dynamical equation for the turbulent kinetic energy and a closure for F_{kin} to compute D_{f}(ϵ). The equivalents of c_{1} and c_{2} 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 doublediffusive convection (Canuto 2011b) or overshooting (Canuto 2011c). The latter paper provides a detailed discussion of the computation of ϵ, which considers Eq. (9) and for nonlocal 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} cm^{2} 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 Astars 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
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 
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 Astar envelope models discussed in Kupka & Montgomery (2002). 

In the text 
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 
Fig. 4. Ratio of τ/τ_{b} as a function of convective stability from a solution of the nonlocal 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 
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 
Fig. 6. Estimate of the anisotropy factor according to Eq. (24) for a 3equation model without limited dissipation lengthscale Λ. The profile of the turbulent kinetic energy of this model is shown in Fig. 1. 

In the text 
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 
Fig. 8. Comparison of the different terms in the TKE equation (Eq. (A.4)) in the Kuhfuß 1equation (panel a), original 3equation (panel b) and improved 3equation (panel c) model. The buoyant driving term, the dissipation term and the nonlocal flux term are shown with a red, black, and blue line here. 

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