Overshooting in simulations of compressible convection

TL;DR

This study uses 3D simulations to analyze how overshooting depth in compressible convection scales with energy flux, revealing a shallow power law dependence and implications for solar convection zone modeling.

Contribution

It provides a detailed analysis of overshooting depth scaling with flux using different heat conductivity profiles and highlights the importance of the Prandtl number and heat conductivity treatment.

Findings

Overshooting depth scales as approximately flux to the power of 0.08 to 0.12.

Smooth heat conductivity profiles lead to shallower overshoot scaling.

Estimated solar overshooting depth is about 0.2 pressure scale heights, higher than helioseismic estimates.

Abstract

(abridged) Context: Convective motions overshooting to regions that are formally convectively stable cause extended mixing. Aims: To determine the scaling of overshooting depth ($d_{\rm os}$) at the base of the convection zone as a function of imposed energy flux ($\mathscr{F}_{\rm n}$) and to estimate the extent of overshooting at the base of the solar convection zone. Methods: Three-dimensional Cartesian simulations of compressible non-rotating convection with unstable and stable layers are used. The simulations use either a fixed heat conduction profile or a temperature and density dependent formulation based on Kramers opacity law. The simulations cover a range of almost four orders of magnitude in the imposed flux. Results: A smooth heat conduction profile (either fixed or through Kramers opacity law) leads to a relatively shallow power law with $d_{\rm os}\propto \mathscr{F}_{\rm n}^{0.08}$ for low $\mathscr{F}_{\rm n}$. A fixed step-profile of the heat conductivity at the bottom of the convection zone leads to a somewhat steeper dependency with $d_{\rm os}\propto \mathscr{F}_{\rm n}^{0.12}$. Experiments with and without subgrid-scale entropy diffusion revealed a strong dependence on the effective Prandtl number which is likely to explain the steep power laws as a function of $\mathscr{F}_{\rm n}$ reported in the literature. Furthermore, changing the heat conductivity artificially below the convection zone is shown to lead to substantial underestimation of overshooting depth. Conclusions: Extrapolating from the results obtained with smooth heat conductivity profiles suggest that the overshooting depth for the solar flux is of the order of $0.2$ pressure scale heights at the base of the convection zone which is two to four times higher than estimates from helioseismology.

Figures

Tables

Table 1. Table 1: Summary of the runs with smoothly varying profiles of K 𝐾 K .

Table 2. Table 2: Summary of the runs with a fixed step profile of K 𝐾 K .

Table 3. Table 3: Summary of the runs with a double-step profile for K 𝐾 K .

Table 4. Table 4: Summary of the runs with varying ν 𝜈 \nu and χ SGS ( 1 ) superscript subscript 𝜒 SGS 1 \chi_{\rm SGS}^{(1)} .

Table 5. Table 5: Power-law exponents and standard mean errors from fitting d os ∝ ℱ n α proportional-to subscript 𝑑 os superscript subscript ℱ n 𝛼 d_{\rm os}\propto\mathscr{F}_{\rm n}^{\alpha} .

Equations

$${\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.$$

$$\frac{D\ln \rho }{Dt}$$

$$\frac{D\mathbf{u}}{Dt}$$

$$T\frac{Ds}{Dt}$$

$${\mathsf{S}}_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})-\frac{1}{3}{\delta }_{ij}\nabla \cdot \mathbf{u}.$$

$${\mathbf{F}}_{\mathrm{rad}}=-K\nabla T,$$

$$K=\frac{16{\sigma }_{\mathrm{SB}}{T}^3}{3\kappa \rho },$$

$$\kappa ={\kappa }_0(\rho /{\rho }_0{)}^a(T/T_0{)}^b,$$

$$K(\rho ,T)=K_0(\rho /{\rho }_0{)}^{-(a+1)}(T/T_0{)}^{3-b}.$$

$${\mathbf{F}}_{\mathrm{SGS}}=-\rho T{\chi }_{\mathrm{SGS}}^{(1)}\nabla s^{\prime },$$

$${\Gamma }_{\mathrm{cool}}=-{\Gamma }_{0}f(z)(T_{\mathrm{cool}}-T),$$

$$\frac{\partial u_x}{\partial z}=\frac{\partial u_y}{\partial z}=u_z=0.$$

$$\frac{\partial T}{\partial z}=-\frac{F_{\mathrm{bot}}}{K_{\mathrm{bot}}},$$

$$[x]=d,[t]=\sqrt{d/g},[\rho ]={\rho }_0,[s]=c_{\mathrm{P}},$$

$${\mathrm{Pr}}_{\mathrm{SGS}}^{(1)}=\frac{\nu }{{\chi }_{\mathrm{SGS}}^{(1)}},$$

$$\mathrm{Pr}(z)=\frac{\nu }{\chi (z)},$$

$$\frac{K_i}{K_j}=\frac{n_i^{\prime }+1}{n_j^{\prime }+1},$$

$${\mathcal{F}}_{\mathrm{n}}=F_{\mathrm{bot}}/{\rho }_{\mathrm{bot}}{c}_{\mathrm{s},\mathrm{bot}}^3,$$

$$\mathrm{Re}=\frac{u_{\mathrm{rms}}}{\nu k_1},{\mathrm{Pe}}_{\mathrm{SGS}}=\frac{u_{\mathrm{rms}}}{{\chi }_{\mathrm{SGS}}^{(1)}{k}_1},$$

$${\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}}},$$

$${\mathrm{Pr}}_{\mathrm{OZ}}^{\mathrm{eff}}=\frac{\nu }{({\chi }_{\mathrm{SGS}}^{(1)}+{\overline{\chi }}_{\mathrm{OZ}})}.$$

$${\mathrm{Ra}}^{\mathrm{Rad}}$$

$${\mathrm{Ra}}^{\mathrm{SGS}}$$

$${\overline{F}}_{\mathrm{rad}}$$

$${\overline{F}}_{\mathrm{enth}}$$

$${\overline{F}}_{\mathrm{kin}}$$

$${\overline{F}}_{\mathrm{visc}}$$

$${\overline{F}}_{\mathrm{cool}}$$

$${\overline{F}}_{\mathrm{conv}}={\overline{F}}_{\mathrm{enth}}+{\overline{F}}_{\mathrm{kin}}.$$

$${\overline{F}}_{\mathrm{rad}}=-\overline{K}\frac{\partial \overline{T}}{\partial z}=\frac{\overline{K}g}{c_{\mathrm{P}}}\frac{\nabla }{{\nabla }_{\mathrm{ad}}},$$

$$F_{\mathrm{tot}}=\frac{\overline{K}g}{c_{\mathrm{P}}}\frac{{\nabla }_{\mathrm{rad}}}{{\nabla }_{\mathrm{ad}}},{\overline{F}}_{\mathrm{rad}}^{\mathrm{ad}}=\frac{\overline{K}g}{c_{\mathrm{P}}},$$

$$\mathrm{Nu}=\frac{F_{\mathrm{tot}}}{{\overline{F}}_{\mathrm{rad}}^{\mathrm{ad}}}=\frac{{\nabla }_{\mathrm{rad}}}{{\nabla }_{\mathrm{ad}}}.$$

$$d_{\mathrm{os}}=⟨z_{\mathrm{CZ}}(t)-z_{\mathrm{OS}}(t){⟩}_t,$$

$${\mathrm{Ra}}^{\mathrm{Flux}}=\mathrm{Nu}{\mathrm{Ra}}^{\mathrm{Rad}}=\frac{F_{\mathrm{tot}}{d}^4}{\rho {\chi }^2\nu }\left(-\frac{1}{c_{\mathrm{P}}}\frac{\mathrm{d}s}{\mathrm{d}z}\right),$$

Full text

*Aims. *We aim to determine the scaling of the overshooting depth ($d_{\rm os}$) at the base of the convection zone as a function of imposed energy flux ($\mathscr{F}_{\rm n}$) and to estimate the extent of overshooting at the base of the solar convection zone.

*Methods. *Three-dimensional Cartesian simulations of hydrodynamic compressible non-rotating convection with unstable and stable layers were used. The simulations used either a fixed heat conduction profile or a temperature- and density-dependent formulation based on Kramers opacity law. The simulations covered a range of almost four orders of magnitude in the imposed flux, and the sub-grid scale diffusivities were varied so as to maintain approximately constant supercriticality at each flux.

*Results. *A smooth heat conduction profile (either fixed or through Kramers opacity law) leads to a relatively shallow power law with $d_{\rm os}\propto\mathscr{F}_{\rm n}^{0.08}$ for low $\mathscr{F}_{\rm n}$. A fixed step-profile of the heat conductivity at the bottom of the convection zone leads to a somewhat steeper dependency on $d_{\rm os}\propto\mathscr{F}_{\rm n}^{0.12}$ in the same regime. Experiments with and without subgrid-scale entropy diffusion revealed a strong dependence on the effective Prandtl number, which is likely to explain the steep power laws as a function of $\mathscr{F}_{\rm n}$ reported in the literature. Furthermore, changing the heat conductivity artificially in the radiative and overshoot layers to speed up thermal saturation is shown to lead to a substantial underestimation of the overshooting depth.

*Conclusions. *Extrapolating from the results obtained with smooth heat conductivity profiles, which are the most realistic set-up we considered, suggest that the overshooting depth for the solar energy flux is about 20% of the pressure scale height at the base of the convection zone. This is two to four times higher than the estimates from helioseismology. However, the current simulations do not include rotation or magnetic fields, which are known to reduce convective overshooting.

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; Roxburgh & Simmons 1993; Hurlburt et al. 1994; 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 convection zone, 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 convection zone. According to mixing length arguments, the enthalpy flux $\overline{\bm{F}}_{\rm enth}=c_{\rm P}\overline{(\rho{\bm{u}})^{\prime}T^{\prime}}$, where $c_{\rm P}$ is the heat capacity at constant pressure, $\rho$ is the density, ${\bm{u}}$ is the convective velocity, and $T$ is the temperature, and the apostrophes (overbar) denote fluctuations (averages) that can be approximated as $\overline{F}_{\rm enth}=\phi\overline{\rho}u^{\prime 3}$ (Brandenburg 2016). Assuming that $\overline{F}_{\rm enth}\approx F_{\odot}$, it is possible to construct a normalised energy flux (e.g. Brandenburg et al. 2005),

[TABLE]

The last approximation is justified by the fact that the factor $\phi$ 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 convection zone (see, e.g. the solar model of Stix 2002). Using solar values at $r_{0}=0.71\leavevmode\nobreak\ R_{\odot}$, where $R_{\odot}=7\cdot 10^{8}$ m is the solar radius, $F_{\odot}=L_{\odot}/(4\pi r_{0}^{2})\approx 6\cdot 10^{25}$ W m*-2*, where $L_{\odot}=3.84\cdot 10^{26}$ W is the solar luminosity, $\rho\approx 200$ kg m*-3*, and $c_{\rm s}\approx 200$ km s*-1*, the normalised flux at the base of the solar CZ is $\mathscr{F}_{\rm n}^{(\odot)}\approx 4\cdot 10^{-11}$. Thus the Mach number in the deep convection zone is about $10^{-4}$. 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 $\mathscr{F}_{\rm n}=5\cdot 10^{-7}$, which is at the limits of current numerical feasibility, and obtained a power law $d_{\rm os}/H_{\rm p}\propto\mathscr{F}_{\rm n}^{0.31}$ for the overshooting depth $d_{\rm os}$. This led Hotta (2017) to estimate that the overshooting at the base of the solar convection zone 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_{\rm os}$ on $\mathscr{F}_{\rm 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 convection zone to dynamically adapt to changes in the thermodynamic state of the system (Käpylä et al. 2019b) 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 Code111https://github.com/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

[TABLE]

are solved, where $D/Dt=\partial/\partial t+{\bm{u}}\cdot\bm{\nabla}$ is the advective derivative, $\rho$ is the density, ${\bm{u}}$ is the velocity, $\bm{g}=-g\hat{\bm{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 $\nu$ is the constant kinematic viscosity. Furthermore, ${\bm{F}}_{\rm rad}$ and ${\bm{F}}_{\rm SGS}$ are the radiative and turbulent subgrid scale (SGS) fluxes, respectively, and $\Gamma_{\rm cool}$ describes cooling at the surface (see below). $\bm{\mathsf{S}}$ is the traceless rate-of-strain tensor with

[TABLE]

An optically thick fully ionised gas is considered, where radiation is modelled through diffusion approximation. The ideal gas equation of state $p=(c_{\rm P}-c_{\rm V})\rho T={\cal R}\rho T$ applies, where ${\cal R}$ is the gas constant, and $c_{\rm V}$ is the specific heat at constant volume. The radiative flux is given by

[TABLE]

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(\rho,T)$. In the latter case, $K$ is given by

[TABLE]

where $\sigma_{\rm SB}$ is the Stefan-Boltzmann constant and $\kappa$ is the opacity. $\kappa$ is assumed to obey a power law

[TABLE]

where $\rho_{0}$ and $T_{0}$ are reference values of density and temperature. Equations (7) and (8) combine into

[TABLE]

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, $\chi=K/(c_{\rm P}\rho)$, additional turbulent SGS diffusivity is used in the entropy equation to keep the simulations numerically feasible. Here the SGS flux is formulated as

[TABLE]

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_{\rm SGS}^{(1)}$ is constant in the whole domain. ${\bm{F}}_{\rm SGS}$ has a negligible contribution to the net horizontally averaged energy flux, such that $\overline{{\bm{F}}}_{\rm SGS}\approx 0$. The current SGS formulation is similar to those used in Käpylä et al. (2007), Brown et al. (2010), and Käpylä et al. (2019b), for example.

The cooling at the surface is described by

[TABLE]

where $\Gamma_{0}$ is a cooling luminosity, $T=e/c_{\rm V}$ is the temperature where $e$ is the internal energy, and where $T_{\rm cool}=T_{\rm 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_{\rm bot}\leq z\leq z_{\rm top}$ is the vertical coordinate, where $z_{\rm bot}/d=-0.45$, $z_{\rm 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_{\rm 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_{\rm 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 P2I), the two lower layers are polytropic with polytropic indices $n_{1}=3.25$ ($z_{\rm bot}/d\leq z/d\leq 0$) and $n_{2}=1.5$ ($0\leq z/d\leq 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_{1}$ 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. 2019b). 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 P3) all three layers are polytropic with indices $(n_{1},n_{2},n_{3})=(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_{2}=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{dg}$.

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

[TABLE]

The temperature gradient at the bottom boundary is set according to

[TABLE]

where $F_{\rm bot}$ is the fixed input flux and $K_{\rm bot}(x,y,z_{\rm bot})$ is the value of the heat conductivity at the bottom of the domain. On the upper boundary a constant temperature $T=T_{\rm 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

[TABLE]

where $\rho_{0}$ is the initial value of the density at $z=z_{\rm top}$. The models with Kramers heat conductivity are fully defined by choosing the value of the kinematic viscosity $\nu$, the gravitational acceleration $g$, the values of $a$, $b$, $K_{0}$, $\rho_{0}$, $T_{0}$ and the SGS Prandtl number

[TABLE]

along with $z$-dependent cooling profile $f(z)$. The values of $K_{0}$, $\rho_{0}$, and $T_{0}$ are subsumed into a new variable $\widetilde{K}_{0}=K_{0}\rho_{0}^{a+1}T_{0}^{b-3}$ that is fixed by assuming the radiative flux at $z_{\rm bot}$ to equal $F_{\rm 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_{\rm p}^{\rm top}/d=\mathcal{R}T_{\rm top}/gd$ sets the initial pressure scale height at the surface, thus determining the initial density stratification. All of the current simulations have $\xi_{0}=0.054$.

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

[TABLE]

where $\chi(z)=K(z)/c_{\rm P}\rho(z)$, sets the relative importance of viscous to temperature diffusion. We note that in general $\Pr=\Pr({\bm{x}},t)$ because $\rho=\rho({\bm{x}},t)$. In cases where $K$ has a piecewise constant profile, it can be represented in terms of polytropic indices $(n^{\prime}_{1},n^{\prime}_{2},n^{\prime}_{3}),$ 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_{2}$ 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

[TABLE]

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

The dimensionless normalised flux is given by

[TABLE]

where $\rho_{\rm bot}$ and $c_{\rm s,bot}$ are the density and the sound speed, respectively, at $z_{\rm bot}$ in the initial non-convecting state.

The input energy flux determines the overall convective velocity realised in the simulations via $u_{\rm conv}\propto\mathscr{F}_{\rm n}^{1/3}$, see Eq. (1) and Käpylä et al. (2019a). Thus if only $\mathscr{F}_{\rm 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 $\mathscr{F}_{\rm n}^{1/3}$. In addition to changing the diffusion coefficients, the cooling luminosity $\Gamma_{0}$ at the surface is also scaled proportionally to $\mathscr{F}_{\rm 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_{\rm bot}=-F_{\rm bot}/(\partial T/\partial z)_{z_{\rm bot}}$. 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 $\nu$ and $\chi_{\rm SGS}^{(1)}$ on $\mathscr{F}_{\rm n}$, however, eliminates these dependences such that the supercriticality and effective Prandtl numbers are nearly constant over the range of $\mathscr{F}_{\rm n}$ considered here. This is discussed in detail in Section 4.1.

In addition to the explicit viscosity, SGS, and radiative diffusion, the advective terms in each of the Equation (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,

[TABLE]

where $u_{\rm rms}$ is the volume-averaged rms velocity, and $k_{1}=2\pi/d$ is an estimate of the largest eddies in the system. Typically $\chi\ll\chi_{\rm SGS}^{(1)}$ in the CZ in most of the current simulations. Furthermore, $\chi$ 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,

[TABLE]

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

[TABLE]

Precise definitions of $z_{\rm CZ}$ and $d_{\rm os}$ are given in Section 3.

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

[TABLE]

Supercriticality of convection is roughly determined by $\mbox{min}({\rm Ra}^{\rm Rad},{\rm Ra}^{\rm 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 (${\rm Ra}^{\rm SGS}=0$) and any change in the flux would be reflected by the then decisive radiative Rayleigh number. For example, a system where ${\rm Pr}$ is fixed ${\rm Ra}^{\rm rad}\propto\mathscr{F}_{\rm n}^{-4/3}$. However, taking this path severely limits the computationally feasible range of $\mathscr{F}_{\rm n}$ because the Reynolds and Péclet number also increase in proportion to $\mathscr{F}_{\rm n}$ and the resolution requirements quickly become prohibitive.

Contributions to the vertical energy flux are

[TABLE]

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:

[TABLE]

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

[TABLE]

where $c_{\rm P}\nabla_{\rm ad}=c_{\rm P}(1-1/\gamma)=\mathcal{R}$, and $g=|\bm{g}|$. Above, $\nabla$ 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

[TABLE]

where $\nabla_{\rm rad}$ is a hypothetical radiative gradient in the absence of convection and $F_{\rm tot}=F_{\rm bot}$. The ratio of the total-to-adiabatic flux is the Nusselt number (e.g. Brandenburg 2016),

[TABLE]

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 ${\rm 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 convection zone 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}_{\rm conv}>0$, whereas in the overshoot zone (OZ) $F_{\rm conv}<0$. This is motivated by the work of Deng & Xiong (2008), who used a similar definition but employed $\overline{F}_{\rm enth}$ instead of $\overline{F}_{\rm conv}$. A definition of the overshooting depth based on $\overline{F}_{\rm 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 convection zone. The bottom of the CZ is denoted by $z_{\rm CZ}$.

The mean overshooting depth $z_{\rm OS}$ is taken to be the position where the horizontally averaged $\overline{F}_{\rm kin}$ drops below 1% of its value at $z_{\rm 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_{\rm CZ}$ is computed from the simulation data using the definition given above. The values of $z_{\rm OS}$ and $z_{\rm CZ}$ are obtained by linear interpolation from the grid points closest to the respective transitions. Furthermore, $z_{\rm CZ}$ and $z_{\rm OS}$ are functions of time. The overshooting depth is defined as

[TABLE]

where $\langle\cdots\rangle_{t}$ denotes a time average over the statistically stationary part of the time series. Error estimates for $d_{\rm 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_{\rm OS}$ , and the buoyancy zone (BZ) is where $\overline{F}_{\rm 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}_{\rm 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_{\rm SGS}^{(1)}=0,$ and in Set DSS1 (DSS2) the SGS diffusion is included with $\Pr_{\rm 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 $\nu$ and $\chi_{\rm 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 P3 set-up, and the remaining sets were initialised with the P2I set-up.

4.1 Basic characterisation of the solutions

Typical flow patterns for Runs R4 and R7 are shown in Figure 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 ${\rm Re}\approx 20\ldots 40$ such as in Run R4 in the left panel of Figure 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 Figure 1 for Run R7.

Figure 2(a) 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_{2}/K_{1}=K_{3}/K_{1}=4/17\approx 0.235$. In Run DS5, $K_{2}/K_{1}=2/85\approx 0.0333$ and the uppermost layer above $z/d=1$ has another constant value of $\overline{K}$ corresponding to $K_{2}/K_{1}=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^{\prime}_{3}=1.5$ is achieved only near the upper surface.

The vertical profiles of $\Pr(z)$ and $\Pr_{\rm SGS}^{(1)}$ are shown in Figure 2(b). The strong variation of $\Pr$ as function of depth is to be contrasted with the constant value of $\Pr_{\rm SGS}^{(1)}$. Although $\Pr\propto\mathscr{F}_{\rm 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 $\mathscr{F}_{\rm n}$ (see Tables 1–4). This is in stark contrast to the solar convection zone, where $\Pr\ll 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 $\Pr_{\rm 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$\mathscr{F}_{\rm n}$ refers to varying the magnitude of $\overline{K}$ proportionally and thus $\chi\propto\mathscr{F}_{\rm n}$. In combination with $\nu\propto\mathscr{F}_{\rm n}^{1/3}$ and the mixing-length estimate $-(H_{\rm p}/c_{\rm P}){\rm d}s/{\rm d}z=\nabla-\nabla_{\rm ad}\propto\mathscr{F}_{\rm n}^{2/3}$ (Vitense 1953), this involves changing ${\rm Ra}^{\rm Rad}\propto\mathscr{F}_{\rm n}^{-2/3}$. This is to be contrasted with the SGS Rayleigh number, where $\chi$ is replaced by $\chi_{\rm SGS}^{(1)}\propto\mathscr{F}_{\rm n}^{1/3}$ , which leads to ${\rm Ra}^{\rm SGS}$ being independent of $\mathscr{F}_{\rm n}$. These scalings are in accordance with the simulations, as can be seen from Figure 3 for Sets Kh, DS, and DSS1. As mentioned earlier, the supercriticality of convection is determined roughly by ${\rm min}({\rm Ra}^{\rm Rad},{\rm Ra}^{\rm SGS})$. In the current simulations the SGS Rayleigh number is almost always smaller. This also means that with ${\rm Ra}^{\rm SGS}\approx$ 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)

[TABLE]

can be seen to vary as ${\rm Ra}^{\rm Flux}\propto\mathscr{F}_{\rm n}^{-2/3}$ given the dependences above. We note that ${\rm Nu}{\rm Ra}^{\rm SGS}$ is independent of $\mathscr{F}_{\rm 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 $\mathscr{F}_{\rm n}^{1/3}$ dependence, with the coefficient of proportionality varying between 1.3 and 1.8. This is in agreement with mixing length 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. 2019a).

Results regarding the energy fluxes and force balance from representative runs from Sets K, S, and DS are shown in Figure 5. The left panels show the contributions to the energy flux and the superadiabatic temperature gradient $\nabla-\nabla_{\rm ad}$. The fluxes for Runs K4 and S4 in Figure 5(a) 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 $\mathscr{F}_{\rm n}$; $1.8\cdot 10^{-5}$ here versus $9\cdot 10^{-6}$ 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_{2}/K_{1}$ 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)

[TABLE]

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 5(b), (d), and (f) show the horizontally averaged vertical total and viscous forces, $\overline{f}_{z}=\overline{\rho Du_{z}/Dt}$, and $\overline{f}_{\rm visc}=2\nu\partial_{i}(\overline{\rho S_{iz}})$, respectively, and the resulting power of the forces ($\overline{P}_{z}=\overline{u_{z}f_{z}}$ and $\overline{P}_{\rm visc}=\overline{u_{z}f_{\rm visc}}$) separately for the upflows ($\uparrow$) and the downflows ($\downarrow$) from Runs K4, S4, and DS5. The force balance in Run K4 (Figure 5b), is very similar to the corresponding Run K in Käpylä et al. (2017), see their Fig. 2(b). 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 convection zone 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 Figure 5(c) 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 $\mathscr{F}_{\rm 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 $\mathscr{F}_{\rm 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 Figure 6(a) where the vertical coordinates of the bottoms of the convection ($z_{\rm CZ}$) and overshoot ($z_{\rm OZ}$) zones are shown as functions of $\mathscr{F}_{\rm n}$. The results for $z_{\rm 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 $\mathscr{F}_{\rm 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 convection zone is generally reduced in comparison to Sets S and Sh. In these runs the depth of the CZ increases as $\mathscr{F}_{\rm n}$ decreases. This is contrasted by the results from Sets P where a Kramers-like, but static, profile of $K$ is used (see Figure 2): here $z_{\rm CZ}$ is practically fixed in the current range of $\mathscr{F}_{\rm n}$.

Figure 6(a) shows that the location of the base of the OZ ($z_{\rm OZ}$) increases monotonically as a function of $\mathscr{F}_{\rm 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 $\mathscr{F}_{\rm n}$ runs, rendering the results of these simulations unusable in the following analysis. The depth of the overshoot zone, $d_{\rm os}$, as a function of $\mathscr{F}_{\rm n}$ is shown in Figure 6b. The results for Sets K, Kh, and P fall almost on top of each other. The data suggests two power laws: $\mathscr{F}_{\rm n}^{0.08}$ for $\mathscr{F}_{\rm n}\lesssim 10^{-5}$ and $\mathscr{F}_{\rm n}^{0.19}$ for $\mathscr{F}_{\rm n}\gtrsim 10^{-5}$; 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_{\rm 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 Section 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 $\mathscr{F}_{\rm n}\approx 10^{-5}$, such that the data points for lower values of $\mathscr{F}_{\rm n}$ lie below those of Set Sh. However, power-law fits for both full and partial ranges are consistent with a $\mathscr{F}_{\rm 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_{\rm os}\propto\mathscr{F}_{\rm n}^{0.27}$), and Sets DSS1 with $d_{\rm os}\propto\mathscr{F}_{\rm n}^{0.07}$ and DSS2 with $d_{\rm os}\propto\mathscr{F}_{\rm 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 ${\rm Re}_{\rm OZ}$ and ${\rm Pe}_{\rm OZ}^{\rm 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 ${\rm Pr}_{\rm SGS}^{(1)}$ increases through Sets DSS1 and DSS2 to DS are indicative that the latter set of runs is sufficiently well resolved numerically. Figure 6(c) 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_{\rm b}$ and values of density and pressure at the bottom of the domain (denoted here as $\rho_{\rm b}$ and $p_{\rm b}$). Thus the normalised flux in both cases is computed from $\mathscr{F}_{\rm n}=F_{\rm b}/F_{\rm n}$ , where $F_{\rm n}=\rho_{\rm b}c_{\rm s}^{3}=\rho_{\rm b}^{-1/2}(\gamma p_{\rm b})^{3/2}$ with $c_{\rm s}^{2}=\gamma p_{\rm b}/\rho_{\rm b}$ and $\gamma=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_{\rm os}\propto\mathscr{F}_{\rm n}^{0.33}$ scaling. On the other hand, Tian et al. (2009) obtained a very shallow dependence with this definition ($d_{\rm os}\propto\mathscr{F}_{\rm n}^{0.05}$) and a much steeper one ($d_{\rm os}\propto\mathscr{F}_{\rm 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 7(a) shows that the effective Péclet numbers in OZ in all of these sets are comparable. However, as is shown in Section 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 $\mathscr{F}_{\rm n}$ , whereas in the latter two, they approach constant values for low $\mathscr{F}_{\rm n}$, see Figs. 3 and 7(b). 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 Figure 8(a). This means that in Set DS the temperature fluctuations become increasingly more important in the enthalpy transport as $\mathscr{F}_{\rm 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 8(b) shows that in Set DSS1 the ratio of the temperature and velocity fluctuations remains practically constant as a function of $\mathscr{F}_{\rm n}$. The break in the power law in Sets K and Kh can be understood similarly by the decreasing effective Prandtl number as $\mathscr{F}_{\rm n}$ increases.

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 $\nu_{\rm S}$ was computed from the flow and entropy diffusion according to $\chi_{\rm S}={\rm Pr}_{\rm S}\nu_{\rm S}$ with a constant SGS Prandtl number of ${\rm Pr}_{\rm S}={\textstyle{1\over 3}}$ was used. While this leads to the same average dependence of the diffusion coefficients as in the present study777This is due to $\nu_{\rm S}=(C_{k}\ell)^{2}\sqrt{\bm{\mathsf{S}}^{2}}\propto u^{\prime}\ell$ with $u^{\prime}\propto\mathscr{F}_{\rm n}^{1/3}$, $C_{k}=$ const. and $\ell\propto\Delta x$ const., and where $\Delta x$ is the grid spacing., 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, $\nu_{\rm sl}\propto u^{\prime}\Delta x$ and $\chi_{\rm sl}\propto s^{\prime}\Delta x$. The velocity (entropy) fluctuations scale like $\mathscr{F}_{\rm n}^{1/3}$ ($\mathscr{F}_{\rm n}^{2/3}$) (cf. Käpylä et al. 2019a) and thus this would lead to $\mathscr{F}_{\rm n}^{-1/3}$ scaling for the effective Prandtl number ${\rm Pr}_{\rm sl}=\nu_{\rm sl}/\chi_{\rm 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_{\rm 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 Figure 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 Figure 8).

Extrapolating from the current data of the Kramers runs (Sets K and Kh) to solar conditions suggests that the overshooting depth for $\mathscr{F}_{\rm n}^{\odot}$ is $\mathcal{O}(0.2H_{\rm 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-${\rm Re}$ branch the grid resolution was kept fixed and the diffusivities were increased (Runs R1 and R2). In the high-${\rm 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 ${\rm Re}={\rm Pe}_{\rm SGS}$ are shown in Figure 9.

The current results suggest that overshooting is roughly constant as a function of ${\rm Re}$. We note, however, that the data points with the highest values of ${\rm Re}$ (Runs R6 and R7) in Figure 9 could not be run sufficiently long to establish that they are truly in a statistically stationary state. Thus the values of $d_{\rm 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 ${\rm Re}$. This, however, is likely because Hotta (2017) modified the heat conductivity in the radiative layer to speed up thermal relaxation (see Section 4.5), and possibly exacerbated by the varying effective Prandtl number in his models (Section 4.2).

The inset of Figure 9 shows $d_{\rm os}/H_{\rm p}$ as a function of ${\rm 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 ${\rm 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 Figure 10. The value of $\nabla-\nabla_{\rm ad}$ in the RZ is close to that of the hydrostatic solution with $\nabla T=\mbox{const.}$ in a polytropic atmosphere with adiabatic index $n=13/4$, that is, $\nabla_{\rm rad}^{\rm(hs)}-\nabla_{\rm ad}=-17/85\approx-0.165$ (see also Käpylä et al. 2019b). The transition from nearly adiabatic to the radiative gradient becomes increasingly sharper as $\mathscr{F}_{\rm 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 $\mathscr{F}_{\rm n}$ , suggesting penetration in the nomenclature of Zahn (1991).

The convection zone is characterised by $\nabla\approx\nabla_{\rm ad}$ , whereas in the radiative zone, $\nabla\approx\nabla_{\rm rad}^{\rm(hs)}$. Thus a rough estimate of the width of the transition between the zones and its dependence on $\mathscr{F}_{\rm n}$ is obtained by computing the vertical derivative of $\nabla-\nabla_{\rm 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 $\delta_{\rm trans}$ from Set K are shown as a function of $\mathscr{F}_{\rm n}$ in Figure 11. The results indicate a power-law $\delta_{\rm trans}\propto\mathscr{F}_{\rm n}^{0.30}$ for $\mathscr{F}_{\rm n}\lesssim 10^{-4}$. For higher values of $\mathscr{F}_{\rm n}$ the approximation $\nabla\approx\nabla_{\rm ad}$ is no longer accurate in the CZ and the results deviate from the general trend. Extrapolating from $\delta_{\rm trans}\approx 0.1H_{\rm p}$ for $\mathscr{F}_{\rm n}=1.9\cdot 10^{-7}$, a value of $\delta_{\rm trans}\approx 7.9\cdot 10^{-3}H_{\rm p}$ is obtained for the solar value of $\mathscr{F}_{\rm n}\approx 4\cdot 10^{-11}$. 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.

4.5 Modified heat conductivity in the radiative zone

The study of Hotta (2017) reached the lowest value of $\mathscr{F}_{\rm n}$ in the literature thus far, with $\mathscr{F}_{\rm n}=5\cdot 10^{-7}$. 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_{\rm tot}$ in the RZ and OZ (Hotta 2017), that is,

[TABLE]

The original and modified profiles of $\overline{K}$ for Runs S7 and S7m are shown in Figure 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.

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 $\mathscr{F}_{\rm 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 Figure 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.

Figure 13 also shows that the modification of the heat conductivity has serious repercussions for the overshooting depth: $d_{\rm os}$ is reduced by roughly 30 per cent from $0.33H_{\rm p}$ in Run K7 to $0.23H_{\rm 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 $\mathscr{F}_{\rm 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_{\rm os}\propto\mathscr{F}_{\rm n}^{0.08}$ dependence for $\mathscr{F}_{\rm n}\lesssim 10^{-5}$. Furthermore, $d_{\rm os}$ is consistent with a constant as a function of the Reynolds and Péclet number in the range ${\rm Re}=9\ldots 523$ in cases where ${\rm Pr}_{\rm SGS}^{(1)}=1$. A somewhat steeper power, $\mathscr{F}_{\rm 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_{\rm os}\propto\mathscr{F}_{\rm n}^{0.27}$) similar to the power laws reported earlier. Otherwise identical runs with relatively weak SGS diffusion with ${\rm Pr}_{\rm SGS}^{(1)}=5$ and 10, on the other hand, produced shallower dependences ($d_{\rm os}\propto\mathscr{F}_{\rm n}^{0.07}$, and $d_{\rm os}\propto\mathscr{F}_{\rm 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 $\mathscr{F}_{\rm 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 convection zone (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 $\mathscr{F}_{\rm 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^{3}$. 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 $\mathcal{O}(0.2),H_{\rm p}$ which is somewhat higher than the canonical estimates of $(0.05\ldots 0.1)H_{\rm 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.

Acknowledgements.

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).