© ESO, 2017

1. Introduction

Convection is a fundamental stellar process which underlies heat transport, mixing, shear, and the dynamo. During many phases of stellar evolution, a star’s interior is characterized by layers that are convectively unstable. Outside of a convection zone are layers that are stable to convection. Convective fluid motions can cross the boundary between a convectively unstable region and a neighboring stable region. When the Péclet number is large, these convective motions slowly erode the density and thermal stratification of the stable region. When this occurs, the process is termed convective penetration. We study convective motions that may contribute to convective penetration below a convection zone, and thus can have an impact on the structure and evolution of a star.

Studies of stellar evolution utilize one-dimensional calculations that evolve physical quantities as a function of the radial position interior to a star. The effects of convection on other physical quantities are typically modeled using stellar mixing length theory, which depends on the local temperature gradient (e.g., Vitense 1953; Böhm-Vitense 1958; Abbett et al. 1997; Trampedach 2010; Brandenburg 2016). The phenomenon of convective penetration has been included in mixing length theory in two competing ways, based on different non-local convection theories. One is an overshooting length, expressed as a fraction of the pressure scale height at the convective boundary, which defines a region outside of the convection zone where the stellar model assumes full convective mixing. The other is a region of enhanced one-dimensional diffusivity (Noels et al. 2010; Zhang 2013), where the intensity of convective mixing changes with radius. Using statistical observations, this work addresses both of these methods and thereby contributes to the 321D link.

Stellar convection is a nonlinear and non-local phenomenon. The extent to which convective motions can extend beyond the theoretical boundary of a convection zone is intuitively linked to the local velocity, density, and temperature of the fluid at the boundary¹. In this work we will refer to convective flow structures that possess radial velocity as plumes, without invoking geometry or physics specific to plume studies. Plumes frequently move a short distance beyond the convective boundary. Less frequently, plumes penetrate the stable zone to a larger extent. These plumes that extend deeper clearly have a greater effect on the star, because they penetrate with larger momentum and affect the stratification deeper inside the stable zone. Convective penetration is therefore naturally a process in which the intermittency is significant. A theoretical description of such flows should be built upon statistics gathered by observing convective penetration over long times. The Multidimensional Stellar Implicit Code (MUSIC), which uses implicit time integration, allows us to cover efficiently the long times necessary for converged statistics on convective penetration. This work thus builds on earlier explorations (e.g., Hurlburt et al. 1986, 1994; Brummell et al. 2002; Rogers & Glatzmaier 2005b; Rogers et al. 2006), which simulate convective overshooting or penetration over shorter times, typically between two and five convective turnover times.

In addition to observing the long-time behavior of convective penetration, a statistical treatment is necessary to model meaningfully instances of large convective penetration. Extreme value theory (EVT) is a statistical theory that has been widely used to predict events that have large impacts. Examples are record rainfall, floods, dangerously high winds, earthquakes, and avalanches. For the case of convective flows penetrating a stable zone, an analysis using extreme value theory is natural.

This work is structured as follows. In Sect. 2 we discuss the young Sun model and the numerical framework of our simulations. In Sect. 3 we examine convective penetration based on simulation data that covers a range of times, and a range of different spatial resolutions and extents. Motivated by these results, in Sect. 4 we introduce the framework of extreme value theory, and compare statistical results for simulations in a truncated domain that model penetration due to local convection, and simulations in a large domain that model penetration due to non-local convection. In Sect. 5 we build on these statistical results to formulate a one-dimensional radial diffusion coefficient that is relevant in the large Péclet number limit. In Sect. 6 we summarize and discuss the implications of these results.

2. Simulations

In this work we examine a model of a prototypical pre-main sequence star, the young Sun, using the numerical set-up benchmarked in Goffrey et al. (2017). In Pratt et al. (2016) the structure and properties of the stellar model of the young Sun are described in detail. Here we briefly summarize the salient features of our simulations. We perform two-dimensional implicit large eddy simulations (LES) of the young Sun using the MUSIC code. The stellar radius of our young Sun model is approximately three times larger than the current Sun, it is one solar mass, and has homogeneous chemical composition. The radial profiles of density and temperature for the young Sun model are typical for a pre-main sequence star that is no longer accreting and is gradually contracting. The luminosity is increasing with the interior radius. Based on the radial entropy profile and an evaluation of the Schwarzschild criterion, a central radiative zone is expected below a large convection zone. This convection zone spans the outer 1.2 × 10¹¹ cm of the total radius of 2.13 × 10¹¹ cm. A layer exists between the convection zone and the radiative zone, where mixing of convective flows into the radiative zone can take place. This layer is the focus of our work, and is located far from the physically complicated near-surface layers. In this layer, the stiffness (e.g., Hurlburt et al. 1994; Brummell et al. 2002; Rogers et al. 2006) is approximately S ≈ 13. Our simulations do not include additional stellar physics such as rotation, a tachocline, chemical mixing, or magnetic fields, that we intentionally omit from our statistical study of convective penetration. This simplifies the dynamics at the boundary between the interior radiative zone and the convection zone, where we observe penetrative convective motions. These additional effects may influence the results we obtain. In three dimensions the shape of plumes may be different from two dimensions (Schmalzl et al. 2004; van der Poel et al. 2015), and this could affect mixing properties. Rotation and shear may tilt plumes, modifying the angle of penetration (e.g. Brun et al. 2017). A magnetic field may shift the Schwarzschild criteria (Gough & Tayler 1966) and fundamentally change the mechanisms of plume braking (Hughes & Proctor 1988). These additional physical effects will be explored in future work.

The MUSIC code solves the inviscid compressible hydrodynamic equations for density ρ, momentum ρu, and internal energy ρe: using a finite-volume method, a MUSCL scheme for interpolation (Van Leer 1977), and a van Leer flux limiter (Van Leer 1974; Roe 1986). Time integration in the MUSIC code is fully implicit, and uses a Jacobian free Newton-Krylov (JFNK) solver (Knoll & Keyes 2004) with a variant of physics-based preconditioning (Viallet et al. 2016). For all simulations examined in this work, an adaptive time-step maintains a general CFL number of ten, while a CFL number based on simple advection is restricted to be 0.5. We find that this produces good convergence of relevant basic quantities, such as the average kinetic energy.

MUSIC simulations are designed to contribute to the 321D link (Arnett 2014; Arnett & Meakin 2009), the effort to improve one-dimensional stellar evolution models by studying critical short phases using stellar hydrodynamics in two and three dimensions. One aspect that makes our simulations relevant to the 321D link is the use of an equation of state and realistic opacities standard in one-dimensional stellar evolution calculations. Opacities are interpolated from the OPAL (Iglesias & Rogers 1996) and Ferguson et al. (2005) tables, which cover a range in temperature suitable for the description of the entire structure of a low-mass star. The compressible hydrodynamic equations, given by Eqs. (1)−(3), are closed by determining the pressure p(ρ,e) and temperature T(ρ,e) from a tabulated equation of state for a solar composition mixture. This equation of state accounts for partial ionization of atomic species by solving the Saha equation, and neglects partial degeneracy of electrons; it is suitable for the description of our solar model at a young age. The initial state for MUSIC simulations is produced using data extracted from a one-dimensional model calculated using the Lyon stellar evolution code (Baraffe & El Eid 1991; Baraffe et al. 1997, 1998), which uses the same opacities and equation of state implemented in MUSIC. In Eq. (2), g is the gravitational acceleration, a spherically-symmetric vector identical to that used in the stellar evolution code, and not evolved in our simulations.

The thermal diffusivity in MUSIC is also realistic for the young Sun model, and has not been enhanced. In Eq. (3)the thermal conductivity is defined using the Stefan-Boltzmann constant σ and the Rosseland mean opacity . An artificially high thermal diffusivity is a key issue that has been found to cause differences between mixing length theory and numerical results (Rempel 2004), and can result in structural changes for a star. Despite these consequences, enhancing the thermal diffusivity has become a conventional step in numerical studies of convective overshooting (Browning et al. 2004; Rogers et al. 2006; Tian et al. 2009; Strugarek et al. 2011), because it allows the star to more rapidly achieve energy balance by shortening the Kelvin-Helmholtz time scale. For the young Sun model considered here, the Kelvin-Helmholtz time scale is several orders of magnitude greater than its convective turnover time. Thus although the young Sun model has not achieved energy balance, during the time that our simulations span, there is no measured change in the thermal profile, and the convection may be considered to be in steady state.

When a realistically small thermal diffusivity is used, a high resolution is typically required to resolve the corresponding small-scale temperature fluctuations. No attempt is made to fully resolve the small scales related to the realistically small thermal diffusivity that we use; we rely on the implicit large eddy simulation paradigm to approximate the effect of this small-scale dissipation. Our use of the realistic thermal diffusivity for the young Sun model allows us to study convective penetration because it produces a large Péclet number; this is discussed fully in Sect. 2.2. In MUSIC, numerical truncation errors contribute to both the viscosity and thermal diffusivity. Because MUSIC simulations also use an explicit thermal diffusivity related to the thermal conductivity in Eq. (3), the Prandtl number is expected to be everywhere less than one. The characteristic length scales, velocities, viscosity, and thermal diffusivity vary throughout the radius of a star. The viscosity and thermal diffusivity may also be dependent on local dynamics, the wavenumber, and the estimation method (Domaradzki et al. 2003; Zhou et al. 2014; Radice et al. 2015; Schranner et al. 2015; Zhou & Thornber 2016). For these reasons, the Rayleigh number and Reynolds number are not specified in a general sense for our simulations.

2.1. Spherical-shell geometry and boundary conditions

Equations (1)−(3)are solved in a two-dimensional spherical shell using spherical coordinates: radius r and colatitude θ. A two-dimensional geometry has been chosen in order to allow us to produce a long time series of data that also has a satisfactory resolution near the bottom of the convection zone. Convective penetration is the result of large-scale convection, and our simulations do not model the significantly three-dimensional effects of turbulence, rotation, or magnetic fields present in a realistic star. Therefore a two-dimensional simulation provides an approximation to three-dimensional stellar convection. It has long been recognized that two-dimensional stellar simulations result in higher velocity magnitudes than three-dimensional simulations (Muthsam et al. 1995; Meakin & Arnett 2007). However, because our simulations do not approach the high velocities realistic to a star, the difference between two- and three-dimensional velocities do not seriously detract from a study of the form of convective penetration.

Fig. 1 Typical time snapshots of radial velocity in simulation YS2 (left), YS5 (middle), and YS6 (right). The zero point is colored black, while inward flows are gray and outward flows are red.

The placement of boundaries and the choice of boundary conditions affect the physical outcome of a hydrodynamic simulation. We use boundary conditions targeted to maintain the original radial profiles of density and temperature of the one-dimensional stellar evolution model of the young Sun. Each simulation volume begins at 20° from the north pole, and ends 20° before the south pole, allowing for a wide convection zone with multiple descending and ascending plumes. Periodicity is imposed on all physical quantities at the boundaries in θ. Radial boundaries require more sophisticated treatment. In velocity, we impose non-penetrative and stress-free boundary conditions on the radial boundaries. The energy flux and luminosity are held constant at the inner radial boundary at the value of the energy flux at that radius in the one-dimensional stellar evolution calculation. Typically we hold the energy flux and surface luminosity constant on the outer radial boundary as well; the surface luminosity of our model for the young Sun is 2.32 times the luminosity of the current Sun. An alternative is to allow the surface to radiate energy with the local surface temperature. In this case the energy flux varies as where σ is the Stefan-Boltzmann constant and T_s(θ,t) is the temperature along the surface. This boundary condition can only be effectively used when the temperature gradient near the surface is sufficiently resolved; otherwise it results in artificially high cooling rates. On the inner radial boundary of the spherical shell, a constant radial derivative on the density is imposed, as discussed in Pratt et al. (2016). At the outer radial boundary a hydrostatic equilibrium boundary condition (Grimm-Strele et al. 2015) that maintains hydrostatic equilibrium by assuming constant internal energy and constant radial acceleration due to gravity in the boundary cells is applied to the density.

Consider the suite of two-dimensional simulations of spherical-shell convection summarized in Table 1. For all of these simulations, the table lists radial grid spacing and the radial extent of the spherical shell in units of the total stellar radius R. Rogers et al. (2006) report using a grid of N_r × N_θ = 2048 × 1500, with 620 grid cells in the radiative zone to study convective overshooting in the current Sun. By comparison, in simulation YS5 we use a grid of N_r × N_θ = 2432 × 2048 with 922 grid cells in the radiative zone to study convective penetration in the young Sun. The large section of the radiative zone simulated at high resolution is useful to precisely resolve convective penetration in the young Sun, which has a larger convection zone than the current Sun.

Table 1 provides the parameters of simulations YS0-6, which consist of three different sets of simulations. Simulations YS0-2 model convection in the local approximation, using a spherical shell of limited radial extent around the penetration layer. These simulations are designed to examine penetration due to convection that arises and is dynamically limited to the layer around the bottom of the convection zone, similar to a box-in-a-star approach. Simulations YS3-5 model convection in the non-local approximation, using a spherical shell that includes nearly the whole star, terminating in the near-surface layers. These simulations are designed to examine penetration due to convection that is dynamically accurate to the large convection zone of the young Sun model, where convective dynamics may interact non-locally across more than half of the stellar radius. Among each of these sets of simulations, three different resolutions are examined for otherwise identical physical set-ups and for boundary conditions that include a constant energy flux at the outer surface. In MUSIC simulations, which do not include an explicit viscosity, higher resolution produces higher characteristic velocities in the convection zone. Thus these different resolution simulations are performed to show the robustness of our statistical results with respect to convective velocities. An additional simulation YS6 is performed in order to isolate the effect of allowing the energy flux at the stellar surface to vary according to the local temperature. A typical snapshot of the radial velocity in each of these simulation volumes is shown in Fig. 1. The effect of the placement of the simulation boundaries and the boundary conditions imposed on the properties of convection was examined in Pratt et al. (2016). In this work we examine exclusively the penetration of plumes below the convection zone.

For each simulation, Table 1 provides a time span in units of the convective turnover time of that simulation. This time span indicates a period after steady convection has developed, during which the global kinetic energy of the system is approximately constant. The average radial structure of the star’s density and temperature does not evolve over this period. All time-averages and other statistics are calculated over this time span.

2.2. Characterization based on the Péclet number

The Péclet number is a nondimensional ratio that measures the relative importance of two physical effects that determine the deceleration of a fluid element entering a layer that is stable to convection: advection and thermal diffusion. Therefore it is standardly used to distinguish between convective overshooting and convective penetration. Low Péclet numbers (Pe ≤ 1) indicate that thermal diffusion is the dominant process and convective motions overshoot the convection zone without significantly changing the stratification of the stable layer; high Péclet numbers indicate penetration can take place, a process where convective fluid motions modify the thermal and density stratification in the stable layer outside the convection zone (e.g., Zahn 1991). The Péclet number is typically defined Here the pressure scale height h_p represents a characteristic length scale, and the root-mean-square velocity v_RMS represents a characteristic convective velocity. The thermal diffusivity α is produced from our tabulated equation of state, using the definition of the thermal conductivity χ in Eq. (3), and the specific heat capacity at constant pressure c_P. Numerical truncation errors due to the finite resolution of a simulation can increase this thermal diffusivity.

In our simulations, radial profiles of temperature and thermal diffusivity are taken directly from the model of the young Sun produced by a one-dimensional stellar evolution calculation. Because of this, and because we use a realistically small thermal diffusivity, the Péclet number is much greater than one in the area of interest at the bottom of the convection zone, and convective penetration is expected. Figure 2 shows the average radial profile of the Péclet number in simulation YS4. Neglecting the contributions to thermal diffusivity from numerical truncation, the Péclet number is – in the layer at the bottom of the convection zone for all of the simulations in Table 1. By perturbing the background state we are able to measure an approximate effective thermal diffusivity. We find that while numerical truncation errors increase the effective thermal diffusivity beyond the value implied by the explicit physical thermal diffusivity, our simulations are still in the Pe ≫ 1 regime. A Péclet number larger than one is expected; Péclet numbers are generally higher in low and medium mass stars, like the young Sun model, than for higher mass stars (e.g. Meynet & Maeder 2000). Péclet numbers are also expected to be much higher in stellar interiors than in the near-surface layers. In addition to these physical points, two-dimensional hydrodynamic simulations are known for producing somewhat higher Péclet numbers because they exhibit higher velocities than three dimensional simulations (Muthsam et al. 1995; Meakin & Arnett 2007).

Fig. 2 Radial profile of the time-averaged and horizontally-volume-averaged Péclet number, as defined by Eq. (4), from simulation YS4. The contribution to the thermal diffusivity by numerical truncations are not included in this profile. The heavy vertical black line marks the boundary between the stable radiative zone and the convection zone, determined from the radial profile of entropy and the Schwarzschild criterion produced by the one-dimensional stellar model.

Based on the Péclet number, we define convective motions that cross beyond the boundary of the convection zone to be convective penetration. Indeed in test simulations where the Kelvin-Helmholtz time scale of the young Sun is artificially lowered by the use of an enhanced thermal diffusivity and corresponding increased luminosity, we observe penetrative motions that change the density and thermal stratification of the star. However in the simulations presented in this work, the Kelvin-Helmholtz time scale is many orders of magnitude larger than the characteristic time scale associated with convective motions². Therefore we do not observe erosion of the stellar structure due to convective penetration; the temperature and density profiles are approximately constant during our simulations. We do study penetrative flows that have the potential to change stellar structure. However our study may also have implications for convective overshooting motions in the low Péclet number regime.

Over the last decade it has become common to refer to the mixing at the boundary of a convection zone as turbulent entrainment (e.g. Meakin & Arnett 2007). We note that convective penetration has long been considered a type of entrainment in which mixing is due to large-scale convective motions (e.g. Tennekes & Driedonks 1981), and that in a physical star, mixing at a boundary may be due to a combination of fluid effects including convection, turbulence, and also shear (e.g. Jonker et al. 2013). Generally in LES of stellar convection, the range of scales that correspond to turbulent mixing include scales considerably smaller than the grid spacing of the simulation. Indeed for all of the simulations considered in this work, flows are expected to be in the laminar regime with moderate Reynolds numbers, and the effect of turbulent entrainment is not targeted.

3. Results: the penetration depth

3.1. Determination of the extent of convective penetration

Historically several different definitions have been used for the depth of convective penetration. Broadly speaking, a physical quantity is sought that either drops to a low level, or changes sign where convective motions cease. To measure the penetration depth using a quantity that drops to a low level, an arbitrary percentage of the quantity is chosen. An example of such a quantity is the kinetic energy density . In several earlier works, the point where convective motions cease has been defined as the point where the kinetic energy density reaches 1% (Brummell et al. 2002), or 5% (Rogers et al. 2006) of its peak value near the bottom of the convection zone.

In our simulations, difficulty with this type of measurement arises after convective plumes interact with the underlying radiative region. Convective penetration and overshooting generate a spectrum of waves in the radiative zone that propagate in the angular direction and feed back on the convective plumes. These waves have been identified as internal gravity waves in many cases (Hurlburt et al. 1986; Andersen 1994; Dintrans et al. 2003; Rogers & Glatzmaier 2005a; Rogers et al. 2006; Lecoanet & Quataert 2013; Shiode et al. 2013; Brun et al. 2013; Alvan et al. 2014, 2015; Pinçon et al. 2016). We observe waves in the radiative zone evolving over tens of convective turnover times, and they spoil an evaluation of a percentage of the peak convection zone kinetic energy density. This problem is clear from the long-time-average of the kinetic energy density shown in Fig. 3 for simulation YS3. In this long-time average, these waves appear as a small peak in the kinetic energy density at a radius 0.35 <r/R< 0.4, precisely in the region where plumes that have a large penetration depth terminate. A second issue arises when we observe steady convection over tens of convective turnover times: the peak value of the kinetic energy density near the bottom of the convection zone evolves. Thus a criterion based on a percentage of the peak value changes during the long time span of our simulations.

Fig. 3 Radial profile of the time-averaged and horizontally volume-averaged kinetic energy, , scaled to its maximum, from simulation YS3. The heavy vertical black line marks the convective boundary.

For these reasons it is advantageous to base the penetration depth on a quantity that changes sign at the point in the radiative zone where convective motions cease. Two such quantities have been proposed (as discussed in, e.g. Hurlburt et al. 1986, 1994; Ziegler & Rüdiger 2003; Rogers et al. 2006; Tian et al. 2009; Chan et al. 2010): the vertical kinetic energy flux and the vertical heat flux. We examine both of these measures, and compare the results. The vertical kinetic energy flux F_k is defined as (6)where ρ is the density, and v is the velocity. The penetration layer can be defined by the region of positive vertical kinetic energy flux. A radial profile of the vertical kinetic energy flux, calculated by averaging both in time and using a volume weighted-average in the horizontal (θ) direction, is shown in Fig. 4 for simulation YS4. Previous works have used this type of averaged profile to determine the penetration depth. In Fig. 4, the vertical kinetic energy flux exhibits a smooth positive peak; the width of the layer where convective mixing takes place is well defined, but small. In any time-snapshot, larger penetration may be identified at particular angles.

Fig. 4 Radial profile of the time-averaged and horizontally volume-averaged vertical kinetic energy flux Fk defined in Eq. (6), and vertical heat flux Fc defined in Eq. (7)for simulation YS4. Both of these fluxes have been normalized to their absolute maximum value. The heavy vertical black line marks the convective boundary.

The second criteria used to measure the penetration depth is the vertical heat flux: (7)where the temperature perturbation δT is defined relative to a time and horizontally volume-averaged temperature profile. Although the time and horizontally volume-averaged temperature profile is allowed to evolve during our simulations, it does not change substantially; this simplifies the calculation of the temperature perturbation. A layer defined by convective penetration is associated with negative vertical heat flux. Figure 4 also shows the time-averaged and horizontally volume-averaged vertical heat flux for simulation YS4. The penetration layer is clearly identifiable as the negative peak of vertical heat flux surrounding the bottom of the convection zone. There is significant discrepancy between the penetration depth defined from the averaged vertical kinetic energy flux and the averaged vertical heat flux. Because of this, we examine both measures in detail. We note that the use of Lagrangian coherent structures (e.g. Hadjighasem et al. 2016) to define the shape of convective plumes offers an alternative way to understand this discrepancy.

Fig. 5 Angular structure of the penetration layer at three arbitrary times spread throughout simulation YS5. The penetration depth in this illustration is determined by zeros of the vertical kinetic energy flux. The boundary between the convection zone and the stable radiative zone is indicated by a solid black line. The vertical axis is in units of the pressure scale height hp at this boundary. A dashed black line indicates the average penetration depth at this time.

Fig. 6 PDF of the penetration position, ro, based on plumes penetrating at all angles and sampled at a fixed time interval throughout the full time span of each simulation indicated in Table 1. a) Penetration depth calculated from the vertical kinetic energy flux for the local convection simulations YS0-2. b) Penetration depth calculated from the vertical kinetic energy flux in the non-local convection simulations YS3-5. c) Penetration depth calculated from the vertical heat flux in the local convection simulations YS0-2. d) Penetration depth calculated from the vertical heat flux in the non-local convection simulations YS3-5. The heavy vertical black line marks the convective boundary.

Horizontal and time averages have been used historically to define the width of the penetration layer. However the full statistics of convective penetration have not been previously examined, nor has the use of a simple average to represent them been specifically justified. Straightforward averages can mask contributions from the intermittent convective plumes that penetrate deeper into the stable radiative zone than the average plume. The distribution of penetration lengths may not be symmetric, or may have other relevant features that are not captured by a mean. Just such intermittent plumes can have a potentially larger impact on the physics of the penetration layer. Indeed Christensen-Dalsgaard et al. (2011) find that consideration of temporal and spatial inhomogeneity of overshooting could lead to a thermal profile that agrees with observations from helioseismology. We therefore approach the full statistics of penetrating plumes with the goal of quantifying the contributions of intermittency.

3.2. Shape of the penetration layer

Using the first zero of the vertical kinetic energy flux and the vertical heat flux in the stable zone as our criteria, we examine the depth of convective penetration for simulations YS0-5, described in Table 1. In hydrodynamic simulations without rotation, the penetration is statistically independent of the colatitude θ. However at any instant in time the depth that plumes penetrate beyond the convection zone varies in θ, defining the structure of the penetration layer. Three typical snapshots that illustrate the shape of plumes penetrating into the radiative zone are shown in Fig. 5. In this figure, the bottom of the convection zone, established by the entropy profile and the Schwarzschild criterion, is shown as a thick horizontal line. The depth of convective penetration at any angular position in the spherical shell, based on the zero point of the vertical kinetic energy flux, is shown as a shaded area below this line. At the angular resolution of simulation YS5, penetrating plumes often fill several grid spaces. A general penetration depth calculated by a straightforward average in angle is indicated by the dashed horizontal line in each case; indeed we find that in Fig. 5a, for example, an average will hide the physically important more deeply penetrating plumes at approximately 40° and approximately 150°. This figure is calculated from the vertical kinetic energy flux; a qualitatively similar picture is produced from the vertical heat flux. We see similar pictures for each of the simulations considered in this work, independent of the resolution or boundary conditions used. Comparing with three-dimensional hydrodynamic simulations during steady convection, we also observe similar pictures, although statistics based on long-times are not available for three-dimensional stellar convection.

To build a statistical model of convective penetration, we first consider the extent that plumes penetrate into the radiative zone at each angle defined in our simulation grid, and we sample the simulation data at a fixed time interval identical for each simulation. This fixed time interval is on the order of τ_conv/ 10³, a value selected to capture the fastest moving plumes entering the penetration layer. We define the penetration position r_o as the position where convective motions cease based on either of our criteria. Figure 6 shows the probability density function (PDF) of this penetration position for each set of simulations. Each PDF is calculated using bins of equal size, which are identical for each set of simulations, and is properly normalized. From these PDFs, two important observations can be made. First, the PDFs possess either a multi-modal shape or a significant shoulder. This identifies two distinct layers within the region below the bottom of the convection zone. Immediately below the convection zone is a shallow layer where weak plumes frequently penetrate the radiative zone. Below this layer is a second layer where stronger plumes penetrate, and more dramatically affect the stable radiative zone. This strong penetration occurs with a typical probability between approximately 1 / 2 and 1 / 7th of the shallow layer. The second important observation arises from the probability of larger penetration: these PDFs have a heavy tail that is non-Gaussian in appearance. Both of these observations suggest that a straightforward average does not accurately describe the data.

Unlike the averaged profiles, the PDFs obtained based on the zeroes of the vertical kinetic energy flux and the zeroes of the vertical heat flux in Fig. 6 yield characteristically similar results. The PDFs calculated from the vertical heat flux are slightly more sharply peaked for small penetration extent than those calculated from the vertical kinetic energy flux. However, the peaks and shoulders are present and their approximate positions are unchanged regardless of the physical criterium used.

Meaningful differences between the PDFs calculated for the local simulations YS0-2 and the non-local simulations YS3-5 are clear from Fig. 6. The PDF of the non-local simulations show that penetration in simulations YS3-5 reaches deeper into the radiative zone. In addition, the multi-model behavior of the PDFs and the heavy tail is more significant. This may be directly related to the higher velocities that arise for identically-resolved simulations when a large portion of the star is simulated, as reported by Pratt et al. (2016). Additional effects, such as the larger range of spatial scales of plumes that appear in the non-local simulations, may also be responsible, but are more difficult to quantify.

4. Extreme value statistics of maximal convective penetration

4.1. Formulation

The heavy tails evident in Fig. 6 motivate a more detailed statistical examination to characterize the depth of convective penetration that is meaningful for modeling stellar evolution. To pursue this, we define a maximal penetration depth to be the lowest position in the radiative zone that is reached at a given time: (8)Here we select the minimum over data at different angles, at the same time. In Fig. 7 the PDF of this maximal penetration depth is compared with the PDF of the penetration depth of all plumes (as shown in Fig. 6) and the time- and horizontally volume-averaged Péclet number. The Péclet number reaches a minimum near the bottom of the convection zone, above r/R = 0.4. Between 0.35 <r/R< 0.4, however, the Péclet number raises again; this is the signature of waves excited by convective penetration, that form in a narrow band in radius below the convection zone in the young Sun. The rms velocity profile and the rms radial velocity profile are impacted by the velocity of the waves in this radial range. The PDF of the maximal penetration depth correlates with the position where the Péclet number has a maximum due to these waves. The quantity r_max thus pinpoints an important physical consequence of convective penetration, while the PDF of the penetration depth for all plumes r_o does not.

Fig. 7 PDF of the convective penetration depth for all plumes ro (dashed line), PDF of the maximal penetration depth at any time rmax (solid line), and the time- and horizontally volume-averaged Péclet number (dotted line) for simulation YS4. Each of these quantities is vertically shifted and scaled for comparison. The convective boundary is delineated by a heavy vertical line.

Fig. 8 Cumulative distribution function F of the maximal penetration length, Δrmax defined in Eq. (9). a) Maximal penetration length calculated from the vertical kinetic energy flux for the local convection simulations YS0-2. b) Maximal penetration length calculated from the vertical kinetic energy flux for the non-local convection simulations YS3-5. c) Maximal penetration length calculated from the vertical heat flux for the local convection simulations YS0-2. d) Maximal penetration length calculated from the vertical heat flux for the non-local convection simulations YS3-5. Data is shown as point symbols, while a fit to the GEVD form is shown as a line.

Because the PDF of the maximal penetration depth targets the waves that are the most obvious dynamical effect of penetration, it is reasonable to suggest that it should determine the relevant penetration length (or overshooting length) for stellar evolution calculations. We therefore define a maximal penetration length to be the difference between the maximal penetration depth and the convective boundary: (9)Here r_B is the radial position of the convective boundary and the maximum is taken over all angles θ at a single time t. Although in this work we treat penetration below a convection zone, the absolute value allows us to refer to penetration above or below a convective boundary. This construction allows us to focus on the larger amount of convective penetration that is likely to influence stellar structure.

The maximal penetration length Δr_max could be used straightforwardly as an overshooting length in simulations of stellar evolution. A simple average of Eq. (9)produces an overshooting length ℓ_ov ~ 0.3h_p in our simulations, where h_p is the pressure scale height at the convective boundary. A calculation of an overshooting length based on a time and horizontally averaged profile of the vertical kinetic energy flux, as in Fig. 4, produces a much smaller overshooting length, ℓ_ov ~ 0.1h_p from our simulations. In recent years stellar calculations have adopted overshooting lengths in a wide range 0.05h_p ≤ ℓ_ov ≤ 0.4h_p (Basu 1997; Chen & Han 2002; Brummell et al. 2002; Rogers et al. 2006; Politano et al. 2010; Liu et al. 2012; Montalbán et al. 2013; Jin et al. 2015), which includes both of these values. Both of these values are consistent with previously reported results for simulations with a similar stiffness to the young Sun (Brummell et al. 2002; Rogers et al. 2006). Those earlier results demonstrate that the overshooting length is dependent on the Prandtl number, the Rayleigh number, and the resolution of the simulation, as well as the stiffness. However our analysis suggests that larger values of the overshooting length than previously predicted would better reproduce the physics of large Péclet number convective penetration in the interior of a young low-mass star.

4.2. Application of extreme value theory

The maximal penetration length in Eq. (9)suggests concepts from extreme value theory (Castillo et al. 2005). In extreme value theory, unlike the central limit theorem that underlies much of probability theory, the interest is not the central values of an underlying distribution, but in accurately modeling the tails. The centerpiece of extreme value theory is the generalized extreme value distribution (GEVD), which can be used to model the probability of maximal events. The cumulative distribution function (CDF) identified with the GEVD has the form: (10)Here κ is commonly called the shape parameter, μ is the location parameter, and λ is the scale parameter (see e.g. Sect. 9.1.1 Eqs. (9.3) and (9.4) of Castillo et al. 2005; Charras-Garrido & Lezaud 2013; Gomes & Guillou 2015). Historically the well-documented Frechét distribution corresponds to a positive shape parameter κ> 0, the Weibull distribution corresponds to κ< 0, and the Gumbel distribution corresponds to κ = 0. The shape parameter is related directly to the heaviness of the tail of the underlying distribution. For this reason κ is also sometimes called the extreme value index (EVI). For the Gumbel case of κ = 0 the GEVD may be expressed in the simpler form: (11)Thus the graph of the natural log of the negative natural log of F will be a straight line if the distribution is of the Gumbel form. It will curve downward if the distribution is of the Weibull form, and upward if the distribution is of the Frechét form. When a strict upper limit exists for the variable considered, the Weibull case is predicted. Our maximal penetration length is defined so that 0 < Δr_max< 0.43R, meaning that convective plumes cannot penetrate beyond the center of the star. Thus a Weibull distribution is expected.

The CDFs of the maximum penetration length for our simulations are shown in Fig. 8. The forms of the CDFs in Fig. 8 indeed appear to have slight downward curvature, consistent with the Weibull form. For all of these CDFs, most data lie at lower values of Δr_max, and thus the distribution is clearest in this region; a lack of data at high values of Δr_max, due to the comparative rarity of highly penetrating plumes, may contribute to the downward curvature in this region. However simulations with different resolution produce remarkably similar CDFs, although different resolutions produce different convective velocities, and different spans of time are examined. Increased resolution of a simulation leads to an improved characterization of the physical flows. In a direct sense, it also leads to a larger sample size over which our maximum is taken to define the length Δr_max. However in a fluid simulation the data sampled is not strictly independent, so that it is not clear to what extent a higher resolution is related to an improvement in the use of extreme value theory, where taking a maximum over a large sample of independent data is assumed. We do find that identical distributions are produced when our data is sampled ten times less frequently in time, on the order of the autocorrelation time of the flow at the bottom of the convection zone.

To precisely determine the parameters of the GEVD for each of our simulations, we use the package evd (Stephenson 2002; Penalva et al. 2013) publicly available for R (The R Project for Statistical Computing)³. Independently, we confirm these parameters using a nonlinear least squares fit to the GEVD. Although in our fluid simulations the data is not strictly independent, as assumed by extreme value theory, we find that the GEVD provides an excellent fit. The shape, location, and scale parameters resulting from our fit are summarized in Table 2. In the table a value of κ = 0 is indicated when the best fit is of the Gumbel form. The value of the shape parameter κ is either small and negative, or zero for all of our simulations. The available time-series of data for our highest resolution simulations is shorter than for lower resolutions. Despite this, a statistically significant trend is observed in which the value of κ is explicitly zero in the high resolution, high velocity simulations YS4 and YS5. Based on this observation, we predict that the Weibull characterization of our distributions should converge toward the simpler Gumbel distribution at high resolution, realistic stellar velocities, and when stellar convection is observed over long times.

4.3. Radiation with the local surface temperature

To understand whether the treatment of the surface produces different results for convective penetration, we perform an additional simulation YS6. This simulation has identical resolution to simulations YS0 (local convection) and YS3 (non-local convection), but in simulation YS6 the near-surface layers are included and the surface of the star is allowed to radiate energy with the local surface temperature. The PDF of the penetration depth of all plumes, r_o, for YS6 is highly similar to simulation YS3; this comparison is shown in Fig. 9. The penetration depth calculated from the vertical kinetic energy flux in Fig. 9a reveals a slightly larger probability of penetration between 0.39 <r< 0.42 than in simulation YS3. The penetration depth calculated from the vertical heat flux in Fig. 9b reveals a slightly larger probability of penetration between 0.42 <r< 0.43 than in simulation YS3. Although each of these simulations is observed over more than a hundred convective turnover times, the small difference in the statistics may still result from the limited length of time that these simulations are observed, particularly because convective penetration is intermittent.

Fig. 9 PDF of the penetration extent of all plumes, ro. a) Penetration depth calculated from the vertical kinetic energy flux. b) Penetration depth calculated from the vertical heat flux. The convective boundary is delineated by a heavy vertical line.

The parameters for the GEVD of the maximal penetration length in simulation YS6 are included in Table 2. For both the vertical kinetic energy flux and the vertical heat flux criteria, the GEVD for simulation YS6 is best fit by a Gumbel distribution. The treatment of the surface in simulation YS6 is physically more accurate than in YS3. That these two different treatments of the surface radiation produce similar results may be due to the fact that both treat the surface only approximately and neglect important physical effects. However the variable surface radiation in YS6 leads to higher velocity throughout the convection zone compared with either simulation YS3 or YS0; this was observed in Pratt et al. (2016). This is a significant observation, because it reinforces the finding that the overshooting length is not solely linked to the velocity of structures entering the penetration layer (Schmitt et al. 1984; Hurlburt et al. 1994; Saikia et al. 2000; Zahn 2000).

5. Determination of a diffusion coefficient for convection

5.1. Background

Convective overshooting and penetration have long been modeled by analyzing the enhanced diffusion in the penetration layer due to intermittent convective flows (e.g. as described by Freytag et al. 2010; Noels et al. 2010; Zhang 2013). Such models are based on the definition of a one-dimensional diffusion coefficient, as in (12)Here A is any scalar quantity, and D is a diffusion coefficient. Based on dimensional arguments, a diffusion coefficient has typically been estimated for astrophysical applications by either a characteristic length scale multiplied by a characteristic velocity scale (e.g. Van Ballegooijen 1982; Andrássy & Spruit 2015), or a characteristic squared velocity multiplied by a characteristic time scale.

Equation (12)typically describes processes like molecular diffusion, and does not explicitly account for additional mixing due to large-scale convective motions, which enhance diffusion. In the setting of steady two-dimensional large-scale convection⁴, the relation between the diffusion coefficient and the Péclet number of a convective flow has been studied in both the small and large Péclet number regimes (Moffatt 1983; Rosenbluth et al. 1987; Shraiman 1987; Haynes & Vanneste 2014). In the small Péclet number regime this relation has been analytically and numerically shown to be: (13)Here D is the enhanced diffusion coefficient and D₀ is the diffusion coefficient based on small-scale diffusive processes, which may be molecular diffusion in a liquid, or small-scale turbulent diffusion in a fluid. The constant c₁ is related to the aspect ratio of the convection rolls. In the large Péclet number regime the corresponding relation is: (14)where c₂ is a constant based on the aspect ratio of the convection cells. Rosenbluth et al. (1987) evaluate c₂ for a wide range of possible convection roll aspect ratios, and find that with only small variation. The enhanced diffusion coefficients in Eqs. (13)and (14)are derived for an ideal setting rather than stellar convection, but are suggestive of diffusive scalings relevant to different Péclet number regimes.

Near the surface of a star, the small Péclet number limit is most appropriate; the Péclet number in simulations of this layer has been reported to be of order one. In this setting, Eq. (13)indicates that an enhancement of the diffusion coefficient dependent on the squared Péclet number is reasonable. If the thermal diffusivity and characteristic length scale are assumed to be approximately constant, the Péclet number is proportional to the velocity. This situation can be related directly to the enhanced diffusivity proposed by Freytag et al. (1996). In their radiation hydrodynamics simulations, Freytag et al. (1996) observe exponential decay of the rms radial velocity beyond the boundary of the convection zone. They use this result to define a piecewise-continuous diffusion coefficient that is constant in the convection zone. Beyond the convection boundary, the piecewise-continuous diffusion coefficient is defined using the Ansatz that it should decay exponentially as the velocity: (15)Here v_B and r_B are the rms radial velocity at the convective boundary, and the position of the convective boundary, respectively. The pressure scale-height h_p at the convective boundary is used to define the characteristic time scale τ = h_p/v_B and the velocity scale-height at the convective boundary h_v = fh_p, where f is a constant of proportionality. The exponentially decaying diffusion coefficient of Eq. (15)has been incorporated in several stellar evolution codes, including MESA (Paxton et al. 2013), ATON 3.1 (Ventura et al. 2008), and GARSTEC (Weiss & Schlattl 2008).

5.2. Diffusion coefficient in the large Péclet number regime

Fig. 10 Time-averaged and horizontally volume-averaged rms radial velocity near the bottom of the convection zone in simulation YS4. The time-averaged and horizontally volume-averaged square-root of the Péclet number is shown, scaled and shifted, for comparison. The heavy vertical black line marks the convective boundary.

In contrast with the early simulations examined by Freytag et al. (1996), in our comparatively high-resolution young Sun simulations the lower convective boundary is in the large Péclet number regime. For this physically and numerically distinct case, we observe that the rms radial velocity does not decay exponentially beyond the convective boundary. The decay of the natural log of rms radial velocity is shown in Fig. 10, and it does not decay linearly with radius in the penetration layer. Motivated by Eq. (14), we examine the square-root of the Péclet number in this figure in addition to the rms radial velocity, and find that predictably it exhibits a trend of decay outside the convection zone similar to the rms radial velocity. The observed departure from exponential decay motivates our proposal of a new form for the diffusion coefficient for stellar interiors in the large Péclet number regime.

Fig. 11 Time and horizontally volume-averaged square-root of the Péclet number compared with the data for the cumulative distribution function F(rmax/R) = 1 − F(Δrmax/R) calculated from the vertical kinetic energy flux. The best fit to this data from simulation YS4 is a Gumbel distribution with the parameters indicated in Table 2. The Péclet number has been normalized and shifted vertically. The dashed line indicates a Péclet number that is also shifted in radius so that the common shape is clear. The heavy vertical black line marks the convective boundary.

Based on our statistical analysis in Sect. 4, we predict that the maximal penetration length should match the profile of the Gumbel distribution: exp( − exp( − x)). Figure 11 shows that the radial profile of the square-root of the Péclet number closely matches the form of the CDF of maximal penetration length. We propose a diffusion coefficient of the form: (16)Here penetration beneath a convection zone is assumed, and Pe_B is a characteristic Péclet number in the convection zone. In constructing this diffusion coefficient, we have used the identity: (17)to relate the cumulative distribution functions of the maximal penetration depth r_max and the maximal penetration length Δr_max. This identity allows us to express the diffusion coefficient in the same form as Freytag et al. (1996), as a difference between the radial coordinate and the position of the convective boundary indicated by r_B. The prefactor D₀Pe_B^{1 / 2} in Eq. (16)is representative of diffusion in the convection zone, and may be estimated by comparison with stellar evolution models for convective diffusion. Using mixing length theory, diffusion in the convection zone is characterized by D_MLT = 1 / 3 L_MLTv_MLT, where the length scale is proportional to the pressure scale height⁵. The values for μ and λ in Eq. (16)may be estimated from simulations YS3-5 in Table 2. In units of the pressure scale height at the base of the convection zone, h_p ≈ 0.178 R, these parameters give a range of We note that μ is approximately equivalent to the maximal penetration length estimated in Sect. 4. As with D_F96, the diffusion coefficient in Eq. (16)should be applied only in the layer of the star affected by convective penetration, so that it does not introduce infinitesimally small amounts of mixing into deeper layers of the star. This proposed diffusion coefficient is shown in Fig. 12, calibrated using the mixing-length-theory diffusion coefficient, D_MLT, used to produce the young Sun model. Both D_MLT, and D_F96 are shown for comparison. The mixing-length-theory diffusion coefficient does not include any diffusion in the penetration layer. A detailed comparison of the exponentially decaying diffusion coefficient D_F96 with our diffusion coefficient in Eq. (16)naturally depends on the parameters that define the exponential decay. These are commonly defined by treating the constant of proportionality between the velocity scale-height and pressure scale-height as a free parameter f ≤ 0.1 (e.g. Angelou et al. 2015). For these small values of f our diffusion coefficient models a higher level of diffusion throughout the penetration layer than the exponentially decaying diffusion coefficient of D_F96. Because of the higher level of diffusion predicted, and the different form of its decay below the convective boundary, it is worth exploring the effect of this new diffusion coefficient on the abundance of light elements like lithium and beryllium in solar-type stars (Boesgaard 1976; Pinsonneault 1994), and to characterize its observational signatures. In simulations that include magnetic fields, this differing amount of diffusion due to convective penetration may affect magnetic flux in the penetration layer, and thus have consequences for the stellar dynamo (e.g. Nordlund et al. 1992; Rüdiger & Brandenburg 1995; Tobias et al. 2001; Rempel 2003; Brummell et al. 2010).

Fig. 12 Diffusion coefficient based on EVT proposed in Eq. (16)compared with a diffusion coefficient defined from mixing length theory DMLT, and the decaying exponential diffusion coefficient DF96 of Freytag et al. (1996) defined in Eq. (15). Here the peak of the mixing length theory diffusion coefficient has been used to calibrate diffusion in the convection zone.

6. Summary and discussion

We have studied flows that penetrate below a convection zone due to large-scale two-dimensional stellar convection in a prototypical young low-mass star. We have analyzed and compared convective penetration in two fundamental types of simulations: (1) those that model local convection by truncating the spherical shell to a minimal radial layer around the convective boundary, and (2) those that model non-local convection by simulating nearly the whole star. For each of these situations we examine stellar convection at three different resolutions, allowing us to observe penetration relating to a range of characteristic convective velocities. To compare with our simulations that hold the energy flux constant at the outer boundary, we also consider a non-local convection simulation where the surface radiation is allowed to vary with the surface temperature. This simulation is characterized by higher velocities in the convection zone. Our simulations do not account for rotation, shear flows, chemical mixing, or magnetic fields. These processes may change the properties of the flow at the convective boundary and alter the quantitative results. Analysis of these effects will be explored in future work.

Historically, a time and horizontally volume-averaged radial profile of the vertical kinetic energy flux has often been used to determine the depth of the penetration layer. Different results are produced from a comparison of the time and horizontally volume-averaged radial profiles of the vertical kinetic energy flux with the vertical heat flux, a second way to define the penetration layer. Here we show that when plumes at all angles and sampled at a fixed time interval are considered, PDFs of convective penetration depth calculated from these two different fluxes are characteristically similar. The PDFs are non-Gaussian, with a heavy tail. Based on these PDFs, we define two distinct layers that form between the convection zone and the stable radiative zone: a shallow layer where convective plumes frequently penetrate, and a deeper layer where convective plumes penetrate intermittently.

Motivated by the non-Gaussian character of convective penetration, we examine the maximal penetration depth at any time. The PDF of this quantity is centered closely around the position where waves are excited by convective penetration in the radiative zone. This suggests that the maximal penetration depth is a physically significant length scale to define the penetration layer. Based on a statistical analysis of our simulation data, we determine that a reasonable estimate of the overshooting length in the young Sun is ℓ_ov ~ 0.3h_p. We observe that the cumulative distribution function of the maximal penetration length can be modeled accurately by a Weibull distribution with a small shape parameter. As both higher resolution and higher velocities are examined, this Weibull distribution appears to converge toward a Gumbel distribution. This approach toward measuring convective penetration is new. Further work is required to demonstrate that it accurately describes convective penetration when rotation, shear flows or magnetic fields are considered. However, it provides a promising avenue to quantify the mixing due to overshooting or penetration in stars. Analysis of the properties of the waves generated by convective penetration into the radiative zone will be explored in future work.

Building on these statistical results, as well as scalings for two-dimensional convection, we propose a new form for the convection-enhanced diffusion coefficient suitable for one-dimensional stellar evolution calculations. Unlike the D_F96 diffusion coefficient (Freytag et al. 1996), our proposed diffusion coefficient is targeted for the large Péclet number flows characteristic of stellar interiors. This represents a step forward for stellar evolution modeling, because the D_F96 diffusion coefficient has been broadly applied to both small and large Péclet number flows, outside of its range of demonstrated physical validity. Follow-up studies to explore the effect of our diffusion coefficient in one-dimensional stellar evolution calculations are now required to analyze possible differences that may result (Baraffe et al., in prep.). Isolating observable signatures to discern the structure of overshooting and penetration layers is currently a challenge, but may also shed light on the differences between these two forms of the diffusion coefficient.

In this work we have focused on identifying the amount of mixing due to convective penetration, which can be modeled by a simple diffusion coefficient. However convective penetration also has a significant impact on the thermal profile of a penetration layer. A complete model of penetration therefore should include a treatment of heat transport relevant for stellar evolution calculations. Such a model is currently being pursued based on the simulations presented.

Acknowledgments

The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework (FP7/2007-2013)/ERC grant agreement No. 320478.

I.B. thanks MOCA at Monash University for warm hospitality.

This work used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National E-Infrastructure.

This work also used the University of Exeter Supercomputers, ZEN and ISCA, funded by the STFC/DiRAC, the Large Facilities Capital Fund of BIS and the University of Exeter.

References

All Tables

All Figures

	Fig. 1 Typical time snapshots of radial velocity in simulation YS2 (left), YS5 (middle), and YS6 (right). The zero point is colored black, while inward flows are gray and outward flows are red.
In the text

Fig. 2 Radial profile of the time-averaged and horizontally-volume-averaged Péclet number, as defined by Eq. (4), from simulation YS4. The contribution to the thermal diffusivity by numerical truncations are not included in this profile. The heavy vertical black line marks the boundary between the stable radiative zone and the convection zone, determined from the radial profile of entropy and the Schwarzschild criterion produced by the one-dimensional stellar model.

In the text

	Fig. 3 Radial profile of the time-averaged and horizontally volume-averaged kinetic energy, , scaled to its maximum, from simulation YS3. The heavy vertical black line marks the convective boundary.
In the text

	Fig. 4 Radial profile of the time-averaged and horizontally volume-averaged vertical kinetic energy flux Fk defined in Eq. (6), and vertical heat flux Fc defined in Eq. (7)for simulation YS4. Both of these fluxes have been normalized to their absolute maximum value. The heavy vertical black line marks the convective boundary.
In the text

Fig. 5 Angular structure of the penetration layer at three arbitrary times spread throughout simulation YS5. The penetration depth in this illustration is determined by zeros of the vertical kinetic energy flux. The boundary between the convection zone and the stable radiative zone is indicated by a solid black line. The vertical axis is in units of the pressure scale height hp at this boundary. A dashed black line indicates the average penetration depth at this time.

In the text

Fig. 6 PDF of the penetration position, ro, based on plumes penetrating at all angles and sampled at a fixed time interval throughout the full time span of each simulation indicated in Table 1. a) Penetration depth calculated from the vertical kinetic energy flux for the local convection simulations YS0-2. b) Penetration depth calculated from the vertical kinetic energy flux in the non-local convection simulations YS3-5. c) Penetration depth calculated from the vertical heat flux in the local convection simulations YS0-2. d) Penetration depth calculated from the vertical heat flux in the non-local convection simulations YS3-5. The heavy vertical black line marks the convective boundary.

In the text

	Fig. 7 PDF of the convective penetration depth for all plumes ro (dashed line), PDF of the maximal penetration depth at any time rmax (solid line), and the time- and horizontally volume-averaged Péclet number (dotted line) for simulation YS4. Each of these quantities is vertically shifted and scaled for comparison. The convective boundary is delineated by a heavy vertical line.
In the text

Fig. 8 Cumulative distribution function F of the maximal penetration length, Δrmax defined in Eq. (9). a) Maximal penetration length calculated from the vertical kinetic energy flux for the local convection simulations YS0-2. b) Maximal penetration length calculated from the vertical kinetic energy flux for the non-local convection simulations YS3-5. c) Maximal penetration length calculated from the vertical heat flux for the local convection simulations YS0-2. d) Maximal penetration length calculated from the vertical heat flux for the non-local convection simulations YS3-5. Data is shown as point symbols, while a fit to the GEVD form is shown as a line.

In the text

	Fig. 9 PDF of the penetration extent of all plumes, ro. a) Penetration depth calculated from the vertical kinetic energy flux. b) Penetration depth calculated from the vertical heat flux. The convective boundary is delineated by a heavy vertical line.
In the text

	Fig. 10 Time-averaged and horizontally volume-averaged rms radial velocity near the bottom of the convection zone in simulation YS4. The time-averaged and horizontally volume-averaged square-root of the Péclet number is shown, scaled and shifted, for comparison. The heavy vertical black line marks the convective boundary.
In the text

Fig. 11 Time and horizontally volume-averaged square-root of the Péclet number compared with the data for the cumulative distribution function F(rmax/R) = 1 − F(Δrmax/R) calculated from the vertical kinetic energy flux. The best fit to this data from simulation YS4 is a Gumbel distribution with the parameters indicated in Table 2. The Péclet number has been normalized and shifted vertically. The dashed line indicates a Péclet number that is also shifted in radius so that the common shape is clear. The heavy vertical black line marks the convective boundary.

In the text

	Fig. 12 Diffusion coefficient based on EVT proposed in Eq. (16)compared with a diffusion coefficient defined from mixing length theory DMLT, and the decaying exponential diffusion coefficient DF96 of Freytag et al. (1996) defined in Eq. (15). Here the peak of the mixing length theory diffusion coefficient has been used to calibrate diffusion in the convection zone.
In the text