The magnitude of viscous dissipation in strongly stratified two-dimensional convection
Laura K. Currie Matthew K. Browning
astro-ph.SR physics.flu-dyn



















Department of Physics and Astronomy, University of Exeter, Stocker Road, EX4 4QL Exeter, UK; lcurrie@astro.ex.ac.uk






















Convection in astrophysical systems must be maintained against dissipation. Although the effects of dissipation are often assumed to be negligible, theory suggests that in strongly stratified convecting fluids, the dissipative heating rate can exceed the luminosity carried by convection. Here we explore this possibility using a series of numerical simulations. We consider two-dimensional numerical models of hydrodynamic convection in a Cartesian layer under the anelastic approximation and demonstrate that the dissipative heating rate can indeed exceed the imposed luminosity. We establish a theoretical expression for the ratio of the dissipative heating rate to the luminosity emerging at the upper boundary, in terms only of the depth of the layer and the thermal scale height. In particular, we show that this ratio is independent of the diffusivities and confirm this with a series of numerical simulations. Our results suggest that dissipative heating may significantly alter the internal dynamics of stars and planets.  
























§ INTRODUCTION

Convection occurs in the interiors of many astrophysical bodies and must be sustained against viscous and ohmic dissipation. This dissipation is often neglected in astrophysical models, e.g., in standard stellar 1D evolution codes <cit.> though its effects have lately been considered in a few specific contexts <cit.>.

Astrophysical convection often occurs over many scale heights. While for incompressible fluids the contribution of dissipative heating to the internal energy budget is negligible <cit.>, <cit.> (hereafter HMW) showed that in strongly stratified systems, it is theoretically possible for the rate of dissipative heating to exceed the luminosity. This was supported numerically by <cit.> for the case of a compressible liquid with infinite Prandtl number, Pr, (the ratio of viscous and thermal diffusivities), appropriate for models of the Earth's interior.

In this study we aim to establish the magnitude of dissipation for conditions more akin to those encountered in stellar interiors. Specifically, we consider dissipation in a stratified gas at finite Pr, and examine how the total heating changes as system parameters are varied. To begin, we briefly review some relevant thermodynamic considerations that underpin our work.
 

 


 §.§ Thermodynamic constraints on dissipative heating

For a volume V of convecting fluid enclosed by a surface S with associated magnetic field 𝐁, in which the normal component of the fluid velocity 𝐮 vanishes on the surface, and either all components of 𝐮, or the tangential stress, also vanish on the surface, local conservation of energy gives that the rate of change of total energy is equal to the sum of the net inward flux of energy and the rate of internal heat generation (e.g., by radioactivity or nuclear reactions). This implies 

    ∂/∂t(ρe+1/2ρu^2.   .+B^2/2μ_0-ρΨ)=-∇·(ρ(e+1/2u^2-Ψ)𝐮.
       .+(𝐄×𝐁)/μ_0+P𝐮-τ·𝐮-k∇T)+H

where ρ is the fluid density, e is the internal energy of the fluid, Ψ is the gravitational potential that satisfies 𝐠=∇Ψ, P is the pressure, τ_ij is the contribution to the total stress tensor from irreversible processes, k is the thermal conductivity, T is the temperature, H is the rate of internal heat generation, and 𝐄×𝐁/μ_0 is the Poynting flux (𝐄 is the electric field and μ_0 is the permeability of free space). 
Integrating (<ref>) over V gives the global relation
 
    ∫_Sk∂T/∂x_i dS_i+∫_VH dV=0,

assuming both a steady state and that the electric current, 𝐣, vanishes everywhere outside V. Equation (<ref>) implies that the net flux out of V is equal to the total rate of internal heating. Viscous and ohmic heating do not contribute to the overall heat flux: dissipative heating terms do not appear in equation (<ref>).

To examine dissipative heating, we consider the internal energy equation:

    ρ(∂e/∂t+(𝐮·∇)e)=∇(k∇T)-P(∇·𝐮)+τ_ij∂u_i/∂x_j+j^2/σ+H

where σ is the conductivity of the fluid.
Integrating over V, and assuming a steady state, (<ref>) becomes

    ∫_V(𝐮·∇)P dV+Φ =0.

Here 

    Φ=∫_Vτ_ij∂u_i/∂x_j+j^2/σ dV

is the total dissipative heating rate including viscous and ohmic heating terms.
Equation (<ref>) implies that the global rate of dissipative heating is cancelled by the work done against the pressure gradient. Equation (<ref>) is only equivalent to HMW's equation (22) when considering an ideal gas (so that αT=1, where α is the coefficient of thermal expansion); however, in arriving at (<ref>), we made no assumption about the fluid being a gas. <cit.> note that this inconsistency arises because HMW assume c_p to be constant in their derivation, which is not valid when α T≠1.

Alternatively, from the first law of thermodynamics, we have

    Tds=de-P/ρ^2dρ

where s is the specific entropy, so (<ref>) can also be written as 

    Φ=∫_VρT(𝐮·∇)s dV=-∫_Vρs(𝐮·∇)T dV

where we have invoked mass continuity in a steady state (∇·(ρ𝐮)=0).
Hence the global dissipation rate can also be thought of as being balanced by the work done against buoyancy <cit.>.

HMW used the entropy equation to derive an upper bound for the dissipative heating rate in a steadily convecting fluid that is valid for any equation of state or stress-strain relationship. 
For the case of convection in a plane layer, that upper bound is 

    Φ/L_u<T_max-T_u/T_u

where L_u is the luminosity at the upper boundary, T_max is the maximum temperature and T_u is the temperature on the upper boundary. 

One consequence of this bound is that, for large enough thermal gradients, the dissipative heating rate may exceed the heat flux through the layer; this is perhaps counter-intuitive, but is thermodynamically permitted, essentially because the dissipative heating remains in the system's internal energy <cit.>.   

The above considerations should hold for both ohmic and viscous dissipation. However, HMW further considered the simple case of viscous heating in a liquid (neglecting magnetism) and showed that the viscous dissipation rate is not only bounded by (<ref>) but that

    E≡Φ/L_u=d/H_T(1-μ/2)

where d is the height of the convective layer, H_T is the (constant) thermal scale height and 0≤μ≤1 is the fraction of internal heat generation.
Interestingly, the theoretical expression (<ref>) is dependent only on the ratio of the layer depth to the thermal scale height and the fraction of internal heat generation.

As expected, (<ref>) implies that the dissipative heating rate is negligible when compared with the heat flux in cases where the Boussinesq approximation is valid (i.e., when the scale heights of the system are large compared to the depth of the motion).
But it follows from (<ref>) that Φ is significant compared to L_u if d is comparable to H_T, i.e., if the system has significant thermal stratification. Stellar convection often lies in this regime, so it is not clear that dissipative heating can be ignored.

This paper explores these theoretical predictions using simulations of stratified convection under conditions akin to those encountered in stellar interiors.   Previous numerical simulations conducted by HMW considered only 2D Boussinesq convection and neglected inertial forces (infinite Pr approximation); later work by <cit.> within the so-called anelastic liquid approximation considered stronger stratifications but likewise assumed a liquid at infinite Pr.  We extend these by considering an ideal gas (so that αT=1) at finite Pr, so inertial effects are important and compressibility is not negligible. 

In section <ref>, we describe the model setup before presenting results from numerical simulations. In section <ref> we offer a discussion of the most significant results that emerge before providing conclusions.




§ SIMULATIONS OF DISSIPATIVE CONVECTION



 §.§ Model setup

We consider a layer of convecting fluid lying between impermeable boundaries at z=0 and z=d. We assume thermodynamic quantities to be comprised of a background, time-independent, reference state and perturbations to this reference state. The reference state is taken to be a polytropic, ideal gas with polytropic index m given by 

    T̅=T_0(1-β z), ρ̅=ρ_0(1-β z)^m, p̅=ℛρ_0T_0(1-β z)^m+1,

where β=g/c_p,0T_0. Here, g is the acceleration due to gravity, c_p is the specific heat capacity at constant pressure, ℛ is the ideal gas constant and a subscript 0 represents the value of that quantity on the bottom boundary. β is equivalent to the inverse temperature scale height and so is a measure of the stratification of the layer, although we shall use the more conventional 

    N_ρ=-mln(1-β d)

to quantify the stratification, with N_ρ the number of density scale heights across the layer. We assume a polytropic, monatomic, adiabatic, ideal gas, therefore m=1.5. Here we consider only the hydrodynamic problem; i.e., all dissipation is viscous.

We use anelastic equations under the Lantz-Braginsky-Roberts (LBR) approximation <cit.>; these are valid when the reference state is nearly adiabatic and when the flows are subsonic <cit.>, as they are here.

The governing equations are then

    ∂𝐮/∂t   +(𝐮·∇)𝐮=-∇p̃+gs/c_p𝐞̂_̂𝐳̂
       +ν[1/ρ̅∂/∂x_j(ρ̅(∂u_i/∂x_j+∂u_j/∂x_i))-2/3ρ̅∂/∂x_i(ρ̅∂u_j/∂x_j)]
 

    ∇·(ρ̅𝐮)=0


    ρ̅T̅(∂s/∂t+(𝐮·∇)s)=∇·(κρ̅T̅∇s)+τ_ij∂u_i/∂x_j+H,

where 𝐮 is the fluid velocity, p̃=p/ρ̅ is a modified pressure and ν is the kinematic viscosity. The specific entropy, s, is related to pressure and density by

    s=c_vlnp-c_plnρ.

We assume the perturbation of the thermodynamic quantities to be small compared with their reference state value. Therefore the entropy is obtained from

    s=c_vp/p̅-c_pρ/ρ̅

and the linearised equation of state is

    p/p̅=T/T̅+ρ/ρ̅.

In (<ref>) κ is the thermal diffusivity and 
    τ_ij=νρ̅(∂u_i/∂x_j+∂u_j/∂x_i-2/3δ_ij∇·𝐮)

is the viscous stress tensor (δ_ij is the Kronecker delta).
Here, we only consider cases with H=0 (i.e., no internal heat generation), and instead impose a flux (F) at the bottom boundary.
Note the LBR approximation diffuses entropy (not temperature); see <cit.> for a discussion of the differences.
We assume a constant ν and κ. 

We solve these equations using the Dedalus pseudo-spectral code <cit.> with fixed flux on the lower boundary and fixed entropy on the upper boundary. 
We assume these boundaries to be impermeable and stress-free. We employ a sin/cosine decomposition in the horizontal, ensuring there is no lateral heat flux.
We employ the semi-implicit Crank-Nicolson Adams-Bashforth numerical scheme and typically use 192 grid points in each direction with dealiasing (so that 128 modes are used). In some cases, 384 (256) grid points (modes) were used to ensure adequate resolution of the solutions.
For simplicity, and to compare our results with those of HMW, we consider 2D solutions so that 𝐮=(u,0,w) and ∂/∂y≡0. This also allows us to reach higher supercriticalities and N_ρ with relative ease.

As we neglect magnetism, the total dissipation rate, Φ, is given by (<ref>) with 𝐣=0 and τ_ij as given by (<ref>).

An appropriate non-dimensionalisation of the system allows the parameter space to be collapsed such that the dimensionless solutions (in particular E) are fully specified by m, N_ρ, Pr, together with F̂_̂0̂= Fd/κc_p,0ρ_0T_0 (a dimensionless measure of the flux applied at the lower boundary) and a flux-based Rayleigh number <cit.>      

    Ra=gd^4F_u/νκ^2ρ_0c_p,0T_0.

The parameters used in our simulations are given in Table <ref>.  

In a steady state, an expression for the luminosity L at each depth z=z' can be obtained by integrating the internal energy equation (<ref>) over the volume contained between the bottom of the layer and the depth z=z':

    L=   FA=∫_V_z'∇·(ρ̅T̅s𝐮) dV+∫_V_z'-∇·(κρ̅T̅∇s) dV
       +∫_V_z'-sρ̅(𝐮·∇)T̅ dV+∫_V_z'-τ_ij∂u_i/∂x_j dV,

where A is the surface area.
The divergence theorem allows the first two integrals to be transformed into surface integrals giving

    L=   FA=∫_S_z'ρ̅T̅sw dS_L_conv=AF_conv+∫_S_z'-κρ̅T̅∂s/∂z dS_L_cond=AF_cond
       +∫_V_z'-sρ̅(𝐮·∇)T̅ dV_L_buoy=A∫_0^z'Q_buoy dz+∫_V_z'-τ_ij∂u_i/∂x_j dV_L_diss=A∫_0^z'Q_diss dz,

where the surface integrals are over the surface at height z=z'.
The first and second terms define the horizontally-averaged heat fluxes associated with convection (F_conv) and conduction (F_cond) respectively, along with associated luminosities.
The third and fourth terms define additional sources of heating and cooling (Q_diss and Q_buoy) associated with viscous dissipation and with work done against the background stratification, respectively. These two terms must cancel in a global sense i.e., when integrating from z=0 to z=d, but they do not necessarily cancel at each layer depth. 

An alternative view of the heat transport may be derived by considering the total energy equation (<ref>), which includes both internal and mechanical energy. 
In a steady state (with entropy diffusion), the local balance gives

    ∇·(ρ̅(e+1/2u^2-Ψ)𝐮+p𝐮-τ·𝐮-κρ̅T∇s)=H

which when integrated over the volume for an ideal gas gives <cit.>

    L=   FA=∫_S_z'ρ̅c_pwT' dS_L_e=AF_e+∫_S_z'-κρ̅T̅∂s/∂z dS_L_cond=AF_cond
       +∫_S_z'1/2ρ̅|u^2|w dS_L_KE=AF_KE+∫_S_z'-(τ_iju_i)·ê_𝐳 dS_L_visc=AF_visc,

defining the horizontally-averaged enthalpy flux (F_e), kinetic energy flux (F_KE) and viscous flux (F_visc).
Note that (<ref>) and (<ref>) are equivalent; whether decomposed in the manner of (<ref>) or the complementary fashion of (<ref>), the transport terms must sum to the total luminosity L.
L_visc represents the total work done by surface forces, whereas L_diss represents only the (negative-definite) portion of this that goes into deforming a fluid parcel and hence into heating. 					      



 §.§ Relations between global dissipation rate and convective flux

For the model described in section <ref>, equation (<ref>) becomes

    Φ=   -∫_Vρ̅s(𝐮·∇)T̅ dV
    
    =   g/c_p,0∫_Vsρ̅w dV=gA/c_p,0∫_0^dF_conv/T̅ dz,

Often it is assumed that in the bulk of the convection zone, the total heat flux is just equal to the convective flux as defined above (i.e., F_conv≈F). We show later that this a poor assumption in strongly stratified cases, but it is reasonable for approximately Boussinesq systems. In the case F_conv≈F, (<ref>) becomes

    Φ=gAF/c_p,0T_0∫_0^d1/1-βz dz=-L_uln(1-βd)

and 

    E=-ln(1-βd)=βd+…≈d/H_T,0.

However, in strongly stratified cases F≈F_conv+F_other where F_other=∫_0^z'(Q_buoy+Q_diss) dz from (<ref>), or alternatively, F_other=F_p+F_KE+F_visc from (<ref>) (the conductive flux is small in the bulk convection zone). Here F_p=1/A∫_S_z'wp dS is the difference between the enthalpy flux F_e and the convective flux F_conv.  Physically, F_other is equivalent to the steady-state transport associated with processes other than the convective flux as defined above. In this case, (<ref>) becomes

    Φ=gAF/c_p,0∫_0^d(1-F_other/F)1/T̅ dz,

where we note that in general F_other is a function of depth and (1-F_other/F)≥1.
A complete theory of convection would specify F_other a priori, and thereby constrain the dissipative heating everywhere.
In the absence of such a theory,  we turn to numerical simulations to determine the magnitude of Φ for strong stratifications.







 §.§ Dissipation in simulations: determined by stratification

We examine the steady-state magnitude of Φ for different values of N_ρ and Ra. 
Figure <ref> shows the ratio of the global dissipation rate to the luminosity through the layer, E=Φ/L_u, for varying stratifications. First, we highlight the difference between simulations in which the dissipative heating terms were included (red squares) and those where they were not (black circles). At weak stratification, 
there is not much difference in the dissipative heating rate between these cases, but differences become apparent as N_ρ is increased. Including the heating terms in a self-consistent calculation leads to a much larger value of E than if Φ is only calculated after the simulation has run (i.e., if heating is not allowed to feedback on the system). When heating terms are included, the global dissipative heating rate exceeds the flux passing through the system (i.e., E>1) when N_ρ>1.22. 

As expected, the expression for E, in the Boussinesq limit, given by (<ref>), is a good approximation to E for small N_ρ, but vastly underestimates E at large N_ρ (see Figure <ref>, dash-dot line).
In the cases where the heating terms are not included, E cannot exceed unity for all N_ρ.
This might have been expected, since in this case none of the dissipated heat is returned to the internal energy of the system; instead, the dissipated energy is simply lost (i.e., energy is not conserved). This has the practical consequence that the flux emerging from the top of the layer is less than that input at the bottom.
In these cases E is very well described by the dashed line which is given by d/H_T,0, the leading order term from the expression for E in (<ref>).

The theoretical upper bound derived by HMW is shown on Figure <ref> by the solid black line. It is clear that all of our cases fit well within this upper bound, even at strong stratifications. This upper bound is equivalent to d/H_T,u in this system, where H_T,u is the value of H_T on the upper boundary.

Cases in which the heating terms were included are well described by 

    E=d/H̃_̃T̃,

where

    H̃_̃T̃ = H_T,0H_T,u/H_T,z^*

is a modified thermal scale height involving H_T at the top, bottom and at a height z^*, defined such that half the fluid (by mass) lies below z^* and half sits above; for a uniform density fluid, z^*=d/2. This expression resembles that originally proposed by HMW, on heuristic grounds, for a gas (E≈d/H_T); in our case H_T is not constant across the layer and we find that the combination H̃_̃T̃ is the appropriate “scale height" instead. Like HMW's suggestion, it depends only on the layer depth and temperature scale heights of the system.

For 2D convection, at Pr=1 and the Ra considered here, the solutions are steady (time-independent) <cit.>; the convection takes the form of a single stationary cell occupying the layer. To assess if the same behaviour occurs for chaotic (time-dependent) solutions, we have included some cases at Pr=10 (orange triangles), since then the flow is unsteady. In the cases included here, this unsteady flow is characterised by the breakup of the single coherent convection cell (seen at Pr=1); these time-dependent solutions seem also to be well described by the line given by (<ref>). This behaviour is sampled in Figure <ref>, Supplementary Material, which shows the velocity and entropy fields in a simulation with Pr=10, N_ρ=1.31, Ra=4.13×10^8 and F̂_0=0.14.
At higher Ra, the solutions transition to turbulence <cit.>. 




 §.§ Dissipation in simulations: independent of diffusivities

The results of section <ref>, specifically equation (<ref>), suggest that the amount of dissipative heating is determined by the stratification, not by other parameters such as Ra. To probe this further, we consider how/if E changes as Ra is varied. Figure <ref> shows the results for three different stratifications. For N_ρ≈0.1, the fluid is close to being Boussinesq and it is clear that E remains constant (and equal to the value given by (<ref>)) for many decades increase in Ra. This result complements that of HMW obtained from Boussinesq simulations at infinite Pr. For increasing N_ρ, we find that for large enough Ra, E approaches the constant given by (<ref>). That E becomes independent of Ra at large enough Ra for all N_ρ was also found by <cit.>, albeit for liquids at infinite Pr. 

Figure <ref> indicates that the solutions have to be sufficiently supercritical in order for the theory to be valid. It also suggests that stronger stratifications require simulations to be more supercritical in order to reach the asymptotic regime. (All the simulations displayed in Figure <ref> approach this asymptotic regime, except possibly the uppermost point at N_ρ=2.8. That simulation has Ra/Ra_c ≈ 9 ×10^5, but it is likely that still higher Ra would yield somewhat greater values of E at this stratification.)







§ DISCUSSION AND CONCLUSION

We have demonstrated explicitly that the amount of dissipative heating in a convective gaseous layer can, for strong stratifications, equal or exceed the luminosity through the layer. 
A principal conclusion is that the ratio of the global viscous heating rate to the emergent luminosity is approximated by a theoretical expression dependent only on the depth of the layer and its thermal scale heights.  This ratio, akin to one originally derived for a simpler system by HMW, is given (for the cases studied here) by (<ref>). Interestingly, this relation does not depend on other parameters such as the Rayleigh number. Our simulations confirm that this expression holds for 2D convection in an anelastic gas, provided the convection is sufficiently supercritical.  This regime is attainable in our 2D simulations, and is surely reached in real astrophysical objects, but may be more challenging to obtain in (for example) 3D global calculations <cit.>.

The dissipative heating appears in the local internal energy (or entropy) equation, in the same way as heating by fusion or radioactive decay.  Where it is large, we therefore expect it will modify the thermal structure, just as including a new source of heating or cooling would have done.  It must be reiterated, though, that in a global sense this heating is balanced by equivalent cooling terms; i.e., L_diss and L_buoy in equation (<ref>) cancel in a global sense; no additional flux emerges from the upper boundary. Stars are not brighter because of viscous dissipation. Locally, however, these terms do not necessarily cancel, as explored in Figure <ref>.  There we show the net heating and cooling at each depth in two simulations; in Figure <ref>a, the fluid is weakly stratified, and in (b) is has a stratification given by N_ρ=2.08. In both cases the sum of the terms must be zero at the top and bottom of the layer, but not in between.
Furthermore, in (a) the terms are small compared to the flux through the layer (typically a few %) but in the strongly stratified case, the local heating and cooling become comparable to the overall luminosity. 
In general, stronger stratifications lead to stronger local heating and cooling in the fluid.


In a steady state the imbalance between this local heating and cooling is equivalent to certain transport terms as discussed in section <ref>; these are assessed for our simulations in figure <ref> where the terms are plotted as luminosities and labelled correspondingly. Turning first to Figure <ref>a, we show the components of the total flux of thermal energy (as described by (<ref>)), namely L_conv, L_cond, L_buoy and L_diss. The conductive flux is small throughout the domain except in thin boundary layers and the dissipative heating (L_diss) is comparable to the convective flux (L_conv) throughout the domain. The sum of the four transport terms is shown as the black line (L) and is constant across the layer depth, indicating thermal balance.
Figure <ref>b assesses the total energy transport using the complementary analysis of (<ref>), using L_KE, L_cond, L_e and L_visc. The primary balance is between the positive L_e and the negative L_KE. Viewed in this way, the viscous flux (L_visc) is small except near the lower boundary, but (as discussed in section <ref>) this does not necessarily mean the effect of viscous dissipation is also small.
In figure <ref>c we highlight the equivalence of some transport terms, by showing the term AF_other together with its different constituent terms from either the total or thermal energy equations. As expected, AF_other is the same in both cases; it is the sum of L_diss and L_buoy, or equivalently, it is the sum of L_p, L_KE and L_visc. That is, changes in the dissipative heating are reflected not just in Q_diss (if analysing internal energy) or F_visc (if analysing total energy); the other transport terms (F_KE, F_p, F_e, F_conv, Q_buoy) also change in response.
To emphasise the importance of dissipative heating in modifying the transport terms, we include in Figure <ref>d, L_KE^nh , L_e^nh , L_cond^nh  and L_visc^nh i.e., the kinetic energy, enthalpy, conductive and viscous fluxes (expressed as luminosities) respectively, in the case where heating terms were not included. It is clear that these are much smaller than in the equivalent simulation with heating (Figure <ref>b), demonstrating explicitly that the inclusion of dissipative heating influences the other transport terms.
In particular, the maximum value of the kinetic energy flux is 3.2 times larger when the heating terms are included. The black line in Figure <ref>d shows that when heating is not included the flux emerging at the upper boundary is smaller than the flux imposed at the lower boundary; in this case it is approximately 27% of L.

The local heating and cooling (or, equivalently, the transport term F_other that must arise from this in a steady state) described above is not included in standard 1D stellar evolution models, and we do not yet know what effects (if any) would arise from its inclusion. 
In some contexts they may be negligible; the total internal energy of a star is enormously greater than its luminosity L⋆, so even internal heating that exceeds L⋆ may not have a noticeable effect on the gross structure. If, however, this heating is concentrated in certain regions (e.g., because of spatially varying conductivity) or occurs in places with lower heat capacity, its impact may be more significant.

If the results explored here also apply to the full 3D problem with rotation and magnetism – which clearly must be checked by future calculation – then the total dissipative heating is determined non-locally, dependent as it is on the total layer depth.  Simple modifications to the mixing-length theory (which is determined locally) may not then suffice to capture it. We have begun to explore these issues by modification of a suitable 1D stellar evolution code, and will report on this in future work.































































We acknowledge support from the European Research Council under ERC grant agreements No. 337705 (CHASM). The simulations here were carried out on the University of Exeter supercomputer, a DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS and the University of Exeter. We also acknowledge PRACE for awarding us access to computational resources Mare Nostrum based in Spain at the Barcelona Supercomputing Center, and Fermi and Marconi based in Italy at Cineca. We thank the referee for a thoughtful review that helped to improve the manuscript.



§ SIMULATION PARAMETERS

CCCCCC
Simulation parameters used in figures <ref>-<ref> 
6
1
0pt

Pr    
N_ρ     Ra     F̂_̂0̂    
E     Figure

1     0.1050     3.83 ×10^5     3.26 ×10^-4     0.0630     2

1     0.1050     3.83 ×10^6     3.26 ×10^-3     0.0662     2

1     0.1050     2.63 ×10^7     2.24 ×10^-2     0.0678     2

1     0.1050     6.13 ×10^7     5.22 ×10^-2     0.0682     2

1     0.1050     3.83 ×10^8     3.26 ×10^-1     0.0689    1-3

1     0.2776     6.58 ×10^7     5.60 ×10^-2     0.1828    1

1     0.3828     8.77 ×10^7     7.47 ×10^-2     0.2557    1

1     0.5819     8.01 ×10^7     1.07 ×10^-3     0.4014    1

1     0.7060     6.65 ×10^4     3.62 ×10^-3     0.4159    2

1     0.7060     6.65 ×10^5     3.62 ×10^-2     0.4594    2

1     0.7060     9.36 ×10^6     1.24 ×10^-4     0.4875    2

1     0.7060     1.05 ×10^8     2.64 ×10^-3     0.5008    2

1     0.7060     2.72 ×10^8     3.62 ×10^-3     0.5038    2

1     0.7060     4.88 ×10^8     4.40 ×10^-3     0.5057    1-2

1     0.7967     1.03 ×10^8     1.37 ×10^-3     0.5770    1

1     0.9887     1.20 ×10^8     1.60 ×10^-3     0.7533    1

1     1.3104     8.45 ×10^7     1.12 ×10^-3     1.0908    1

1     2.0846     1.33 ×10^5     7.24 ×10^-3     1.5830    2

1     2.0846     1.33 ×10^6     2.68 ×10^-3     1.8726    2

1     2.0846     1.63 ×10^7     2.17 ×10^-4     2.0882    2

1     2.0846     5.44 ×10^7     7.24 ×10^-4     2.1656    1-4

1     2.7938     1.23 ×10^8     1.63 ×10^-3     3.5951    1
     
10     0.1050     2.63 ×10^7     1.43 ×10^-1     0.0668   1

10     0.2776     1.15 ×10^8     9.78 ×10^-3     0.1822   1

10     0.3828     3.25 ×10^6     1.59 ×10^-1     0.2454   1

10     0.9887     1.37 ×10^9     4.66 ×10^-1     0.7413   1

10     1.3104     4.13 ×10^8     1.40 ×10^-1     1.0594   1




































































[Alboussiere & Ricard(2013)]AlboussiereRicard2013 Alboussiere, T., & Ricard, Y. 2013, JFM, 725, R1
[Alboussiere & Ricard(2014)]AlboussiereRicard2014 Alboussiere, T., & Ricard, Y. 2014, JFM, 751, 749

[Aubert et al.(2017)]Aubertetal2017 Aubert, J., Gastine, T., & Fournier, A. 2017, JFM, 813, 558
[Backus(1975)]Backus1975 Backus, G. E. 1975, PNAS, 72, 1555

[Batygin & Stevenson(2010)]BatyginStevenson2010 Batygin, K., & Stevenson, D., J. 2010, ApJL, 714, L238

[Braginsky & Roberts(1995)]BraginskyRoberts1995 Braginsky, S. I., & Roberts, P. H. 1995, GApFD, 79, 1
[Browning et al.(2016)]Browningetal2016 Browning, M. K., Weber, M. A.,  Chabrier, G., & Massey, A. P. 2016, , 812, 189

[Burns et al.(in prep.)]dedalus Burns, K. J., Vasil, G. M., Oishi, J. S., Lecoanet, D., Brown, B. P. & Quataert, E. http://dedalus-project.org (In preparation)


[Chabrier & Baraffe(1997)]ChabrierBaraffe1997 Chabrier, G., & Baraffe, I. 1997, A&A, 327, 1039

[Duarte et al.(2016)]Duarteetal2016 Duarte, L. D., Wicht, J., Browning, M. K., & Gastine, T. 2016. MNRAS, 456, 1708

[Featherstone & Hindman(2016)]FeatherstoneHindman2016 Featherstone, N. A., & Hindman, B. W. 2016, , 818, 32

[Gough(1969)]Gough1969 Gough, D. O. 1969, JAtS, 26, 448

[Hewitt et al.(1975)]Hewittetal1975 Hewitt, J. M., McKenzie, D. P., & Weiss, N. O. 1975, JFM, 68, 721

[Jarvis & McKenzie(1980)]JarvisMcKenzie1980 Jarvis, G. T., & McKenzie, D. P. 1980, JFM, 96, 515
[Jones & Kuzanyan(2009)]JonesKuzanyan2009 Jones, C. A., & Kuzanyan, K. M. 2009, Icar, 204, 227

[Kundu(1990)]Kundu1990 Kundu, P. K. 1990, Fluid Mechanics (Academic Press)

[Lantz(1992)]Lantz1992 Lantz, S. R. 1992, Ph.D. Thesis, Cornell University.
[Lantz & Fan(1999)]LantzFan1999 Lantz, S. R., & Fan, Y. 1999, ApJS, 121, 247
[Lecoanet et al.(2014)]Lecoanetetal2014 Lecoanet, D., Brown, B. P., Zweibel, E. G., Burns, K. J., Oishi, J. S., & Vasil, G. M. 2014, , 797, 94

[Ogura & Phillips(1962)]OguraPhillips1962 Ogura, Y. & Phillips, N. A. 1962, JAtS, 19, 173

[Paxton et al.(2011)]Paxtonetal2011 Paxton, B., Bildsten, L., Dotter, A., Herwig, F., Lesaffre, P., & Timmes, F. 2011, ApJS, 192, 3

[Rogers et al.(2003)]Rogersetal2003 Rogers, T. M., Glatzmaier, G. A., & Woosley, S. E. 2003, PhRvE, 67, 026315

[Viallet et al.(2013)]Vialletetal2013 Viallet, M., Meakin, C., Arnett, D., & Mocák, M. 2013, , 769, 1

[Vincent & Yuen(1999)]VincentYuen1999 Vincent, A. P., & Yuen, D. A. 1999, PhRvE, 60, 2957
















