© ESO, 2011

1. Introduction

The first theoretical calculations of the Cepheids instability strip done in the 60’s assumed that convection was steady with respect to oscillations. Unfortunately, this “frozen-in convection” approximation led to a cooler red edge than the observed one, because the strong coupling between convection and pulsations occurring in cool Cepheids was ignored (e.g. Baker & Kippenhahn 1965). Then, following the pioneering works of Unno (1967) and Gough (1977), several time-dependent convection (TDC) models have been developed to investigate the influence of the convection on the pulsational stability (e.g. Stellingwerf 1982; Kuhfuß 1986; Gehmeyr & Winkler 1992a). The last up-to-date TDC models actually succeed in predicting both a red edge close to the observed one and realistic luminosity curves (e.g. Yecko et al. 1998; Bono et al. 1999; Feuchtinger 1999; Kolláth et al. 2002).

However, all of these TDC models suffer from a common weakness due to the numerous free parameters, usually known as α coefficients, which describe the turbulent convection (e.g., the eight dimensionless parameters in Smolec & Moskalik 2010). These parameters are either obtained from a fit to the observations or hardly constrained by theory when taking the asymptotic limit of stationary convection (Vitense 1953; Böhm-Vitense 1958). A parametric study carried out by Yecko et al. (1998) emphasises the intrinsic degeneracy of TDC models, because similar instability strips have been obtained with different sets of parameters.

Direct numerical simulations (hereafter DNS) are able to constrain these TDC models because they fully account for the non-linearities involved in the convection-pulsation coupling. These nonlinear simulations are challenging as they require both large-amplitude oscillations and convective motions. In our last 2-D simulations, we improved the way acoustic waves are generated by reproducing the self-consistent excitation operating in variable stars, known as the κ-mechanism (Gastine & Dintrans 2011, hereafter GD2011). The resulting coupling of acoustic modes with convection is therefore more consistent because the mode amplitude is not imposed artificially. We showed in GD2011 that the convective plumes may either quench the radial oscillations or coexist with the acoustic modes, depending mainly on the density contrast of the equilibrium model.

The purpose of this letter is to compare the fully nonlinear results with two main TDC models: (i) the first one refers to an initial formulation of Stellingwerf (1982) that was used and improved by Bono & Stellingwerf (1994) and Bono et al. (1999); (ii) the other one was developed by Kuhfuß (1986) and Gehmeyr & Winkler (1992a) and implemented both in the Vienna hydrocode (e.g. Wuchterl & Feuchtinger 1998; Feuchtinger 1999) and the Florida-Budapest one (e.g. Yecko et al. 1998; Kolláth et al. 2002). In these models, a single equation for the turbulent kinetic energy ℰ_t is added to the classical mean-field equations and the main second-order correlations, such as the convective flux, are expressed as a function of ℰ_t alone (see e.g. Baker 1987; Gehmeyr & Winkler 1992b; Buchler 2000, 2009).

We investigate here a particular simulation where the oscillations strongly modulate the convective flux in more detail. The nonlinear results are first compared at each snapshot with the TDC recipes by computing χ²-statistics of the relevant α coefficients. Second, the mean values of α are used to compare the optimal TDC fluxes with the DNS ones. The formulation of Stellingwerf is found to be closer to the nonlinear result than the Kuhfuß one: (i) the temporal variation of the α coefficient is weaker; (ii) the mean convective flux is closer to the DNS one, especially in the description of the overshooting area.

2. The hydrodynamic model

We consider a local 2-D box of size L_x × L_z, filled by a perfect monatomic gas and centred on both sides of an ionisation region. Its associated opacity bump is shaped by a hollow in the temperature-dependent radiative conductivity profile K(T). Such a configuration can lead to unstable acoustic modes due to the driving term ∇·K(T)∇T in the energy equation; i.e., this is the κ-mechanism (Gastine & Dintrans 2008a,b). Furthermore, this hollow in K(T) is deep enough that the equilibrium temperature gradient is locally superadiabatic and convective motions develop here according to Schwarzschild’s criterion.

The governing hydrodynamic equations are written in non-dimensional form by choosing L_z as the length scale and $\sqrt{{\mathit{c}}_{\mathit{p}}{\mathit{T}}_{\mathrm{top}}}$ as the velocity one, hence the time scale ${\mathit{L}}_{\mathit{z}}\mathit{/}\sqrt{{\mathit{c}}_{\mathit{p}}{\mathit{T}}_{\mathrm{top}}}$ (with c_p the specific heat and T_top the surface temperature). The resulting nonlinear set of equations is advanced in time with the high-order, finite-difference pencil code¹, which is fully explicit except for the radiative diffusion term that is solved implicitly thanks to a parallel alternate direction implicit (ADI) solver. With the chosen units, the simulation box spans about 10% of the star radius on both sides of the ionisation region, while the time-step is about 1 minute, such that a simulation typically spans over 4500 days (see GD2011).

3. Results

The time evolution of the convective flux obtained in fully nonlinear 2-D simulations is compared with the following TDC expressions developed by Stellingwerf (1982) and Kuhfuß (1986):

$\{\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}{\mathrm{ℱ}}_{\mathrm{St}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{=}{\mathit{\alpha }}_{\mathrm{St}}\frac{\mathit{A}}{\mathit{B}}{\mathrm{ℰ}}_{\mathrm{t}}\mathrm{sign}\mathrm{(}\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}\mathrm{)}\sqrt{\mathrm{|}\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}\mathrm{|}}\mathit{,}\begin{array}{}\end{array}\\ \begin{array}{}\end{array}{\mathrm{ℱ}}_{\mathrm{Ku}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{=}{\mathit{\alpha }}_{\mathrm{Ku}}\mathit{A}\sqrt{{\mathrm{ℰ}}_{\mathrm{t}}}\left(\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}\right)\mathit{,}\end{array}$(1)where ∇ = dlnT/dlnp, ∇_ad = 1 − c_v/c_p and

${\mathrm{ℰ}}_{\mathrm{t}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{=}⟨\frac{{{\mathit{u}}_{\mathit{z}}^{\mathrm{\prime }}}^{\mathrm{2}}}{\mathrm{2}}⟩\mathit{,}\mathit{A}\mathrm{=}{\mathit{c}}_{\mathit{p}}⟨\mathit{\rho }⟩⟨\mathit{T}⟩\mathrm{and}\mathit{B}\mathrm{=}\sqrt{{\mathit{c}}_{\mathit{p}}⟨\mathit{T}⟩{\mathrm{\nabla }}_{\mathrm{ad}}}\mathit{,}$(2)with p and ρ the pressure and density, respectively, and the brackets denote a horizontal average. Two free, dimensionless parameters α_St and α_Ku are introduced that the simulations allow to constrain. We first focus on a 2-D DNS that is similar to the G8 simulation in GD2011, which is a simulation in which the total kinetic energy is almost entirely contained in acoustic modes excited by the κ-mechanism (80%), with the rest in convective plumes (20%). The convective motions do not affect the growth and the nonlinear saturation of the unstable radial acoustic mode much. We can therefore expect that the heat transport is strongly modulated by wave motions, which is an ideal frame for testing the accuracy of time-dependent convection models.

3.1. Temporal modulation of the convective flux in the DNS

Fig. 1 a) Mean vertical profiles of radiative ℱrad (solid blue line), turbulent enthalpy ℱconv (dashed green line), kinetic ℱkin (dotted black line), modes ℱmod (dot-dashed magenta line), and total ℱtot (dotted red line) fluxes, normalised to the bottom flux Fbot. b) Temporal power spectrum for the convective flux alone.

Fig. 2 Evolution of the convective flux in a (t,  z) plane: a) in the DNS; b) with Stellingwerf’s formalism ℱSt; c) with Kuhfuß’s formalism ℱKu (for αSt = αKu = 1). The vertical extent is centred on both sides of the convection zone to emphasise its oscillations. These snapshots were computed after 1800 periods of oscillations, i.e. well after the nonlinear saturation was achieved.

The imposed bottom flux F_bot is transported mainly through the computational domain by the radiative ℱ_rad, enthalpy, and kinetic ℱ_kin fluxes. Following GD2011, the enthalpy flux is divided into the classical convective flux ℱ_conv and the contribution coming from the (unstable) fundamental mode (hereafter ℱ_mod) with

$\{\begin{array}{c}\\ {\mathrm{ℱ}}_{\mathrm{conv}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}& \mathrm{=}{\mathit{c}}_{\mathit{p}}⟨\mathit{\rho }{\mathit{u}}_{\mathit{z}}{\mathit{T}}^{\mathrm{\prime }}⟩\mathit{,} & {\mathrm{ℱ}}_{\mathrm{mod}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{=}{\mathit{c}}_{\mathit{p}}⟨\mathit{\rho }{\mathit{u}}_{\mathit{z}}⟩\mathit{\theta ,}\\ \\ {\mathrm{ℱ}}_{\mathrm{rad}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}& \mathrm{=}\mathrm{-}⟨\mathit{K}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{\nabla }\mathit{T}⟩\mathit{,} & {\mathrm{ℱ}}_{\mathrm{kin}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{⟨}\mathit{\rho }{\mathit{u}}_{\mathit{z}}{\mathit{u}}^{\mathrm{2}}\mathrm{⟩}\mathit{,}\end{array}$(3)where the primed quantities denote the fluctuations about the horizontal average, and θ is the temperature eigenfunction of the fundamental mode. The resulting time-averaged and normalised fluxes are given in Fig. 1a.

The bulk of the total flux is transported by the radiative flux, except in the convective zone where ℱ_conv/F_bot ≃ 20%, while the kinetic flux is negligible (ℱ_kin/F_bot ≤ 1%). The value of ℱ_mod is hardly as large as ℱ_conv in the convective zone. This quantity is a good signature of the amplitude of the acoustic modes – the higher ℱ_mod, the larger the radial oscillations (GD2011) – therefore confirming the efficiency of the κ-mechanism in this simulation.

The signature of the temporal modulation of the convective flux is extracted from its Fourier spectrum in time; that is, we first compute $\begin{array}{}{}_{\mathrm{􏽢}}\\ {\mathrm{ℱ}}_{\mathrm{conv}}\end{array}\mathrm{\left(}\mathit{z,\omega }\mathrm{\right)}$, with ω the frequency and then integrate over the vertical direction z to get the power spectrum $\begin{array}{}{}_{\mathrm{􏽢}}\\ {\mathrm{ℱ}}_{\mathrm{conv}}\end{array}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ (Fig. 1b). Several discrete peaks appear about given frequencies, of which the physical nature is emphasised after superimposing the linear acoustic eigenfrequencies (the vertical dashed blue lines). The match with the fundamental mode frequency ω₀₀ = 3.85 is perfect, while the weak-amplitude secondary peaks instead correspond to harmonics of ω₀₀ (i.e. 2ω₀₀,3ω₀₀,...   , shaped by downward-directed vertical grey arrows in Fig. 1b) than to the acoustic overtones ω₀₁,  ω₀₂,...    It means that the amplitude of the fundamental acoustic mode is large enough to generate several harmonics through a nonlinear cascade, and this again indicates both the robustness of the κ-driving and the relevance of this kind of convection-pulsation simulation to check the TDC recipes.

3.2. DNS vs. TDC models: convective patterns

We first compare the time evolution of the horizontally averaged convective flux obtained in the DNS (Eq. (3)) with its TDC counterparts ℱ_St and ℱ_Ku (Eq. (1)). As we are interested here in the qualitative agreement between the convective patterns in the plane (t,  z), we simply assume α_St = α_Ku = 1.

The three resulting patterns are displayed in Fig. 2 for a time interval spanning about four periods of the fundamental acoustic mode. The black areas denote positive values for the convective flux and therefore delimit the convective zone. An oscillatory behaviour is clearly visible in each panel, with a period that looks similar to the one of the fundamental acoustic mode; that is, almost 4 oscillation cycles are depicted. This is consistent with the large peak shown around ω₀₀ in the power spectrum in Fig. 1b.

The two TDC patterns also display good agreement with the nonlinear simulation regarding the overshooting phenomenon. We indeed recover the same dark red structures below the convective zone, and they correspond to the downdrafts entering the radiative zone (where this penetration is associated with a negative convective flux). Nevertheless, we note that the overshooting seems to be more vigorous in Kuhfuß’s model than in Stellingwerf’s, because these dark red structures fill a larger surface in Fig. 2c in the bottom radiative zone than in Fig. 2b.

3.3. DNS vs. TDC models: statistics of coefficients α

Fig. 3 Time evolution around the mean value of coefficients αSt (upper panel) and αKu (bottom panel). Horizontal grey spans mark the limits of the associated relative standard deviation.

A one-to-one comparison between the convective flux in the DNS and its TDC predictions requires finding the optimal values of α_St and α_Ku. This is done by performing several χ²-tests at different snapshots in the simulation to track the variations in coefficients α. The resulting fluctuations of α_St and α_Ku around their mean values are shown in Fig. 3.

The dispersion of the Stellingwerf coefficient α_St (upper panel) is weaker than the Kuhfuß one (lower panel), since its values are almost within a 5% range around the mean $\overline{{\mathit{\alpha }}_{\mathrm{St}}}\mathrm{=}\mathrm{0.462}$. In contrast, several outliers with values greater than 10% are found in the evolution of α_Ku which then appears more chaotic around the mean value $\overline{{\mathit{\alpha }}_{\mathrm{Ku}}}\mathrm{=}\mathrm{0.076}$. As a consequence, the relative standard deviation (depicted in gray in Fig. 3) is weaker in Stellingwerf’s case than in Kuhfuß’s.

3.4. DNS vs. TDC models: convective fluxes with optimal α’s

We recall that TDC models assume that the coefficients α entering Eq. (1) are constant. The optimal value of these coefficients has been deduced by adjusting the TDC recipes with the instantaneous convective flux throughout the nonlinear simulation. The final test, given in Fig. 4, then consists in comparing the mean nonlinear convective flux taken over the entire simulation and its best TDC approximations built from these optimal α values.

This figure emphasises that Stellingwerf’s formulation gives a better agreement than Kuhfuß’s one. Indeed, the Kuhfuß model overestimates the overshooting as the (negative) convective flux remains non-negligible until the bottom of the radiative zone. In contrast, the Stellingwerf profile accounts for the local penetration of convective plumes better and shows the same exponential-like decay in the negative convective flux when sinking in the radiative zone (e.g. Dintrans 2009). However, the two models are fairly similar in the bulk of the convective zone where convection is fully developed. One also notes that they both predict a negative flux at the top of the convective zone, which is an upper overshooting of convection motions near the surface that is not observed in the nonlinear simulation.

4. Conclusion

Fig. 4 Mean convective flux in the DNS (solid black line), compared with the best TDC predictions based on the optimal α values that came out of the χ2-test in Sect. 3.3 (Stellingwerf’s model is a dashed blue line and Kuhfuß’s a dotted green line).

The main weakness of all theories of time-dependent convection (TDC) lies in the large number of free parameters. This is particularly awkward when convection is strongly coupled with pulsations, such as near the red border of the Cepheid instability strip where similar results are obtained with different sets of parameters (Yecko et al. 1998). Additional constraints must be found to reduce the intrinsic degeneracy of models. But the modelling of the interplay between the turbulent convection and oscillations is really a difficult task owing to the strong nonlinearities involved in this coupling.

One solution to tackle this problem and to bring new constraints consists in performing fully nonlinear simulations, and this is the path we chose in this work following our first study in GD2011. Two widely used TDC theories, namely the Stellingwerf (1982) and Kuhfuß (1986) ones, are compared with results coming from 2-D nonlinear simulations of compressible convection in which strong acoustic oscillations are self-sustained by the κ-mechanism. The heat transport is then modulated by the fundamental acoustic mode such that this kind of simulation is relevant for investigating the convection-pulsation coupling (Figs. 1, 2).

Focusing on the two TDC formulae for the convective flux, we computed the evolution of free coefficients “α” from a χ²-test applied to the fully nonlinear results (Fig. 3). A large temporal variability is found in both cases, and it weakens the robustness of the TDC assumption α = const. Moreover, the mean values $\overline{\mathit{\alpha }}$ are not universal. By applying the same method to other simulations performed in GD2011 (Table 1), we indeed do not recover the same $\overline{\mathit{\alpha }}$ with a relative standard deviation of about 12% (Stellingwerf) and 18% (Kuhfuß).

Within these limits, Stellingwerf’s formulation is found to agree with the nonlinear simulations better than does Kuhfuß’s: (i) the final mean convective flux ℱ_St is closer to its DNS counterpart with a much better estimation of the bottom overshooting; (ii) the temporal dispersion of the α_St coefficient is weaker; that is, the stability is enhanced (Fig. 4). This result means that the time-dependent convective flux scales better with a law ${\mathrm{ℱ}}_{\mathrm{conv}}\mathrm{\propto }{\mathrm{ℰ}}_{\mathrm{t}}\sqrt{\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}}$ than ${\mathrm{ℱ}}_{\mathrm{conv}}\mathrm{\propto }\sqrt{{\mathrm{ℰ}}_{\mathrm{t}}}\mathrm{(}\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}\mathrm{)}$, and the Stellingwerf formulation may probably be preferred in the 1-D hydrocode used in, e.g., the topical Cepheids models. However, this study emphasises that both formalisms lead to an artificial overshooting at the top of the convection zone. We also noted that this TDC test involves the exact value of the turbulent kinetic energy ℰ_t provided by the nonlinear simulation, and not by the dedicated 1-D TDC equation, which is inherently less accurate. As a consequence, the obtained profiles are certainly the best we can expect from these TDC recipes.

In this work, the temporal modulation of the convective flux was ensured by the acoustic modes excited by κ-mechanism. An interesting prospect could be the case of a modulation based on the internal gravity waves excited by convection itself. Indeed, convection can excite gravity waves in variable stars, either by convective elements penetrating stably stratified regions, such as in solar-type stars (e.g. Dintrans et al. 2005); or through the so-called “convective blocking” mechanism, as in white dwarfs or Gamma Doradus stars (e.g. Pesnell 1987). In both cases, the convective flux is ipso facto modulated by gravity waves, and it may be interesting in that respect to also check the accuracy of time-dependent convection models.

Online material

Appendix A: Setup of the simulations

A.1. The equilibrium model

As in GD2011, our system represents a zoom around an ionisation region. Since we are computing local simulations, the vertical gravity g =  −ge_z and the kinematic viscosity ν are assumed to be constant. Following our purely radiative model of the κ-mechanism (Gastine & Dintrans 2008a,b), the ionisation region is represented by a temperature-dependent radiative conductivity profile that mimics an opacity bump:

$\mathit{K}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathit{K}}_{\max }\left[\mathrm{1}\mathrm{+}\mathrm{𝒜}\frac{\mathrm{-}\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{+}\arctan \mathrm{\left(}\mathit{\sigma }{\mathit{T}}^{\mathrm{+}}{\mathit{T}}^{\mathrm{-}}\mathrm{\right)}}{\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{+}\arctan \mathrm{\left(}\mathit{\sigma }{\mathit{e}}^{\mathrm{2}}\mathrm{\right)}}\right]\mathit{,}$(A.1)

with

$\mathrm{𝒜}\mathrm{=}\frac{{\mathit{K}}_{\max }\mathrm{-}{\mathit{K}}_{\min }}{{\mathit{K}}_{\max }}\mathit{,} {\mathit{T}}^{\mathrm{\pm }}\mathrm{=}\mathit{T}\mathrm{-}{\mathit{T}}_{\mathrm{bump}}\mathrm{\pm }\mathit{e,}$(A.2)where T_bump is the position of the hollow in temperature and σ, e, and denote its slope, width, and relative amplitude, respectively. We assume both radiative and hydrostatic equilibria; that is,

$\{\hspace{0.17em}\begin{array}{}\begin{array}{}\end{array}\frac{\mathrm{d}{\mathit{p}}_{\mathrm{0}}}{\mathrm{d}\mathit{z}}\mathrm{=}\mathrm{-}{\mathit{\rho }}_{\mathrm{0}}\mathit{g,}\begin{array}{}\end{array}\\ \begin{array}{}\end{array}\frac{\mathrm{d}{\mathit{T}}_{\mathrm{0}}}{\mathrm{d}\mathit{z}}\mathrm{=}\mathrm{-}\frac{{\mathit{F}}_{\mathrm{bot}}}{{\mathit{K}}_{\mathrm{0}}\mathrm{\left(}{\mathit{T}}_{\mathrm{0}}\mathrm{\right)}}\mathit{,}\end{array}$(A.3)

where F_bot is the imposed bottom flux. Following GD2011, we chose L_z as the length scale, i.e.  [x]  = L_z, top density ρ_top and top temperature T_top as density and temperature scales, respectively. The velocity scale is then $\sqrt{{\mathit{c}}_{\mathit{p}}{\mathit{T}}_{\mathrm{top}}}$, while time is given in units of $\mathrm{\left[}\mathit{t}\mathrm{\right]}\mathrm{=}{\mathit{L}}_{\mathit{z}}\mathit{/}\sqrt{{\mathit{c}}_{\mathit{p}}{\mathit{T}}_{\mathrm{top}}}$.

Table A.1 then summarises the parameters of the numerical simulations presented in this study in these dimensionless units. The penultimate column of this table contains the value of the frequency ω₀₀ of the fundamental unstable radial mode excited by the κ-mechanism, which lies between 3 and 4 for every DNS. The last column gives the value of the Rayleigh number, which quantifies the strength of the convective motions. It is given by

$\mathrm{Ra}\mathrm{=}\frac{\mathit{g}{\mathit{L}}_{\mathrm{conv}}^{\mathrm{4}}}{\mathit{\nu \chi }{\mathit{c}}_{\mathit{p}}}\left|\frac{\mathrm{d}\mathit{s}}{\mathrm{d}\mathit{z}}\right|\mathit{,}$(A.4)

where L_conv is the width of the convective zone, χ = K₀/ρ₀c_p the radiative diffusivity, and s the entropy.

A.2. The nonlinear equations

With the parallel version of the alternate direction implicit (ADI) solver for the radiative diffusion implemented in the pencil code (see GD2011), we advance the following hydrodynamic equations in time:

$\{\begin{array}{c}\\ \frac{\mathrm{D}\ln \mathit{\rho }}{\mathrm{D}\mathit{t}}& \mathrm{=}\mathrm{-}\mathrm{div}{u}\mathit{,}\\ \\ \frac{\mathrm{D}{u}}{\mathrm{D}\mathit{t}}& \mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho }}\mathbf{\nabla }\mathit{p}\mathrm{+}{g}\mathrm{+}\mathrm{2}\mathit{\nu }\left(\mathbf{\nabla }\mathrm{·}\mathbf{S}\mathrm{+}\mathbf{S}\mathrm{·}\mathbf{\nabla }\ln \mathit{\rho }\right)\mathit{,}\\ \\ \frac{\mathrm{D}\mathit{T}}{\mathrm{D}\mathit{t}}& \mathrm{=}\frac{\mathrm{1}}{\mathit{\rho }{\mathit{c}}_{\mathit{v}}}\mathrm{div}\mathit{K}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathbf{\nabla }\mathit{T}\mathrm{-}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\mathit{T}\mathrm{div}{u}\mathrm{+}\mathrm{2}\mathit{\rho \nu }\mathbf{S}\mathrm{2}\mathit{,}\end{array}$(A.5)

where ρ,  u, and T denote density, velocity, and temperature, respectively, while K(T) is given by Eq. (A.1). The operator D/Dt = ∂/∂t + u·∇ is the usual total derivative, while S is the (traceless) rate-of-strain tensor given by

${\mathrm{S}}_{\mathit{ij}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathit{\partial }{\mathit{u}}_{\mathit{i}}}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}\mathrm{+}\frac{\mathit{\partial }{\mathit{u}}_{\mathit{j}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\delta }}_{\mathit{ij}}\mathrm{div}{u}\right)\mathit{.}$(A.6)

Finally, we impose the condition that all fields are periodic in the horizontal direction, while stress-free boundary conditions (i.e., u_z = 0 and du_x/dz = 0) are assumed for the velocity in the vertical one. Concerning the temperature, a perfect conductor at the bottom (i.e., flux imposed) and a perfect insulator at the top (i.e., temperature imposed) are applied.

To ensure that both the nonlinear saturation and thermal relaxation are achieved, the simulations were computed over very long times, typically t ≳ 3000. As the eigenfrequency of the unstable acoustic mode ω₀₀ ∈  [3 − 4]  (see Table A.1), this corresponds approximately to 1800 periods of oscillation.

Acknowledgments

This work was granted access to the HPC resources of CALMIP under the allocation 2010-P1021 (http://www.calmip.cict.fr).

References

All Tables

All Figures

	Fig. 1 a) Mean vertical profiles of radiative ℱrad (solid blue line), turbulent enthalpy ℱconv (dashed green line), kinetic ℱkin (dotted black line), modes ℱmod (dot-dashed magenta line), and total ℱtot (dotted red line) fluxes, normalised to the bottom flux Fbot. b) Temporal power spectrum for the convective flux alone.
In the text

	Fig. 2 Evolution of the convective flux in a (t,  z) plane: a) in the DNS; b) with Stellingwerf’s formalism ℱSt; c) with Kuhfuß’s formalism ℱKu (for αSt = αKu = 1). The vertical extent is centred on both sides of the convection zone to emphasise its oscillations. These snapshots were computed after 1800 periods of oscillations, i.e. well after the nonlinear saturation was achieved.
In the text

	Fig. 3 Time evolution around the mean value of coefficients αSt (upper panel) and αKu (bottom panel). Horizontal grey spans mark the limits of the associated relative standard deviation.
In the text

	Fig. 4 Mean convective flux in the DNS (solid black line), compared with the best TDC predictions based on the optimal α values that came out of the χ2-test in Sect. 3.3 (Stellingwerf’s model is a dashed blue line and Kuhfuß’s a dotted green line).
In the text