Modelling turbulent stellar convection zones: sub-grid scales effects

Modelling turbulent stellar convection zones: sub-grid scales effects

A. Strugarek P. Beaudoin A. S. Brun P. Charbonneau S. Mathis P. K. Smolarkiewicz Département de physique, Université de Montréal, C.P. 6128 Succ. Centre-Ville, Montréal, QC H3C-3J7, Canada Laboratoire AIM Paris-Saclay, CEA/DSM Université Paris-Diderot CNRS, IFRU/SAp, F-91191 Gif-sur-Yvette, France. European Centre for Medium-Range Weather Forecasts, Reading RG2 9AX, UK

The impressive development of global numerical simulations of turbulent stellar interiors unveiled a variety of possible differential rotation (solar or anti-solar), meridional circulation (single or multi-cellular), and dynamo states (stable large scale toroidal field or periodically reversing magnetic fields). Various numerical schemes, based on the so-called anelastic set of equations, were used to obtain these results. It appears today mandatory to assess their robustness with respect to the details of the numerics, and in particular to the treatment of turbulent sub-grid scales. We report on an ongoing comparison between two global models, the ASH and EULAG codes. In EULAG the sub-grid scales are treated implicitly by the numerical scheme, while in ASH their effect is generally modelled by using enhanced dissipation coefficients. We characterize the sub-grid scales effect in a turbulent convection simulation with EULAG. We assess their effect at each resolved scale with a detailed energy budget. We derive equivalent eddy-diffusion coefficients and use the derived diffusivities in twin ASH numerical simulations. We find a good agreement between the large-scale flows developing in the two codes in the hydrodynamic regime, which encourages further investigation in the magnetohydrodynamic regime for various dynamo solutions.

convection – turbulence – dynamo – stars: interiors – stars: kinematics and dynamics
journal: Advances in Space Research

1 Introduction

Cool stars are known to possess a substantial convection zone forming their outer layer. The convective motions participate in the self-organization of the interior of the star and in particular in the sustainment of a large scale differential rotation (e.g. Brun and Toomre 2002; Featherstone and Miesch 2015, and references therein), and a potentially cyclic magnetism (e.g. Ghizaru et al. 2010; Käpylä et al. 2012; Augustson et al. 2015; Brun et al. 2015, and references therein). Understanding the properties of solar and stellar convection is hence of prior importance to understand the joint evolution of the stellar rotation rate and large-scale magnetic fields along the Main Sequence.

Thanks to helioseismology (Christensen-Dalsgaard et al. 1991; Basu 1997), we know that the solar stratification is close to adiabatic down to 0.71 solar radii, leaving little doubt that convective motions extend down to such depths. However, the amplitude of these deep convective flows in the Sun is today the subject of a stimulating controversy. Using time-distance helioseismology, Hanasoge et al. (2010, 2012) obtained surprisingly low upper-limits ( cm/s at depth Mm) for the large-scale convection motions in the solar convection zone. It is a possibility that solar convection models over-estimate by a few orders of magnitude the deep convective flow. For instance, an inadequate cut in the turbulent scales due to finite numerical resolution alters the spectral repartition of energy and hence may lead to incorrect estimates of the large-scale convective flows. Another related possibility is that convection may actually be driven by the strong cooling layer at the solar surface (Spruit 1997), through the so-called entropy rain of small turbulent scales. On the other hand, more recent observational results using ring-diagram analysis were recently obtained by Greer et al. (2015), who found convective flows two orders of magnitude faster than Hanasoge et al. (2010) at depth Mm. More cross-comparisons between helioseismology techniques used to estimate convective flows in the Sun are today needed to carefully pin down those observational constraints.

There is yet a somewhat more fundamental issue regarding our understanding of stellar convective turbulence, disregarding temporarily the exact mechanism exciting it. Stellar interiors are mainly composed of stratified, fully ionized hydrogen which can be modelled as a one-fluid plasma under the magnetohydrodynamic (MHD) approximation (which is essentially composed of a form of the Navier-Stokes equations combined with the heat transport equation and coupled to an induction equation). The microscopic dissipation capabilities of the stellar plasma are relatively low: in the Sun, the microscopic viscosity in the convection zone typically lies in the cm/s range (see, e.g., Miesch 2005). Considering the lowest large-scale convective velocities estimates from Hanasoge et al. (2012) ( cm/s), and a typical solar convection zone depth , the typical Reynolds number of the large-scale flows in the solar convection zone is at the very least . Such a tremendously high Reynolds number is, for the time being, unfortunately inaccessible to both numerical simulations and laboratory experiments. As a result, the spectral repartition of energy (large-scales vs small-scales, direct vs inverse cascades of energy) in the solar convection zone is today unknown. The microscopic Prandtl number in the solar convection zone () is furthermore very challenging to reach with current numerical simulations, since it implies a difference of at least three orders of magnitude for viscous and heat dissipation time-scales. Finally, stellar convection zones are also magnetized and are thought to generally support dynamo action, making their theoretical and numerical modelling an outstanding challenge.

In order to approach the extreme parameters of solar (and more generally stellar) convection, which are unreachable with present computational power, various numerical simulation techniques have been developed. Two complementary paths can generally be followed. First, numerical simulations can be designed on localized, small portions of the solar convection zone (e.g. Rempel and Cheung 2014; Kitiashvili et al. 2015, and references therein). Even with this approach, solar parameters are extremely hard to achieve. Furthermore, the problem is truncated at the largest scales (due to the small box extent) and hence cannot address the important issue of large-scale convective motions and the sustainment of differential rotation or large-scale magnetism. A second path may thus be followed, where the global convection zone is modelled but the small-scales are parametrized. Both paths follow the general approach of large-eddy simulations (LES), on which we will focus in this work.

The simplest parametrization of sub-grid scales (SGS) consists in modelling their effects as enhanced viscosity and heat diffusivity. This approach has the advantage of giving the modeller full control on the sub-grid scales model, but lacks a physical justification in the context of solar (and stellar) convection. In contrast to classical turbulence, an inertial range for convective turbulence can be defined as the range of scales where the non-linear, local advective energy transport balances the buoyancy source of convection or the turbulent pressure gradient (see, e.g., Bolgiano 1959; Rincon 2006). If the smallest resolved scales of the LES model lie within this modified inertial range of the turbulent spectrum, the so-called dynamic Smagorinsky procedure (see Smagorinsky 1963; Germano et al. 1991) can be used to mimic a given self-similar spectrum for the smallest resolved scales of the model (see Nelson et al. 2013 for an implementation of such a method in the context of solar convection). Finally, another approach has been pursued with the so-called implicit-LES (ILES) methods. Indeed, no physically-rooted SGS model of the full MHD equations exists today (for recent advances in this direction, see Yokoi 2013; Chernyshov et al. 2014). A pragmatic approach can hence consist in minimizing (e.g. down to the numerical stability limit) the effect of the unresolved scales on the scales resolved by the model. The advantage of this method is that it ensures that the effect of sub-grid scales is minimized for a given grid size. However, the SGS model is fully subsumed by the numerical method, and its role in the development of the large scales is, unlike for explicit LES, non-trivial to estimate a posteriori.

The aim of this work is to compare comprehensively the ILES and LES modelling techniques for an idealized turbulent convection zone. We design a turbulent convection zone simulation based on the anelastic benchmark of Jones et al. (2011). We use the EULAG code (Prusa et al. 2008) to compute ILES based on the MPDATA algorithm (see, e.g. Smolarkiewicz and Charbonneau 2013), and use the ASH code (Clune et al. 1999; Miesch et al. 2000; Brun et al. 2004) to compute the LES counterparts. In a pioneering work, Elliott and Smolarkiewicz (2002) first attempted to reconcile empirically simulations of convective turbulence carried with LES and ILES. Here, we quantify the dissipation of the ILES with an original analysis of the energy transfers in spectral space. We show that the effect of the SGS modelling in the ILES can be interpreted in terms of enhanced viscosity and diffusivity operators similar to LES.

The paper is organized as follows. The set of anelastic equations used to describe convective turbulence in stars is given in Section 2, along with a description of the two numerical codes (EULAG and ASH) used in this work. A fiducial convective turbulence simulation in spherical geometry with the EULAG code is presented in Section 3. In Section 4, we develop an original method based on energy transfers in the spectral spherical harmonics space to estimate effective dissipation coefficients. These coefficients are then used in an ASH LES that satisfyingly mimics the EULAG ILES (Section 4.3). We finally conclude and summarize our results in Section  5.

2 Modelling stellar deep convection

2.1 Formulation of the anelastic equations

The equations solved by ASH and EULAG are the Lantz-Braginsky-Roberts (LBR, see Lantz and Fan 1999; Braginsky and Roberts 1995) —or, equivalently, the Lipps-Hemler (Lipps and Hemler 1982, 1985)— set of anelastic equations (see Vasil et al. 2013 for a recent review on sets of anelastic equations). The perturbed equations are written with respect to an ambient state (hereafter denoted with the subscript ) that theoretically may differ from the background state (hereafter denoted with bars) around which the generic anelastic equations were derived (both states will be detailed in Section 2.2). The background state is supposed to be isentropic, and both the ambient and background states are supposed to satisfy the hydrostatic equilibrium. The anelastic equations written in the stellar rotating frame are


where the perturbed quantities are denoted without prime for the sake of simplicity, and is the material derivative. We recall that we use standard notation for the basic fluid quantities, i.e. is the fluid velocity, its density, its pressure, and its specific entropy. The dissipative terms, when present, are defined by


where is the strain rate tensor and the identity tensor. Note that here we choose not to consider viscous heating, as in ILES with EULAG there is no equivalent in the implicit treatment of sub-grid scales. This choice has the potential drawback of not formally conserving the total energy in the system. In addition, we use a standard perfect gas equation of state which is linearized around the background state.

In the preceding equations, convection is forced by the conjunction of the advection of the unstable ambient entropy profile , and a Newtonian cooling term with a characteristic timescale (for details, see Prusa et al. 2008; Smolarkiewicz and Charbonneau 2013). The Newtonian cooling damps entropy perturbations over the timescale which is always chosen to exceed the convective overturning time. This ensures that on long time-scales, the model mimics a stellar convection zone remaining in thermal equilibrium (e.g. Cossette et al. 2016). We have modified the equations solved in ASH to include this Newtonian cooling term in order to force convection in exactly the same way in both codes. As a result, the same set of anelastic equations are solved in both cases, with the exception of explicit dissipation operators in ASH.

The anelastic equations can equivalently be specified in terms of potential temperature , which is related to the specific entropy through


Equations (2) and (3) can be written in terms of potential temperature through


Note that the ambient temperature is formally defined similarly to , i.e. . When deriving Equation (9) from (3) we have assumed that , which is a reasonable assumption given that the ambient entropy profile differs only by a small amount from the background entropy profile (see hereafter in Section 2.2).

2.2 Background and ambient states

Our numerical setup closely follows the anelastic benchmark of Jones et al. (2011). We consider a spherical shell of aspect ratio , and we note . We assume a gravity profile , for which the anelastic equations admit an equilibrium (denoted with bars) polytropic solution (see, e.g., Jones et al. 2011)


where is the polytropic index, are the density, pressure and temperature at the bottom of the domain, and the constants and are given by


where is the number of density scale heights in the layer. The background entropy profile is given by


where the standard adiabatic exponent for a perfect gas is . We choose a polytropic exponent to naturally ensure an isentropic background state.

The ambient state needs to be specified only in terms of entropy and potential temperature. The entropy jump throughout the domain, , is used to define the ambient entropy profile by


which we recall is related to the aforementioned ambient potential temperature profile by .

It should be noted that the background and ambient states used in this work differ from the ones used in past published ASH and EULAG simulations. We implemented those profiles in EULAG to be able to compare easily our results with ASH simulations. Furthermore, in this work only the convective layer is modelled, with no underlying stable layer.

2.3 Numerical methods

We use two codes based on different numerical methods.

The Eulerian-Lagrangian (EULAG) code is designed to use either Eulerian (flux form) or semi-Lagrangian (advective form) integration schemes (see Prusa et al. 2008; Smolarkiewicz and Charbonneau 2013). In the case presented here, Equations (6) and (8) are written as a set of Eulerian conservation laws and projected on a geospherical coordinate system (Prusa and Smolarkiewicz 2003). EULAG solves the evolution equations using MPDATA (multidimensional positive definite advection transport algorithm), which belongs to the class of nonoscillatory Lax-Wendroff schemes (Smolarkiewicz 2006), and is more specifically a second-order-accurate nonoscillatory forward-in-time template. Since all dissipation is delegated to MPDATA, this provides an implicit turbulence model (Domaradzki et al. 2003). Implicit dissipation diminishes if explicit dissipation is introduced in the model, providing seamless transition between ILES and LES (see Margolin et al. 2006, and references therein). In EULAG, all linear forcing terms are integrated in time using a second-order Crank-Nicholson scheme.

The Anelastic Spherical Harmonics (ASH) code is a pseudo-spectral code (see Boyd 1989; Glatzmaier 1984; Clune et al. 1999) based on a spherical harmonics decomposition, which avoids the classical issues related to the convergence of meridians at the poles of a sphere. As in EULAG, the linear terms of the anelastic equations are treated with an implicit Crank-Nicholson scheme of order 2. An Adams-Bashford scheme is used for the non-linear terms. The latter are evaluated in physical space, making the numerical method overall pseudo-spectral. In the radial direction, variables can be described either via a Chebyshev decomposition, or via a finite difference method (Alvan et al. 2014). We use the latter here with a fourth-order finite difference scheme.

Parameter Value
[] 1
[] 1
[10 rad s] 6.488
[dyne-cm g] 6.673 10
[erg g K] 3.4 10
[g cm] 0.2
Polytropic exponent 3/2
[s] 5.184 10
[erg g K] 2
Table 1: Fiducial stellar convection zone parameters. The solar radius is cm, and the solar mass g.

We consider in both codes stress-free, impermeable boundaries at the top and bottom of the domain such that


The entropy gradient is set to zero on the upper and lower boundaries. The simulations presented in this work are initialized with random, small-amplitude perturbations around the background profiles defined in Section 2.2.

3 Fiducial stellar convection zone with EULAG

Figure 1: Top panel: Background density (left vertical axis) and temperature (right vertical axis) profiles as a function of spherical radius. Bottom panel: Background and ambient potential temperature and entropy profiles. The left vertical axis corresponds to the potential temperature, and the right axis to the entropy.

We consider a fiducial stellar convection zone computed with the EULAG code. The parameters of our fiducial case are indicated in Table 1. We consider a convection zone with a solar-like aspect ratio and which rotates times faster than the Sun, which ensures that the Rossby number (where is the average vorticity in the middle of the convection zone) remains sufficiently smaller than one, and consequently that the model is strongly influenced by rotation. The background density and temperature profiles are shown in the upper panel of Figure 1. A moderate density contrast (1.5 density scale-heights) is adopted to limit the size of the smallest turbulent scales near the top of the domain. In the lower panel we show the background and ambient potential temperature (and equivalently entropy) profiles. A moderate entropy constrast over the convective shell is chosen to ensure the model is above the critical onset of convection, but remains in an accessible turbulent regime.

The numerical method in EULAG implicitly supplies the dissipation needed to maintain numerical stability, which in turns varies with the grid resolution. We consider two physically identical cases, labelled E1 and E2, in which the grid resolution is respectively = 51 64 128, and 101 128 256, with respective time steps of and seconds. By increasing the resolution (in both time and space) by a factor of 2, the implicit dissipation of the numerical scheme of EULAG is expected to be reduced by a factor of 4 (see, e.g., Prusa et al. 2008).

Figure 2: Left panels: Mollweide projections in the middle of the convection zone of the radial velocity (upper panels) and entropy perturbations (lower panels) for cases E1 (left) and E2 (right). The positive and perturbations are denoted in red, and the negative values in blue. Right panels: Differential rotation () profile (averaged over time and longitude) on the meridional plane for cases E1 (left) and E2 (right). The differential rotation profile is normalized to the stellar rotation rate . Rotation faster than is shown in red, slower in blue, and co-rotation appears in white.

We display in Figure 2 the radial velocity and entropy perturbations (left panels) in both cases on a Mollweide projection in the middle of the convection zone. As expected, smaller structures are observed when the resolution is increased. In both cases, the so-called banana cells clearly appear on the equator. Interestingly, the convective luminosity almost does not change between the two cases and peaks at about 111The solar luminosity is erg/s in the middle of the convection zone. The convective kinetic energy only increases by from E1 to E2, and the temperature perturbations compensate this variation which leads to a very similar convective luminosity. The main difference between the two cases hence lies in the kinetic energy of the differential rotation, which results from the complex interplay of the turbulent scales (through Reynolds stresses) involved in the simulation. The differential rotation (right panels) exhibits a solar-like (fast equator, slow poles) pattern in both cases. Case E2 possesses a significantly stronger pole-equator constrast () compared to E1 (). The iso- contours are in both cases mostly parallel to the rotation axis, and an almost co-rotating stripe is observed at mid-to high latitude, indicating that the simulation is indeed in a regime strongly influenced by rotation (low Rossby number).

4 Estimation of the sub-grid scale modelling effects in EULAG

We now quantitatively estimate the effect of the numerical treatment of sub-grid scales in EULAG. We base our analysis on the estimation of energy transfers in the turbulent convection zone simulations. We derive the spectral analysis method in Section  4.1 and detail the results in Section  4.2.

Figure 3: Axisymmetric (left panels) and non-axisymmetric (right panels) kinetic energy (upper panels) and entropy (lower panels) spectra averaged over 40 rotation periods. Case E1 spectra are shown in red and case E2 in blue. The spectra are summed over the spherical harmonics degree and displayed as a function of (see text).

4.1 Spectral Analysis: Method

We analyze the results from the EULAG code with the help of a spectral analysis of energy transfers between scales in the spherical harmonics space. This method was originally developed in Strugarek et al. (2013) in the context of characterizing the transfer of magnetic energy between scales in dynamo simulations. It was more recently adapted to the EULAG code to study the self-consistent development of MHD instabilities in the solar tachocline (Lawson et al. 2015). We extend it here to study kinetic energy and entropy spectral transfers (we refer the reader to Strugarek et al. 2013, and to A for more details about this spectral analysis method). The scalar and vectorial quantities are respectively projected on the standard and vectorial spherical harmonics bases (Rieutord 1987)




where are the Legendre polynomials.

Figure 4: Kinetic energy (upper panels) and entropy (lower panels) balance Equations (21)-(22) for cases E1 (left panels) and E2 (right panels), in the middle of the convection zone. For each equation the various terms are labeled with a coloured symbol, and the total of all the contributions is shown by the gray transparent stars.

We define the kinetic energy and entropy spectra by


where and the exponent denotes the complex conjugate. The subscript corresponds to the sum over the subset of spherical harmonics coefficients with fixed and restricted to a chosen ensemble . In this work we will consider the two ensembles and . We display those axisymmetric and non-axisymmetric spectra at mid-depth for cases E1 and E2 in Figure 3. The axisymmetric kinetic energy spectrum (upper left panel) is dominated by the differential rotation that populates the odd components of the spectrum. The entropy axisymmetric spectrum (lower left panel) is conversely dominated by even components, which correspond to the mean pole-equator temperature contrast that establishes itself in the simulations. The non-axisymmetric kinetic energy spectrum (upper right panel) exhibits a peak around in case E1 ( in case E2) which corresponds to the dominant convective structure that can be observed in the Mollweide projection in Figure 2. As expected, the non-axisymmetric kinetic energy spectrum at scales is shifted to smaller scales (higher ’s) when the resolution is increased, but the overall shape of the non-axisymmetric kinetic energy spectrum remains the same. The non-axisymmetric entropy spectrum peaks at in case E1 and in case E2. When the resolution is doubled, the entropy spectrum is not shifted to smaller scales but the negative slope extends down to the smallest scales resolved in the domain.

Using the anelastic Equations (2) and (3) with no explicit dissipation, we obtain the following evolution equations


In the kinetic energy Equation (21), the various terms in the right hand side correspond to contributions from the pressure gradient, buoyancy, Coriolis force and non-linear advection of momentum. Note that the Coriolis contribution vanishes when summed over all scales as it should, but is able to spectrally redistribute energy among neighbour shells (see also Augier and Lindborg 2013 and A.2.2). In the entropy Equation (22), they correspond to the ambient state advection, Newtonian cooling, and non-linear advection of the entropy perturbations. The detailed expressions of the different terms can be found in A.2 and A.3. We show in Figure 4 the various terms of Equations (21) and (22) for cases E1 and E2 at mid-depth (we focus here solely on the non-axisymmetric spectra). The different contributions are averaged over more than 40 stellar rotations after the simulations have reached a steady-state (i.e. after the total kinetic energy is stabilized). In the kinetic energy balance (upper panels), buoyancy (red squares) is clearly the source of energy for almost all turbulent scales and is generally opposed by the pressure gradient (green triangles). The non-linear advection contribution (magenta diamonds) changes sign at the peak of the spectrum because of a direct energy cascade: scales larger than the maximum energy scale lose energy to the smaller scales. The Coriolis force (blue circles) efficiently acts on the largest scales due to the small latitudinal extent of the higher modes, as expected.

The entropy balance is shown in the lower panels. The entropy perturbations draw energy from the ambient entropy profile at all scales (blue circles) and are stabilized by the non-linear advection (red squares). The Newtonian cooling term (green triangles) contributes very marginally to the entropy balance since we chose a long cooling timescale (see Table 1). We finally note that no cascading process is observed in the entropy transfers in these numerical experiments: the non-linear advection is negative for all the resolved scales.

The total right hand side of the kinetic energy and entropy evolution equations —shown as gray stars in Figure 4— should cancel out in steady-state. In Figure 4, we recall that the various contributions were averaged over a time period of 40 stellar rotations during which the total energy in the system has stabilized. Hence, the system is in steady-state and the imbalance (gray stars) can only be attributed to the effect of the numerical scheme that we do not represent in these plots and that tends to dissipate energy in EULAG. We now propose a possible interpretation of this imbalance to obtain quantitative estimates of the effect of the sub-grid scales modelling in EULAG.

4.2 Effective dissipation coefficients

In EULAG, the implicit treatment of sub-grid scales results in an additional contribution to the kinetic energy and entropy evolution equations that is a priori unknown. In the context of 3D homogeneous turbulence in a cartesian box, Domaradzki et al. (2003) showed that the implicit sub-grid scales treatment of MPDATA (the numerical algorithm behind EULAG) mimics qualitatively an eddy viscosity. We build here on this idea and attempt to match the unknown additional terms of Equations (21) and (22) to an eddy viscosity and an eddy thermal dissipation. Such dissipation terms can be written as


where and respectively depend on an eddy viscosity and an eddy thermal dissipation coefficient (Equations 4 and 5). Using the results from EULAG simulation, we calculate Equations (23) and (24) neglecting at this point the dependency upon radius and scale of the unknown dissipation coefficients and (set to arbitrary values , ). Our procedure to obtain effective dissipation coefficient is as follows. We average in time the right hand side of Equation (21) (the left hand side is vanishingly small because the system has reached a statistical steady-state), and divide it with the analogously averaged Equation (23). The resulting ratio gives us naturally . The exact same procedure is applied to Equations (22) and (24) to obtain . We display the resulting eddy-diffusion coefficients and in the left panels of Figure 5 for case E1.

Our numerical procedure embeds numerous numerical approximations that need to be acknowledged. First, we have no formal proof that the numerical sub-grid scales effects can be matched to an eddy-type dissipation in the case of turbulent convection simulations with EULAG. Second, EULAG is formulated on a regular cartesian grid which is mapped to the spherical geometry. Here, we carry our analysis on a spherical harmonics decomposition that is obtained through an interpolation of EULAG’s results on Gauss-Legendre collocation points in latitude. Hence, some numerical errors can appear in our analysis since this does not reflect directly how Equations (1)-(3) are solved in EULAG. Third, the spherical harmonics decomposition is cut at the maximum accessible spherical harmonics degree , corresponding to the smallest grid size in latitude in EULAG. We repeated our analysis by de-aliasing the spherical harmonics decomposition up to , which did not significantly change our results, giving us confidence that the second and third error sources are not significantly affecting our study. Finally, because we cannot expect the sub-grid scale model to behave entirely as an eddy-type dissipation operator, we consider a criterion to assess the robustness of the effective eddy-dissipation coefficients we obtained. For each value of (and ), we calculate the standard deviation along the time window over which we averaged the various contributions shown in Figure 4. If the standard deviation is of the order of the deduced eddy-dissipation coefficient, we do not confidently rely on its value. In the left panels of Figure 5 we darken the pixels of the colormaps accordingly, and only consider the bright pixels in the following analysis.

Figure 5: Left panels: Effective visocity (upper panel), heat dissipation coefficient (middle panel) and Prandtl number (lower panel), as a function of spherical radius and spherical harmonics degree . The colormap in each panel is displayed on a logarithmic scale. In each panel, gray pixels label non-robust effective dissipation coefficients (see text). Note that in the lower panel, white areas correspond to an effective Prandtl number of unity. Right panels: Averaged, smoothed profiles (blue lines) as a function of spherical radius (left) and spherical harmonics degree (right). The gray area corresponds to the standard deviation of the effective dissipation coefficients around the average value. The red lines denote the fitted dissipation coefficients that are used in ASH simulations in Section  4.3. In the lower panel, the Prandtl number of unity is labeled by the horizontal dashed black line.

The robust eddy-diffusion coefficients show characteristic patterns in and , which are shown in the right panels of Figure 5. The eddy viscosity is maximized near the radial boundaries (which are impenetrable), and increases with . Its average value lies between cm/s and cm/s which agrees with empirically deduced implicit eddy-diffusion coefficients obtained with mean-field models representing EULAG simulations (see, e.g., Simard et al. 2016). The average profile in and are shown in blue, and the standard deviation along and is denoted by the gray areas.

The eddy heat dissipation coefficient shows a reversed profile. It is maximized at the center of the convection zone and decreases strongly (by almost a factor of 8) close to the radial boundaries. In the robust part of the profile () it decreases with .

The lower panels show the effective Prandtl number . It is generally assumed in ILES that the Prandtl number is close to . Here we find that it is indeed very close to unity in the middle of the convection zone ( corresponds to white in the colormap), but rapidly increases by almost an order of magnitude near the radial boundaries. The Prandtl number also slightly increases with , due to the simultaneous increase of and decrease of .

We performed the same analysis for case E2, which is not shown here. The radial profile of the dissipation coefficients is very similar to the one shown in Figure 5, and the profile is simply shifted to higher ’s thanks to the finer resolution. The average value of and is decreased by factor between and , which is compatible with the naive expectation of doubling the overall resolution (in both time and space) of the simulation, leading to an implicit dissipation reduced by a factor of 4.

4.3 Comparison with the ASH code

In order to further validate the deduced dissipation coefficients, we now use them in LES done with the ASH code. The ASH code allows to specify dissipation coefficient depending on both and , which enables to fully take into account the dissipation coefficients inverted in Section  4.2.

We fit the dissipation coefficients with the following formulation


The radial shape of the effective dissipation coefficients is fitted with a standard polynomial. The dependency against is fitted only in the range because the matched and are not statistically significant at low (see the discussion in previous section). We linearly fit both of them for the sake of simplicity.

At large scales, dissipation is small in EULAG and our method is not able to match EULAG’s dissipation to standard enhanced dissipation coefficients. In ASH, we hence arbitrarily set them to the small values cm/s and cm/s using a smooth hyperbolic tangent (red curves in Figure 5). Furthermore, an ASH simulation using these fitted coefficients is not stable when using the same grid resolution as in EULAG. Hence, we double the resolution of case E1 in ASH and arbitrarily increase the value of the fitted coefficients to cm/s and cm/s at small scales (large ), using again a smooth hyperbolic tangent. The fitted profiles used in ASH are shown in red in the middle and right panels of Figure 5, and can be written as (here, for the viscosity)


where the subscripts and respectively refer to the small and large scale branches of the profile, and the generic step function is defined by


where is centering parameter of the hyperbolic tangent.

The ASH results are sensitive to the arbitrary choices of effective dissipation coefficients at large and small scales. In order to illustrate their sensitivity, we show three different cases obtained with ASH when the transition at large scales is slightly shifted (see the solid, dashed and dotted red lines in Figure 5). We label the cases ’A’ where is the centering scale of the hyperbolic tangent at large scales. At small scales, the centering scale is set to for all cases. Figure 6 shows the resulting differential rotation, its radial profiles at various latitudes, and the radial velocity at mid depth, correspondingly for the three cases together with the case E1.

Figure 6: Upper panels Differential rotation for the three ASH cases (see red curves in Figure 5) and case E1. As in Figure 2, blue and red respectively denote slower and faster rotation than the stellar rotation rate shown in white. Middle panels Radial profile of the differential rotation at five different latitudes in the northern hemisphere. The profile at the equator is in blue, and at latitude in red (see color bar). Lower panels Radial velocity colormap at the middle of the convection zone for the ASH cases and case E1. Again, red denotes outward motions and blue downward motions.

The three cases exhibit a differential rotation pattern which qualitatively agrees with case E1. Overall, case A20 matches E1 best, albeit its differential rotation is stronger, and steeper at the equator (blue line). However, A20 better matches E1 at mid latitudes (yellow to red lines) where A14 and A17 evince too strong differential rotation.

The convective patterns (lower panels) change significantly between the three cases. In case A14, convection is very weak on an equatorial band and is concentrated at higher latitudes. In case A17 only an “active nest of convection” (see Brown et al. 2008) is observed on the equator, and in case A20 the convective patterns qualitatively match E1 with slightly stronger convective amplitudes. We recall that the only difference between cases A14 and A20 is the location of the transition to low dissipation coefficients at large scales, i.e. in case A20 a larger range of large scales is evolved with low dissipation coefficients. We recall that the kinetic energy spectrum peaks in between and in the ASH cases. The non-linear balance saturating the peak of the turbulent convective spectrum is thus altered from left to right in Figure 6, with higher viscous and heat dissipation in case A14 than in case A20 at the peak of the turbulent spectrum. In case A14, the radial shear of the differential rotation near the equatorial plane is strong enough to weaken significantly convection, whereas in case A20 the banana cells still live on.

In spite of these differences, the three ASH cases exhibit a similar convective heat transport (not shown) peaked at the center of the convection zone with an equivalent convective luminosity in cases A17 and A20, and in case A14. All three are close to the E1 case, for which the convective transport reaches (see Section 3).

Figure 7: Total kinetic energy (top panel) and entropy (bottom panel) spectra for cases A14 (dashed cyan), A17 (dotted magenta), A20 (plain blue) and E1 (plain red).

Interestingly, the large-scale non-axisymmetric entropy spectrum (Figure 7) is significantly larger in ASH cases compared to EULAG, while the kinetic energy spectrum at those scales remain the same. The slope of the kinetic energy spectrum at mid-scale in EULAG is well recovered by case A20, with a departure near the transition to high dissipation around . The change in the spectra from case A14 to case A20 denotes a non-trivial influence of scale-dependent dissipation coefficients on the correlations between the large-scale heat and momentum fluctuation. This may have important theoretical implications for the amplitude of the large-scale convective fluctuations in the Sun (Hanasoge et al. 2012; Greer et al. 2015; Lord et al. 2014), and will be studied in near future.

Finally, we recall that the viscous heating has been omitted in the governing equations, as in ILES it is generally not accounted for by numerics supplying small-scale diffusion, and hence should not be formally introduced in this ILES-LES comparison. Nevertheless, omitting the viscous heating could formally lead to net energy loss in simulations, we hence intend to explore this effect in a subsequent study.

Numerous numerical approximations have been made in (i) projecting the results from EULAG in spectral space, (ii) matching the imbalance in Equations (21) and (22) to standard laplacian dissipation terms, and (iii) using the matched dissipations coefficients in a LES with the ASH code, introducing arbitrary dissipation coefficients at very large and very small scales. In spite of these approximations, we managed to simulate with ASH (case A20), using fitted explicit dissipation coefficients, a convective state producing a large scale differential rotation that compares adequately with the results obtained with EULAG. This comparison gives us confidence in the methodology we developed to estimate the dissipation properties of ILES with EULAG, and a posteriori suggests that, at least in the bulk part of the turbulent spectrum, EULAG’s implicit dissipation could be interpreted in terms of standard explicit dissipation.

5 Conclusions

In this work we have studied the dissipation properties of implicit large-eddy simulations of an idealized turbulent stellar convection zone under the influence of rotation with the EULAG code. By considering twin simulations where the grid resolution was doubled, we showed that the kinetic energy spectrum of turbulent convection peaks at smaller scales in the most refined model, as expected. The latter model develops a stronger differential rotation, which results from the complex interplay between the turbulent scales, which are themselves affected by the implicit dissipation of EULAG.

In order to characterize this implicit dissipation, we developed a spectral method to a posteriori quantify the effect of the sub-grid scale modelling of EULAG. By evaluating balance equations for kinetic energy and entropy when the system has reached a steady-state, we were able to isolate the implicit dissipation contribution. The sub-grid scale modelling was shown to match quantitatively well a standard laplacian-like operator from medium to small scales. At large scales, the implicit dissipation introduced by EULAG is very weak and could not be matched to such classical formulation.

For a grid resolution comparable to previously published results with EULAG in the context of solar dynamo (Ghizaru et al. 2010; Racine et al. 2011; Beaudoin et al. 2013; Passos and Charbonneau 2014; Lawson et al. 2015), we find the effective viscosity and thermal diffusivities are of the order of cm/s. However, we recall that the effective dissipation does not match to an equivalent enhanced dissipation coefficient for the large scales feature such as differential rotation, for which we showed the effective dissipation to be a least an order of magnitude lower.

In order to further test our estimates of the dissipation coefficients, we used them in a series of LES with the ASH code. Note however that we were compelled to select arbitrary dissipation coefficients at very small scales (below the EULAG grid scale) and very large scales (where the implicit dissipation in EULAG is small and does not match a standard laplacian operator). We showed that our arbitrary choice of dissipation coefficients at large scale could significantly affect the large-scale flows in ASH. However, by choosing adequately these dissipation coefficients, we were able to reproduce the large scale differential rotation of EULAG simulation. Our results thus indicate that results from ILES and LES could be reconciled.

The spectral analysis developed in this paper is generic and can be easily extended to the MHD regime (Strugarek et al. 2013). A natural extension of this work is to perform such joint simulations of cyclic dynamos, in order to isolate to what extent the particular treatment of sub-grid scales chosen in EULAG simulations impacts the existence and characteristics of such solutions, as well as their dependency to rotation. We intend to explore this aspect in a future publication.


Constructive comments from two anonymous referees helped to improve the presentation. The authors thank J.F. Cossette for numerous discussions on convection modelling with EULAG, and N. Wedi and S. Malardel for discussions about the projection of the Coriolis force in spectral energy budgets. A. Strugarek is a National Postdoctoral Fellow at the Canadian Institute of Theoretical Astrophysics. The authors acknowledge support from Canada’s Natural Sciences and Engineering Research Council. P. K. Smolarkiewicz is supported by funding received from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2012/ERC Grant agreement no. 320375). This work was also supported by the INSU/PNST, the ANR 2011 Blanc Toupies, and the ERC grant STARS2 207430. S. Mathis acknowledges funding by the European Research Council through ERC grant SPIRE 647383. We acknowledge access to supercomputers through GENCI (project 1623), Prace (8th call), and ComputeCanada infrastructures.



Appendix A Spectral transfers equations

a.1 Vectorial spherical harmonics basis

a.1.1 The standard vectorial basis

We define from Rieutord (1987); Mathis and Zahn (2005):


where defines the spherical basis and are the spherical harmonics defined by


where are the associated Legendre polynomials. The basis (29) have the following properties :


where the solid angle, means complex conjugate and is the Kronecker symbol. We also have:


and all the other scalar cross products are .

a.1.2 Scalar fields identities

Defining , we get:


where .

a.1.3 Vectorial fields identities

For a vector

we obtain:


a.1.4 An alternative vectorial basis

The vectorial spherical harmonics basis defined in appendix A.1 is very efficient to calculate scalar products or linear differential operator on vectors. Nevertheless, it is quite hard to use it to express vectorial products. Instead we define the following basis (e.g., see Varshalovich et al. 1988):




is the -j Wigner coefficient linked to Clebsch-Gordan coefficients, and the vectors are


where defines the cartesian basis. Note that the equivalent of the conjugation rule (33) is then


a.1.5 Vectorial product

As previously, we decompose a vector on this basis:

Evaluating the vectorial product of two vectors , one gets:




with being the -j Wigner coefficient.

a.1.6 Scalar product

We decompose a vector on this basis in the following way:

Evaluating the scalar product of two vectors , one gets:


where the sum symbol is the same as in Equation (43) and


with being here the -j Wigner coefficient.

a.1.7 Basis change relations

For a vector decomposed in the following manner:

we have the two following relations to change from one basis to the other:


a.2 Kinetic Energy

We start from the momentum equation:


We define the kinetic energy density spectrum by . We multiply Equation (48) by and integrate it over the spherical shell to obtain:


where the various terms are given by


The pressure gradient contribution is easily calculated with Equation (34) and the buoyancy contribution is also straightforward to evaluate. We now detail the three remaining contributions.

a.2.1 Advection

We know that


As a result, one can simply use the standard formulae to evaluate the scalar product (A.1.6) and the vectorial product (A.1.5) part of the preceding equation to compute the full advection term 53.

a.2.2 Coriolis force

We recall the reader that the Coriolis force contribution vanishes for the total energy. It is able though to redistribute energy spectrally among scales, as we show now. By definition we have

Writing , we note that:

The scalar product with is then trivial to compute. This Coriolis contribution is equivalent to Equations (24-25) in Augier and Lindborg (2013), where a spectral kinetic energy equation is derived in the context of General Circulation Models.

a.2.3 Viscous tensor

From the definition

and the decomposition

one may rewrite the components of the symmetric tensor in the following fashion: