© ESO 2019

1. Introduction

Convective mixing in stars has important consequences, for example, in early and late phases of stellar evolution and for the diffusion of light elements. Furthermore, stellar differential rotation (e.g. Rüdiger 1989) and dynamos (e.g. Moffatt 1978; Krause & Rädler 1980; Brandenburg & Subramanian 2005) owe their existence to turbulent fluid motions. The efficiency of mixing in convective and radiative layers in stars differs greatly. The former are vigorously mixed on a timescale much shorter than the evolutionary timescale of the star. Thus it is of great interest to be able to predict where effective convective mixing occurs. The greatest uncertainty in this respect is the amount of overshooting from convection zones (CZ) to adjacent radiative layers.

Stellar structure and evolution models most often apply some variant of the mixing length (ML) model of Vitense (1953) to describe convection. These models are completely local and do not allow overshooting. Non-local extensions to ML models (e.g. Shaviv & Salpeter 1973; Schmitt et al. 1984; Skaley & Stix 1991) yield estimates of overshooting, but the validity of the ML approach has been questioned (e.g. Renzini 1987). More advanced closures of convection based on Reynolds stress (e.g. Xiong 1985; Deng et al. 2006; Garaud et al. 2010; Canuto 2011) are physically more consistent but challenging to implement (see, however, Zhang et al. 2012; Zhang 2013). Furthermore, testing and validation of the Reynolds stress models, for example, by comparison to three-dimensional numerical simulations, is still in its infancy (e.g. Kupka 1999; Snellman et al. 2015; Cai 2018).

A seemingly attractive option to study overshooting is to solve the governing equations directly by means of three-dimensional simulations. Numerical simulations of convection have been used to estimate the overshooting depth in numerous studies (e.g. Hurlburt et al. 1986, 1994; Roxburgh & Simmons 1993; Singh et al. 1995, 1998; Saikia et al. 2000; Brummell et al. 2002; Ziegler & Rüdiger 2003; Rogers et al. 2006; Tian et al. 2009; Pratt et al. 2017; Brun et al. 2017; Hotta 2017; Korre et al. 2019). Early studies indicated overshooting of the order of a pressure scale height at the base of the CZ, which is an order of magnitude more than typical estimates from helioseismology (e.g. Basu 1997). The difference between these studies and the Sun is that the energy flux imposed in the simulations is typically much greater than the corresponding solar flux. This leads to higher convective velocities in the simulations and to an overestimation of the overshooting. Scaling laws, based on the relation of convective velocities with the energy flux, arise in analytic models of overshooting (e.g. Zahn 1991; Rempel 2004) and predict a reduction of the overshooting depth as a function of decreasing flux. The primary aim of the current study is to establish these relations from carefully controlled numerical experiments.

Compressible simulations with realistic solar energy flux are hampered by the disparity of the acoustic and dynamical timescales, or by a low Mach number, in the deep parts of the CZ. According to ML arguments, the enthalpy flux ${\overline{\mathit{F}}}_{\mathrm{enth}}=c_{\mathrm{P}}\overline{(\rho \mathit{u})\prime T\prime }$, where c_P is the heat capacity at constant pressure, ρ is the density, u is the convective velocity, and T is the temperature, and the apostrophes (overbar) denote fluctuations (averages) that can be approximated as ${\overline{F}}_{\mathrm{enth}}=\varphi \overline{\rho }{u}^{\prime 3}$ (Brandenburg 2016). Assuming that ${\overline{F}}_{\mathrm{enth}}\approx F_{\odot }$, it is possible to construct a normalised energy flux (e.g. Brandenburg et al. 2005),

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{n}}^{(\odot )}=F_{\odot }/\rho c_{\mathrm{s}}^3\approx \varphi \rho u{\prime }^3/\rho c_{\mathrm{s}}^3\approx {\mathrm{Ma}}^3.\end{array}$$(1)

The last approximation is justified by the fact that the factor ϕ has been reported to be in the range 4 − 20 (Brandenburg et al. 2005; Brandenburg 2016), whereas convection carries only a fraction of the total flux in the lower part of the solar CZ (see, e.g. the solar model of Stix 2002). Using solar values at r₀ = 0.71 R_⊙, where R_⊙ = 7 × 10⁸ m is the solar radius, $F_{\odot }=L_{\odot }/(4\pi r_0^2)\approx 6\times {10}^{25}$ W m⁻², where L_⊙ = 3.84 × 10²⁶ W is the solar luminosity, ρ ≈ 200 kg m⁻³, and c_s ≈ 200 km s⁻¹, the normalised flux at the base of the solar CZ is ${\mathcal{F}}_{\mathrm{n}}^{(\odot )}\approx 4\times {10}^{-11}$. Thus the Mach number in the deep CZ is about 10⁻⁴. This leads to a short time-step due to the high sound speed and to prohibitively long integration times (e.g. Kupka & Muthsam 2017). Although anelastic methods (e.g Gough 1969) bypass the acoustic time-step problem at solar luminosity, any simulation attempting to do this self-consistently would need to be run for a Kelvin–Helmholtz time to achieve thermal saturation. This is not feasible for solar parameters without resorting, for example, to arbitrarily changing the heat conductivity in the radiative layer (e.g. Brun et al. 2017) with current and any foreseeable supercomputers (e.g. Kupka & Muthsam 2017).

Recently, Hotta (2017) presented results from numerical simulations of fully compressible convection where the overshooting depth was computed from a range of two orders of magnitude in the input flux. He reached values of ℱ_n = 5 × 10⁻⁷, which is at the limits of current numerical feasibility, and obtained a power law $d_{\mathrm{os}}/H_{\mathrm{p}}\propto {\mathcal{F}}_{\mathrm{n}}^{0.31}$ for the overshooting depth d_os. This led Hotta (2017) to estimate that the overshooting at the base of the solar CZ is about 0.4% of the pressure scale height, or roughly 200 km. Earlier numerical studies (e.g. Singh et al. 1998; Tian et al. 2009) have reported results that suggest a similar steep dependency of d_os on ℱ_n.

Here these studies are revisited by a set-up where the heat conductivity is self-consistently computed using the Kramers opacity law. This set-up allows the depth of the CZ to dynamically adapt to changes in the thermodynamic state of the system (Käpylä et al. 2019a) and to produce a smooth transition between convective and radiative layers (Brandenburg et al. 2000; Käpylä et al. 2017). Furthermore, a significantly broader range of imposed flux values is covered than in any of the previous studies. Moreover, particular care is taken to isolate the effect of the input flux by performing models where the supercriticality of convection, degree of turbulence, and effective thermal Prandlt number are approximately constant as the flux varies. An effort is made to link to the earlier studies of Singh et al. (1998) and Tian et al. (2009) by targeted sets of simulations probing the influence of subgrid scale entropy diffusion on convection and the resulting overshooting depth. Finally, a critical assessment of some of the modeling choices of Hotta (2017) is presented.

The PENCIL CODE¹ was used to produce the simulations. At the core of the code is a switchable finite-difference solver for partial differential equations that can be used to study a wide selection of physical problems. In the current study a third-order Runge–Kutta time-stepping method and centred sixth-order finite differences for spatial derivatives are used (cf. Brandenburg 2003).

2. The model

The Kramers set-ups used in the current study are similar to those in Käpylä et al. (2017). The equations for compressible hydrodynamics

$$\begin{array}{cc}& \frac{\mathrm{D}ln\rho }{\mathrm{D}t}=-\mathbf{\nabla }\mathbf{·}\mathit{u},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\mathrm{D}\mathit{u}}{\mathrm{D}t}=\mathit{g}-\frac{1}{\rho }(\mathbf{\nabla }p-\mathbf{\nabla }\mathbf{·}2\nu \rho \mathsf{S}),\hfill \end{array}$$(3)

$$\begin{array}{cc}& T\frac{\mathrm{D}s}{\mathrm{D}t}=-\frac{1}{\rho }\left[\mathbf{\nabla }\mathbf{·}({\mathit{F}}_{\mathrm{rad}}+{\mathit{F}}_{\mathrm{SGS}})\right]+2\nu {\mathsf{S}}^2+{\mathrm{\Gamma }}_{\mathrm{cool}},\hfill \end{array}$$(4)

are solved, where D/Dt = ∂/∂t + u ⋅ ∇ is the advective derivative, ρ is the density, u is the velocity, $\mathit{g}=-g{\widehat{\mathit{e}}}_z$ is the acceleration due to gravity with g >  0, p is the pressure, T is the temperature, s is the specific entropy, and ν is the constant kinematic viscosity. Furthermore, F_rad and F_SGS are the radiative and turbulent subgrid scale (SGS) fluxes, respectively, and Γ_cool describes cooling at the surface (see below). [BOLD]S is the traceless rate-of-strain tensor with

$$\begin{array}{c}\hfill {\mathsf{S}}_{\mathit{ij}}=\frac{1}{2}(u_{i,j}+u_{j,i})-\frac{1}{3}{\delta }_{\mathit{ij}}\mathbf{\nabla }·\mathit{u}.\end{array}$$(5)

An optically thick fully ionised gas is considered, where radiation is modelled through diffusion approximation. The ideal gas equation of state p = (c_P − c_V)ρT = ℛρT applies, where ℛ is the gas constant, and c_V is the specific heat at constant volume. The radiative flux is given by

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathrm{rad}}=-K\mathbf{\nabla }T,\end{array}$$(6)

where K is the radiative heat conductivity. Two qualitatively different heat conductivity prescriptions are considered, where K either has a fixed profile K(z) or is a function of density and temperature, K(ρ, T). In the latter case, K is given by

$$\begin{array}{c}\hfill K=\frac{16{\sigma }_{\mathrm{SB}}{T}^3}{3\kappa \rho },\end{array}$$(7)

where σ_SB is the Stefan-Boltzmann constant and κ is the opacity. κ is assumed to obey a power law

$$\begin{array}{c}\hfill \kappa ={\kappa }_0{(\rho /{\rho }_0)}^a{(T/T_0)}^b,\end{array}$$(8)

where ρ₀ and T₀ are reference values of density and temperature. Equations (7) and (8) combine into

$$\begin{array}{c}\hfill K(\rho ,T)=K_0{(\rho /{\rho }_0)}^{-(a+1)}{(T/T_0)}^{3-b}.\end{array}$$(9)

Here a = 1 and b = −7/2 are used, which correspond to the Kramers opacity law for free-free and bound-free transitions (Weiss et al. 2004). Heat conductivity consistent with the Kramers law was first used in convection simulations by Brandenburg et al. (2000).

Owing to the strong depth dependence of the radiative diffusivity, χ = K/(c_Pρ), additional turbulent SGS diffusivity is used in the entropy equation to keep the simulations numerically feasible. Here the SGS flux is formulated as

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathrm{SGS}}=-\rho T{\chi }_{\mathrm{SGS}}^{(1)}\mathbf{\nabla }s\prime ,\end{array}$$(10)

where $s\prime =s-\overline{s}$ is the fluctuation of the specific entropy. The overbar indicates horizontal averaging here and in what follows. The coefficient ${\chi }_{\mathrm{SGS}}^{(1)}$ is constant in the whole domain. F_SGS has a negligible contribution to the net horizontally averaged energy flux, such that ${\overline{\mathit{F}}}_{\mathrm{SGS}}\approx 0$. The current SGS formulation is similar to those used in Käpylä et al. (2007, 2019a), and Brown et al. (2010), for example.

The cooling at the surface is described by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{cool}}=-{\mathrm{\Gamma }}_{0}f(z)(T_{\mathrm{cool}}-T),\end{array}$$(11)

where Γ₀ is a cooling luminosity, T = e/c_V is the temperature where e is the internal energy, and where T_cool = T_top is a reference temperature corresponding to the fixed value at the top boundary.

2.1. Geometry, and initial and boundary conditions

The computational domain is a rectangular box where z_bot ≤ z ≤ z_top is the vertical coordinate, where z_bot/d = −0.45, z_top/d = 1.05, and where d is the depth of the initially isentropic layer (see below). In a few runs the domain extends to deeper layers such that z_bot/d = −0.75 to accommodate deeper overshooting. The horizontal coordinates x and y run from −2d to 2d. The horizontal size of the box is thus L_H/d = 4, and the vertical extent L_z/d is either 1.5 or 1.8.

The initial stratification consists of three layers. Two configurations of the three-layer set-up are considered here: in the first set-up (hereafter P²I), the two lower layers are polytropic with polytropic indices n₁ = 3.25 (z_bot/d ≤ z/d ≤ 0) and n₂ = 1.5 (0 ≤ z/d ≤ 1). The uppermost layer above z/d = 1 is initially isothermal. This layer mimics a photosphere where radiative cooling is efficient. The choice of n₁ is motivated by fact that in the special case where the temperature gradient in the corresponding hydrostatic state is constant, the solution is a polytrope with index 13/4; see Barekat & Brandenburg (2014) and Appendix A of Brandenburg (2016). Assuming that an extended stable layer forms at the bottom of the domain, its stratification is close to the hydrostatic solution (see, e.g. Käpylä et al. 2019a). This, however, can only be confirmed a posteriori because the depth of the convective layer is not pre-determined in the cases where the Kramers opacity law is used. In the second set-up (hereafter P³) all three layers are polytropic with indices (n₁, n₂, n₃)=(2, 1.5, 1.5). In these cases the radiative diffusion in the uppermost layer is enhanced and no explicit cooling is applied. This configuration is the same as in Singh et al. (1998) and was chosen to accommodate comparisons with that study. The choice of n₂ = 1.5 for the thermal stratification in the middle layer comes from the expectation that the convectively unstable layer is nearly isentropic in the final statistically saturated state. The initial velocity follows a Gaussian-noise distribution with an amplitude of about ${10}^{-4}\sqrt{\mathit{dg}}$.

The horizontal boundaries are periodic, and on the vertical boundaries, impenetrable and stress-free boundary conditions are imposed for the flow such that

$$\begin{array}{c}\hfill \frac{\partial u_x}{\partial z}=\frac{\partial u_y}{\partial z}=u_z=0.\end{array}$$(12)

The temperature gradient at the bottom boundary is set according to

$$\begin{array}{c}\hfill \frac{\partial T}{\partial z}=-\frac{F_{\mathrm{bot}}}{K_{\mathrm{bot}}},\end{array}$$(13)

where F_bot is the fixed input flux and K_bot(x, y, z_bot) is the value of the heat conductivity at the bottom of the domain. On the upper boundary a constant temperature T = T_top, coinciding with the initial value, is assumed.

2.2. Units, control parameters, and simulation strategy

The units of length, time, density, and entropy are given by

$$\begin{array}{c}\hfill [x]=d,[t]=\sqrt{d/g},[\rho ]={\rho }_0,[s]=c_{\mathrm{P}},\end{array}$$(14)

where ρ₀ is the initial value of the density at z = z_top. The models with Kramers heat conductivity are fully defined by choosing the value of the kinematic viscosity ν, the gravitational acceleration g, the values of a, b, K₀, ρ₀, T₀ and the SGS Prandtl number

$$\begin{array}{c}\hfill {\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}=\frac{\nu }{{\chi }_{\mathrm{SGS}}^{(1)}},\end{array}$$(15)

along with z-dependent cooling profile f(z). The values of K₀, ρ₀, and T₀ are subsumed into a new variable ${\stackrel{\sim }{K}}_0=K_0{\rho }_0^{a+1}{T}_0^{b-3}$ that is fixed by assuming the radiative flux at z_bot to equal F_bot in the initial state. The profile f(z) = 1 above z/d = 1 and f(z) = 0 below z/d = 1, connecting smoothly over the interface over a width of 0.025d. Furthermore, ${\xi }_0=H_{\mathrm{p}}^{\mathrm{top}}/d=\mathcal{R}{T}_{\mathrm{top}}/gd$ sets the initial pressure scale height at the surface, thus determining the initial density stratification. All of the current simulations have ξ₀ = 0.054.

In runs where a fixed profile of heat conductivity is used, the profile K(z) is needed instead of specifying ${\stackrel{\sim }{K}}_0$. In these cases the value of K at z = z_bot is fixed similarly as was done for ${\stackrel{\sim }{K}}_0$ in the Kramers cases. In these cases the initial profile of the Prandtl number based on the radiative heat conductivity

$$\begin{array}{c}\hfill \mathrm{Pr}(z)=\frac{\nu }{\chi (z)},\end{array}$$(16)

where χ(z)=K(z)/c_Pρ(z), sets the relative importance of viscous to temperature diffusion. We note that in general Pr = Pr(x, t) because ρ = ρ(x, t). In cases where K has a piecewise constant profile, it can be represented in terms of polytropic indices $(n_1^{\prime },n_2^{\prime },n_3^{\prime })$, which refer to a corresponding non-convective hydrostatic solution. Starting the simulations from such solutions for the thermodynamic quantities is, however, impractical especially if the value n₂ is far away from the final convective state, which is always close to the adiabatic value of 1.5. Thus these indices typically differ from those used to initialise the thermal variables (see also Brandenburg et al. 2005). The heat conductivities in the different layers are connected through

$$\begin{array}{c}\hfill \frac{K_i}{K_j}=\frac{n_i^{\prime }+1}{n_j^{\prime }+1},\end{array}$$(17)

where i and j refer to any of the three layers.

The dimensionless normalised flux is given by

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{n}}=F_{\mathrm{bot}}/{\rho }_{\mathrm{bot}}{c}_{\mathrm{s},\mathrm{bot}}^3,\end{array}$$(18)

where ρ_bot and c_s, bot are the density and the sound speed, respectively, at z_bot in the initial non-convecting state.

The input energy flux determines the overall convective velocity realised in the simulations via $u_{\mathrm{conv}}\propto {\mathcal{F}}_{\mathrm{n}}^{1/3}$, see Eq. (1) and Käpylä et al. (2019b). Thus if only ℱ_n were changed, the relative importance of the diffusion coefficients, measured by Reynolds and Péclet numbers, would change as well. To eliminate these dependences, and to be able to concentrate solely on the effects of varying input flux, the viscosity and SGS entropy diffusion are scaled proportional to ${\mathcal{F}}_{\mathrm{n}}^{1/3}$. In addition to changing the diffusion coefficients, the cooling luminosity Γ₀ at the surface is also scaled proportionally to ℱ_n. The input flux itself is varied by changing the overall magnitude of K such that the value at the bottom is given by K_bot = −F_bot/(∂T/∂z)_zbot. Furthermore, the degree of supercriticality of convection, measured by a Rayleigh number, and the Prandtl number, describing the ratio of viscosity to thermal diffusion, will affect the properties of convection and overshooting if they are allowed to vary. The current choice of the dependency of ν and ${\chi }_{\mathrm{SGS}}^{(1)}$ on ℱ_n, however, eliminates these dependences such that the supercriticality and effective Prandtl numbers are nearly constant over the range of ℱ_n considered here. This is discussed in detail in Sect. 4.1.

In addition to the explicit viscosity, SGS, and radiative diffusion, the advective terms in each of the Eqs. (2)–(4) are written in terms of a fifth-order upwinding derivative with a hyperdiffusive sixth-order correction where the diffusion coefficient depends locally on the flow, see Appendix B of Dobler et al. (2006).

2.3. Diagnostics quantities

The following quantities are outcomes of the simulations that can only be determined a posteriori. These include the global Reynolds and SGS Péclet numbers,

$$\begin{array}{c}\hfill \mathrm{Re}=\frac{u_{\mathrm{rms}}}{\nu k_1},{\mathrm{Pe}}_{\mathrm{SGS}}=\frac{u_{\mathrm{rms}}}{{\chi }_{\mathrm{SGS}}^{(1)}{k}_1},\end{array}$$(19)

where u_rms is the volume-averaged rms velocity, and k₁ = 2π/d is an estimate of the largest eddies in the system. Typically $\chi \ll {\chi }_{\mathrm{SGS}}^{(1)}$ in the CZ in most of the current simulations. Furthermore, χ has a strong depth dependence due to the density stratification. Thus it is useful to define Reynolds and effective Péclet numbers separately for the overshoot zone,

$$\begin{array}{c}\hfill {\mathrm{Re}}_{\mathrm{OZ}}=\frac{u_{\mathrm{rms}}^{\mathrm{OZ}}}{\nu k_{\mathrm{OZ}}},{\mathrm{Pe}}_{\mathrm{OZ}}^{\mathrm{eff}}=\frac{u_{\mathrm{rms}}^{\mathrm{OZ}}}{({\chi }_{\mathrm{SGS}}^{(1)}+{\overline{\chi }}_{\mathrm{OZ}})k_{\mathrm{OZ}}},\end{array}$$(20)

where all quantities are taken from the base of the CZ, z = z_CZ, with ${\overline{\chi }}_{\mathrm{OZ}}={\overline{K}}_{\mathrm{OZ}}/(c_{\mathrm{P}}{\overline{\rho }}_{\mathrm{OZ}})$ being the mean (horizontally averaged) radiative diffusivity, and where k_OZ = 2π/d_os is a wavenumber based on the depth of the overshoot layer d_os. Similarly, an effective Prandtl number can be defined as

$$\begin{array}{c}\hfill {\mathrm{Pr}}_{\mathrm{OZ}}^{\mathrm{eff}}=\frac{\nu }{({\chi }_{\mathrm{SGS}}^{(1)}+{\overline{\chi }}_{\mathrm{OZ}})}·\end{array}$$(21)

Precise definitions of z_CZ and d_os are given in Sect. 3.

To assess the level of supercriticality of convection, radiative and SGS Rayleigh numbers are defined as

$$\begin{array}{cc}\hfill {\mathrm{Ra}}^{\mathrm{Rad}}& =\frac{g{d}^4}{\nu \chi }(-\frac{1}{c_{\mathrm{P}}}\frac{\mathrm{d}s}{\mathrm{d}z}),\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill {\mathrm{Ra}}^{\mathrm{SGS}}& =\frac{g{d}^4}{\nu {\chi }_{\mathrm{SGS}}^{(1)}}(-\frac{1}{c_{\mathrm{P}}}\frac{\mathrm{d}s}{\mathrm{d}z}).\hfill \end{array}$$(23)

Supercriticality of convection is roughly determined by min(Ra^Rad, Ra^SGS). Both quantities vary as functions of height and are quoted near the surface at z/d = 0.85 for all models. Conventionally, the Rayleigh number in the hydrostatic non-convecting state is one of the control parameters. This is still true for the cases with fixed K profile in the current study, but in the runs with Kramers conductivity, the convectively unstable layer in the hydrostatic case is very thin and confined to the near-surface layers (Brandenburg 2016). Thus the Rayleigh numbers are quoted from the thermally saturated and statistically stationary states.

The path in parameter space taken here is artificial in that in real fluids, there is no SGS diffusivity (Ra^SGS = 0) and any change in the flux would be reflected by the then decisive radiative Rayleigh number. For example, a system where Pr is fixed ${\mathrm{Ra}}^{\mathrm{rad}}\propto {\mathcal{F}}_{\mathrm{n}}^{-4/3}$. However, taking this path severely limits the computationally feasible range of ℱ_n because the Reynolds and Péclet number also increase in proportion to ℱ_n and the resolution requirements quickly become prohibitive.

Contributions to the vertical energy flux are

$$\begin{array}{cc}& {\overline{F}}_{\mathrm{rad}}=-\overline{K}\frac{\partial \overline{T}}{\partial z},\hfill \end{array}$$(24)

$$\begin{array}{cc}& {\overline{F}}_{\mathrm{enth}}=c_{\mathrm{P}}\overline{(\rho u_z)\prime T\prime },\hfill \end{array}$$(25)

$$\begin{array}{cc}& {\overline{F}}_{\mathrm{kin}}=\frac{1}{2}\overline{\rho {\mathit{u}}^2{u}_z\prime },\hfill \end{array}$$(26)

$$\begin{array}{cc}& {\overline{F}}_{\mathrm{visc}}=-2\nu \overline{\rho u_i{\mathsf{S}}_{\mathit{iz}}}\hfill \end{array}$$(27)

$$\begin{array}{cc}& {\overline{F}}_{\mathrm{cool}}={\displaystyle {\int }_{z_{\mathrm{bot}}}^{z_{\mathrm{top}}}{\mathrm{\Gamma }}_{\mathrm{cool}}\mathrm{d}z.}\hfill \end{array}$$(28)

Here the primes denote fluctuations and overbars horizontal averages. The total convected flux (Cattaneo et al. 1991) is the sum of the enthalpy and kinetic energy fluxes:

$$\begin{array}{c}\hfill {\overline{F}}_{\mathrm{conv}}={\overline{F}}_{\mathrm{enth}}+{\overline{F}}_{\mathrm{kin}}.\end{array}$$(29)

The radiative flux can be written in terms of the mean double-logarithmic temperature gradient $\mathrm{\nabla }=\partial ln\overline{T}/\partial ln\overline{p}$ as

$$\begin{array}{c}\hfill {\overline{F}}_{\mathrm{rad}}=-\overline{K}\frac{\partial \overline{T}}{\partial z}=\frac{\overline{K}g}{c_{\mathrm{P}}}\frac{\mathrm{\nabla }}{{\mathrm{\nabla }}_{\mathrm{ad}}},\end{array}$$(30)

where c_P∇_ad = c_P(1 − 1/γ)=ℛ, and g = |g|. Above, ∇ is the actual temperature gradient realised in the system. Now it is possible to define the total flux and the flux transported by adiabatic stratification as

$$\begin{array}{c}\hfill F_{\mathrm{tot}}=\frac{\overline{K}g}{c_{\mathrm{P}}}\frac{{\mathrm{\nabla }}_{\mathrm{rad}}}{{\mathrm{\nabla }}_{\mathrm{ad}}},{\overline{F}}_{\mathrm{rad}}^{\mathrm{ad}}=\frac{\overline{K}g}{c_{\mathrm{P}}},\end{array}$$(31)

where ∇_rad is a hypothetical radiative gradient in the absence of convection and F_tot = F_bot. The ratio of the total-to-adiabatic flux is the Nusselt number (e.g. Brandenburg 2016),

$$\begin{array}{c}\hfill \mathrm{Nu}=\frac{F_{\mathrm{tot}}}{{\overline{F}}_{\mathrm{rad}}^{\mathrm{ad}}}=\frac{{\mathrm{\nabla }}_{\mathrm{rad}}}{{\mathrm{\nabla }}_{\mathrm{ad}}}·\end{array}$$(32)

In the current set-up the Nusselt number is fixed from the outset for cases with a static profile of $\overline{K}$. This is because the total flux imposed at the lower boundary is proportional to $\overline{K}$. In the cases where Kramers conductivity is used, $\overline{K}$ can evolve but the change in the Nusselt number is not very large; see, for example, Table 1. The behaviour of Nu is different in cases where fixed temperature is imposed at both boundaries (e.g. Yadav et al. 2016; Gastine et al. 2016) because the total flux is no longer in direct proportion with the heat conductivity.

3. Definitions of CZ and overshooting

To characterise the different layers, the nomenclature introduced in Käpylä et al. (2017) is used, although with somewhat differing definitions. The CZ is defined to be the part of the domain where ${\overline{F}}_{\mathrm{conv}}>0$, whereas in the overshoot zone (OZ) F_conv <  0. This is motivated by the work of Deng & Xiong (2008), who used a similar definition but employed ${\overline{F}}_{\mathrm{enth}}$ instead of ${\overline{F}}_{\mathrm{conv}}$. A definition of the overshooting depth based on ${\overline{F}}_{\mathrm{enth}}$ was also used in Käpylä et al. (2017). However, because the kinetic energy flux carries a substantial fraction of the energy, it is natural to include it in the definition of the CZ. The bottom of the CZ is denoted by z_CZ.

The mean overshooting depth z_OS is taken to be the position where the horizontally averaged ${\overline{F}}_{\mathrm{kin}}$ drops below 1% of its value at z_CZ. Various earlier studies have used a similar definition (Hurlburt et al. 1986, 1994; Singh et al. 1995, 1998; Brummell et al. 2002). In most of the previous studies the location of the bottom of the CZ is assumed to be fixed by the initial conditions. Here this assumption is relaxed and z_CZ is computed from the simulation data using the definition given above. The values of z_OS and z_CZ are obtained by linear interpolation from the grid points closest to the respective transitions. Furthermore, z_CZ and z_OS are functions of time. The overshooting depth is defined as

$$\begin{array}{c}\hfill d_{\mathrm{os}}={⟨z_{\mathrm{CZ}}(t)-z_{\mathrm{OS}}(t)⟩}_t,\end{array}$$(33)

where ⟨⋯⟩_t denotes a time average over the statistically stationary part of the time series. Error estimates for d_os are obtained by dividing the time series into three equally long parts and considering their largest deviation from the time average over the whole time series as the error.

The radiative zone (RZ) is defined as the region below z_OS, and the buoyancy zone (BZ) is where ${\overline{F}}_{\mathrm{conv}}>0$ and ${\partial }_z\overline{s}<0$. Finally, the Deardorff zone (DZ) is characterised by a formally stable stratification with a positive vertical mean entropy gradient (${\partial }_z\overline{s}>0$) and ${\overline{F}}_{\mathrm{conv}}>0$; see, Tremblay et al. (2015). In this layer, the convective energy transport is dominated by a non-local non-gradient contribution to the enthalpy flux introduced by Deardorff (1961, 1966); see also Brandenburg (2016) and Käpylä et al. (2017). Such layers have been reported by various authors from simulations (e.g. Chan & Gigas 1992; Roxburgh & Simmons 1993; Tremblay et al. 2015; Korre et al. 2017; Bekki et al. 2017; Karak et al. 2018; Nelson et al. 2018). Furthermore, the union of OZ, DZ, and BZ is referred to as the mixed zone (MZ).

4. Results

Four main sets of simulations, denoted as K, P, S, and DS were conducted, see Tables 1–3. The sets are named after their heat conductivity prescriptions: in Set K, the Kramers law is used to compute the heat conductivity. In Set P, a static profile corresponding to the heat conductivity computed according to the Kramers law in the initial state of the simulation is used. Set S employs a static step profile of heat conductivity, K = K(z). These three sets correspond to Runs K, P, and S of Käpylä et al. (2017). Additionally, Sets Kh and Sh are the otherwise the same as Sets K and S, but lower viscosity and SGS entropy diffusion were used. These runs were branched off from corresponding thermally saturated snapshots from runs in Sets K and S, respectively. In Sets DS, DSS1, and DSS2 a double-step profile for K, similar to that in Singh et al. (1998), was used. In Set DS ${\chi }_{\mathrm{SGS}}^{(1)}=0,$ and in Set DSS1 (DSS2) the SGS diffusion is included with ${\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}=5(10)$. Runs DS5 and DS5h were the same except for the grid resolution to study numerical convergence of the results. Finally, in Set R the explicit diffusivities ν and ${\chi }_{\mathrm{SGS}}^{(1)}$ were varied while all other parameters were kept fixed using Run K5 as a progenitor, see Table 4. Sets DS, DSS1, and DSS2 used the P³ set-up, and the remaining sets were initialised with the P²I set-up.

4.1. Basic characterisation of the solutions

Typical flow patterns for Runs R4 and R7 are shown in Fig. 1. The flow structure observed in various studies of stratified non-rotating convection is recovered with connected downflows near the surface merging into isolated plumes at larger depths (e.g. Stein & Nordlund 1989, 1998). The downflows are surrounded by broader upflows. In the majority of the current simulations the flows are at best mildly turbulent with Re ≈ 20…40 such as in Run R4 in the left panel of Fig. 1. However, the qualitative large-scale structure of convection does not change at higher resolutions and Rayleigh and Reynolds numbers, see the right panel of Fig. 1 for Run R7.

Fig. 1. Normalised vertical velocity ${\stackrel{\sim }{u}}_z=u_z/{(gd)}^{1/2}$ (colours) and streamlines of the flow from Runs R4 (left) and R7 (right).

Figure 2a shows the profiles of $\overline{K}(z)$ from representative runs in each of the main sets K, P, S, and DS. The Kramers run K4 and the corresponding fixed profile model P4 have a smoothly varying profile as a function of height with very small values of $\overline{K}$ near the surface. In the piecewise polytropic run S4, the profile of $\overline{K}$ is characterised by constant values in the upper (z >  0) and lower (z <  0) layers (e.g. Hurlburt et al. 1986; Nordlund et al. 1992), which can be characterised by the ratios K₂/K₁ = K₃/K₁ = 4/17 ≈ 0.235. In Run DS5, K₂/K₁ = 2/85 ≈ 0.0333 and the uppermost layer above z/d = 1 has another constant value of $\overline{K}$ corresponding to K₂/K₁ = 5/6. However, the transitions of $\overline{K}$ between the layers are smoothed over a distance of 0.05d, due to which the value of $\overline{K}$ corresponding to $n_3^{\prime }=1.5$ is achieved only near the upper surface.

Fig. 2. Horizontally averaged and normalised heat conductivity $\stackrel{\sim }{K}=\overline{K}(z)/K_{\mathrm{bot}}$ (a) and Prandtl number Pr(z) (b) from Runs K4, P4, S4, and DS5. The dashed black (red) line in the lower panel indicates ${\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}=1$ (${\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}=5$).

The vertical profiles of Pr(z) and ${\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}$ are shown in Fig. 2b. The strong variation of Pr as function of depth is to be contrasted with the constant value of ${\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}$. Although Pr ∝ ℱ_n, it is still typically greater than unity at z/d = 0 in almost all cases, with the exception of a few runs with the highest values of ℱ_n (see Tables 1–4). This is in stark contrast to the solar CZ, where Pr ≪ 1 everywhere (e.g. Ossendrijver 2003). Numerical simulations with Prandtl numbers far different from unity are challenging numerically and render parameter studies infeasible. Thus an enhanced SGS diffusivity with ${\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}$ of the order of unity was applied in most of the current simulations. In most of the current simulations the effective Prandtl number is dominated by the SGS diffusion.

Changing the input energy flux ℱ_n refers to varying the magnitude of $\overline{K}$ proportionally and thus χ ∝ ℱ_n. In combination with $\nu \propto {\mathcal{F}}_{\mathrm{n}}^{1/3}$ and the ML estimate −(H_p/c_P)ds/dz = $\mathrm{\nabla }-{\mathrm{\nabla }}_{\mathrm{ad}}\propto {\mathcal{F}}_{\mathrm{n}}^{2/3}$ (Vitense 1953), this involves changing ${\mathrm{Ra}}^{\mathrm{Rad}}\propto {\mathcal{F}}_{\mathrm{n}}^{-2/3}$. This is to be contrasted with the SGS Rayleigh number, where χ is replaced by ${\chi }_{\mathrm{SGS}}^{(1)}\propto {\mathcal{F}}_{\mathrm{n}}^{1/3}$, which leads to Ra^SGS being independent of ℱ_n. These scalings are in accordance with the simulations, as can be seen from Fig. 3 for Sets Kh, DS, and DSS1. As mentioned earlier, the supercriticality of convection is determined roughly by min(Ra^Rad, Ra^SGS). In the current simulations the SGS Rayleigh number is almost always smaller. This also means that with Ra^SGS≈ constant, the supercriticality of convection is also constant, and it is eliminated as an influence on the overshooting depth. For completeness, the flux-based Rayleigh number (e.g. Brun & Browning 2017)

$$\begin{array}{c}\hfill {\mathrm{Ra}}^{\mathrm{Flux}}=\mathrm{Nu}{\mathrm{Ra}}^{\mathrm{Rad}}=\frac{F_{\mathrm{tot}}{d}^4}{\rho {\chi }^2\nu }(-\frac{1}{c_{\mathrm{P}}}\frac{\mathrm{d}s}{\mathrm{d}z}),\end{array}$$(34)

Fig. 3. Rayleigh numbers RaSGS (black) and RaRad (red) as functions of ℱn from Set Kh. The purple (magenta) line shows RaRad (RaSGS) from Set DS (DSS). The dotted lines show scalings expected from ML arguments.

can be seen to vary as ${\mathrm{Ra}}^{\mathrm{Flux}}\propto {\mathcal{F}}_{\mathrm{n}}^{-2/3}$ given the dependences above. We note that NuRa^SGS is independent of ℱ_n.

Figure 4 shows the dependence of the rms velocity within the CZ from all simulations except for those in Set R. The simulation results are close to a ${\mathcal{F}}_{\mathrm{n}}^{1/3}$ dependence, with the coefficient of proportionality varying between 1.3 and 1.8. This is in agreement with ML estimates (e.g. Vitense 1953; Brandenburg et al. 2005; Brandenburg 2016) and earlier numerical findings (Brandenburg et al. 2005; Karak et al. 2015; Käpylä et al. 2019b).

Fig. 4. Normalised rms velocity ${\stackrel{\sim }{u}}_{\mathrm{rms}}=u_{\mathrm{rms}}/\sqrt{\mathit{gd}}$ in the CZ as a function of ℱn for the simulation sets indicated by the legend. The dotted line is proportional to ${\mathcal{F}}_{\mathrm{n}}^{1/3}$. The inset shows ${\stackrel{\sim }{u}}_{\mathrm{rms}}{\mathcal{F}}_{\mathrm{n}}^{-1/3}$ for the same runs. The purple diamond denotes Run DS5h.

Results regarding the energy fluxes and force balance from representative runs from Sets K, S, and DS are shown in Fig. 5. The left panels show the contributions to the energy flux and the superadiabatic temperature gradient ∇ − ∇_ad. The fluxes for Runs K4 and S4 in Fig. 5a and c are qualitatively similar to those of Runs K and S of Käpylä et al. (2017), respectively. The main difference to the latter is the treatment of the near surface layers; cooling layer in the present runs as opposed to an imposed entropy gradient in Käpylä et al. (2017), and the somewhat different values of ℱ_n; 1.8 × 10⁻⁵ here versus 9 × 10⁻⁶ in Käpylä et al. (2017). The subadiabatic Deardorff zone encompasses roughly a quarter of the MZ in Run K4, whereas in Run S4, the DZ is almost absent. The runs in Set DS differ from those in Set S in that convection transports almost all of the energy due to the lower K₂/K₁ ratio. Assuming that the temperature gradient in the final statistically saturated state is nearly adiabatic, the fraction of convective transport can be estimated from (cf. Brandenburg et al. 2005)

$$\begin{array}{c}\hfill \frac{F_{\mathrm{conv}}}{F_{\mathrm{tot}}}\approx 1-\frac{{\mathrm{\nabla }}_{\mathrm{ad}}}{{\mathrm{\nabla }}_{\mathrm{rad}}}=1-{\mathrm{Nu}}^{-1}=1-{\mathrm{\nabla }}_{\mathrm{ad}}(n_2\prime +1).\end{array}$$(35)

Fig. 5. Panels a, c, and e: total (black dash-dotted lines), convective (black), enthalpy (blue), kinetic energy (light purple), radiative (dark purple), and cooling (green) fluxes and the superadiabatic temperature gradient (red) from Runs K4, S4, and DS5, respectively. Panels b, d, and f: total averaged vertical forces (solid lines) and the power of the forces (dashed) on the upflows (red) and downflows (blue) from the same runs. The dotted lines in these panels show the corresponding viscous force and its power. The shaded areas indicate the BZ (darkest), DZ, and OZ (lightest) and the thick black line at the horizontal axis denotes the MZ.

According to this expression, the convective flux transports 60% (96%) in the runs in Set S (DS). This is confirmed by the numerical results.

Figures 5b, d, and f show the horizontally averaged vertical total and viscous forces, ${\overline{f}}_z=\overline{\rho \mathrm{D}{u}_z/\mathrm{D}t}$, and ${\overline{f}}_{\mathrm{visc}}=2\nu {\partial }_i(\overline{\rho S_{\mathit{iz}}})$, respectively, and the resulting power of the forces (${\overline{P}}_z=\overline{u_z{f}_z}$ and ${\overline{P}}_{\mathrm{visc}}=\overline{u_z{f}_{\mathrm{visc}}}$) separately for the upflows (↑) and the downflows (↓) from Runs K4, S4, and DS5. The force balance in Run K4 (Fig. 5b), is very similar to the corresponding Run K in Käpylä et al. (2017), see their Fig. 2b. The downflows appear to adhere to the Schwarzschild criterion such that they are accelerated in the BZ and decelerated in the layers below. This is contrasted by the upflows that are accelerated in the stably stratified OZ and DZ and in the lower part of the BZ. As demonstrated in Käpylä et al. (2017), the upflows are not driven by the convective instability but are a result of matter displaced by the deeply penetrating downflows. In Run S4 the downflows are decelerated already in the lower part of the BZ whereas the force on the upflows is qualitatively similar to that in Run K4. The force balance in Run DS5 is qualitatively similar to that in Run S4 although the magnitude and details of the quantities differ. Interestingly the sign of the total force in the stably stratified layer near the surface is not reversed. The shallowness of the layer is likely contributing to this. Another aspect is the (true) overshooting from the CZ below. It can, however, be concluded that displacement of the matter due to the downflows is driving the upflows in the OZ and DZ also in Runs S4 and DS5. However, the superadiabatic temperature gradient has a local maximum at the bottom of the BZ in these cases, see Figs. 5c and e. Thus it is possible that the convective instability is contributing more to the upward acceleration in these cases. The viscous force is small in all cases and has a noticeable effect only in the near-surface layers above z/d = 0.75.

4.2. Dependence of overshooting on input flux

Reaching the solar value of ℱ_n is currently not possible due to the prohibitive time-step constraint and the long thermal adjustment time involved (e.g. Brandenburg et al. 2000; Kupka & Muthsam 2017). However, it is reasonable to assume that the overshooting depth scales with a power law as a function of ℱ_n (e.g. Schmitt et al. 1984; Zahn 1991). Thus it is in principle possible to estimate the extent of overshooting in the Sun provided that a sufficiently broad range of higher flux values are probed and their results are extrapolated to the solar case. A few such studies can be found in the literature (e.g. Singh et al. 1998; Tian et al. 2009; Hotta 2017).

One of the most restrictive modelling choices in the past has been the use of a static heat conduction profile that effectively enforces the layer structure of the simulation. This can be seen from Fig. 6a where the vertical coordinates of the bottoms of the convection (z_CZ) and overshoot (z_OZ) zones are shown as functions of ℱ_n. The results for z_CZ from Sets S and Sh show that the interface between the CZ and OZ stays at the initial position at z = 0 for all values of ℱ_n. The runs in Sets DS, DSS1, and DSS2 behave similarly, although the bottom of the CZ is shifted downward from its initial position. For Sets K and Kh the depth of the CZ is generally reduced in comparison to Sets S and Sh. In these runs the depth of the CZ increases as ℱ_n decreases. This is contrasted by the results from Sets P where a Kramers-like, but static, profile of K is used (see Fig. 2): here z_CZ is practically fixed in the current range of ℱ_n.

Fig. 6. Panel a: vertical (z) coordinates of the bottom of the CZ (zCZ, solid lines) and OZ (zOZ, dashed). The dotted red line indicates the bottom of the domain. Panel b: overshooting depth dos normalised by the pressure scale height Hp as a function of ℱn for Sets K (black), Kh (grey), P (blue), S (red), Sh (green), DS (purple), DSS1 (cyan), and DSS2 (orange). The purple diamond denotes Run DS5h. The dotted lines show approximate power laws from fits to simulation data; see Table 5. Panel c: comparison of Sets DS, DSS1, and DSS2 with the studies of Singh et al. (1998) (red) and Tian et al. (2009) (blue).

Figure 6a shows that the location of the base of the OZ (z_OZ) increases monotonically as a function of ℱ_n in all sets except for K and Kh. In Set DSS1 (DSS2) the overshooting extends to the lower boundary of the domain in the three (two) highest ℱ_n runs, rendering the results of these simulations unusable in the following analysis. The depth of the overshoot zone, d_os, as a function of ℱ_n is shown in Fig. 6b. The results for Sets K, Kh, and P fall almost on top of each other. The data suggests two power laws: ${\mathcal{F}}_{\mathrm{n}}^{0.08}$ for ℱ_n ≲ 10⁻⁵ and ${\mathcal{F}}_{\mathrm{n}}^{0.19}$ for ℱ_n ≳ 10⁻⁵; see Table 5. These results suggest that the overshooting depth is relatively insensitive to the choice of the heat conduction scheme if the profile of K at the base of the CZ is smooth. Furthermore, d_os is consistently greater in Set Kh with higher Reynolds and SGS Péclet numbers than in the corresponding runs in Set K. This is because the current simulations have relatively modest Reynolds and Péclet numbers that are not in a fully turbulent regime. This aspect is studied in more detail in Sect. 4.3. At first glance, the data of Set Sh appear to be more consistent with a single power, and in Set S there appears to be a break around ℱ_n ≈ 10⁻⁵, such that the data points for lower values of ℱ_n lie below those of Set Sh. However, power-law fits for both full and partial ranges are consistent with a ${\mathcal{F}}_{\mathrm{n}}^{0.12}$ scaling within the error estimates; see Table 5.

A similar set-up as in Singh et al. (1998) and Tian et al. (2009) is adopted in Sets DS, DSS1, and DSS2. The results of Set DS show a steep dependence of the overshooting on the input flux ($d_{\mathrm{os}}\propto {\mathcal{F}}_{\mathrm{n}}^{0.27}$), and Sets DSS1 with $d_{\mathrm{os}}\propto {\mathcal{F}}_{\mathrm{n}}^{0.07}$ and DSS2 with $d_{\mathrm{os}}\propto {\mathcal{F}}_{\mathrm{n}}^{0.12}$ are more in line with the other sets of simulations. The overshooting depth for the higher resolution Run DS5h is statistically consistent with that of Run DS5, although the overall velocities in these runs are somewhat lower. This is in particular reflected by the values of Re_OZ and ${\mathrm{Pe}}_{\mathrm{OZ}}^{\mathrm{eff}}$ (seventh and eighth columns in Table 3). The fact that the results do not change drastically with resolution and that there is an apparently continuous transition as ${\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}$ increases through Sets DSS1 and DSS2 to DS are indicative that the latter set of runs is sufficiently well resolved numerically. Figure 6c shows a comparison of the overshooting depths from Sets DS, DSS1, and DSS2 with the studies of Singh et al. (1998, their Table 1) and Tian et al. (2009, their Table 4). Both studies list the input fluxes F_b and values of density and pressure at the bottom of the domain (denoted here as ρ_b and p_b). Thus the normalised flux in both cases is computed from ℱ_n = F_b/F_n, where $F_{\text{n}}={\rho }_{\text{b}}{c}_{\text{s}}^3={\rho }_{\text{b}}^{-1/2}{(\gamma p_{\text{b}})}^{3/2}$ with $c_{\text{s}}^2=\gamma p_{\text{b}}/{\rho }_{\text{b}}$ and γ = 5/3 has been assumed. The results of Singh et al. (1998) were obtained by a similar definition as used in the present study based on kinetic energy flux falling below a threshold fraction of its value at the base of the CZ. Their results are consistent with a $d_{\mathrm{os}}\propto {\mathcal{F}}_{\mathrm{n}}^{0.33}$ scaling. On the other hand, Tian et al. (2009) obtained a very shallow dependence with this definition ($d_{\mathrm{os}}\propto {\mathcal{F}}_{\mathrm{n}}^{0.05}$) and a much steeper one ($d_{\mathrm{os}}\propto {\mathcal{F}}_{\mathrm{n}}^{0.42}$) when using a criterion based on the enthalpy flux falling below a threshold value of its absolute maximum in the OZ. It is unclear why the results from the two definitions used by Tian et al. (2009) deviate. Both definitions were tested with the current data and the results were found to be in fair agreement.

The drastic change from a steep power law in Set DS to the shallower power laws in Sets DSS1 and DSS2 is related to the difference in the diffusion of temperature fluctuations. Figure 7a shows that the effective Péclet numbers in OZ in all of these sets are comparable. However, as is shown in Sect. 4.3, the overshooting depth is insensitive to the Péclet number in this parameter range provided that the effective Prandtl number is not varied at the same time. The only difference between Sets DS and DSS1 and DSS2 is that in Set DS, the SGS entropy diffusion is omitted and thus the effective Rayleigh and Prandtl numbers vary as a function of ℱ_n, whereas in the latter two, they approach constant values for low ℱ_n, see Figs. 3 and 7b. Changing the effective Prandtl number leads to a dramatic change in the way convection transports energy in that the sound speed (temperature) fluctuations are enhanced over the velocity fluctuations with increasing Prandtl number, see Fig. 8a. This means that in Set DS the temperature fluctuations become increasingly more important in the enthalpy transport as ℱ_n (Pr) decreases (increases). This is immediately reflected in the overshooting depth as a smaller velocity fluctuation is required to carry the same flux. Figure 8b shows that in Set DSS1 the ratio of the temperature and velocity fluctuations remains practically constant as a function of ℱ_n. The break in the power law in Sets K and Kh can be understood similarly by the decreasing effective Prandtl number as ℱ_n increases.

Fig. 7. Effective Péclet (panel a) and Prandtl (panel b) numbers at the base of the OZ. The vertical dotted line shows the approximate position of the break in the power laws in the overshooting depth in Fig. 6b.

Fig. 8. Ratio of the rms vertical velocity and sound speed fluctuations from the runs in Set DS (panel a) and the runs in Set DSS (panel b).

To connect the current findings to earlier studies, it is necessary to study the Prandtl number regimes explored by Singh et al. (1998), Tian et al. (2009), and Hotta (2017). In the two former studies a Smagorinsky viscosity ν_S was computed from the flow and entropy diffusion according to χ_S = Pr_Sν_S with a constant SGS Prandtl number of ${\mathrm{Pr}}_{\mathrm{S}}=\frac{1}{3}$ was used. While this leads to the same average dependence of the diffusion coefficients as in the present study², the SGS entropy diffusion in these studies was set explicitly to zero in radiative layers. This means that while the effective Prandtl number in the CZ is fixed, it increases in the overshoot layer as the flux is decreased. Thus the sensitivity to the Prandtl number discussed above is likely to explain the steep power laws found by Singh et al. (1998) and Tian et al. (2009). On the other hand, Hotta (2017) used a slope-limited diffusion method where the effective diffusion coefficients are also likely to be proportional to the gradients at small (grid) scales, that is, ν_sl ∝ u′Δx and χ_sl ∝ s′Δx. The velocity (entropy) fluctuations scale like ${\mathcal{F}}_{\mathrm{n}}^{1/3}$ (${\mathcal{F}}_{\mathrm{n}}^{2/3}$) (cf. Käpylä et al. 2019b) and thus this would lead to ${\mathcal{F}}_{\mathrm{n}}^{-1/3}$ scaling for the effective Prandtl number Pr_sl = ν_sl/χ_sl. Whether this simplistic picture of the effective Prandtl number with slope-limited diffusion when applied separately for the momentum and entropy equations is correct remains to be tested with numerical experiments. However, if it is true, then the effective Prandtl number is increasing with decreasing flux and is likely to contribute to the steep power law reported by Hotta (2017).

Furthermore, the studies of Singh et al. (1998), Tian et al. (2009), and Hotta (2017) all considered the bottom of the CZ to be fixed and given by the initial non-convecting state. Although all of these studies used a fixed profile for the heat conductivity, which fixes z_CZ, it does not necessarily stay at the same position as in the initial state; compare, for example, the solid blue and red and green curves in Fig. 6a. This is particularly clear for Sets DS, DSS1, and DSS2 which are similar to the set-up of Singh et al. (1998). However, it is hard to assess whether such a systematic error is present in the results of Singh et al. (1998). In the simulations of Hotta (2017) a similar issue is also possible but it appears that this effect may be small (e.g. his Fig. 8).

Extrapolating from the current data of the Kramers runs (Sets K and Kh) to solar conditions suggests that the overshooting depth for ${\mathcal{F}}_{\text{n}}^{\odot }$ is 𝒪(0.2H_p). This is in better agreement with the constraints from helioseismology (e.g. Basu 1997) than the estimates using the steeper power laws of the earlier numerical studies (e.g. Hotta 2017). However, this estimate should be considered as an upper limit because rotation and magnetic fields, both of which reduce overshooting (e.g. Brummell et al. 2002; Ziegler & Rüdiger 2003; Käpylä et al. 2004), were omitted in the present study. The current estimate is also of the same order of magnitude as in early analytic models of overshooting (van Ballegooijen 1982; Schmitt et al. 1984; Pidatella & Stix 1986) and in two-dimensional anelastic convection models with solar-like parameters (Rogers et al. 2006).

4.3. Dependence on Reynolds and Péclet numbers

In an earlier study, Hotta (2017) concluded that the overshooting depth depends strongly on the resolution of the simulations. The numerical models of Hotta (2017) use a numerical diffusion scheme based on slope limiters where the effective Reynolds and Péclet numbers depend on the grid spacing. Here the explicit viscosity and entropy diffusion are varied to study this effect. Run K5 is taken as a reference run, and the diffusion coefficients were varied within current computational limits in Set R. Run K5 is referred to as Run R3 in Set R. The simulation strategy was such that two branches of runs were performed by taking a thermally saturated snapshot of Run K5 as a basis. In the low-Re branch the grid resolution was kept fixed and the diffusivities were increased (Runs R1 and R2). In the high-Re branch (Runs R4-7) the diffusivities were decreased, and if necessary, a snapshot from a previous simulation was re-meshed to a higher grid resolution (Runs R5 and R6). These two branches were ran consecutively such that the previous runs were first run to a thermally saturated state before changing the diffusivities for the next run to avoid long transients. The results for the normalised overshooting depth as a function of Re = Pe_SGS are shown in Fig. 9.

Fig. 9. Overshooting depth normalised by the pressure scale height at zCZ as a function of Re for Set R. The inset shows dos/HP as a function of Ma = urms/(gd)1/2.

The current results suggest that overshooting is roughly constant as a function of Re. We note, however, that the data points with the highest values of Re (Runs R6 and R7) in Fig. 9 could not be run sufficiently long to establish that they are truly in a statistically stationary state. Thus the values of d_os from these runs should be considered as upper limits. In any case, these results are at odds with those obtained by Hotta (2017), who found a steeply declining trend as a function of Re. This, however, is likely because Hotta (2017) modified the heat conductivity in the radiative layer to speed up thermal relaxation (see Sect. 4.5), and possibly exacerbated by the varying effective Prandtl number in his models (Sect. 4.2).

The inset of Fig. 9 shows d_os/H_p as a function of Ma which quantifies the overall magnitude of the convective velocity. The current data suggest that the overshooting depth is independent of the overall velocity. This is at odds with Zahn (1991), for instance, who derived a Ma^3/2 dependence. However, the range of values explored here is too narrow to draw definite conclusions.

4.4. Transition from the nearly adiabatic to the radiative zone

The superadiabatic temperature gradient from the runs in Set K is shown in Fig. 10. The value of ∇ − ∇_ad in the RZ is close to that of the hydrostatic solution with ∇T = const. in a polytropic atmosphere with adiabatic index n = 13/4, that is, ${\mathrm{\nabla }}_{\mathrm{rad}}^{(\mathrm{hs})}-{\mathrm{\nabla }}_{\mathrm{ad}}$=−17/85 ≈ −0.165 (see also Käpylä et al. 2019a). The transition from nearly adiabatic to the radiative gradient becomes increasingly sharper as ℱ_n decreases (e.g. Käpylä et al. 2007). Furthermore, the temperature gradient in the upper part of the OZ also approaches adiabatic as a function of ℱ_n, suggesting penetration in the nomenclature of Zahn (1991).

Fig. 10. Superadiabatic temperature gradient at the bottom of the CZ in Set K. The normalised energy flux in each run is indicated in the legend. The vertical dotted line indicates the position of the bottom of the CZ in the initial state, and ${\mathrm{\nabla }}_{\mathrm{rad}}^{(\mathrm{hs})}$ corresponds to the hydrostatic solution in the case where ∇T = const.

The CZ is characterised by ∇ ≈ ∇_ad, whereas in the radiative zone, $\mathrm{\nabla }\approx {\mathrm{\nabla }}_{\mathrm{rad}}^{(\mathrm{hs})}$. Thus a rough estimate of the width of the transition between the zones and its dependence on ℱ_n is obtained by computing the vertical derivative of ∇ − ∇_ad, taking its maximum, and computing where the derivative drops below a fixed fraction of the maximum. Here the threshold was set at half the maximum value. The results for the computed depth of the transition δ_trans from Set K are shown as a function of ℱ_n in Fig. 11. The results indicate a power-law ${\delta }_{\mathrm{trans}}\propto {\mathcal{F}}_{\mathrm{n}}^{0.30}$ for ℱ_n ≲ 10⁻⁴. For higher values of ℱ_n the approximation ∇ ≈ ∇_ad is no longer accurate in the CZ and the results deviate from the general trend. Extrapolating from δ_trans ≈ 0.1H_p for ℱ_n = 1.9 × 10⁻⁷, a value of δ_trans ≈ 7.9 × 10⁻³H_p is obtained for the solar value of ℱ_n ≈ 4 × 10⁻¹¹. This corresponds to roughly 400 km and indicates a sharp transition between overshoot and radiative zones. This is similar to what Schmitt et al. (1984) found based on a non-local ML model.

Fig. 11. Depth of the transition δtrans from nearly adiabatic to radiative zones as a function of input flux ℱn normalised by the pressure scale height at zCZ from Set K.

4.5. Modified heat conductivity in the radiative zone

The study of Hotta (2017) reached the lowest value of ℱ_n in the literature thus far, with ℱ_n = 5 × 10⁻⁷. However, these results were obtained by modifying the heat conductivity in the radiative and overshoot layers while the simulations were running. This was done to achieve a statistically stationary state without having to run a full thermal diffusion time. This procedure is sometimes used in anelastic simulations in order to avoid having to run a prohibitively long Kelvin–Helmholtz time (e.g. Brun et al. 2011, 2017). While this procedure can potentially shorten the time to saturation considerably, it can have serious repercussions for overshooting. To demonstrate this, a new run was branched off from Run S7 in which the profile of $\overline{K}(z)$ in the OZ and RZ was modified. The procedure entails computing the energy flux from the non-relaxed run and modifying $\overline{K}$ such that the sum of all the fluxes matches F_tot in the RZ and OZ (Hotta 2017), that is,

$$\begin{array}{c}\hfill \overline{K}\prime =\frac{{\overline{F}}_{\mathrm{enth}}+{\overline{F}}_{\mathrm{kin}}+{\overline{F}}_{\mathrm{visc}}-F_{\mathrm{tot}}}{{\partial }_z\overline{T}}·\end{array}$$(36)

The original and modified profiles of $\overline{K}$ for Runs S7 and S7m are shown in Fig. 12. It is important to note that this procedure alters a crucial system parameter of the model and that the modification can be applied at an arbitrary time in the non-relaxed phase of the simulation.

Fig. 12. Profiles of K, normalised by Kbot, from Runs S7 (black) and S7m (red). The bottom of the initially unstable layer is indicated by the vertical dotted line at z = 0.

What happens in practice is that in the early phases of the simulations the cooling at the surface, in combination with weak radiative diffusion in the CZ, drives efficient convection that overshoots significantly into the radiative layer. This leads to a nearly adiabatic temperature gradient and to a reduced radiative flux in the upper part of the RZ. If ℱ_n is low, the radiative flux is not replenished rapidly enough and the initially vigorous convection cannot be maintained. This leads to a long period of slow evolution in which part of the heat coming from below is deposited in the RZ, OZ, and DZ such that the temperature gradient gradually steepens there to ultimately allow for the total flux to be transmitted. In the standard scenario (Run S7), the heat conductivity is fixed and the temperature gradient steepens, which means that the upper part of the RZ becomes less stiff, allowing relatively deep overshooting, see Fig. 13. The situation is exactly the opposite in Run S7m: the temperature gradient remains shallow and the stratification is significantly stiffer in comparison to the case where $\overline{K}$ was not altered.

Fig. 13. Time-averaged total (black dash-dotted), convective (black), enthalpy (blue), radiative (dark purple), kinetic energy (light purple), and cooling (green) fluxes from Runs S7 (solid) and S7m (dashed). The vertical dotted lines indicate the bottoms of the buoyancy (BZ) and Deardorff zones (DZ). The bottom of the overshoot zone (OZ) for Run S7 (S7m) is denoted by the dashed (dot-dashed) vertical line.

Figure 13 also shows that the modification of the heat conductivity has serious repercussions for the overshooting depth: d_os is reduced by roughly 30 per cent from 0.33H_p in Run K7 to 0.23H_p in Run K7m. This result demonstrates that changing $\overline{K}$ during the run leads to a substantial underestimation of the overshooting depth. This is particularly relevant for the higher resolution cases where the modification of $\overline{K}$ presumably has to be made at an earlier stage. This might explain the strongly decreasing overshooting depth as a function of ℱ_n in the study of Hotta (2017).

5. Conclusions

The scaling of convective overshooting at the base of the CZ was studied as a function of the imposed energy flux, Reynolds number, and different heat conduction profiles and prescriptions. Using heat conductivity based on Kramers opacity, or a similar smoothly varying, but fixed, profile leads to a $d_{\mathrm{os}}\propto {\mathcal{F}}_{\mathrm{n}}^{0.08}$ dependence for ℱ_n ≲ 10⁻⁵. Furthermore, d_os is consistent with a constant as a function of the Reynolds and Péclet number in the range Re = 9…523 in cases where ${\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}=1$. A somewhat steeper power, ${\mathcal{F}}_{\mathrm{n}}^{0.12}$ was found in cases with a fixed step profile for the heat conduction. These results thus indicate a much milder dependence on the imposed energy flux than previous studies in the literature (Singh et al. 1998; Hotta 2017). Numerical experiments with set-ups where the SGS diffusion was turned off led to a steep power law ($d_{\mathrm{os}}\propto {\mathcal{F}}_{\mathrm{n}}^{0.27}$) similar to the power laws reported earlier. Otherwise identical runs with relatively weak SGS diffusion with ${\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}=5$ and 10, on the other hand, produced shallower dependences ($d_{\mathrm{os}}\propto {\mathcal{F}}_{\mathrm{n}}^{0.07}$, and $d_{\mathrm{os}}\propto {\mathcal{F}}_{\mathrm{n}}^{0.12}$, respectively). The cause for a steep power law in the case without SGS entropy diffusion is that the effective Prandtl number increases proportional to ℱ_n, causing the temperature fluctuations to increase. In such cases a smaller velocity fluctuation is needed to carry the same flux, which leads to reduced overshooting. This is the most likely cause for the steep power laws reported by Singh et al. (1998) and Tian et al. (2009), where the effective Prandtl number in the overshoot layer is likely much larger than unity. A similar argument can tentatively be made regarding the slope-limited diffusion used by Hotta (2017), but this should be tested with further experiments.

Furthermore, the current results indicate that modifying the heat conductivity in the layers below the CZ (e.g. Brun et al. 2017; Hotta 2017) leads to a substantial underestimation of the overshooting depth. The only way to extract reliable scaling of the overshooting depth as a function of ℱ_n currently is to run the simulations self-consistently to a thermally relaxed state without modifying the system parameters such as the heat conductivity. The present study also demonstrates the limits of this approach in that the runs with the lowest input flux require integration times of the order of several months even at a relatively low resolution of 288³. A more promising alternative to speeding up thermal saturation is to alter the thermodynamic quantities instead (e.g. Hurlburt et al. 1986; Anders et al. 2018). However, even this method has its limitations, and the applicability of this approach, for example, to rotating convection in spherical shells remains to be demonstrated.

The current results suggest that the overshooting depth in the Sun would be of the order of 𝒪(0.2),H_p which is somewhat higher than the canonical estimates of (0.05…0.1)H_p from helioseismology. Furthermore, the transition from overshoot to radiative zone is expected to be abrupt and occur over a depth of roughly 400 km. However, the present models lack rotation and magnetic fields, which have a significant effect on the convective flows in the deep parts of the solar CZ. Possibly the biggest caveat is the unrealistically large SGS Prandtl number we used in the current study. Set-ups where these constraints are relaxed will be explored in future publications.

Acknowledgments

The anonymous referee is acknowledged for the constructive comments on the manuscript. Axel Brandenburg is acknowledged for his insightful comments on the manuscript. The simulations were performed using the supercomputers hosted by CSC – IT Center for Science Ltd. in Espoo, Finland, who are administered by the Finnish Ministry of Education and the Gauss Center for Supercomputing for the Large-Scale computing project “Cracking the Convective Conundrum” in the Leibniz Supercomputing Centre’s SuperMUC supercomputer in Garching, Germany. This work was supported in part by the Deutsche Forschungsgemeinschaft Heisenberg programme (grant No. KA 4825/1-1) and the Academy of Finland ReSoLVE Centre of Excellence (grant No. 307411).

References

All Tables

All Figures

	Fig. 1. Normalised vertical velocity ${\stackrel{\sim }{u}}_z=u_z/{(gd)}^{1/2}$ (colours) and streamlines of the flow from Runs R4 (left) and R7 (right).
In the text

	Fig. 2. Horizontally averaged and normalised heat conductivity $\stackrel{\sim }{K}=\overline{K}(z)/K_{\mathrm{bot}}$ (a) and Prandtl number Pr(z) (b) from Runs K4, P4, S4, and DS5. The dashed black (red) line in the lower panel indicates ${\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}=1$ (${\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}=5$).
In the text

	Fig. 3. Rayleigh numbers RaSGS (black) and RaRad (red) as functions of ℱn from Set Kh. The purple (magenta) line shows RaRad (RaSGS) from Set DS (DSS). The dotted lines show scalings expected from ML arguments.
In the text

Fig. 4. Normalised rms velocity ${\stackrel{\sim }{u}}_{\mathrm{rms}}=u_{\mathrm{rms}}/\sqrt{\mathit{gd}}$ in the CZ as a function of ℱn for the simulation sets indicated by the legend. The dotted line is proportional to ${\mathcal{F}}_{\mathrm{n}}^{1/3}$. The inset shows ${\stackrel{\sim }{u}}_{\mathrm{rms}}{\mathcal{F}}_{\mathrm{n}}^{-1/3}$ for the same runs. The purple diamond denotes Run DS5h.

In the text

Fig. 5. Panels a, c, and e: total (black dash-dotted lines), convective (black), enthalpy (blue), kinetic energy (light purple), radiative (dark purple), and cooling (green) fluxes and the superadiabatic temperature gradient (red) from Runs K4, S4, and DS5, respectively. Panels b, d, and f: total averaged vertical forces (solid lines) and the power of the forces (dashed) on the upflows (red) and downflows (blue) from the same runs. The dotted lines in these panels show the corresponding viscous force and its power. The shaded areas indicate the BZ (darkest), DZ, and OZ (lightest) and the thick black line at the horizontal axis denotes the MZ.

In the text

Fig. 6. Panel a: vertical (z) coordinates of the bottom of the CZ (zCZ, solid lines) and OZ (zOZ, dashed). The dotted red line indicates the bottom of the domain. Panel b: overshooting depth dos normalised by the pressure scale height Hp as a function of ℱn for Sets K (black), Kh (grey), P (blue), S (red), Sh (green), DS (purple), DSS1 (cyan), and DSS2 (orange). The purple diamond denotes Run DS5h. The dotted lines show approximate power laws from fits to simulation data; see Table 5. Panel c: comparison of Sets DS, DSS1, and DSS2 with the studies of Singh et al. (1998) (red) and Tian et al. (2009) (blue).

In the text

	Fig. 7. Effective Péclet (panel a) and Prandtl (panel b) numbers at the base of the OZ. The vertical dotted line shows the approximate position of the break in the power laws in the overshooting depth in Fig. 6b.
In the text

	Fig. 8. Ratio of the rms vertical velocity and sound speed fluctuations from the runs in Set DS (panel a) and the runs in Set DSS (panel b).
In the text

	Fig. 9. Overshooting depth normalised by the pressure scale height at zCZ as a function of Re for Set R. The inset shows dos/HP as a function of Ma = urms/(gd)1/2.
In the text

	Fig. 10. Superadiabatic temperature gradient at the bottom of the CZ in Set K. The normalised energy flux in each run is indicated in the legend. The vertical dotted line indicates the position of the bottom of the CZ in the initial state, and ${\mathrm{\nabla }}_{\mathrm{rad}}^{(\mathrm{hs})}$ corresponds to the hydrostatic solution in the case where ∇T = const.
In the text

	Fig. 11. Depth of the transition δtrans from nearly adiabatic to radiative zones as a function of input flux ℱn normalised by the pressure scale height at zCZ from Set K.
In the text

	Fig. 12. Profiles of K, normalised by Kbot, from Runs S7 (black) and S7m (red). The bottom of the initially unstable layer is indicated by the vertical dotted line at z = 0.
In the text

Fig. 13. Time-averaged total (black dash-dotted), convective (black), enthalpy (blue), radiative (dark purple), kinetic energy (light purple), and cooling (green) fluxes from Runs S7 (solid) and S7m (dashed). The vertical dotted lines indicate the bottoms of the buoyancy (BZ) and Deardorff zones (DZ). The bottom of the overshoot zone (OZ) for Run S7 (S7m) is denoted by the dashed (dot-dashed) vertical line.

In the text