Internally heated and fully compressible convection: flow morphology and scaling laws

In stars and planets natural processes heat convective flows in the bulk of a convective region rather than at hard boundaries. By characterizing how convective dynamics are determined by the strength of an internal heating source we can gain insight into the processes driving astrophysical convection. Internally heated convection has been studied extensively in incompressible fluids, but the effects of stratification and compressibility have not been examined in detail. In this work, we study fully compressible convection driven by a spatially uniform heating source in 2D and 3D Cartesian, hydrodynamic simulations. We use a fixed temperature upper boundary condition which results in a system that is internally heated in the bulk and cooled at the top. We find that the flow speed, as measured by the Mach number, and turbulence, as measured by the Reynolds number, can be independently controlled by separately varying the characteristic temperature gradient from internal heating and the diffusivities. 2D simulations at a fixed Mach number (flow speed) demonstrate consistent power at low wavenumber as diffusivities are decreased. We observe convection where the velocity distribution is skewed towards cold, fast downflows, and that the flow speed is related to the length scale and entropy gradient of the upper boundary where the downflows are driven. We additionally find a heat transport scaling law which is consistent with prior incompressible work.


I. INTRODUCTION
Turbulent convective flows are critical to many astrophysical phenomena.In planets and brown dwarfs, convection drives observable chemical disequilibrium [1][2][3] and cloud formation [4,5].In stars, convection generates differential rotation [6,7] and magnetic dynamos [8,9].While convection is common, it remains poorly understood, as is indicated by the solar convective conundrum [10,11], where reported helioseismic observations of convective flows differ by orders of magnitude from one another and from theory and simulations [12][13][14][15].This disagreement calls for a re-examination of the fundamental processes driving solar convection.In order to understand emergent convectivelydriven phenomena in stars and planets, we must develop robust models of the convection driving them.
Astrophysical convection is often driven by internal heating.In stars, convection is caused by changes in opacity and nuclear fusion, both of which enter the equations of motion as spatially distributed heating processes [e.g., 16].In gas giant planets, convection is driven by changes in radiative conductivity and the latent heat of formation [17].The evolution time of stars is much longer than the dynamical times for convection, and so a non-uniform conductivity can be treated like a stationary-in-time source of internal heating.Simulations which attempt to model astrophysical systems include many physical processes, one of which is internal heating [18][19][20][21][22][23].While many complex simulations drive convection with internal heating, these studies usually examine how the heating strength affects processes including boundary mixing and waves in adjacent stable regions, and there has been less work studying how heating strength affects the fundamental flow measures in fully compressible, stratified fluids [24][25][26][27].It is important to understand how the input parameters of internally heated convection map to important measures such as Mach number number so that we can produce suites of simulations where Mach number is held constant at a known value for a particular system.
In many experiments and numerical simulations, convection is driven at hot or cold boundaries, and the resulting boundary layers throttle the heat transport [28,29].Some recent simulations and experiments of incompressible fluids instead drive convection using internal heating in an effort to reduce boundary effects.Numerical simulations from Goluskin and Spiegel [30] and Goluskin and van der Poel [31] measured heat transport scaling laws in internally heated Boussinesq convection, and these findings were generalized to find mathematically robust limits to heat transport scalings in Goluskin [32] and diffusion-based bounds in Wang et al. [33].Miquel et al. [34] found that internally heated convection with a non-uniform heat source can, surprisingly, produce heat transport scaling in excess of the limit set by mixing length theory.Several laboratory experiments have studied internally heated convection using a dye to absorb light thereby heating the fluid and have found that heat transport scaling depends on the concentration of the dye and thus the length scale over which the fluid is heated [35,36].Earlier internal heating experiments used electrical currents for bulk heating, and incorporated a cooled upper boundary [37].The approach of using internal heating and cooling has also found its way into some studies of astrophysical convection driven through heating and cooling layers since it produces diffusion free scalings, which are relevant in astrophysical settings where bulk heating processes drive convective flows without hard boundaries [38,39].
In this work we study fully compressible convection which is internally heated in the bulk and cooled at the top boundary.This is accomplished using a constant and uniform source term of internal energy acting upon an adiabatically stratified, ideal gas atmosphere.We enforce a fixed temperature upper boundary, giving rise to boundary cooling at the top of the domain.This setup has been studied in the case of Boussinesq convection in the uniform heating case of Kazemi et al. [40].This case is also analogous to the IH3 configuration described in Goluskin [32].This setup resembles convection near the photosphere of lower main sequence and red giant branch stars, where rapid photospheric cooling can be approximated as boundary cooling.
We characterize internally heated fully compressible convection, and in particular how quantities like the Mach number (Ma) and Reynolds number (Re) vary with our chosen input parameters.Mach number may be known for a particular system, and in this case it is important to have a suite of simulations with this particular Mach number and different Reynolds numbers.Mach number is measure of how important compressibility is to the dynamics, and thus when alternative equation sets such as the anelastic or Boussinesq can be used.In this work we determine the input parameters that allow us to change Mach number and Reynolds number independently.
We demonstrate that changing the magnitude of the heating term affects the Mach number, and thus the importance of compressibility.With this setup, flow speed as measured by the Mach number, and turbulence as measured by the Reynolds number, can be controlled independently in 2D simulations and are nearly independent in 3D.We observe convective flows dominated by cool downflows produced at the upper boundary layer.At fixed heating magnitude, we find that heat transport and the Reynolds number follow expected scaling laws from Boussinesq convection.With independent control of Mach number and Reynolds number we can produce a series of simulations with consistent behavior in convective power spectra at low wavenumber.

II. METHODS
We study the evolution of an ideal gas according to the fully compressible Navier-Stokes equations.We follow previous works [e.g., 41,42] and use a temperature and log-density formulation of the equations.We modify the equations by adding a heating term, Q, to the internal energy equation.These fully compressible equations are: where ρ is the density, u = ux + v ŷ + wẑ is the velocity vector, T is temperature, µ is the dynamic viscosity, κ is the thermal diffusion, and Q is an internal heating source term.We take µ, κ, and Q to be constant in time and spatially uniform, which is assumed in the expressions of equations 2-3.We use an ideal gas equation of state for pressure P = R m ρT , where R is the gas constant and m is the mean molecular weight.c V and c P are the specific heat at constant volume and pressure respectively, and their ratio is the adiabatic exponent γ = c P c V .The viscous stress tensor σ and viscous heating term Φ are respectively defined as where I ij is the identity tensor and is the rate-of-strain tensor.We decompose our thermodynamics into an adiabatic and hydrostatically equilibrated portion (which we denote with subscripts 0) and fluctuations around that background (which we denote with subscripts 1).We decompose our fully non-linear thermodynamics about that background state, ln ρ = ln ρ 0 + ln ρ 1 and T = T 0 + T 1 .The temperature gradient of the hydrostatic background is set to the adiabatic temperature gradient, ∇T 0 = ∇T ad = g/c P where g = −gẑ is the gravitational acceleration vector which is taken to be constant and uniform.The domain spans z ∈ [0, L z ] and is a horizontally periodic box with aspect ratio Γ = L x /L z = 4.For 3D simulations we set L x = L y .At the top and bottom boundary we enforce no-slip boundary conditions, with mixed thermal boundary conditions, where we fix the temperature gradient at the bottom and fix the temperature at the top of the domain, We nondimensionalize our equations using values of ρ and T taken at top of the domain to specify our reference temperature, T c , and density ρ c such that T 0 (z = L z ) = 1 and ln ρ 0 (z = L z ) = 0. We choose a nondimensional time scale by setting the ideal gas constant R/m = 1.This choice ensures that the isothermal sound speed, c 2 s = ∂P/∂ρ = R/mT = 1 when evaluated at the top of the domain.We choose our characteristic length scale to be L c = T c /∇T ad .We specify the depth of our domain so that the adiabatic reference state spans N ρ = 3 density scale heights according to the formula: We set the ratio of viscous to thermal diffusivities with the Prandtl number Pr = µc P /κ = 1.
To control the simulations we specify the strength of the source term in the energy equation with a nondimensional heating strength parameter, H, such that where H is constant.H measures a characteristic temperature gradient scale from the heat source, H = ∇T 1 (z = L z )/∇T a d.Since we are considering a constant and uniform heat source, the flux in excess of that carried by the adiabat linearly increases with height and the slope of the flux is determined by H.By separating Q into a part proportional to κ and a part proportional to H we ensure that the hydrostatic thermal equilibrium remains constant as we change the degree of turbulence via changes to the diffusivities.By integration we can produce a constant temperature scale for this static solution The quantity ∆T static is analagous to the quantity ∆ in equation 1.5 of Goluskin [32].Using H as an input parameter removes the scaling factor of kappa from the heating function, which ensures that the Rayleigh number only depends on diffusivities when H is fixed.We specify a Rayleigh number at the top of the domain to control the degree of turbulence in a simulation.The Rayleigh number, Ra, takes the form H is a fundamental input parameter, and the Ra here is directly comparable to the one used to study Boussinesq, internally heated convection in Goluskin and Spiegel [30], Goluskin [32], and Kazemi et al. [40].We perform a linear stability analysis using eigentools [43] to determine the value of Ra crit ; see appendix A for details on how we perform this analysis.We find that the critical Ra is near constant at Ra crit ≈ 50 for values of H < 1.At higher values of H the critical Rayleigh number increases with H.There are 8 nondimensional numbers which control this experiment, Rayleigh number, H, Prandtl number, γ, N ρ , 2 aspect ratios (for x and y), and ∇T (z = 0)/∇T ad which we have set to 1 implicitly with our setup.We derive the number of nondimensional parameters in a similar manner to Graham [44], however we have two additional nondimesional parameters from the addition of a 3rd dimension and our internal heat source.The non-dimensionalized equations are For our numerical experiment we vary Ra and H and hold all other nondimensional inputs constant.With this choice of input parameters, varying Ra at constant H changes the value of κ and µ, while changing H at constant Ra changes the magnitude of the source term while also changing κ and µ to ensure that the ratio of the buoyancy timescale to the diffusive timescale remains constant.By using Ra and H as our fundamental input parameters, we are able to largely separate the posterior Reynolds and Mach number of our simulations.We run 2D simulations across 9 values of H and 6 values of Ra/Ra crit .These 53 simulations are shown in parameter space in the upper left panel of Fig. 1.We run select 3D simulations to validate the behavior of the 2D simulations; these simulations are denoted as boxes, ×, and + to indicate different values of H.
We decompose equations 14 -16 into linear and non-linear terms using the adiabatic background state and then use an implicit-explicit timestepper to implicitly timestep the linear terms and explicitly timestep the nonlinear terms.Our timestep is therefore determined by the CFL constraint of the nonlinear advective term rather than the linear sound waves or diffusion terms, enabling us to simulate low Mach flows efficiently.We timestep the decomposed version of equations 14-16 with the Dedalus v21 pseudospectral solver [45] at resolutions up to 2048 × 1024 spectral coefficients for 2D simulations and up to 512 × 512 × 256 for 3D simulations.We use an implicit-explicit, four-stage, third-order Runge-Kutta RK443 timestepper [46] with a CFL safety factor of 0.75 for 2D simulations and the two-step semi-implicit backwards differentiation SBDF2 timestepper [47] with a CFL safety factor of 0.3 for 3D simulations.All fields are represented as spectral expansions of N z Chebyshev coefficients in the vertical (z) basis, and N x , N y Fourier coefficients in the horizontal (x and y) bases thus making our domain horizontally periodic.To avoid aliasing errors, we use the 3/2-dealiasing rule in all directions.We initialize the simulations with random noise temperature perturbations with magnitude 10 −3 Q 1/3 sin(πz/L z ) added to the initial temperature profile.We define a heating timescale from mixing length theory t h = Q −1/3 to set the output cadence for the simulations.The code used in this paper can be found at https://github.com/whitney-powers/internally_heated_fc_convection.
Simulations whose initial temperature profiles are the adiabatic state T 0 experience a long thermal relaxation time while a superadiabatic thermal boundary layer develops at the upper boundary.To reduce the time spent bringing simulations into thermal equilibrium, we construct an initial stratification characterized by a temperature profile which is adiabatic in the bulk but has a superadiabatic upper boundary layer.We will refer to these initial conditions as "fast-IC".The width of the boundary layer is predicted from a Nusselt number scaling law extrapolated from low-Ra, low resolution simulations which used an adiabatic polytrope for their initial stratification.We solve a nonlinear boundary value problem to find the density stratification for a specified temperature profile which satisfies hydrostatic equilibrium while containing the same total mass as the adiabatic initial conditions.Example stratifications for adiabatic and fast initial conditions are shown in the right panel of Fig. 1.
The "Fast-IC" approach offers significant improvement in the run time of 2D simulations.We discuss the convergence and performance properties of "Fast-IC" in more detail in section III C. We use fast-IC for 2D cases with Ra/Ra crit > 2 × 10 5 and for all 3D cases except the H = 10.667.We calibrated the fast-IC by measuring and extrapolating a Nu (Nusselt number) vs Ra power law (Eqn.28)from 2D simulations with Ra/Ra crit ≤ 2 × 10 5 .For each value of H we fit a scaling law for Nu of form Nu = A * Ra α , where α ≈ 1/5 and the prefactor A ≈ 1.The exact values of A and α are listed in Appendix B. This extrapolation works well for 2D, but for 3D runs the 2D calibration produces initial conditions with too large of a thermal boundary layer implying a subtly different Nu vs. Ra law in 2D and 3D.This discrepancy is especially pronounced at high H where these initial conditions produce a highly unstable region which leads to non-physical negative temperature perturbations which cause a timestepping instability.For the 3D, Ra = 2 × 10 4 Ra crit , H = 10.667case we reverted to the standard adiabatic polytrope for the initial conditions.Nu vs. Ra in the fast-IC approach should be properly calibrated from 3D results to improve performance.For the details of generating "fast-IC" initial conditions see Appendix B.

A. Dynamics
In Fig. 2 we visualize the entropy fluctuations about the adiabatic entropy in four 2D and two 3D simulations.Since entropy in the interior is constant (see Fig. 4) we estimate this adiabatic entropy by taking a time averaged value at the mid z-plane.We observe a cold (blue) upper boundary layer where heat is transported by conduction.In both 2D and 3D systems the dynamics are dominated by convective downflow lanes in blue while convective upflows in red are weaker in magnitude.We use an asymmetric colorbar so that features can be resolved, but note that features in red are at much lower magnitude than features in blue.We do not observe heated upflow plumes at the lower boundary, and the boundary itself is not heated as the superadiabatic flux is zero.Flows at the upper boundary appear similar to standard boundary driven flows with cold plumes launching from a cooled layer.The behavior of the hot and cold plumes is indicative of the internal heating.Downflows cool at the upper boundary and as they fall they are heated, becoming buoyantly neutral (white), and hot upflows gradually heat as they rise, becoming a more intense red.We observe flow morphologies which are similarly dominated by cold downflows in 2D and 3D simulations.Our high Rayleigh number 2D simulations have characteristic features of 2D turbulence, such as vortices in the hot upflow plumes and the formation of large flywheel structures, albeit these flywheels are weaker than is typically observed in boundary driven convection [42].We observe multiple downflow plumes in highly turbulent states in contrast to the single plume flywheel structure previously observed in boundary driven convection.
In Fig.In Fig. 2 d we show turbulent 2D dynamics with a higher value of H = 16.5.The dynamics appear similarly turbulent between Fig. 2 c and d; both have wrapped vorticies and similar levels of small scale detail.This indicates that H has little effect on the Reynolds number of 2D flows.Where these two simulations differ is in the magnitude of fluctuations.The fluctuations in panel d are an order of magnitude larger than in panel c as indicated by the colorbar limits (see the ± values left of "2D").So varying Ra and H control mostly independent properties of the dynamics.
In the right column we show the effect of increasing Rayleigh number in 3D.The simulation in Fig.The asymmetry between upflows and downflows can be observed more quantitatively through probability distribution functions (PDFs) of the entropy fluctuations and vertical velocity as shown in Fig. 3.We plot PDFs of a 2D (green line) and 3D (orange line) simulation at Ra/Ra crit = 2 × 10 5 and H = 0.75 which correlates to panels b and f of Fig. 2. The PDF of entropy fluctuations about the adiabatic value(left panel) are nearly identical in the 2D and 3D case.The peak of the entropy distribution is found near zero (dashed grey line), which is the adiabatic entropy.The distribution is skewed to the left, which indicates that the entropy fluctuation of the cold downflows are larger in magnitude on average than the hot upflows.In contrast the PDFs of vertical velocity differ substantially between 2D and 3D.The asymmetry between upflows and downflows is visible in the vertical velocity distribution, where the upflows on the right side of the dashed grey line (zero velocity) are closer on average to zero than the downflows on the left hand side.The 2D and 3D cases have distinct velocity distributions.The 3D case is a narrower distribution overall with a kurtosis of 7.52 as compared to 4.54 for the 2D case, indicating that the velocities are smaller overall, and the skew towards negative velocity is more pronounced, with skewness of -1.74 as compared to -0.60 for the 2D case.This likely is related to the unique dynamics created with flows restricted to a plane such as the formation of long-lived flywheel structures in 2D, but not in 3D.

B. Structure
The downflow-dominated dynamics can be explained by horizontally averaged profiles of the energy fluxes and the entropy gradient, which are shown in Fig. 4. In the left panel, we show the vertical component of the horizontally and time averaged fluxes.We denote a horizontally averaged quantity with a bar, ie a = 1 LxLy a dx dy where L y and ∂y are dropped for 2D simulations.To derive the fluxes for this system, we first define a total energy equation by summing ρu dotted into Eqn. 2 and ρc V multiplied into 3.After invoking Eqn. 1 and some manipulation we find ∂ ∂t This equation is in conservation form with the sum of the kinetic energy, potential energy, and internal energies in the time derivative, the sum of the fluxes in the divergence term, and our source term Q on the right hand side.We define the conductive, convective, and total fluxes as The components of the convective flux are the enthalpy flux, kinetic energy flux, and potential energy flux: We further define the adiabatic flux, which is the flux conducted along the adiabatic temperature gradient, and the flux imposed by our heat source which is defined such that our fixed flux lower boundary is satisfied as Our systems are in thermal equilibrium, denoted by the fact that on average the total flux through the system matches F imposed .In the interior, conduction transports a flux equivalent to the adiabatic flux, and convection efficiently transports the remainder of the flux.The kinetic energy flux (dashed pink line) is small and directed downward, thus the enthalpy flux is large and positive.In the upper boundary layer convective motions come to a stop to satisfy the impenetrable boundary conditions and the dominant form of heat transport transitions to conduction giving rise to a superadiabatic upper boundary layer (blue layer visible in Fig. 2 at top boundary).These same features are seen in the right panel of Fig. 4 where we plot the specific entropy gradient ∇s, where The entropy gradient follows the adiabatic value in the bulk but becomes strongly superadiabatic near the upper boundary.Strong downflow plumes launch from this superadiabatic boundary layer.Note the lack of a lower boundary layer, which is consistent with the absence of strong upflow plumes in Fig. 2. The size of the upper boundary layer can be measured by the integral of the entropy gradient, ∆s = ∂s∂z = s(z = L z ) − s(z = 0), which scales with both H (which varies the magnitude of the boundary layer), and Ra (which varies the width of the boundary layer).We will discuss the scaling of ∆s with our control parameters as well as its impact on the Mach number in section III D.

C. Scalings
We calculate time and volume averaged Mach number, Reynolds number, and Nusselt numbers, which we respectively define as We display time traces for three 2D simulations at H = 0.75 at Ra/Racrit = 1.9 × 10 3 , 1.9 × 10 5 and 1.9 × 10 7 .The Ra/Racrit = 1.9 × 10 7 simulation uses fast-IC, and the Ra/Racrit = 1.9 × 10 3 simulation uses an adiabatic polytrope for its initial condition.We show both adiabatic (orange line) and fast initial conditions (brown line) for the simulation at Ra/Racrit = 1.9 × 10 5 in the Reynolds number plot.
The simulations using fast-IC reach a thermally converged state quickly.At high Ra, the instantaneous value of Nu becomes chaotic, so we also plot the rolling time average over 40 heating timescales which is denoted by the darker colored trace.
Here, ⟨a⟩ = 1 LxLyLz a dx dy dz represents a volume-average (where we drop L y and ∂y for 2D simulations), and we only use the vertical components of the fluxes (defined in equations 18 -26) to calculate Nu.Our Nusselt number takes the same form as the compressible Nusselt number used in Graham [44], Hurlburt et al. [48] and Anders and Brown [42].This Nusselt number is the ratio of the total flux in excess of the adiabatic flux divided by the conductive component of the flux in excess of the adiabat.As we approach the Boussinesq limit, F ad approaches zero, and we recover the standard Boussinesq Nusselt number.
The time evolution of Re and Nu are displayed in Fig. 5.We show time traces from three 2D simulations at H = 0.75 with different values of Rayleigh number.The highest Rayleigh number run (purple line) uses Fast-IC, and it reaches a thermally equilibriated solution at approximately t = 1000 whereas the lower Rayleigh number run at Ra/Ra crit = 2 × 10 5 (orange line) takes until approximately t = 2500 despite having a shorter diffusive timescale.These simulations converge to a statistically similar state.We find that the simulation shown in orange in Fig. 5 reaches Re = 582 with the standard initial conditions and Re = 589 with "Fast-IC".The standard deviation of the fluctuations in Re are σ = 14, so these two simulations reach values of Re that are well within one standard deviation.While fluctuations of Re are relatively small compared to the mean, the fluctuations of Nu are large and grow significantly with Rayleigh number, indicating that nonlinearity becomes more important.In this highly nonlinear state, such as the orange and purple trace in Fig. 5, minima in the Nusselt number are caused when cold fluid is swept upwards and hot fluid is swept downwards.Because of the high Peclet number (Pe ≈ 5 × 10 3 ), these hot and cold fluid parcels retain their thermal signature for several domain crossings, as is seen in Refs.[49][50][51].
In the upper left panel of Fig. 6, we show the scaling of average Ma against H.We find a scaling that is roughly consistent with a H 1/2 power law at H ≲ 1 and a weaker scaling at large H.This transition at H ∼ 1 is where ∇T 1 (z = L z ) ∼ ∇T ad The scaling deviates from the H 1/2 power law more rapidly at lower Ra.While we do not have results from higher values of H, and so we cannot be conclusive, it appears that in the high H, higher Mach number limit, Mach number may become constant as H increases.In the upper right panel, we plot Re vs. H.Re is independent of H when H ≲ 6. Re increases with H as we transition to higher Mach number flows at H ≳ 6.While we start to see some scaling of Re with H, the Reynolds number still scales like Ra 1/2 in this regime.In the lower left panel we show the scaling of average Ma against supercriticality.We see that our 2D simulations have a roughly constant Mach number for a fixed value of H, however Mach number decreases with Ra for our 3D simulations.In the lower right panel we plot average Reynolds number against supercriticality.Our simulations roughly match a Re ∝ Ra 1/2 power law which is the expected behavior.Re scales weakly with H, so we achieve more turbulent flows at fixed supercriticality.
We perform a least-squares fit for Re and Ma for 2D simulations with H < 6.We exclude values at higher H due to the lack of clean power law behavior in Ma as can be seen in Fig. 6 (e.g., top left panel).We find The Reynolds number scales as Importantly for 2D simulations, we find little Re dependence on H, and little Ma dependance on Ra.This validates our experimental design.We performed a least squares fit for the scaling or Re and Ma against Ra for the suite of 3D simulation with constant H = 0.75.We find that while the 2D simulations produce roughly constant Ma for a given value of H, Ma decreases with Ra in 3D as Ma ∝ Ra −0.10 .Additionally we find that the 3D simulations achieve a weaker Re scaling of Re ∝ Ra 0.38 than our 2D simulations.The Nusselt number is a measure of how efficiently convection transports heat.It is expected that the Nusselt number in astrophysical convection follows a diffusion-free power-law scaling which is referred to as the "ultimate regime" [52].We measure a Nusselt number scaling law using our input Rayleigh number, and given this choice we expect the "ultimate regime" to be a Nu ∝ Ra 1/3 power law.We measure the scaling of the Nusselt number with Ra/Ra crit and show these results in Fig. 7.We find Nu ∝ (Ra/Ra crit ) 0.21 for 2D and Nu ∝ (Ra/Ra crit ) 0.20 for 3D, both are very close to a 1/5 power law.As we transition to a higher Mach number regime at the upper range of H we find that the Nusselt number increases with H and that a 1/5 power law does not fit our data as effectively at high H where we observe a steeper trend at high Ra.This suggests that compressible effects modify the boundary layer formation and scaling for high Mach number convection.

D. Effect of Upper Boundary Layer
To study the origin of our observed scaling of Mach number with H we measure the change in specific entropy, s, across the domain.Since the entropy profile is adiabatic everywhere except the upper boundary, this domain-average measures the properties of the upper boundary layer.∆s is a function of both H and Ra, where increasing H increases the magnitude the boundary layer, and increasing Ra decreases its width.In the low-Mach number regime where the anelastic approximation is valid, a balance between buoyancy and nonlinear inertia takes the form [53] where ℓ is a dominant flow length scale.This suggests that we should find Ma ∼ (∆s) 1/2 .In Fig. 8 we plot the behavior of ∆s against our control parameters Ra (left panel) and H (middle panel).Increasing Rayleigh number causes the superadiabatic upper boundary to be thinner, causing a negative scaling, while increasing H causes the amplitude of the superadiabatic boundary to increase, casuing a positive scaling.The scaling of ∆s is similar for both 2D simulations (circles) and 3D simulations (squares, ×, and + symbols).3D simulations have slightly smaller boundary layers and thus smaller values ∆s, but the scaling exponent remains the same.This indicates that the size of the boundary layer scales in a consistent manner for both 2D and 3D simulations.
In the right panel of Fig. 8 we plot Ma against ∆s.We find that our simulations are broadly consistent with a Ma ∝ (∆s) 1/2 power law (dashed black line) suggesting that this is the dominant nonlinear trend driving our dynamics.Our 3D simulations are well described by a (∆s) 1/2 power law, but the 2D simulations have more deviation from the trend line indicating the presence of a second order effect from our control parameters on the Mach number of the flows.To better show this effect we plot only the simulations with H = 0.75 in the inlay.We see that the 3D simulations all fall on the (∆s) 1/2 trend line, however the 2D simulations depart from this trend line where we find that higher Ra simulations produce higher Ma flows than predicted from the anelastic theory.This discrepancy between 2D and 3D simulations cannot be explained by the size of the boundary layer.Instead the way in which the dynamics respond to the boundary layer changes between the 2D and 3D case.Since all the 3D simulations fall on a single trend line, it should be possible to design 3D simulations that reach constant Ma against variations in supercriticality as long as the change in entropy is held constant.

E. Power Spectra
We calculate the convective power spectra at the mid-z plane of our 2D simulations at H = 0.75.To do so, We output the interpolated velocity at z = L z /2, and compute the discrete Fourier transform.We define the Fourier transform of the velocity as Nx kxn (32) We define the power spectral density as the product of the Fourier transform of velocity and its complex conjugate normalized by the horizontal resolution squared, P(k , where * denotes the complex conjugate and N x is the horizontal resolution. For 3D simulations we perform a 2D Fourier transform of the velocity at z = L z /2 for simulations at H = 0.75.We define the 2D discrete Fourier transform as Ra/Ra crit FIG. 8.In the left panel we display a scatter plot of ∆s against Ra/Racrit.Colors correspond to values of H, and squares, x, and + symbols correspond to 3D runs.Circles denote 2D runs, with large circles indicating a 2D run with a corresponding 3D run at the same control parameters.In the middle panel we show a scatter plot of ∆s against H, with the same marker styling as in the left panel, but with colors now corresponding to Rayleigh number.In these two panels we show that the scaling of ∆s with our control parameters is consistent between 2D and 3D simulations and the small differences that are seen are not from a change in the power law exponent.In the right panel we show how Mach number scales with ∆s.We find a scaling consistent with Ma ∝ (∆s) 1/2 for 3D simulations and for 2D simulations at a fixed value of Ra.However, when we increase Ra in 2D simulations ∆s decreases.We show an inlay plot with just the simulations with H = 0.75.We see that for 3D, all simulations fall on the Ma ∝ (∆s) 1/2 trend line, but the 2D simulations do not, with increased Ra leading to smaller ∆S and slightly decreased Ma.
We define the power spectral density as ) and then transform the power spectra density from k x and k y space to k θ and k h space where k h = k 2 x + k 2 y is the horizontal wavenumber and k θ = arctan(k y /k x ) we then interpolate and average over k θ to produce a horizontal wavenumber power spectral density P(k h ).
For 2D simulations at low wavenumber we find that the power spectrum magnitude changes little as Ra is increased.We expect this behavior in 2D since Mach number and Reynolds number are independent of each other under our choice of parameterization.The spectra decay with roughly a k −3 power law, which is associated with a two dimensional Ra/Ra crit FIG. 9. (left panel) Convective power spectra of 2D simulations accumulated over 40 heating timescales.We plot the power spectra of the mid-z plane velocity vector for a set of simulations with constant H.We find consistent power at low wavenumber for these simulations.The less turbulent simulations enter the viscous regime at lower kx than the more turbulent simulations.We plot canonical k −5/3 and k −3 power laws, and note that the k −3 power law better matches the energy cascade in our simulations.The sharp vertical cutoff of the spectra is found at the highest possible wavenumber given the simulation's resolution.(right panel) Convective power spectra of 3D simulations accumulated over 2 heating timescales.We plot power spectra of the mid-z plane velocity vector for a suite of simulations at constant H.The low wavenumber power decreases with increasing Ra.
turbulent cascade of vorticity from low to high wavenumber [54,55].The characteristic scale of our internal heat source is the size of the domain, so it makes sense that we find the k −3 forward vorticity cascade as we are driving convection over a large length scale.At high k x we see the viscous cutoff wavenumber increase with Ra, but the nature of the spectra at lower k x remains unchanged.This suggests that it may be possible to extrapolate the large-scale dynamics present in our simulations to higher Rayleigh number cases.However, the large scales of our 2D simulations are dominated by flywheel structures, which are not found in astrophysical convection, limiting the utility of this behavior.We find that the spectra of our 3D simulations does not have the same consistent behavior as seen in 2D.As Rayleigh number increases, the magnitude of the low wavenumber modes decreases.We expect this behavior from the 3D Mach number scaling since the magnitude of the power spectrum depends on the Mach number.This unfortunately means that less turbulent simulations cannot be easily used to predict low wavenumber flows in a more turbulent system.

F. Comparison to Other Work
It is useful to compare our results to simulations of fully compressible boundary driven convection, as well as incompressible internally heated convection.The former allows us to interpret the effect of different thermal driving models, while the latter allows us to study how compressibility affects simulations with otherwise similar thermal driving.Our dynamics are distinct from the characteristic behaviour of boundary driven convection.We observe dynamics with strong downflows and weak, diffuse upflows.In our 3D simulations downflows form individual plumes and are not pushed into downflow lane structures.This is counter to what is seen in a boundary driven system [e.g., 42] The downflow dominated dynamics we observe are reminiscent of the convective dynamics observed by Käpylä et al. [56], who studied a depth-dependent internal heating, despite the fact that we study a constant heat source.
We measure scalings of Reynolds number, Mach number, and Nusselt number.Our Reynolds number scalings of Re ∝ Ra 0.45 for 2D and Re ∝ Ra 0.38 for 3D are similar to those in Anders and Brown [42] (AB17) who found a Re ∝ Ra 3/4 power law for 2D low Mach simulations and a Re ∝ Ra 1/2 power law for 3D simulations and Ma ∼ 1 2D simulations.However our scaling exponent in 3D is smaller.This similarity is unsurprising since our Rayleigh number scales viscosity and diffusivity in the same basic manner.Our Mach number scaling of Ma ∝ Ra −0.02 for 2D and Ma ∝ Ra −0.10 for 3D however is in stark contrast with the scalings found by AB17, where 2D simulations followed a Ma ∝ Ra 1/4 power law and 3D simulations had a constant Mach number with a fixed superadiabatic excess (ϵ) which is comparable to H. Like AB17 we find stronger positive scalings in 2D than in 3D.It should also be noted that AB17 used fixed temperature boundary conditions at both boundaries whereas we use mixed flux-temperature boundary conditions.
Our Nusselt number scaling law of Nu ∝ Ra 1/5 is comparable to the Nusselt number scalings of the boundary driven fully compressble convection simulations of AB17 in their high-Ra 2D simulations, however their low-Ra 2D simulations follow a 1/3 power law, and their 3D simulations follow a 2/7 power law.The 2/7 power law was also observed with fully compressible boundary driven convection by [49] and a steeper scaling exponent of 0.38 was observed by [50].We can also compare our Nusselt number scaling to incompressible internally heated convection studies.We find the same Nusselt number scaling as [30,31], however their simulations used fixed temperature at the upper and lower boundaries.Therefore, the uniform heating case of [40] is the more analogous incompressible study for our setup.Note that [31,40] report a scaling law based on a "diagnostic" Rayleigh number which is equivalent to Ra/Nu.We have converted their results to an input Rayleigh number scaling to be consistent with our choice of Ra and Nu.Kazemi et al. [40] found a Nu ∝ Ra 1/4 scaling, which is close to our scaling law, but not identical.It is striking that we find a similar Nusselt number scaling to [40], and a near identical scaling to [30,31] despite the different boundary conditions.

IV. CONCLUSIONS
In this work we simulated fully compressible, internally heated convection.We used 2D and 3D simulations and we varied the heating, H, and Rayleigh number, Ra.We observed convection dominated by cold downflows.We did not observe the formation of hot upflow plumes at the lower boundary.This agrees with our expectations as the lower boundary is thermodynamically stable.The dynamics we observe appear similar to those produced by Käpylä et al. [56], Cossette and Rast [57], where the downflows are buoyantly driven and the upflows are pressure driven and heated as they rise.
We found that the Mach and Reynolds numbers could be independently controlled in 2D by using H and Rayleigh number as input parameters.Varying Ra produces a Re ∝ Ra 1/2 and a weak Ma scaling.Varying H produces a Ma ∝ H 1/2 scaling and a weak Re scaling.In 3D we find a weaker scaling for Reynolds number of Re ∼ Ra 1/3 , and a negative Ra scaling for Mach number of Ma ∼ Ra −1/10 .We found that the entropy jump in the upper boundary layer determines the average Mach number of the flow, and thus the magnitude of the kinetic energy power spectrum.While the change in entropy across the boundary layer scales similarly in 2D and 3D simulations, the impact of the boundary layer on the Mach number differs, with 2D simulations having an additional effect tied to the Rayleigh number.This likely reflects flywheel modes since these make up the large scale motions of our 2D simulations.
More robust mechanisms for controlling convective driving, such as an internal cooling layer which does not vary in size between experiments, could improve our ability to separate the resultant Re and Ma of the flows in 3D, and thus the behavior of power spectra.Internal heating and cooling layers have been used to reduce the effect of the boundaries on heat transport in simulations [38,40] and experiments [35,36] of Boussinesq convection.Since we find consistent heat transport scaling with the analogous internally heated Boussinesq setup of Goluskin & Spiegel [30], we expect that the use of a cooling layer may have similar effects in the fully compressible regime.Critically, however, the relationship between a cooling layer width and heating magnitude may differ.Therefore further study on the shape of a non-uniform heating function is needed.
We additionally find heat transport scaling laws which are consistent with internally heated Boussinesq convection [30,31,40].This suggests that certain properties of incompressible flows also hold for stratified flows with low Mach number and that studies of incompressible flows [e.g., 30,32,36,40] are useful for understanding astrophysical convection.However, this similarity may be limited to the low to intermediate Mach number regime which we studied in this work (Ma ≤ 0.2).
We also found that velocity power spectra of 2D simulations were invariant at low wavenumber while turbulence (Re) increased when the Mach number was held constant.This behavior may be useful for studying large scale flows since their properties do not change significantly as Ra increases, thus allowing us to study them with low resolution, low Ra simulations.However, with 3D simulations, the power at low wavenumber decreases as Ra increases.We may be able to achieve consistent low wavenumber behavior in 3D by using a non-uniform heat source with heating and cooling layers, as the properties of the flow would not be set by a thermal diffusion dominated upper boundary.
Stars and gas giant planets do not have boundaries which act like hard walls, so the no-slip boundaries we used are likely a poor model for astrophysical flows and stress-free boundary conditions are likely more appropriate.We choose no-slip boundary conditions to prevent the onset of parasitic mean shear flows in 2D simulations.Mean shear flows reduce thermal transport by convection and suppress convective flows with predator-prey like dynamics [58,59].Wider aspect ratios have been shown to suppress shear flows in 2D simulations with stress-free boundaries Wang et al. [60], however this approach requires higher resolutions and longer convergence times.Shearing states typically do not emerge in 3D simulations, so future work relying exclusively on 3D simulations could utilize stress-free boundary conditions without this concern.
In this work we studied flows at Mach numbers between 10 −2 and 0.2.We find that the behavior of the Reynolds number and Nusselt number in our highest Mach number simulations deviated from the behavior at lower Mach number.Scaling laws for high Mach number flows could be inconsistent with the heat transport scaling from internally heated Boussinesq convection, and future studies should investigate this regime of flow.Future work should also study how internal heating interacts with other complications such as rotation [61,62] and chemical mixing [3].

Appendix B: Fast Initial Conditions
To speed up the evolution to a thermally equilibrated, nonlinear convective state, we initialize our simulation with an approximation of the time averaged steady state (see Fig. 4).We refer to this method as "Fast-IC".Convection efficiently mixes the temperature gradient to the adiabatic value in the interior.To satisfy thermal equilibrium and our boundary conditions, the time-averaged superadiabatic temperature gradient at the upper boundary must be −H.We impose that the half-width of the thermal boundary layer is set by the Nusselt number as This scaling is chosen to approximate the width of a thermally equilibriated IVP solution.Our approximation of the steady state solution is where erf is the error function.To ensure that the new initial conditions are in hydrostatic equilibrium and have the same mass as our original initial conditions, we solve a nonlinear boundary value problem to find ln ρ with equation A2 and integral constraint A3 while taking eqn.B2 as ∂T /∂z.T eq is obtained by integrating B2 subject to boundary condition A5.After running our lower resolution 2D simulations with the adiabatic polytrope for initial conditions we fit a power law Nu = A(Ra/Ra crit ) α for each value of H with A and α listed in Table I.We use this powerlaw to find the value of Nu in Eqn.B1 and then use the nonlinear boundary value problem to produce initial conditions for high Rayleigh number 2D runs and most 3D runs.We show an example of the specified shape for T and the resulting solution for ρ in Fig 1.
This method offers considerable reduction in CPU time required for running turbulent simulations in 2D.In Fig. 11 we show a comparison of the time-evolution of volume-averaged measures between two simulations with Ra/Ra crit = 2 × 10 5 and H = 0.75, one using "Fast-IC", the other using adiabatic initial conditions.We plot (left) the Reynolds number, a measure of the dynamics, and (right) ∆s, a measure of the structure.We calculate the time required for the Reynolds number to equilibrate to within 3% of its final value.The dynamics convergence times are t d fast = 320 for "Fast-IC" and t d ad = 3800 for adiabatic initial conditions.We calculate a 3% convergence time for ∆s as well.We find that with "Fast-IC" the ∆s begins within 3% of the final value, and the structure convergence time for the adiabatic initial conditions is t s ad = 4500.By predicting the final structure with "Fast-IC", ∆s starts near its final value, whereas the adiabatic initial conditions slowly approach the converged state."Fast-IC" starts with u = 0, so the dynamics still take time to converge, however the time required is significantly shorter than with the adiabatic initial conditions.
In the left column of Fig. 12 we show the dynamics at t d fast for both "Fast-IC" and adiabatic initial conditions.With the adiabatic initial conditions in panel a the simulation has reached convective onset with cold plumes forming  at the top of the domain.The thermal boundary layer is shallow so the entropy perturbation of the plumes is small, therefore they loose their cool signature due to the heating before they transit the domain and halt partway down.
As the thermal structure evolves the plumes reach deeper into the domain, eventually reaching the bottom.With "Fast-IC" in panel c we skip this slow evolution and convective plumes immediately cross the entire domain after convective onset and the dynamics quickly reach a converged state.We show probablity density functions (PDFs) of vertical velocity, w, for both cases in panel e.The two PDFs are distinct; the adiabatic-IC case has a narrower distribution and smaller average velocities, where the PDF of the Fast-IC case has a wider distribution which is close to the final velocity distribution seen in panel f.In the right column we show dynamics for the adiabatic-IC case (panel b) and the fast-IC case (panel d) as well as PDFs for both cases (panel f) at t = 6000.At this time the dynamics and structure have converged for both initial condition choices.The dynamics not only look similar between panels b and d, but the distribution of velocities have converged in panel f.The "Fast-IC" case reached a final distribution with a mean of -0.013, standard deviation of 0.048, skewness of -4.24, and kurtosis of 21.74.The adiabatic initial conditions reached a final state with a mean of -0.012, standard deviation of 0.048, skewness of -4.10, and kurtosis of 20.37.Racrit and H = 0.75 with adiabatic initial conditions (panel a) and "Fast-IC" (panel c).The run with adiabatic initial conditions has reached convective onset, however cold downflow plumes stall out part way through the atmosphere, becoming more buoyant (redder) than the nearby fluid.As time progresses convective downflows reach further into the atmosphere.In the run with "Fast-IC" the downflow plumes cross the domain immediatly after convective onset and the dynamics quickly reach a converged state.The colorbar is calculated from the adiabatic initial conditions case in the same manner as in Fig. 2 and is used in all dynamics panels.In panel e we show probability density functions (PDFs) for entropy fluctuations averaged over 12 heating timescales for both adiabatic-IC and fast-IC simulations at t = t d fast .The adiabatic-IC run has not converged at this time as evidenced by the discrepancy betweeen the PDFs.In the right column (panels b, d, and f) we show dynamics and PDFs at t = 6000, which is after dynamics and structure have converged for both choices of initial conditions.The dynamics of the adiabatic-IC simulation (Panel b) are comparable to the dynamics of the Fast-IC run (Panel d).In panel f we show the PDFs for both runs where the distributions have converged, indicating that the dynamics have reached a statistically similar final state.

Appendix C: Table of Simulations
In table III we display all simulations used in this study where we show control parameters Ra/Ra crit and H, horizontal resolution n horiz and vertical resolution n z , dimensionality, heating timescale t h eat, simulation run time t sim , averaging window t avg , and time and volume averaged Reynolds number, Mach number, Nusslet Number, and change in entropy across the domain.All timescales are expressed in nondimensionalized time units.
FIG. 1. (Left panel) Symbols indicate the simulations conducted in this work in H -Ra/Racrit parameter space.2D simulations are indicated with a solid circle.3D simulations are indicated with the square outline, ×, and + symbols.The marker color indicates the value of Ra and H for each point.The left side of each 2D point indicates the value of Ra (upper color bar) and the right side indicates the value of H (right color bar).3D points are colored by the value of Ra.While these colors are redundant with the position on the plot, these colors are used consistently throughout the rest of the figures to refer to Ra/Racrit and H. (Right panel) Initial conditions for T and ρ where dashed lines correspond to an adiabatic polytrope stratification and solid lines show the "Fast-IC" initial conditions, where we approximate the temperature structure of the final state and solve a boundary value problem for the density profile which satisfies hydrostatic equilibrium and mass conservation, see Appendix B. The "Fast-IC" profile shown is for Ra/Racrit = 2 × 10 5 and H = 10.667.
2 a-c, we show how dynamics change with increasing supercriticality in 2D simulations at low H.The simulation in Fig. 2 a is laminar.As the Rayleigh number increases in Fig. 2 b and c we see more complex turbulent dynamics emerge.These include long-lived wrapped vortices, where hot fluid is encircled by cold fluid indicating high Péclet number flows.
2 e shares the same values of Ra and H as Fig. 2 a, and Fig. 2 f shares the same values as Fig. 2 b.As 3D flows become more turbulent we see more convection cells at the upper boundary and a transition from coherent cell structure at the midplane to flows dominated by discrete downflow plumes distributed randomly across the midplane.The behavior of the cool downflows in Fig. 2 f are consistent with weak upflows, as indenpendent downflow plumes are not pushed into coherent downflow lanes.

FIG. 2 .
FIG.2.We show entropy fluctuations with the adiabatic entropy subtracted.For each panel we set the colorbar limits at the 3rd and 97th percentile values of entropy from that simulation.These minimum and maximum values are displayed above each panel.Note the substantial imbalance between the extent of the minimum and maximum values.We show 2D dynamics in the left column, and 3D dynamics in the right column.In panels a), b), and c) we show 2D dynamics at H = 0.75 and with Ra/Racrit = 2 × 10 3 , 2 × 10 5 , and 2 × 10 7 , respectively.As Rayleigh number increases, the flows become more turbulent.In panel d) we show a simulation with the same Rayleigh number as panel c), but with stronger internal heating at H = 16.5.In panel e) we display the 3D simulation with the same parameters as panel a).and in panel f) we display the 3D simulation with the same parameters as panel b).An animation of this figure can be found at https://vimeo.com/793788958

1 FIG. 3 .
FIG.3.Time-averaged probability distribution functions (PDFs) of entropy fluctuations with the adiabatic entropy subtracted (left panel) and vertical velocity (right panel) are shown for 2D (green line) and 3D (orange line) simulations.Both simulations use Ra/Racrit = 2 × 10 5 and H = 0.75.The PDFs are averaged over 12 heating timescales.The entropy PDF is skewed towards negative values and is nearly identical between the 2D and 3D case.In contrast, the velocity PDFs exhibit different behavior in 2D and 3D.The 3D case has a narrower velocity distribution with more pronounced skew towards downward flows.

FIG. 4 .
FIG. 4. (Left panel) Time-and horizontally-averaged vertical energy flux profiles for the 3D simulation with Ra/Racrit = 2×10 4 and H = 10.667plotted against height.The imposed flux from the source term is shown with the solid black line and the adiabatic flux (defined in equation 25) with the dashed black line.The time-averaged total flux F total is shown with the orange line which converges to imposed heat flux F imposed .The green line shows the conductive flux F cond which is adiabatic at the lower boundary and interior.The conductive flux becomes large in the superadiabatic upper boundary layer.We show the convective flux Fconv with the purple line.We also display the kinetic energy flux, a component of the convective flux, with a dashed pink line.Due to the asymmetry introduced by stratification the KE flux is negative.In the right panel we show the entropy gradient.The entropy gradient is zero everywhere except the superadiabatic upper boundary.The integral of the entropy gradient, denoted by the shaded grey region, is the change in entropy, ∆s, which is important in predicting the Mach number of the flows as we discuss in section III D.

3 NuFIG. 5 .
FIG.5.Time evolution of Reynolds number (left panel) and Nusselt number (right panel).We display time traces for three 2D simulations at H = 0.75 at Ra/Racrit = 1.9 × 10 3 , 1.9 × 10 5 and 1.9 × 10 7 .The Ra/Racrit = 1.9 × 10 7 simulation uses fast-IC, and the Ra/Racrit = 1.9 × 10 3 simulation uses an adiabatic polytrope for its initial condition.We show both adiabatic (orange line) and fast initial conditions (brown line) for the simulation at Ra/Racrit = 1.9 × 10 5 in the Reynolds number plot.The simulations using fast-IC reach a thermally converged state quickly.At high Ra, the instantaneous value of Nu becomes chaotic, so we also plot the rolling time average over 40 heating timescales which is denoted by the darker colored trace.

FIG. 11 .
FIG.11.An example of time traces of Reynolds number (left panel) and ∆s (right panel) for simulations with Ra = 2×105 Racrit and H = 0.75 with "Fast-IC" (orange) and adiabatic initial conditions(green).Grey horizontal lines show the time-and volumeaveraged values reported in Appendix C. Arrows show the time where a rolling average over 500 data points has converged to within 3% of the the accepted value.

FIG. 12 .
FIG.12.(Panels a and c) An example of entropy fluctuations at t = t d fast for simulations with Ra = 2 × 10 5 Racrit and H = 0.75 with adiabatic initial conditions (panel a) and "Fast-IC" (panel c).The run with adiabatic initial conditions has reached convective onset, however cold downflow plumes stall out part way through the atmosphere, becoming more buoyant (redder) than the nearby fluid.As time progresses convective downflows reach further into the atmosphere.In the run with "Fast-IC" the downflow plumes cross the domain immediatly after convective onset and the dynamics quickly reach a converged state.The colorbar is calculated from the adiabatic initial conditions case in the same manner as in Fig.2and is used in all dynamics panels.In panel e we show probability density functions (PDFs) for entropy fluctuations averaged over 12 heating timescales for both adiabatic-IC and fast-IC simulations at t = t d fast .The adiabatic-IC run has not converged at this time as evidenced by the discrepancy betweeen the PDFs.In the right column (panels b, d, and f) we show dynamics and PDFs at t = 6000, which is after dynamics and structure have converged for both choices of initial conditions.The dynamics of the adiabatic-IC simulation (Panel b) are comparable to the dynamics of the Fast-IC run (Panel d).In panel f we show the PDFs for both runs where the distributions have converged, indicating that the dynamics have reached a statistically similar final state.
[30]7.We show Nusselt number against Rayleigh number for our suite of simulations.We find a Nu ∝ (Ra/Racrit) 1/5 scaling, which is consistent with the incompressible internally heated simulations from Goluskin and Spiegel[30].
FIG.10.Left panel: Example of a thermal equilibrium solution with H = 10.667 to be used as input to the linear stability analysis.We plot temperature, T (green) and density, ρ (orange), against height, z.Middle panel: Example of growth rates with H = 10.667plotted against horizontal wavenumber on the x-axis and Rayleigh number on the y-axis.Black dots indicate the interpolated value of Ra where the growth rate is zero for the given wavenumber.Right panel: Values of Racrit plotted against H.We find that Racrit is roughly constant below H ∼ 1, and increases for H ≳ 1.