The Spectral Amplitude of Stellar Convection and its Scaling in the High-Rayleigh-Number Regime
aa r X i v : . [ a s t r o - ph . S R ] N ov The Spectral Amplitude of Stellar Convection and its Scaling inthe High-Rayleigh-Number Regime
Nicholas A. FeatherstoneDepartment of Applied Mathematics, University of Colorado, Boulder, CO 80309-0526 [email protected]
Bradley W. HindmanJILA & Department of Astrophysical and Planetary Sciences, University of Colorado,Boulder, CO 80309-0440Received ; accepted 2 –
ABSTRACT
Convection plays a central role in the dynamics of any stellar interior, andyet its operation remains largely-hidden from direct observation. As a result,much of our understanding concerning stellar convection necessarily derives fromtheoretical and computational models. The Sun is, however, exceptional in thatregard. The wealth of observational data afforded by its proximity provides aunique testbed for comparing convection models against observations. Whensuch comparisons are carried out, surprising inconsistencies between those mod-els and observations become apparent. Both photospheric and helioseismic mea-surements suggest that convection simulations may overestimate convective flowspeeds on large spatial scales. Moreover, many solar convection simulations havedifficulty reproducing the observed solar differential rotation due to this appar-ent overestimation. We present a series of 3-dimensional (3-D) stellar convectionsimulations designed to examine how the amplitude and spectral distributionof convective flows are established within a star’s interior. While these simula-tions are non-magnetic and non-rotating in nature, they demonstrate two robustphenomena. When run with sufficiently high Rayleigh number, the integratedkinetic energy of the convection becomes effectively independent of thermal diffu-sion, but the spectral distribution of that kinetic energy remains sensitive to bothof these quantities. A simulation that has converged to a diffusion-independentvalue of kinetic energy will divide that energy between spatial scales such thatlow-wavenumber power is overestimated, and high-wavenumber power is under-estimated relative to a comparable system possessing higher Rayleigh number.We discuss the implications of these results in light of the current inconsistenciesbetween models and observations. 3 –
Subject headings:
Stars: kinematics and dynamics, Sun: helioseismology, Sun:interior, Sun: magnetic fields, Stars: interior, Stars: fundamental parameters
1. Introduction
Convection is a pervasive phenomenon within stars and serves as the principle transportmechanism by which fusion-generated energy is transmitted across significant portions ofthe stellar interior. This transport may occur primarily within a convective core, as in thecase of massive stars, or within an outer convective envelope, as in cool stars such as theSun. As that convection takes place within a rotating, electrically conductive plasma, itundoubtedly plays an active role in the generation of stellar magnetism. More specifically,rotation imbues convective motions with helicity which may serve to amplify magnetic fieldsthrough the so-called α -effect (e.g., Moffat 1978). Further, the propensity of convectiveflows to redistribute angular momentum—thereby driving a differential rotation—may alsorender convective motion central to the generation of coherent, global-scale magnetic fieldsthrough the Ω-effect (i.e., large-scale rotational shearing of the field).Convection in the stellar context is thought to be highly turbulent; typical estimatesfor the Reynolds and Rayleigh numbers are 10 and 10 respectively (e.g., Ossendrijver2003). These extreme parameters make stellar convection challenging to study, whetherexperimentally or theoretically. This challenge is compounded by the fact that suchconvection occurs within a rotating, magnetized environment.Investigations into stellar convection are most readily initiated within the context ofthe Sun. The wealth of observational data afforded by the proximity of that star make itan ideal test bed against which predictions of stellar convection models may be tested andcalibrated before extending those models to other stars. Despite the abundance of solardata, the convective motions central to its dynamo remain largely hidden from view inall but the shallowest layers of the convection zone. Whatever its detailed structure, thatconvection must satisfy two robust observational constraints: it must transport one solarluminosity’s worth of energy across the convection zone, and it must efficiently redistribute 5 –angular momentum in order to sustain the Sun’s helioseismically inferred differentialrotation.Exploring the magnetohydrodynamical balances at work in stellar interiors has largelybeen the purview of three-dimensional (3-D), nonlinear, numerical convection simulations.In the case of the Sun, these simulations are often run by assigning solar values to mostmodel variables, notably the rotation rate and the luminosity, while admitting diffusivitiesthat are computationally tractable but also necessarily orders of magnitude larger thantheir intended microscopic counterparts (e.g., Brun et al. 2004; Ghizaru et al. 2010; Hottaet al. 2014). It is now becoming increasingly evident that such quasi-laminar simulationshave succeeded in maintaining a solar-like differential rotation and energy flux in part dueto the out-sized role played by diffusion in many of those models.As advances in computational power have enabled the use of lower and ostensibly morestellar-like diffusivities, inconsistencies between models and observations have begun toarise. Notable among these are the so-called anti-solar states of differential rotation whichexhibit an equator of retrograde rotation and polar regions of prograde rotation. Such statesdevelop when the characteristic convective turnover time drops below some critical fractionof the rotation period (e.g., Gilman 1977; Glatzmaier & Gilman 1982). This transitionoccurs naturally in solar simulations as the level of diffusion is decreased (Featherstone &Miesch 2015), and it tends to arise even in those models that do not employ an explicit,Fickian diffusivity (Hotta et al. 2014). While almost certainly relevant within a broaderstellar context (see e.g., K˝ov´ari et al. 2007), these anti-solar states stand at odds with thewell-established internal rotation profile of the Sun (e.g., Thompson et al. 2003).Further inconsistencies arise when the distribution of velocity power among differentspatial scales is examined in detail. Photospheric observations and helioseismic observationsof the deeper interior both suggest that numerical simulations may over-estimate the 6 –amplitude of the solar convection on large spatial scales (Hanasoge et al. 2012; Lord etal. 2014). The helioseismic consensus concerning deep flows remains murky, however.Recent ring-analysis measurements of subsurface flows throughout the near-surface sheerlayer exhibit good agreement with convection models (Greer et al. 2015). While thequestion of subsurface flow speeds remains an open question, we note that the time-distancemeasurements of Hanasoge et al. (2012), if confirmed, are particularly troubling as theysuggest that convection on spatial scales larger than 70 Mm can be at most 5–6 m s − .This stands in contrast to theoretical estimates based on gyroscopic pumping which suggestthat convective flows must be at least 30 m s − on those spatial scales (Miesch et al. 2012).Weaker flows would have difficulty generating Reynolds stresses which are strong enough tomaintain the solar differential rotation.These theoretical and observational disagreements suggest that something fundamental,related to the transport of heat and angular momentum by magnetized convection, may bemissing in stellar convection models. Until these issues can be resolved in the context of theSun, for which we have such an abundance of observational data, the relationship of any stellar convection model to an actual star remains highly questionable.We present a set of numerical experiments designed to examine how the structure andamplitude of the convective flow field is determined within a stellar interior. Specifically,we ask, “How is the amplitude and spectral-power distribution of stratified convectionrelated to the stellar luminosity and level of diffusion present in the system.” For thispaper, we address that question without considering the effects of rotation and magnetismwhich, while admittedly important, would substantially extend the parameter space of thisalready challenging study. We choose instead to focus our initial study on how densitystratification, luminosity, and diffusivity relate to the amplitude and spatial power spectrumof the resulting convection. Our results demonstrate that the convective kinetic energy is 7 –related to these three parameters through an unambiguous scaling relationship. In addition,we demonstrate that models run in the typical stellar-convection fashion may naturallyoverestimate low-wavenumber power in the convective power spectrum.This paper is organized as follows. In § §
3, we discuss our numerical model andthe parameter space explored in our numerical experiments. We then discuss our results in §
4, followed by a discussion of their implications in §
2. Numerical Model
This study is based around a series of 3-D, nonlinear convection models run usingthe
Rayleigh convection code, developed as a community code for the ComputationalInfrastructure for Geodynamics (CIG).
Rayleigh simulates convection using the spectraltransform approach described in Glatzmaier (1984). We employ a spherical geometry andrepresent the horizontal variation of all variables along spherical surfaces using sphericalharmonics Y mℓ ( θ, φ ). Here, ℓ is the spherical harmonic degree, and m is the azimuthalmode order. In the radial direction, we employ a Chebyshev collocation method, expandingall variables using Chebyshev polynomials T n ( r ), where n is the degree of the Chebyshevpolynomial. Derivatives are calculated accurately within the spectral representation byexploiting the properties of these two sets of basis functions. We work on a de-aliasedgrid, such that the number of collocation points in each dimension is larger than thecorresponding number of spectral modes by a factor of 3/2.As our approach is pseudo-spectral in nature, nonlinear terms are calculated inphysical space after first transforming the relevant variables from the spectral configuration.Time-integration is carried out in the spectral configuration using a hybrid implicit-explicitapproach. A Crank-Nicolson method is used for linear terms, and an Adams-Bashforth 8 –approach is employed for the nonlinear terms. Both components of the time-stepping possess2nd-order accuracy. Algorithmically, it is the parallelization of Rayleigh , which parallelizesefficiently on up to O(10 ) cores, that differs significantly from that of Glatzmaier (1984).We have verified our particular implementation of that approach against two establishedbenchmark tests, and we provide results from those tests in Appendix A.Our study is concerned with convection as it manifests deep within stellar interiors,far removed from the photospheric surface where radiative processes may contributeconsiderably to the energetics of the convection. In such a region of the star, plasmamotions are subsonic and perturbations to thermodynamic variables are small comparedto their mean, horizontally-averaged values. Under such conditions, the anelasticapproximation, which we choose to employ, provides a convenient means of describing thesystem’s thermodynamics (Gough 1969; Gilman & Glatzmaier 1981).When using this approximation, thermodynamic variables are linearized about aspherically symmetric, time-independent reference state with density ¯ ρ , pressure ¯ P ,temperature ¯ T , and specific entropy ¯ S . Fluctuations about this state are denoted, withoutoverbars ( ρ , P , T , and S ). A further consequence of the anelastic approximation is that themass flux is solenoidal, reducing the continuity equation to ∇ · ( ¯ ρ v ) = 0 , (1)where v = ( v r , v θ , v φ ) is the velocity vector expressed in spherical coordinates. The lack ofany time derivative in Equation (1) means that sounds waves are naturally filtered out as aconsequence of this approach. The divergence-free constraint for the mass flux is enforcedby projecting v onto poloidal and toroidal streamfunctions ( W and Z respectively), suchthat ¯ ρ v = ∇ × ∇ × ( W e r ) + ∇ × ( Z e r ) . (2)The unit vector in the radial direction is indicated by e r . The momentum equation is given 9 –by ¯ ρ D v Dt = − ¯ ρ ∇ P ¯ ρ − ¯ ρSc p g − ∇ · D , (3)where g is the gravitational acceleration, and c p is the specific heat at constant pressure. Inwriting the momentum equation this way, we have employed the so-called Lantz-Braginsky-Roberts approximation, which is exact for adiabatic references states such as those employedin this study (Lantz 1992; Braginsky & Roberts 1995). The viscous stress tensor D is givenby D ij = − ρν (cid:20) e ij −
13 ( ∇ · v ) δ ij (cid:21) , (4)where e ij is the strain rate tensor. The kinematic viscosity is denoted by ν , and δ ij is theKronecker delta. Our thermal energy equation is given by¯ ρ ¯ T DSDt = ∇ · [ κ ¯ ρ ¯ T ∇ S ] + 2 ¯ ρν (cid:20) e ij e ij −
13 ( ∇ · v ) (cid:21) + Q, (5)where the thermal diffusivity is denoted by κ . Sources of internal heating and cooling areencapsulated by the functional form of Q . A linearized equation of state closes our set ofequations and is given by ρ ¯ ρ = P ¯ P − T ¯ T = Pγ ¯ P − Sc p , (6)assuming the ideal gas law ¯ P = R ¯ ρ ¯ T . (7)The specific heat at constant pressure is denoted by c p , R is the gas constant, and γ is theadiabatic index of the gas.
3. The Numerical Experiment
We have constructed a series of 63 model stellar convection zones designed to explorehow the convective kinetic energy depends on three model parameters: the degree of density 10 –
TABLE 1Variable Input ParametersParam Value M i × g ρ i × − g cm − c p × n N ρ r i × cm r o × cm Note. — Polytropic reference state parameters for all cases. The reference state prescription differsbetween simulations in the value of N ρ chosen only. stratification, level of thermal diffusion, and the stellar luminosity. A detailed listing of allmodel parameters may be found in Appendix B.Each of our models is constructed using a polytropic background state following Joneset al. (2011). This approach has the advantage that the thermodynamic background maybe specified analytically, making it easily reproducible. Our background may be completelyspecified by seven numbers: the inner radius of the shell r i , the outer radius of the shell r o , the polytropic index n , the number of density scale heights occurring within the shell N ρ , the mass interior to the shell M i , the density at the inner boundary ρ i , and the specificheat c p . With the exception of N ρ , these values are the same for all experiments reportedon here, and we list them in Table 1. All background states are adiabatically stratified,possessing a polytropic index n of 1.5 and an adiabatic index γ of 5/3. Further details ofthe polytropic background are provided in Appendix C for completeness. 11 – sun D en s i t y ( g c m - ) N ρ = 4N ρ = 3N ρ = 2N ρ = 1 0.75 0.80 0.85 0.90r / R sun π r F / L s t a r a b Fig. 1.—
Polytropic reference states used in this study. ( a ) Radial variation of density corresponding toeach value of N ρ used in our simulations (indicated by line coloring). The N ρ = 3 cases closely resemble theModel-S density stratification (black diamonds; Christensen-Dalsgaard et al. 1996). ( b ) Radial variation ofthermal energy flux that must be transported by convection and thermal conduction for each value of N ρ (coloring as in panel a ). Fluxes have been normalized by the stellar flux. The flux corresponding to ModelS is indicated by the dashed black line. We have chosen to examine a range of density stratifications ( N ρ ) for this study. Oneconsequence of this approach is that shells with a higher degree of density stratificationpossess a lower mass than their more weakly stratified counterparts. We note that whenconstructed using the parameters from Table 1, the thermodynamic background stateclosely resembles that of the Sun for the value N ρ = 3. We plot the density variation withradius for each of our four chosen values of N ρ , along with that from a standard solar model(Model S; Christensen-Dalsgaard et al. 1996), in Figure 1 a .For all simulations, we have adopted impenetrable and stress-free boundary conditionssuch that v r ( r = r i , r o ) = ∂ ( v θ /r ) ∂r (cid:12)(cid:12)(cid:12)(cid:12) r = r i ,r o = ∂ ( v φ /r ) ∂r (cid:12)(cid:12)(cid:12)(cid:12) r = r i ,r o = 0 . (8)The radial entropy gradient is forced to vanish at the lower boundary of the convectionzone, and the entropy perturbations are forced to vanish at the upper boundary, with ∂S∂r (cid:12)(cid:12)(cid:12)(cid:12) r = r i = 0 , S ( r o ) = 0 . (9) 12 –Thus, there is no diffusive entropy flux across the lower boundary. Instead, heat enters thesystem through internal deposition by Q , which drops to zero at the upper boundary. In allsimulations, we adopt a functional form of Q that depends only on the background pressureprofile such that Q ( r, θ, φ ) = α ( ¯ P ( r ) − ¯ P ( r o )) . (10)The normalization constant α is chosen so that L ⋆ = 4 π Z rori Q ( r ) r dr, (11)where L ⋆ is the stellar luminosity. The thermal energy flux F ( r ) that convection andconduction must transport across a spherical surface at radius r is then given by F ( r ) = 1 r Z rri Q ( x ) x dx. (12)The value of L ⋆ is varied between simulations, as is the value of N ρ , and that of the diffusioncoefficients ν and κ . In all instances, we have adopted a Prandlt number P r = ν/κ of unity.The functional form of Q means that as the value of N ρ is increased, the heatingbecomes more focused near the lower boundary. This behavior is illustrated in Figure 1 b ,where we have plotted F ( r ), normalized by the stellar flux at each radius. Our prescriptionfor the internal heating in the N ρ = 3 cases resembles that derived using solar temperatureand opacity values tabulated in Model S (dashed line). We note that small differences,comparable to those observed between reference states with differing values of N ρ , do existbetween the N ρ = 3 case and the solar model. Those differences arise primarily from thefact that we have chosen the base of our model convection zones to lie at a radius of 5 × cm, somewhat higher than that of the Sun. The results discussed in § κ , and the chosenluminosity. Specifically, the time-averaged value of the entropy gradient must equilibrate to ∂S∂r (cid:12)(cid:12)(cid:12)(cid:12) r = r o = L ⋆ πr o κρ ( r o ) T ( r o ) (13) All simulations were run in a manner analogous to one commonly used in solarconvection studies (e.g., Brun et al. 2004; Ghizaru et al. 2010; Hotta et al. 2014,Featherstone & Miesch 2015). Namely, we have run these simulations using dimensional parameters for the equation coefficients, and we assign “stellar” values to those parameters,with the notable exception of the diffusivities. The values chosen for our diffusioncoefficients are instead motivated purely by matters of computational feasibility; thus theyare necessarily orders of magnitude larger than the microscopic values. It is customary tosuggest that this prescription for the diffusivity constitutes a form of subgrid-scale model,representing the mixing of heat and momentum by small-scale eddy motions that cannot besufficiently resolved on the coarse computational grid. We make no such assertion here, andremark simply that our diffusivities are rather large when compared to stellar values. Thediffusion coefficients ν and κ are taken to be constant values within each of our simulations.In particular, they possess no variation in radius.Across our set of numerical convection experiments, we have varied the degree ofthermal stratification (characterized by N ρ ), the luminosity of the star L ⋆ , and the thermaldiffusivity κ . This three-dimensional parameter space may be collapsed by considering anappropriate non-dimensionalization of the system. As we have chosen a P r of unity andhave neglected the effects of rotation and magnetism, all runs may be well-characterized bya single non-dimensional number, the Rayleigh number. 14 –For Boussinesq convection, the Rayleigh number is typically defined in terms of thetemperature, but in our anelastic formulation, the thermal properties of the fluid are insteadconveyed by entropy. Our simulations also differ from the classical Rayleigh-B´enard setup,which possesses fixed-temperature boundaries, in that we have effectively imposed a flux byspecifying Q and the entropy boundary conditions of Equation (9). It is thus appropriatefor our system to define a flux Rayleigh number in terms of entropy. Throughout this paper,the term Rayleigh number will be used to refer to a flux
Rayleigh number Ra F , which wedefine as Ra F = e g e F H c p e ρ e T νκ . (14)We choose to use volume averages for all quantities possessing a tilde in the definitionof Ra F , making it a bulk Rayleigh number. Here H is a typical length scale, and wechoose the shell depth r o − r i as its value. Our numerical experiments span the range10 ≤ Ra F ≤ × . Each simulation was initialized using a small random thermalperturbation, evolved until the kinetic energy reached a statistically steady state, andfurther evolved for at least one thermal diffusion time further following the onset of thatstatistically steady state.
4. Survey of Results4.1. Kinetic-Energy Scaling and Spectral Distribution
We begin our discussion of the convective energetics by looking at the integratedkinetic energy KE , defined as KE = 12 Z V ¯ ρ ( r ) | v ( r, θ, φ ) | dV. (15)The relationship of KE to thermal diffusivity is illustrated in Figure 2 a . There, we haveplotted KE vs. κ − , using colored symbols to distinguish between the level of stratification 15 – -13 -12 κ -1 (cm -2 s)10 K i ne t i c E ne r g y ( e r g ) N ρ = 4N ρ = 3N ρ = 2N ρ = 1 10 Ra F -1 N on - d i m en s i ona l KE a b low-Ra F regime high-Ra F regime Fig. 2.—
Dimensional and non-dimensional view of kinetic energy scaling. ( a ) Dimensional kinetic energy KE vs. κ − for all cases run with L ⋆ = L ⊙ . Colored symbols indicate the level of density stratification.At sufficiently low κ , KE approaches a constant value within each N ρ series. ( b ) Non-dimensional kineticenergy d KE vs. flux Rayleigh number Ra F for all cases (regardless of the value of L ⋆ ). Colored symbolsindicate N ρ as in panel a . When viewed non-dimensionally, all N ρ series collapse onto a single curve andapproach an asymptotic scaling law with d KE ∝ Ra . F (dashed line). The dilineation between high- andlow- Ra F regimes is indicated by the dash-dotted line. employed in each model. Each data point represents an average of 3000 days of simulationtime. This is approximately one thermal diffusion time across the shell for those cases with κ ≥ × cm s − . Only two of our models possess a lower value of κ (5 × cm s − ), but this averaging interval still encompasses several tens of convective overturningtimes. We have also omitted those cases run with a non-solar value of L ⋆ in this figure sothat the general trend may be clearly seen. For each set of curves, the integrated KE risesuntil some critical value of κ is reached. Below this critical value of κ , KE asymptotes toa nearly constant value that is essentially independent of the level of diffusion κ . Thesesimulations have thus reached a state of so-called turbulent free-fall (e.g., Spiegel 1971),which we discuss further in § κ data in Figure 2 a yields therelationship KE ∝ κ − . ± . , (16)which is only weakly dependent on κ . 16 –The point of transition to this asymptotic regime depends on N ρ . More weakly statifiedcases transition to the asymptotic regime at a lower value of κ than their more stronglystratified counterparts. In addition, the saturation level of KE may be weakly dependenton N ρ , as evinced by the N ρ = 1 curve (red). It is possible, however, that those weaklystratified cases are simply converging more slowly to the same value as the N ρ > N ρ > KE .The similarity of behavior in KE between systems with different N ρ is made morereadily apparent by examining how KE depends on the Rayleigh number Ra F . In so doing,we choose a non-dimensional measure of the kinetic energy d KE such that d KE ≡ H κ M KE. (17)The nondimensionalization has been carried out using the depth of the domain, the thermaldiffusion timescale, and the mass M contained within the spherical shell (which is a functionof N ρ ). In Figure 2 b , we have plotted d KE vs. Ra F . When plotted in this fashion, all resultsfall along a similar curve that describes a power-law relationship beween KE and Ra F .Performing a least-squares fit to the high- Ra F portion of this plot ( Ra F ≥ ), we findthat d KE ∝ Ra . ± . F . (18)The measured value of the exponent is very close to . This value indicates that beyond Ra F of 10 , which we designate as the “high- Ra F ” region of this curve, diffusion is nolonger playing a substantial a role in determining the convective amplitudes.While achieving a truly “stellar” value for the Rayleigh number is clearly impossiblein a modern convection simulation, these results nevertheless imply a clear prerequisite forcomparing any convection simulation against stellar observations. Namely, a simulation’sparameters should place it well within the high- Ra F region of Figure 2 b beforeany comparison is made against actual data.
17 –
Upper Convection Zone V e l o c i t y P o w e r ( a ) Middle Convection Zone ( b ) Lower Convection Zone N ρ = ( c ) V e l o c i t y P o w e r ( d ) ( e ) N ρ = ( f ) V e l o c i t y P o w e r ( g ) ( h ) N ρ = ( i ) l V e l o c i t y P o w e r High Ra F Low Ra F ( j ) l ( k ) l N ρ = ( l ) Fig. 3.—
Time-averaged velocity power spectra (dimensional) for all cases run using L ⋆ = L ⊙ . Each rowcorresponds to a single value of N ρ (indicated on the right), and each column corresponds to a single depthwithin the convective shell. Within each panel, spectra for all cases at that depth and N ρ are displayed.Low- Ra F cases are indicated in blue tones, and high- Ra F cases in red tones. The highest- Ra F case withineach series is further indicated by the dashed, dark red line. As Ra F is increased, power at low ℓ -valuesincreases initially. At high Ra F , it decreases as high-wavenumber power is generated at the expense oflow-wavenumber power.
18 –This result has interesting consequences for the relative spectral distribution of kineticenergy between high- and low- Ra F systems. In Figure 3, we have plotted the kinetic-energyspectra for those cases in our study that possess the solar luminosity. Spectra are plottedat three depths in each simulation (1 depth per column). These depths are taken nearthe upper boundary, the mid-shell, and near the lower boundary. Different values of N ρ correspond to different rows. Within each panel, all spectra for high- Ra F runs at that N ρ have been plotted. Low values of Ra F , indicated by dark blue tones, have been plottedfor the N ρ = 1 cases in order to illustrate the trend for that region of parameter space.They have been largely omitted for other values of N ρ to enhance legibility of the high- Ra F behavior. The highest-Rayleigh-number spectrum in each panel is indicated by a dashed,dark red line. Each spectrum represents a time-average over approximately 1,000 days ofsimulated time, corresponding to several tens of convective overturnings.A systematic trend is visible across all cases plotted here. As Ra F increases froman initially very low value, the integral of the kinetic-energy spectrum increases, as doesthe point-wise amplitude at roughly all wavenumbers. This trend is commensurate withthe low- Ra F /high- κ behavior depicted in Figures 2 a, b . At sufficiently high Rayleighnumber, however, power in the high- ℓ portion of the spectrum continues to increase, butthe low-wavenumber end begins to diminish in power. The spectrum with the highest Ra F in each panel (red dashed line) thus has less power at large scales than its low- Ra F counterparts. The “low-wavenumber” end of the spectrum appears to occur in theapproximate range 0 ≤ ℓ ≤
10 for most simulations in this study.This phenomenon arises naturally from the fact that at sufficiently high Rayleighnumber (and sufficiently low κ ), the integrated kinetic energy becomes essentiallyindependent of the degree of diffusion. At the same time, as κ is decreased and Ra F isincreased, the flow field becomes more turbulent and develops structure on increasingly 19 –smaller scales. Because the integral of the power spectrum is conserved, high-wavenumberpower can only come at the expense of low-wavenumber power. This generaltrend is independent of depth and appears so be somewhat more pronounced for systemspossessing larger values of N ρ .Before proceeding, we wish to emphasize that the results encapsulated by Figures 2and 3 constitute the main points of this paper. In summary:1. When a global stellar convection simulation is not run with sufficiently high Rayleighnumber, its resulting convective flow speeds will depend on the (unphysically high)value of κ , in addition to stellar parameters such as mass and luminosity. A simulationthat is not in the high- Ra F regime of Figure 2 b can be expected to yield rms convective amplitudes which are systematically lower than those expected in an actualstar.2. For sufficiently turbulent (high- Ra F ) convection simulations, this problem does notarise. The time-averaged, integrated kinetic energy approaches a constant value,independent of the level of diffusion employed.3. Even in this asymptotic regime, models that correctly capture the integrated kineticenergy may naturally overestimate low-wavenumber power because their Rayleighnumber is still orders of magnitude smaller than that of an actual star. Thus,simulated convective features that derive from large-scale motions may bear littlerelation to their intended stellar counterparts.Changes in Rayleigh number also impact the spatial distribution of kinetic energy.This is most clear when considering the variation of radial velocity with depth. We haveplotted the normalized rms vertical velocity V r in Figure 4 a . We define this quantity as 20 – Radius (cm)0.00.20.40.60.81.01.2 r m s V r ( c m s - ) Ra F = 2.8 x 10 Ra F = 5.5 x 10 Ra F = 2.0 x 10 Ra F δ V N ρ = 1 N ρ = 2 N ρ = 3 N ρ = 4 a b Fig. 4.—
Near-surface and deep-convection-zone flow speed difference as a function of Rayleigh numberand N ρ . ( a ) Normalized rms radial velocity V r as a function of depth for N ρ = 4 and three values of Ra F (indicated by colored lines). The difference in convective flow speed between mid and upper convection zonebecomes more pronounced with increasing Ra F . ( b ) Relative speed difference δV for all cases plotted as afunction of Rayleigh number and N ρ (indicated by the colored symbols). δV increases with both Ra F and N ρ but is most sensitive to N ρ . V r ( r ) ≡ u rms ( r )max u rms , (19)where u rms is the rms vertical velocity at each radius. The mean in the rms is spatial andtemporal; it has been taken over spherical shells and over 1,000 days of simulation time.As the Rayleigh number is increased, the variation in radial velocity across the domainincreases and tends to peak near the upper boundary. We suggest that this behavior isthe likely result of plumes being launched from the upper boundary layer whose thicknessand gradient are both changing with Rayleigh number. At high Ra F , the boundary layeris characterized by a sharp entropy gradient, and is thus amenable to the generation ofrapid, small-scale downflow plumes. The location and width of the velocity peak are thusan indirect measure of the effective penetration depth of these plumes.This effect is most pronounced at high Ra F and high N ρ , and we may quantify it by 21 –examing the relative velocity jump δV , which we define as δV ≡ max u rms − u rms ( r mid ) u rms ( r mid ) , (20)where r mid is the radial depth halfway between r i and r o . The relative velocity jumpmeasured for each case is plotted in Figure 4 b . The trend present in panel a of that figure ismore clearly illustrated here. δV increases with Ra F for all values of N ρ studied, reachingthe highest values for the N ρ = 4 cases. We note that it is due to the logarithmic axes that δV appears to be asymptoting to a constant value within each N ρ series. While the growthof δV is slowing at high Ra F , this quantity is still clearly increasing within each series.We note further that we have not succeeded in identifying a clear scaling law linking thebehavior of δV to the value of Ra F . Nevertheless, the apparent trend suggests that a starsuch as the Sun, whose convection zone is characterized by much higher values of both N ρ and Ra F , should possess a δV greater than 4. This behavior has implications for the solarnear-surface shear layer, as we discuss in § Rotation and magnetism aside, the convective dynamics realized in our simulationswill likely differ significantly from those found within a star due to both the large valuesof diffusivity we employ and the conductive boundary layer that develops near the top ofthe domain. We have seen that diffusion plays a minimal role in the determination of theconvective kinetic energy, and we now seek to quantify its role in the energy transport acrossthe layer. The thermo-mechanical transport of energy in our system may be characterizedby three radial energy fluxes, namely the enthalpy flux F e , the kinetic energy flux F KE , andthe conductive flux F c , which we define as 22 – F e = ¯ ρc p h v r T i , (21) F KE = 12 ¯ ρ (cid:10) v r | v | (cid:11) , (22)and F c = κ ¯ ρ ¯ T (cid:28) ∂S∂r (cid:29) (23)respectively. We consider averages of these quantities taken over 1,000 days of simulationtime and over spherical surfaces, denoting that average by angular brackets. Thus F e , F KE ,and F c are a function of radius only. At each radius, the sum of these three fluxes mustequate to the net flux F , established by the internal heating Q (Equation 12), such that F ( r ) = F e ( r ) + F KE ( r ) + F c ( r ) . (24)We may now quantify the contribution of conduction by considering the fractional convectiveflux f conv , which we define as f conv ≡ R V F e + F KE dV R V F dV = 1 − R V F cond dV R V F dV . (25)A value of zero for f conv indicates a lack of convective heat transport. A value of f conv near unity indicates that convection plays a dominant role over thermal conduction. Thevariation of time-averaged f conv with Rayleigh number is plotted for all simulations inFigure 5 a . This quantity is clearly asymptoting toward unity as Ra F is increased andreaches a value of 0.958 for our highest Ra F , N ρ = 4 run. Thermal conduction is thusplaying a minimal role in the average transport of thermal energy at high Rayleigh number.The role of conduction may be further characterized by considering the extent of thethermal boundary layer. We define the thermal boundary layer width in terms of the 23 – Ra F f c on v N ρ = 4N ρ = 3N ρ = 2N ρ = 1 10 Ra F B ounda r y La y e r W i d t h ( Mm ) Ra BL B ounda r y La y e r W i d t h ( Mm ) a b c Fig. 5.—
Thermal boundary layer and energy transport scaling with Rayleigh number. ( a ) Fractionalconvective flux f conv vs. Ra F . f conv approaches unity as Ra F is increased, indicating that conduction isplaying a minimal role in the energy transport. ( b ) Boundary layer width plotted vs. Ra F and ( c ) theboundary-layer dependent Rayleigh number Ra BL . The boundary layer width follows a clear scaling law of w BL ∝ Ra − / BL (dashed line, panel c ). time-averaged mean entropy, such that w BL = 1max h S ( r ) i Z h S ( r ) i dr. (26)The variation of w BL with Rayleigh number is illustrated in Figure 5 b .The boundary layer width decreases systematically with increasing Ra F and exhibits aclear dependence on the degree of density stratification N ρ . At comparable Ra F , simulationswith low values of N ρ possess a wider thermal boundary layer than those with higher valuesof N ρ . This difference may be largely understood by noting that the thermal boundarylayer is an upper-convection zone phenomenon and is better described by properties localto that region.We may define an alternative formulation of the Rayleigh number, Ra BL , bysubstituting point-wise values taken at the upper boundary for the volume-averagedquantities in Equation (14). Using this definition of the Rayleigh number, the multiplecurves apparent in Figure 5 b nearly collapse onto a single curve, as depicted in Figure 5 c . 24 –A fit to the scaling exponent yields the (unsurprising) relationship w BL ∝ Ra − . BL (27)In other words, the boundary layer width scales very nearly as √ κ . We note that for ourhighest Rayleigh number run, the boundary layer width is roughly 3 Mm, or about 2%of the shell depth. Thus, while our thermal boundary layers lack the radiative processesat work in the Sun, their physical extent is confined to a similarly small region of theconvective domain.
5. Perspectives and Conclusions
Our results hold interesting consequences for several aspects of convection zonedynamics. Many of these results will undoubtedly be sensitive to the effects of rotation andmagnetism which were omitted in this study. Nevertheless, we discuss some implications ofthis work and remark that we are now pursuing a complementary set of simulations thatincorporate those effects.
The kinetic-energy scaling observed in these simulations may be understood in terms ofa free-fall argument (e.g., Spiegal 1971). We outline such an argument as follows, assumingthat
P r = 1 (as is true for the models presented here) and that the kinetic energy derivesfrom the potential energy
P E BL associated with the boundary layer.Our conductive boundary layer possesses an entropy gradient given by Equation (13)and a width that is roughly proportional to Ra − / BL (i.e., √ κ ). Thus, the entropy contrast 25 –∆ S across that layer should scale as∆ S ≈ ∂S∂r (cid:12)(cid:12)(cid:12)(cid:12) r = r o w BL ∝ L ⋆ Ra − / BL πr o κ ¯ ρ ( r o ) ¯ T ( r o ) . (28)The boundary layer possesses an entropy deficit with respect to the adiabatic interior ofroughly ∆ S/
2. This entropy deficit translates into an overdensity within that region, withan associated mass M BL given by M BL ≈ ¯ ρ ( r o ) c p ∆ S × πr o w BL ∝ L ⋆ Ra − / BL κc p ¯ T ( r o ) (29)The potential energy, and thus the kinetic energy associated with the convection then scalesas KE ≈ P E BL ≈ M BL g H ∝ (cid:18) L ⋆ gc p ¯ T ( r o ) (cid:19) / (cid:18) πr o ¯ ρ ( r o ) H (cid:19) / , (30)which is independent of the thermal diffusivity owing to the two factors of √ κ arising fromthe boundary layer thickness. As we observe only a very weak dependence on diffusion inour KE scaling (see Equation 16), we suggest this is the dominant balance struck in ourhigh- Ra F models. This scaling appears as a common feature of high- Ra convection in thelaboratory setup (see e.g., Ahlers et al. 2009 and references therein). Moreover, Boussinesqstudies carried out using spherical geometry and symmetric, fixed temperature boundaryconditions (in lieu of internal heating) have demonstrated a similar scaling relationship(Gastine et al. 2015).The scaling properties of our boundary layer hold interesting consequences forextrapolation to the Sun. Perhaps it is most intuitive to think that the most realisticmodels of solar convection would be achieved by simulations with the same diffusivities orRayleigh number as the Sun. However, the radiative cooling which characterizes the upperboundary layer in the Sun differs markedly from the thermal conduction employed in ourmodels. As a result, extrapolating to the solar microscopic diffusivities may not actuallyreproduce the properties of solar convection. The convection in both the Sun and in our 26 –simulations is largely driven by intense buoyancy fluctuations within the thermal boundarylayer. Since, the structure of that convection is likely tied to the boundary layer’s width, itmight be best to extrapolate our models to a common boundary-layer thickness.This exercise requires far less extrapolation in parameter space because the Sun’sthermal boundary layer is only somewhat thinner than those achieved here. Our mostturbulent case possesses an Ra F of roughly 3 × and a boundary layer width of about3 Mm. Were we to seek a boundary layer width of 1 Mm, consistent with the size of atypical solar granule, we would need to decrease κ for that case by another factor of nineand increase Ra F by a factor of about 300. While this would certainly be challengingto accomplish using modern computational resources, it is not impossible. Alternatively,convective shells embodying a higher degree of density stratification might be examined tosimilarly achieve a thinner boundary layer. In either case, fully spectral simulations such asthese are, somewhat surprisingly, capable of resolving many of the relevant spatial scales inthe solar convection zone.If we were instead to follow our initial instinct and extrapolate our simulations to aRayleigh number more characteristic of stellar convection zones (e.g., 10 ), the boundarylayer would become unphysically thin. The resulting convection would likely be over-drivenat small scales. It is thus possible that the convective velocity spectra achieved in thesesimulations is more similar to that occuring in an actual star than in a hypothetical“stellar”-Rayleigh-number case. This subtlety should be kept in mind when extrapolatingresults such as these to stellar-like parameter regimes. Finally, we note that the measuredprefactor for our scaling is almost certainly determined by the conductive physics at workin the boundary layer, and that the implementation of radiative cooling will likely modifythe amplitude of those curves illustrated in Figures 2 a , b . 27 – While our simulations were non-rotating, it is nonetheless interesting to speculate onhow the behavior of the kinetic energy spectrum might impact the differential rotation. Themost reliable result from helioseismology with regard to the dynamics of the deep convectionzone continues to be the Sun’s internal-rotation profile obtained through global inversionsof p -mode frequency splittings (e.g., Thompson et al. 2003). Such inversions indicate thatthe surface differential rotation possesses a variation in angular velocity (decreasing by30% from equator to pole) that persists throughout the convection zone with little radialvariation.Because of this robust observational constraint, the onset of anti-solar differentialrotation in high Ra F , rotating convection is troubling. This behavior has been studiedin a variety of stellar and astrophysical contexts, and it has been unambiguously linkedto the degree of rotational influence felt by the convection (Gilman 1977; Glatzmaier &Gilman 1982; Gastine et al. 2013; Guerrero et al. 2013; Gastine et al. 2014; K¨apyl¨a etal. 2014; Featherstone & Miesch 2015). The rotational constraint felt by the convectionis characterized by the Rossby number, Ro , which expresses the ratio of a rotation periodto a convective time scale. It is typically defined as Ro = u/ H , where u is a typicalconvective flow speed, Ω is the rotation rate of the star, and H is some typical length scale(typically taken as the depth of the convective layer).A transition occurs between solar and anti-solar behavior when the convective androtational timescales become comparable, and Ro approaches unity. Alternatively, one maydefine a local Rossby number Ro ℓ = u rms / Ω ℓ c , where ℓ c is a characteristic length scaleassociated with the convection (not to be confused with the spherical harmonic degree ℓ ). Such a definition is equally useful for characterizing the transition between solar andanti-solar states (Gastine et al. 2014). The trend evident in the spectra of Figure 3 suggests 28 –that ℓ c will continue to decrease as Ra F is increased. At the same time, our results indicatethat u rms will remain fixed. If this trend continues, Ro ℓ will continue to increase as Ra F isincreased, and anti-solar states, which naturally arise at high values of Ro ℓ , may becomemore prevalent.We note that magnetism will almost certainly modify this discussion. In fact, severalstudies indicate that the presence of a Lorentz force, which tends to inhibit the mixing ofmomentum by convection, increases the transitional value of Ro in rotating systems (e.g.,Fan & Fang 2014; Karak et al. 2015; Simitev et al. 2015). How those magnetic effectsrelate to the spectral redistribution of energy discussed here remains to be examined. Our results may also hold interesting consequences for the study of the solar near-surface shear layer. That region of the convection zone, which spans the outer 5% of theSun by radius, is characterized by a diminishment of the angular velocity established withinthe bulk of the convection zone (e.g., Howe 2009). Precisely how this region of shear isestablished is not yet completely understood, though it is likely determined by a transitionfrom deep convection, which is thought to be rotationally-constrained, to near-surfaceconvection, which is only weakly influenced by rotation (Miesch & Hindman 2011; Hottaet al. 2014). A number of numerical convection studies suggest that density stratificationalso plays an important role in the development of such a region of shear (K¨apyl¨a et al.2011; Gastine et al. 2013; Hotta et al. 2014). One reason for this behavior may be thetendency of a background density stratification to accentuate differences in the speed ofnear-surface flows relative to deep flows, as shown in Figure 4 b . Higher values of Ra F alsomake this difference more pronounced (Figure 4 a ), though density stratification seems to bethe dominant effect in the parameter space explored here. Our results suggest that studies 29 –designed to systematically explore the development of a near-surface shear layer shouldfocus on the high- N ρ and high- Ra F regime. The results from this study indicate that care must be taken in the interpretationof simulated convective flows and, in particular, when carrying out comparisons betweenthose flows and their solar/stellar counterparts. The numerical results discussed hereillustrate one means by which stellar convection simulations may naturally overestimatethe low-wavenumber power in a convective velocity spectrum. Namely, at sufficiently highRayleigh number Ra F , the integrated, dimensional kinetic energy saturates at a constantvalue. That value is independent of the level of thermal diffusion. At the same time, asdiffusion is decreased and Ra F is increased, the convection naturally develops smaller-scalestructure and a corresponding increase in high-wavenumber power. As the integratedkinetic energy remains approximately constant, that high-wavenumber power must come atthe expense of low-wavenumber power.Simulations that are not run with a sufficiently high value of Ra F will thus possessspectra that naturally disagree with the low-wavenumber measurements accessible tohelioseismology. This general trend may contribute to the disagreement between simulationsand the deep time-distance measurements discussed in Hanasoge et al. (2012). It is alsopossible that this trend may help explain the disagreement between photospheric convectionsimulations and the observed photospheric power spectra as discussed in Lord et al. (2014).While such photospheric convection models incorporate radiative-transfer effects that werenot considered in this study, their convection is similarly constrained to transport a solarflux, and so it may develop a similar phenomenology. Those simulations also employ openboundary conditions and Cartesian geometry, however, and thus a more thorough study 30 –of this effect needs to be carried out within that experimental setup before definitiveconclusions may be drawn.What is clear from this work is that stellar convection simulations must be run in thehigh-Rayleigh-number regime if they are to correctly capture even the grossest propertyof a star’s convection, namely the integrated kinetic energy. Unfortunately, even in thehigh-Rayleigh-number regime, it remains unclear how well such simulations can address thespectral distribution of that energy, particularly if they do not properly capture the scaleof convective driving. For the parameter space explored here, the spectral range ℓ ≤ κ ≤ × cm s − ) are required to reach the high-Rayleigh number regime indentified here. However,virtually all stellar convection simulations incorporate rotation, and many now incorporatemagnetism. Our simulations incorporated neither. The inclusion of those effects will almostcertainly increase the critical Rayleigh number of the system (e.g., Chandrasekhar 1961),and so too the lower limit of the high-Rayleigh-number regime. It is thus timely to examinehow magnetism and rotation modify these findings, and then assess where current modelslie within the framework discussed here. In the interim, the utility of detailed spectralcomparisons between observations and simulations will remain questionable.We wish to thank Mark Miesch and Juri Toomre for several useful discussions overthe years that helped to frame the central question posed in this work. We further wishto thank Keith Julien for his many thoughtful comments that substantially improved thepresentation of these results. Finally, we wish to thank the staff at NASA Pleiades, NASADiscover, and the Argonne Leadership Computing Facility (ALCF; particularly Wei Jiang),without whom these simulations would not have been possible. 31 –This work was supported by NASA grants NNX09AB04G, NNX14AC05G,NNX11AJ36G, and NNX14AG05G. Featherstone and the development of Rayleigh werefurther supported by the Computational Infrastructure for Geodynamics (CIG), whichis supported by the National Science Foundation award NSF-094946. This work usedcomputational resources provided through NASA HEC support on the Pleiades and Discoversupercomputers. Further support was provided by an award of computer time on the Mirasupercomputer through the Innovative and Novel Computational Impact on Theory andExperiment (INCITE) program. This research used resources of the ALCF, which is a DOEOffice of Science User Facility supported under Contract DE-AC02-06CH11357.
A. Algorithmic Accuracy Considerations
The accuracy with which
Rayleigh solves the system of equations enumerated in § Rayleigh in Tables A.1 and A.2 below. Forreference, we also provide representative values from those reported in Jones et al. (2011).Specifically, we provide the values achieved by the Glatzmaier code, along with relativedifferences between the Glatzmaier results and the
Rayleigh results.
Rayleigh’s resultstypically agree with the accepted benchmark results to within one tenth of one percent. 32 –
TABLE A.1Anelastic MHD Benchmark ResultsObservable Rayleigh Glatzmaier % DifferenceKinetic Energy (erg) 5.56967 × × -1.09510 × − Zonal Kinetic Energy (erg) 6.37973 × × -1.97462 × − Meridional Kinetic Energy (erg) 1.49799 × × -1.73536 × − Entropy (erg g − K − ) 7.9449 × × -3.7759 × − u φ (cm s − ) 6.8638 × × -5.6355 × − Drift Frequency (rad s − ) 3.10509 × − × − -9.6615 × − Note. — Hydrodynamic benchmark results for
Rayleigh . The simulation was run with a resolution of128 × ×
384 ( N r × N θ × N φ ), corresponding to 96 Chebyshev modes and a maximum spherical harmonicdegree of 127. All quantities reported were averaged over 10,000 time steps following equilibration of thesteady-state solution. The time step size was 30 seconds, and the total evolution time was 4.8 × seconds.
33 –
TABLE A.2Anelastic MHD Benchmark ResultsObservable Rayleigh Glatzmaier % DifferenceKinetic Energy (erg) 8.03633 × × × − Zonal KE (erg) 1.15320 × × × − Meridional KE (erg) 1.01585 × × -1.96876 × − Magnetic Energy (erg) 6.13362 × × × − Zonal ME (erg) 4.62063 × × × − Meridional ME (erg) 3.24922 × × -1.53881 × − Entropy (erg g − K − ) 6.0886 × × -4.4340 × − u φ (cm s − ) -2.9374 × -2.9422 × -1.6314 × − B θ (G) 2.7187 × × -3.8473 × − Drift Frequency (rad s − ) 4.3075 × − × − -2.3215 × − Note. — Hydrodynamic benchmark results for
Rayleigh . The simulation was run with a resolution of128 × ×
384 ( N r × N θ × N φ ), corresponding to 96 Chebyshev modes and a maximum spherical harmonicdegree of 127. All quantities reported were averaged over 10,000 time steps following equilibration of thesteady-state solution. The time step size was 200 seconds, and the total evolution time was 5 × seconds.
34 –
B. Summary of Model Parameters and Results
In the four tables that follow, we report on the input and output parameters fromour set of 63 simulations. A separate table of parameters is presented for each value of N ρ used in this study. Each simulation is specified by three physical input parameters:the thermal diffusivity κ , the luminosity of the simulation L ⋆ , and the bulk flux Rayleighnumber Ra F as defined by Equation (14). We report the thermal diffusivity κ in units of10 cm s − , denoting it by κ in these tables. Similarly, we report luminosities that havebeen normalized by the solar luminosity L ⊙ = 3 . × erg s − . Additionally, we reporton the resolution used in each simulation, listing the maximum spherical harmonic andChebyshev polynomial degrees employed ( ℓ max and n max respectively). All variables havebeen de-aliased using the 2/3 rule. The number of collocation points in radius N R is givenby N R = ( n max + 1). Similarly, the number of collocation points in theta N θ is given by N θ = ( ℓ max + 1). The number of azimuthal points N φ is always twice N θ .We also report on the measured outputs for each simulation. These are the dimensionalkinetic energy KE , the non-dimensional kinetic energy d KE , the relative velocity difference δV , the fractional convective flux f conv , and the boundary layer width w BL . Definitions forthese values are provided in Equations (15), (17), (20), (25), and (26) respectively. Finally,while not discussed in this paper, we report on two Reynolds numbers measured for eachsimulation. The first is the bulk Reynolds number, which we denote by Re and define as Re = qg | v | Hν , (B1)where H and the tilde retain the same meaning as in the rest of this paper, representing theshell depth and a volume average respectively. The density stratification employed in thesesimulations means that the velocity amplitude varies substantially throughout the shell. As 35 –such, we also report the peak Reynolds number Re peak , which we define as Re peak = max v rms ( r ) Hν , (B2)where v rms is the rms velocity amplitude calculated at each depth. Both measures of thevelocity were averaged in time over 1,000 days. 36 – Table B.1: N ρ = 1 Simulation ParametersInput Parameters Measured Output κ L ⋆ /L ⊙ Ra F n max ℓ max KE d KE δV f conv w BL Re Re peak (10 erg) (Mm)0.5 1.00 3 . ×
85 1023 26.99 36,649.8 0.082 0.934 6.47 271.7 355.11 1.00 4 . ×
85 511 25.41 8,626.2 0.056 0.894 9.32 130.5 167.02 1.00 5 . ×
85 255 23.86 2,025.1 0.027 0.830 14.09 63.3 77.74 1.00 6 . ×
85 255 19.82 420.5 0.013 0.725 20.00 29.0 36.38 1.00 7 . ×
85 255 10.66 56.5 0.012 0.533 26.02 10.8 15.29 1.00 5 . ×
85 127 8.77 36.7 0.013 0.469 29.13 8.7 12.310 1.00 4 . ×
85 127 6.95 23.6 0.015 0.379 33.91 6.9 9.911 1.00 3 . ×
85 127 5.12 14.4 0.013 0.345 35.35 5.4 7.912 1.00 2 . ×
85 127 3.62 8.5 0.014 0.249 39.42 4.2 6.213 1.00 1 . ×
85 127 2.21 4.4 0.015 0.166 42.37 3.0 4.3
37 –
Table B.2: N ρ = 2 Simulation ParametersInput Parameters Measured Output κ L ⋆ /L ⊙ Ra F n max ℓ max KE d KE δV f conv w BL Re Re peak (10 erg) (Mm)1 1.00 1 . ×
85 511 34.67 16,095.4 0.416 0.932 6.80 189.4 243.22 1.00 1 . ×
85 255 33.08 3,839.6 0.310 0.887 9.66 91.9 120.04 1.00 1 . ×
85 255 30.92 897.0 0.179 0.807 13.86 44.8 60.88 1.00 2 . ×
85 255 23.94 173.6 0.084 0.689 19.65 20.5 30.19 1.00 1 . ×
85 85 20.28 116.3 0.091 0.644 19.70 16.7 24.210 1.00 1 . ×
85 85 19.41 90.1 0.073 0.620 21.62 15.0 23.511 1.00 9 . ×
85 85 17.25 66.2 0.102 0.565 21.62 12.5 18.512 1.00 7 . ×
85 85 15.53 50.1 0.099 0.526 22.72 10.9 16.313 1.00 5 . ×
85 85 12.71 34.9 0.087 0.500 24.15 9.3 15.314 1.00 4 . ×
85 85 10.84 25.7 0.099 0.422 25.97 7.8 12.415 1.00 3 . ×
85 85 8.99 18.6 0.100 0.382 27.17 6.9 10.916 1.00 3 . ×
85 85 7.06 12.8 0.106 0.330 28.79 5.5 9.317 1.00 2 . ×
85 85 5.64 9.1 0.118 0.290 29.75 4.7 8.218 1.00 2 . ×
85 85 3.33 4.8 0.115 0.194 32.23 3.5 5.919 1.00 1 . ×
85 85 2.36 3.0 0.111 0.139 33.50 2.7 4.620 1.00 1 . ×
85 85 1.16 1.3 0.104 0.074 34.89 1.8 3.1
38 –
Table B.3: N ρ = 3 Simulation ParametersInput Parameters Measured Output κ L ⋆ /L ⊙ Ra F n max ℓ max KE d KE δV f conv w BL Re Re peak (10 erg) (Mm)1 1.00 2 . ×
85 1023 36.22 20,057.7 1.323 0.948 4.96 229.8 320.42 1.00 2 . ×
85 511 34.73 4,807.6 1.094 0.912 6.93 114.1 166.54 1.00 3 . ×
85 263 33.42 1,156.7 0.799 0.847 9.83 56.3 85.66 1.00 9 . ×
85 127 32.92 506.4 0.617 0.791 12.18 37.5 58.28 1.00 4 . ×
85 127 31.62 273.5 0.452 0.736 14.09 27.5 42.412 1.00 1 . ×
85 127 22.08 84.9 0.275 0.623 16.03 15.4 24.016 1.00 5 . ×
42 127 12.02 26.0 0.266 0.442 19.06 9.0 15.624 1.00 1 . ×
42 63 2.05 2.0 0.198 0.122 23.77 2.5 4.50.5 0.19 3 . ×
170 85 12.27 27,173.6 1.353 0.953 4.65 265.3 363.51 0.29 6 . ×
85 1023 15.37 8,512.9 1.189 0.929 6.03 151.4 215.82 0.44 1 . ×
85 511 19.68 2,723.7 1.011 0.890 7.93 86.5 129.84 0.66 2 . ×
85 263 25.39 878.7 0.710 0.830 10.60 49.0 74.46 0.84 8 . ×
85 127 28.96 445.4 0.562 0.781 12.70 35.3 55.01 0.12 2 . ×
85 1023 8.71 4,820.8 1.096 0.912 6.94 113.8 166.12 0.25 6 . ×
85 511 13.44 1,860.7 0.898 0.873 8.75 71.2 105.94 0.50 1 . ×
85 263 21.00 726.7 0.713 0.815 11.09 44.3 68.56 0.75 7 . ×
85 127 27.12 417.1 0.545 0.773 12.85 34.2 52.7
39 –
Table B.4: N ρ = 4 Simulation ParametersInput Parameters Measured Output κ L ⋆ /L ⊙ Ra F n max ℓ max KE d KE δV f conv w BL Re Re peak (10 erg) (Mm)1 1.00 2 . ×
170 1023 34.38 20,892.4 3.216 0.958 3.37 254.0 388.42 1.00 3 . ×
85 511 32.21 4,892.5 2.691 0.925 4.66 124.2 193.04 1.00 4 . ×
85 511 30.91 1,173.8 1.961 0.864 6.63 61.5 98.38 1.00 5 . ×
85 255 29.49 280.0 1.145 0.755 9.42 29.9 48.59 1.00 3 . ×
85 127 29.04 217.9 1.060 0.732 9.99 26.5 43.510 1.00 2 . ×
85 127 28.72 174.5 0.901 0.712 10.54 23.5 38.311 1.00 2 . ×
85 127 27.03 135.7 0.800 0.688 11.00 20.9 34.212 1.00 1 . ×
85 127 25.64 108.2 0.697 0.661 11.47 18.8 31.113 1.00 1 . ×
85 127 22.53 81.0 0.479 0.624 11.99 16.4 27.614 1.00 1 . ×
85 127 20.44 63.4 0.410 0.594 12.32 14.6 24.915 1.00 8 . ×
85 127 18.45 49.8 0.409 0.557 12.57 13.0 22.216 1.00 6 . ×
85 127 16.32 38.7 0.383 0.522 12.81 11.6 20.017 1.00 5 . ×
85 85 14.18 29.8 0.380 0.486 13.04 10.2 18.018 1.00 4 . ×
85 85 12.38 23.2 0.410 0.436 13.37 9.2 16.319 1.00 4 . ×
85 85 10.40 17.5 0.379 0.393 13.62 7.9 14.320 1.00 3 . ×
85 85 7.21 10.9 0.177 0.336 14.09 6.4 12.521 1.00 3 . ×
85 85 5.49 7.6 0.186 0.270 14.50 5.3 10.422 1.00 2 . ×
85 85 3.69 4.6 0.196 0.196 15.01 4.2 8.323 1.00 2 . ×
85 85 3.10 3.6 0.498 0.180 15.03 3.8 7.224 1.00 2 . ×
85 85 2.96 3.1 0.475 0.166 15.15 3.4 6.5
40 –
C. Polytropic Background State
For these experiments, we employ a polytropic background state as formulated inthe anelastic benchmark paper of Jones et al. (2011). We assume that the gravitationalacceleration g ( r ) varies as GM i /r within the shell, where M i is the mass interior to thebase of the convection zone and G is the gravitational constant. The anelastic equationsthen admit an adiabatically stratified, polytropic atmosphere as the reference state: ρ = ρ i (cid:18) ζζ i (cid:19) n , T = T i ζζ i , p = p i (cid:18) ζζ i (cid:19) n +1 , (C1)where the subscript “ i ” denotes the value of a quantity at the inner boundary, and n is thepolytropic index. The radial variation of the reference state is captured by the variable ζ ,defined as ζ = c + c Hr , (C2)where H = r o − r i is the depth of the convection zone. The constants c and c are given by c = 2 ζ o − β − − β , c = (1 + β )(1 − ζ o )(1 − β ) , (C3)with ζ o = β + 1 β exp( N ρ /n ) + 1 , ζ i = 1 + β − ζ o β . (C4)Here ζ i and ζ o are the values of ζ on the inner and outer boundaries, β = r i /r o , and N ρ isthe number of density scale heights across the shell. 41 – REFERENCES
Ahlers, G., Grossmann, S. & Lohse, D., 2009, RvMP, 81, 503Braginsky, S.I., Roberts, P.H., 1995, GAFD, 79, 1Brun, A.S. & Toomre, J., 2002, ApJ, 570, 865Brun, A.S., Miesch, M.S. & Toomre, J., 2004, ApJ, 614, 1073Chandrasekhar, S., 1961, Hydrodynamic and Hydromagnetic Stability (Oxford: Clarendon)Charbonneau, P., 2010, LRSP, 7, 3Christensen-Dalsgaard, J., et al. 1996, Science, 272, 1286Fan, Y. & Fang, F., 2014, ApJ, 789, 35Featherstone, N.A. & Miesch, M.S., 2015, ApJ, 804, 67Gastine, T., Wicht, J. & Aurnou, J.M., 2013, Icar 225, 156Gastine, T., Yadav, R.K., Morin, J., Reiners, A. & Wicht, J., 2014, MNRAS, 438, L76Gastine, T., Wicht, J. & Aurnou, J.M., JFM, 778, 721Ghizaru, M., Charbonneau, P. & Smolarkiewicz, P.K., 2010, ApJL, 715, L133Gilman, P.A., 1977, GApFD, 8, 93Gilman, P.A. & Glatzmaier, G.A., 1981, ApJS, 45, 335Glatzmaier, G.A., 1984, J. Comp. Phys. 55, 461Glatzmaier, G.A. & Gilman, P.A., 1982, ApJS, 256, 316Gough, D.O., 1969, JAtS, 26, 448 42 –Greer, B.J., Hindman, B.W., Featherstone, N.A. & Toomre, J., 2015, ApJL, 803, 17Guerrero, G., Smolarkiewicz, P.K., Kosovichev, A.G. & Mansour, N.N., 2013, ApJ, 779, 176Hanasoge, S.M., Duvall, T.L., Jr. & Sreenivasan, K.R., 2012, PNAS, 109, 11928Hotta, H., Rempel, M. & Yokoyama, T., 2015, 798, 51Howe, R., 2009, LRSP, 6, 1Jones, C.A., Boronski, P., Brun, A.S., Glatzmaier, G.A., Gastine, T. Miesch, M.S. & Wicht,J., 2011, Icarus, 216,120Jouve, L, & Brun, A.S., 2007, A&A, 474, 239K¨apyl¨a, P., Mantere, M. & Brandenburg, A., 2011, AN, 332, 883K¨apyl¨a, P., Mantere, M. & Brandenburg, A., 2014, A&A, 570, 43Karak, B.B., K¨apyl¨a, P.J., K¨apyl¨a, M.J., Brandenburg, A., Olspert, N. & Pelt, J., 2015,A&A, 576, 26K˝ov´ari, Zs., Bartus, J., Strassmeier, K.G., Vida, K., ˘Svanda, M. & Ol´ah, K. , 2007, A&A,474, 165Lantz, S.R., 1992,