DYNAMICS OF TURBULENT CONVECTION AND CONVECTIVE OVERSHOOT IN A MODERATE-MASS STAR

Until recently only one-dimensional (1D) stellar models based on the mixing-length theory (MLT) were available. However, rapidly growing computational capabilities have enabled three-dimensional (3D) stellar dynamics simulations on both global (e.g Charbonneau 2013; Guerrero et al. 2013; Featherstone & Miesch 2015) and local scales (e.g., Trampedach 2010; Beeck et al. 2012; Kitiashvili et al. 2012). The global simulations aim to model the dynamics of the entire star, including differential rotation and dynamo activity.

The overshoot region at the bottom of the convection zone plays a particularly important role in stellar dynamics and magnetism. Helioseismology has found that this region coincides with the tachocline, a zone of strong rotational shear (Kosovichev 1996; Elliott & Gough 1999) that is considered to be a primary source of the solar dynamo. Our basic knowledge about solar tachocline properties comes from helioseismic inversions of p-mode frequencies of medium angular degree (e.g., Kosovichev 1997), which find a pronounced peak in the sound-speed deviations from the standard evolutionary model at the bottom of the convection zone (Christensen-Dalsgaard et al. 1996) as well as deviations in other properties, such as density, entropy gradient, and the adiabatic exponent (Kosovichev 1999). Attempts to explain these deviations in terms of rotationally induced mixing, variations of opacity, and 1D "semi-empirical" overshooting models have not been successful (for recent results and references, see the paper of Christensen-Dalsgaard et al. 2011).

Convective overshooting has been intensively studied by using two-dimensional (2D) and 3D numerical simulations. In particular, Brummell et al. (2002) performed a detailed parametric study for an idealized case of two polytropic layers and found that the overshooting region is mostly subadiabatic with a smooth transition between the convective and radiative zones. Rogers et al. (2006) performed simulations of convective overshoot in a 2D model of the solar equatorial plane using an anelastic approximation and found that convective penetration leads to a slightly extended, mildly subadiabatic temperature gradient beneath the convection zone. A slightly higher temperature is maintained in the overshoot region by adiabatic heating from overshooting plumes, suggesting that this enhanced temperature may partially account for the sound-speed discrepancy between the standard solar model and helioseismology. The anelastic approximation has been recently used in 3D simulations of Brun et al. (2011) and Guerrero et al. (2013), who studied the dynamical coupling and formation of a tachocline in the overshoot region. However, the anelastic simulations do not include the surface layers where the convection is driven by radiative cooling and is strongly compressible. Currently, fully compressible radiative hydrodynamic simulations of the whole solar convection zone cannot be performed because of its large depth and long relaxation time. However, as suggested by Freytag et al. (1996), such simulations can be performed for more massive F- and A-type stars with much shallower outer convection zones and could shed more light on structure and properties of the stellar and solar overshooting zones.

In this Letter, we present results of initial 3D radiative hydrodynamic simulations for a star with mass 1.47 M_⊙, including the whole convection zone, the overshoot region, and the outer part of the radiative zone. We discuss the physical properties of the turbulent convection and the structure and dynamics of the overshoot region. Comparison with the corresponding 1D standard stellar model reveals a striking similarity to the deviations of the solar stratification from the standard solar model previously found using helioseismology. While this comparison is only qualitative, it reveals important features of convective overshooting that were missing in previous studies.

The stellar convection simulations are performed using the 3D radiative MHD code "StellarBox" (Wray et al. 2015) that was specially developed for realistic modeling of solar and stellar convection. The simulations are large-eddy simulations (LES) that use the Smagorinsky model (Smagorinsky 1963; Moin et al. 1991) to describe subgrid-scale turbulence (which also supplies diffusion necessary to stabilize the numerical scheme for solving the PDEs). Radiative transfer between fluid elements is calculated using a 3D, multi-spectral-bin, long-characteristic method, in which the ray-tracing transport is calculated with the Feautrier (1964) method for 14 rays. The opacity wavelength-dependent code provided by the Opacity Project (Seaton 1995; Badnell et al. 2005; http://cdsweb.u-strasbg.fr/topbase/TheOP.html) and the OPAL equation-of-state tables were used assuming a solar-type composition: X = 0.702, Y = 0.28, Z = 0.018 (see Iglesias & Rogers 1996; Rogers et al. 1996 for details).

The initial conditions are calculated using the stellar structure model obtained with the 1D stellar evolution code CESAM (Morel 1997; Morel & Lebreton 2008) using the same composition and OPAL tables and a standard value of the mixing-length parameter α =1.8. The turbulent pressure and overshoot parameters are set to zero. With increasing stellar mass the thickness of the convection zone decreases. This allows us to investigate the dynamics of the entire convection zone of a moderate-mass star. For this study, we selected a star with the following properties: $M=1.47\quad {M}_{\odot }$, ${T}_{\mathrm{eff}}=7063$ K, g = 10^4.279 cm s⁻², L/L_⊙ = 4.731, $R=1.456\quad {R}_{\odot }$. The depth of the convection zone is z_cz =  −28.5 Mm. The 1D stellar model was extended by adding an isothermal 1 Mm atmospheric layer. This allowed us to accurately model the radiative cooling of the surface layers by keeping the optical depth at the upper boundary at about 10⁻², despite the expansion of the upper convection zone caused by the turbulent dynamics.

The simulation domain used in the current simulations is 102.4 × 102.4 × 51 Mm in size (1024 × 1024 × 512 grid points). The computational domain covers 1 Mm of the atmosphere and 50.5 Mm of the stellar interior. This includes the entire convection zone and the upper part of the radiation zone. The grid spacing is 100 km in the horizontal directions and variable in the vertical direction, from 28.6 km at the photosphere and above to 183 km at the bottom boundary of the domain. The boundary conditions are periodic in the horizontal directions. The bottom boundary is closed for flows and open for the radiative energy flux. To suppress vorticity generation near the bottom boundary, we increased the dissipation by imposing a viscous layer (occupying four vertical grid points) near the bottom boundary. The upper boundary condition is formulated in terms of the Riemann characteristics and is open for mass, momentum, and energy fluxes and also for the radiation flux. For other computational details and performance of the "StellarBox," see Wray et al. (2015).

The stellar photosphere is covered by cells, that, like granules on the Sun, are characterized by upflows in the middle and narrow downflows in the surrounding lanes (Figure 1(a) and corresponding animation). Large-scale fluctuations in this animation reflect propagation of acoustic waves excited by the turbulent convection (e.g., Stein & Nordlund 2001; Jacoutot et al. 2008; Kitiashvili et al. 2011a, 2011b). The size of stellar granules varies from ∼2 Mm to ∼12 Mm, and clusters of several granules (resembling mesogranulation) are evident. Most granules of an intermediate size are stretched by larger-scale upflows. The multi-scale property of the convection vanishes in subsurface layers, where the scale of the convective eddies increases with depth (Figures 1(a)–(f)) and reaches a size of 40–45 Mm at about 25 Mm depth (near the bottom of the convection zone). Deeper down, the influence of the overshoot region becomes more and more noticeable as prominent fine structures of upflows are formed around edges of strong downflows. Figure 1(g) illustrates a horizontal slice of the overshoot region at 30 Mm depth, where high-speed downflows (Figure 2(a)) penetrate through the whole convection zone and hit the radiative zone, causing splashes in the horizontal directions and generating numerous small-scale turbulent helical structures and upflows around the plumes. We define the overshoot region as a zone of penetration of convective flows beneath the adiabatic boundary at z_cz (see Sections 4 and 5 for more detail). Below the overshoot region, at a depth of 40 Mm, in the radiative zone, which is a convectively stable layer the velocity perturbations quickly decrease, and the dynamics is dominated by internal gravity waves (Figure 1(h)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our simulation results show that narrow high-speed downdrafts (Figure 2(a)) driven by surface radiative cooling can often penetrate through the whole convection zone, but they quickly stop in the overshoot layer. The downflows are accelerated up to 21 km s⁻¹ and initiate strong perturbations in the radiative zone. Such perturbations are clearly visible in the density, temperature, and sound-speed fluctuations (Figures 2(b)–(c)). The downdrafts are associated with positive density perturbations and negative temperature perturbations. However, the penetrating downflows locally heat the overshoot region. The downdraft impacts also lead to the excitation of internal gravity waves (g-modes) in the radiative interior.

To summarize the dynamical structuring from the stellar photosphere down to the radiative zone, it is convenient to consider the mean vertical profiles of the fluctuations for various turbulence quantities and the convective and kinetic energy flux (Figure 3). The rms velocity profiles show the relative strength of the vertical and horizontal flows. The horizontal flows gradually decrease in strength from the photosphere to the bottom of the convection zone. In the overshoot region, the amplitude of the horizontal velocities increases slightly due to splashing of downdrafts.

Zoom In Zoom Out Reset image size

The rms vertical velocity profile shows a strong deviation from the horizontal flows near the stellar surface layer. In particular, an additional broad bump at about 5 Mm below the surface corresponds to one of the characteristic scales of the granulation layer. In deeper layers, the vertical velocity gradually decreases. It is interesting to note that, in the range of depths of 6–10 Mm, the horizontal and vertical flows behave very similarly. Below 10 Mm, the fluctuations of the vertical velocities are stronger than the horizontal ones, probably reflecting the influence of the penetrating high-speed downdrafts. In the overshoot layer, the vertical velocity sharply decreases.

The distribution of temperature fluctuations with depth (Figure 3(b), black curve) shows the strongest variations near the surface layers and in the overshoot region. Despite the similarity of the fluctuation amplitudes, their nature is different. The near-surface fluctuations (peaked at ∼5 Mm) are mostly related to strong radiative cooling in the intergranular lanes and are associated with downflows that determine the granulation scales. A sharp increase in temperature fluctuations between −10 and −5 Mm corresponds to the He ii ionization zone. In this region the rms of vertical velocity (Figure 3(a)) and the convection energy flux (Figure 3(f)) also increase, indicating an enhancement of turbulent convection in the ionization zone, perhaps confirming theoretical predictions of Rast & Toomre (1993).

The strong fluctuations in temperature at z ≈ −34 Mm, which are associated with the overshoot region, reflect the strong local heating in these layers due to plasma compression by downdrafts that penetrate through the convection zone. The rms sound-speed perturbations (Figure 3(b), blue curve) have a similar shape, but their amplitude in the overshoot region is much smaller than in the subsurface layers. The mean profiles of the enstrophy and kinetic helicity density (Figure 3(c)) gradually decrease with depth through the convection zone, and then sharply decrease at the bottom of the overshoot region. These quantities are important for the dynamo action.

The vertical profiles for the rms fluctuations of density and gas pressure (Figure 3(d)) are qualitatively different. The density fluctuations show variations that are qualitatively similar to the temperature fluctuations, whereas the gas pressure fluctuations gradually increase from the surface to about 20 Mm below the surface, and then remain almost constant in the deeper layers. However, the total turbulent pressure, ${P}_{{\rm{turb}}}$, displays a prominent peak in the upper half of the overshoot region (Figure 3(e)), while the turbulent pressure of vertical flows, ${P}_{\mathrm{turb},z}$, gradually decreases. The convective enthalpy flux, ${F}_{{\rm{conv}}}$, calculated following the formulation of Nordlund & Stein (2001, Equations (22)–(23)), is negative in the overshoot region as suggested by non-local mixing-length theories (e.g., Shaviv & Salpeter 1973).

In 1D stellar evolution models, the convection zone structure is calculated using MLT (Böhm-Vitense 1958), which describes convection as buoyant bubbles traveling in mixing-length steps from deep convective layers to the surface. However, numerical simulations have shown that a significant role in solar and stellar convection is played by high-speed downdrafts driven by surface cooling. Therefore, it is not surprising that there are significant differences between the MLT-based models and the corresponding mean properties of the 3D simulations of turbulent convection.

In the case of a star with mass 1.47 M_⊙, we can investigate the differences between the 1D MLT model and 3D LES simulations through the whole convection zone. Figure 4 shows the mean profiles of temperature T, adiabatic exponent $\gamma \ ={(\partial \mathrm{log}P/\partial \mathrm{log}\rho )}_{{\rm{ad}}}$, temperature gradient ${\rm{\nabla }}=d\mathrm{log}T/d\mathrm{log}P$, and the Ledoux parameter of convective stability, ${A}^{*}=1/\gamma (d\mathrm{log}P/d\mathrm{log}r)-d\mathrm{log}\rho /d\mathrm{log}r$, where P is the gas pressure, ρ the mass density, and r the radius, for the 1D model and 3D simulations.

Zoom In Zoom Out Reset image size

In the 3D stellar simulations, the convection zone expands in both directions: about 800 km up into the atmosphere with a corresponding increase of the stellar photospheric radius. Such an increase of the stellar radius due to turbulent dynamics may be relevant to resolving the discrepancy between the solar seismic and photospheric radii (Schou et al. 1997). In the 3D simulations, the mean profile of the temperature is slightly higher than in the 1D model, but the deviation is strongest in the near-surface layers where the effects of radiative losses and turbulent dynamics are strong. At about 10 Mm, where the He ii ionization zone starts, the temperature discrepancies are decreasing and become smallest at the lower boundary of the zone. In the deeper layers, the difference in the temperature profiles increases again and is highest in the overshoot region. Our simulation results show significant changes in the value of γ for the hydrogen and helium ionization zones (Figure 4(b)): decreasing for the hydrogen and increasing for the He ii ionization zones. For comparison, we plot the values of γ (thin curve) calculated from the horizontally averaged values of the density and internal energy of the 3D model. The differences of the curve from the 3D model (thick curve) in the near-surface layers indicate that the nonlinear dynamical fluctuations caused by turbulence and waves play an important role in the thermodynamic properties. These deviations from the 1D models require further investigation and have to be taken into account when estimating helium abundances by helioseismology and asteroseismology.

The region of convective overshoot is defined in terms of the temperature gradient as a transition zone between the adiabatic and radiative gradients (e.g., Christensen-Dalsgaard et al. 2011). It is indicated by two vertical lines in Figures 4(c)–(d), from −28.5 Mm to −36.5 Mm in depth, and is approximately 1/2 of the pressure scale height. This closely corresponds to the region of negative enthalpy flux (Figure 3(h)), which is a distinct feature of convective overshooting. The simulation results clearly show that the overshooting region consists of two parts: the upper part represents an almost adiabatic extension of the convection zone, while the lower part shows a significant peak of increased subadiabaticity. However, most previous "semi-empirical" models (see Christensen-Dalsgaard et al. 2011 for a review) described the overshoot region as an extension of the depth of the convection zone with a smooth monotonic transition between the adiabatic and radiative temperature (or entropy) gradients. The numerical simulations presented in the previous section show that the braking of penetrating convective plumes results in heating of the overshoot region. This leads to a sharp decrease of the temperature gradient in the lower part of the overshoot region. This effect has to be taken into account in the construction of 1D models of convective overshoot based on modification of the temperature gradient in the standard solar and stellar models.

The sound-speed variations are of particular interest because asteroseismic estimations based on p-mode frequencies are mostly sensitive to the sound-speed profile. Usually, the sound-speed profile calculated from a 1D stellar model is used as a reference; therefore, it is important to investigate the deviations in sound speed between the 3D simulation and 1D model. Figure 5(a) shows the relative sound-speed difference. This profile reveals strong deviations from the 1D model. These deviations are qualitatively similar to the findings of global helioseismology analysis reproduced in Figure 5(e) (e.g., Kosovichev et al. 1997). In particular, the 3D simulations reproduce the characteristic "bump" at the bottom of the convection zone (z ≈ −33 Mm in Figure 5(a)). This bump corresponds to the overshoot layer, which is about 8 Mm thick. Our simulation results show a qualitative agreement for this bump's shape with helioseismology results, but the amplitude of the variations in the stellar simulations is much larger than in the solar measurements, obviously, because the convection is much more vigorous in the F-type star than in the Sun. Nevertheless, the surprising similarity supports an idea about the influence of the convective overshoot in the solar sound-speed inversions. In other panels of Figure 5, we illustrate the qualitative agreement for the relative variations of density, the parameter of convective stability, and for the adiabatic exponent. The good correspondence leads us to the suggestion that the structure of the overshooting region in the Sun is similar to the overshoot structure obtained in our simulations.

Zoom In Zoom Out Reset image size

We have performed, for the first time, detailed 3D fully compressible radiative hydrodynamic simulations of the convection zone dynamics of a moderate-mass F-type star through its whole depth and into the transition layer between the convection and radiative zones ("overshoot zone"). A characteristic feature of the convection zone is a co-existence of several granulation scales at the photosphere. The convection pattern represents well-defined clusters of several granules. In deeper layers of the convection zone, the scale of the convection patterns gradually increases (Figure 1). The downflows in the intergranular lanes of the bigger granules penetrate through the whole convection zone, reaching velocities of more than 20 km s⁻¹, close to the local speed of sound (Figure 2). These downdrafts penetrate into the radiative zone, form an overshoot layer, and cause local heating, thus increasing the sound speed and initiating strong density variations that can be a source of internal gravity oscillations (g-modes).

Comparison of the initial 1D model calculated from the stellar evolution code (CESAM) with the corresponding mean profiles from the 3D numerical simulation shows that in the 3D simulations the stellar radius is increased by about 800 km, and the convection zone is expanded in depth due to the development of an overshoot region about 8 Mm (1/2 of the pressure scale height) thick (Figure 4). The overshoot region is characterized by a negative enthalpy flux and, most importantly, has a two-layer structure: the upper part represents an almost adiabatic extension of the convection zone, while the lower part shows a significant peak of increased subadiabaticity, caused by heating of the overshooting region due to braking of penetrating convective plumes. In addition, the simulations reveal a layer of enhanced turbulence, associated with the ionization zones, which affects variations of the adiabatic exponent.

The deviations of the mean properties between the 3D and 1D models are remarkably similar to the corresponding deviations between helioseismology inversions and the standard solar model (Figure 5), indicating that, perhaps, the helioseismology results including the characteristic sound-speed bump may be explained by dynamical effects of solar convection and overshooting. The differences in the radius between the 3D and 1D models may be partially explained by adjusting the mixing-length parameter, but may also be related to the observed difference between the photospheric and seismic radii of the Sun (Schou et al. 1997). Further understanding and modeling of these discrepancies are important for correct interpretation of asteroseismology and helioseismology observational data and better understanding of the turbulent dynamics of stellar convection zones.