Convective heat transport in stratified atmospheres at low and high Mach number

We study fully compressible convection in the context of plane-parallel, polytropically stratified atmospheres. We perform a suite of 2D and 3D simulations in which we vary the initial superadiabaticity ($\epsilon$) and the Rayleigh number (Ra) while fixing the initial density stratification, aspect ratio, and Prandtl number. The evolved heat transport, quantified by the Nusselt number (Nu), follows scaling relationships similar to those found in the well-studied, incompressible Rayleigh-B\'{e}nard problem. This scaling holds up in both 2D and 3D and is not appreciably affected by the magnitude of $\epsilon$.


I. INTRODUCTION
Convection transports energy in stellar and planetary atmospheres where flows are compressible and feel the atmospheric stratification. This stratification is significant in regions such as the convective envelope of the Sun, which spans 14 density scale heights. In the bulk of these systems, particularly in the deep interior, flows are at very low Mach number (Ma). Unfortunately, numerical constraints have restricted most studies of compressible convection to high Ma. These prior studies [1][2][3][4][5][6] have provided insight into the nature of convection in the low temperature, high Ma region near the Sun's surface. Few fundamental properties of low Ma compressible convection, such as the scaling of convective heat transport, are known.
In the widely-studied Rayleigh-Bénard problem of incompressible Boussinesq convection (RBC), a sufficiently negative temperature gradient causes convective instability. In the evolved solution, upflows and downflows are symmetrical, the temperature in the interior becomes isothermal, and the conductive flux (∝ ∇T ) approaches zero there. For compressible convection in a stratified atmosphere, a negative entropy gradient causes convective instability. Early numerical experiments of moderate-to-high Ma compressible convection in two [1][2][3][4] and three [5][6][7] dimensions revealed a different evolved state from RBC. Downflow lanes become fast and narrow, and upflow lanes turn into broad, slow upwellings. Furthermore, the entropy gradient is negated by convection in the interior, so a significant temperature gradient and conductive flux can persist despite efficient convection.
In RBC, there exist two primary dynamical control parameters: the Rayleigh number (Ra, the ratio of buoyant driving to diffusive damping) and the Prandtl number (Pr, the ratio of viscous to thermal diffusivity). These numbers control two useful measures of turbulence in the evolved solution: the Reynolds number (Re, the strength of advection to viscous diffusion) and the Peclet number (Pe, advection vs. thermal diffusion). In stratified atmospheres, the magnitude of the unstable entropy gradient joins Ra and Pr as a third important and independent control parameter. This superadiabatic excess, , sets the scale of the atmospheric entropy gradient [1]. We find here that primarily controls the Ma of the evolved solution.
Here we study the behavior of convective heat transport, quantified by the Nusselt number (Nu), in plane-parallel, two-and three-dimensional, polytropically stratified atmospheres. We vary and Ra while holding Pr, aspect ratio, boundary conditions, and initial atmospheric stratification constant. We also examine the behavior of flow properties, as quantified by Ma and Re. We find here that the scaling of Nu in stratified, compressible convection is similar to that in Rayleigh-Bénard convection, and that this scaling is not appreciably changed by the magnitude of the superadiabaticity.

II. EXPERIMENT
We examine a monatomic ideal gas with an adiabatic index of γ = 5/3 whose equation of state is P = RρT . This is consistent with the approach used in earlier work [1][2][3][4][5][6][7] and is the simplest stratified extension of RBC. The atmospheres studied here are initially polytropically stratified, where m is the polytropic index and L z is the depth of the atmosphere. The polytropic index is set by the superadiabatic excess, = m ad − m, where m ad = (γ − 1) −1 is the adiabatic value of m. The height coordinate, z, increases upwards in the range [0, L z ]. Subscript 0 indicates initial conditions and subscript t indicates values at z = L z . Stratified atmospheres have a fourth non-dimensional parameter, the number of density scale heights, n ρ = ln [ρ 0 (z = 0)/ρ t ]. We specify the depth of the atmosphere, L z = e nρ/m − 1, by choosing the initial value of n ρ . Throughout this work we set n ρ = 3. Satisfying hydrostatic equilibrium sets the value of gravity, g = RT t (m + 1), which is constant with depth. We study atmospheres with aspect ratios of 4 where both the x and y coordinates have the range [0, 4L z ]. In our 2D cases, we only consider x and z. These domains are nondimensionalized by setting R = T t = ρ t = 1 at z = L z . By this choice, the non-dimensional length scale is the inverse temperature gradient scale and the timescale is the isothermal sound crossing time, τ I , of this unit length. Meaningful convective dynamics occur on timescales of the atmospheric buoyancy time, t b = τ I L z m c P /g n ρ , where c P = Rγ/(γ − 1) = 2.5 is the specific heat at constant pressure.
At fixed n ρ , convective dynamics are controlled by as well as the atmospheric diffusivities. At a fixed value of , the diffusivities are set by the Rayleigh number (Ra) and the Prandtl number (Pr), where ∆S 0 = ln(1 + L z ) = n ρ /m is the initial specific entropy difference across the domain, χ is the thermal diffusivity, and ν is the kinematic viscosity. Throughout this work we specify that Pr = 1 and is depth invariant. The initial thermal conductivity, κ 0 = χρ 0 , is constant with depth, such that (1) is in thermal equilibrium (∇·[κ 0 ∇T 0 ] = 0). By these choices, ν(z) ≡ χ(z) ≡ χ t /ρ 0 . This formulation sets Ra at the bottom of the domain greater than Ra t by a factor of e 2nρ . Henceforth when we specify Ra we are referring to Ra t . The full values of κ = ρχ and µ = ρν (the dynamic viscosity) change as the density profile evolves. The diffusivities scale as χ t , ν t ∝ gL 3 z (∆S 0 /c P )/Ra t . Defining the thermal diffusion timescale as t χ ≡ τ I L 2 z /χ, the ratio of t χ to the buoyancy time is We carry out two experiments in this study. In the first, we fix and increase Ra, thus increasing the ratio in (3). In the second, we fix Ra and vary , scaling the dynamical timescales (t b , t χ ) as −1/2 relative to the speed of sound; we see this reflected in the evolved Mach number scaling ( Fig. 1). We use ln ρ and T as our thermodynamic variables and solve the Fully Compressible Navier-Stokes equations, with the viscous stress tensor given by where δ ij is the Kronecker delta. Taking an inner product of (5) with ρu and adding it to ρc V ×(6) reveals the full energy equation, where F conv ≡ F enth +F KE +F PE +F visc is the convective flux and F cond = −κ∇T is the conductive flux. The individual contributions to F conv are the enthalpy flux, F enth ≡ ρu(c V T + P/ρ); the kinetic energy flux, F KE ≡ ρ|u| 2 u/2; the potential energy flux, F PE ≡ ρuφ (with φ ≡ −gz); and the viscous flux, F visc ≡ −ρνu ·σ. Understanding how each of these fluxes interact is crucial in characterizing convective heat transport. We utilize the Dedalus 1 pseudospectral framework [8] to time-evolve (4)-(6) using an implicit-explicit (IMEX), third-order, four-step Runge-Kutta timestepping scheme RK443 [9]. Thermodynamic variables are decomposed such that T = T 0 + T 1 and ln ρ = (ln ρ) 0 + (ln ρ) 1 , and the velocity is u = wẑ + ux + vŷ. In our 2D runs, v = 0. Subscript 0 variables, set by (1), have no time derivative and vary only in z. Variables are time-evolved on a dealiased Chebyshev (vertical) and Fourier (horizontal, periodic) domain in which the physical grid dimensions are 3/2 the size of the coefficient grid. Domain sizes range from 64x256 coefficients at the lowest values of Ra to 1024x4096 coefficients at Ra > 10 7 in 2D, and from 64x128 2 to 256x512 2 in 3D. By using IMEX timestepping, we implicitly step the stiff linear acoustic wave contribution and are able to efficiently study flows at high (∼ 1) and low (∼ 10 −4 ) Ma. Our equations take the form of the FC equations in [10], extended to include ν and χ which vary with depth, and we follow the approach there. This IMEX approach has been successfully tested against a nonlinear benchmark of the compressible Kelvin-Helmholtz instability [11]. We impose impenetrable, stress free, fixed temperature boundary conditions at the top and bottom of the domain, with w = ∂ z u = T 1 = 0 at z = {0, L z }. T 1 is initially filled with random white noise whose magnitude is infinitesimal compared to T 0 . We filter this noise spectrum in coefficient space, such that only the lower 25% of the coefficients have power. All reported results are taken from time averages over many t b beginning {100, 40}t b after the start of our {2D, 3D} simulations to ensure our results are not biased by the convective transient.

A. Evolved fluid numbers & flow morphology
We measure the adiabatic Mach number (Ma = |u|/ √ γT ), and find that it is a strong function of and a weak function of Ra. In our 2D runs, when Ma < 1, we observe a scaling law of Ma(Ra, ) ∝ ∼ 1/2 Ra 1/4 . This relation breaks down as the mean Ma approaches 1 (Fig. 1). This transition occurs near Ra/Ra crit ≈ {10 2 , 10 3 } for = {1, 0.5}. In our limited 3D runs, Ma appears to be a function of alone, with Ma ∝ ∼ 1/2 , so at high Ra, Ma 3D < Ma 2D . We conjecture that the scaling of Ma with Ra in the 2D runs is due to the formation of coherent high-velocity "spinners," which form between upflow and downflow lanes. These structures, which are reminiscent of flywheel modes in RBC, do not appear in our 3D runs at these parameters [12,13]. Simulations in the range of Ra/Ra crit > 10 3 at = 10 −4 exhibited "windy" states of convection, in which a large-scale horizontal shearing flow replaced the more standard upflow/downflow morphology of convection. Similar states have been studied in RBC [14]. These runs are represented in Figs. 1, 3, & 4 as hatched points, and while this phenomenon does not appear to greatly modify the scaling of fluid properties measured in this work, these states warrant further investigation.
In 2D, low Ma flows (e.g., = 10 −4 ) display the classic narrow downflow and broad upflow lanes of stratified convection (Fig. 2a). At high Ma (e.g., = 0.5, Ra/Ra crit 10 3 ), bulk thermodynamic structures are similar but shock systems form in the upper atmosphere near downflow lanes (Fig. 2b&c), as reported previously [4,15]. At large Ra, the diffusion timescale becomes long (3), and thermodynamic structures form small eddies which traverse the domain repeatedly before diffusing (Fig. 2c). As evidenced by the colorbar scalings, the amplitudes of thermodynamic fluctuations scale with .
In 3D, the same upflow/downflow asymmetry is seen, but other aspects of the flow are distinctly different. Fig.  2d-f show select snapshots of a 3D simulation with the same input parameters as the 2D case in Fig. 2a. In 2D, large-scale, coherent spinners dominate the flow, leading to a single upflow and downflow. New downflowing plumes at the upper boundary are efficiently swept into the large coherent structure (near x ∼ 1 and spanning the vertical domain). The behavior of downflows in 3D is strikingly different (Fig. 2d). In 3D, many individual plumes detach from the upper boundary, but do not organize into a single dominant downflow in the same fashion. Horizontal cuts near the top of the domain (Fig. 2e) reveal a network of narrow downflow lanes surrounding broad upflows. Stronger clusters of downflows near the surface are linked to sheets of low entropy at the midplane of the domain (Fig. 2f). As the flows evolve in time, new downflows appear at the top of the domain in the middle of upflows and join the surrounding downflow network, causing the convective structures to fragment. There is no preferred orientation in the newly forming downflows, and the convective flow field constantly evolves, which appears to prevent the occurrence of either spinners or windy states.
The efficiency of convection is quantified by the Nusselt number (Nu). Nu is well-defined in RBC as the total flux normalized by the steady-state conductive flux [16,17]. In stratified convection Nu is more difficult to define, and we use a modified version of a traditional stratified Nusselt number [1,3], where F conv,z and F cond,z are the z-components of F conv and F cond , and are volume averages. F A ≡ − κ ∂ z T ad is the conductive flux of the proper corresponding adiabatic atmosphere. For a compressible, ideal gas in hydrostatic equilibrium, ∂ z T ad ≡ −g/c P [18]. It is important to measure the evolved value of κ = ρχ , which is nearly κ 0 for small but changes appreciably for large values of . In incompressible Boussinesq convection, where ∇S = 0 only when ∇T = 0, and where κ is constant with depth and time, this definition reduces to the traditionally defined Nusselt number [16,17]. The variation of Nu with Ra is shown in Fig. 3a. We find that the Nu depends primarily on Ra, not on , except where dynamical regimes change. In 2D and at low to moderate Ra, Nu ∝ ∼ Ra 1/3 regardless of , reminiscent of scaling laws in Rayleigh-Bénard boundary layer theory [19][20][21]. As the flow becomes supersonic, Nu ∝ ∼ Ra 1/5 . It is also important to note that, in 2D, the value of Nu is heavily dependent upon the specific thermodynamic structures of the solution, and slight changes in Ra can result in a simulation converging to one solution or another. Select simulations were run at higher aspect ratios (8 and 16), and similar flow morphologies were obtained, suggesting that these states are not highly sensitive to aspect ratio. In our limited 3D runs, it appears that Nu ∝ ∼ Ra 2/7 , a classic scaling law seen in RBC [16].
The rms Reynolds number (Re = |u|L z /ν) and Peclet number (Pe = Pr Re) compare the importance of advection to diffusion in the evolved convective state. For Pr = 1, Pe = Re. Our choice of {ν, χ} ∝ ρ −1 0 drastically changes the value of Re between the top and bottom of the atmosphere. We report values of Re at the midplane (z = L z /2) of the atmosphere in Fig. 3b. We find that Re depends largely on Ra, but not , except when the flow regime changes. In 2D Re ∝ ∼ Ra 3/4 at low Ra. When the 2D flows become supersonic, Re ∝ ∼ Ra 1/2 , as expected from (3). In our limited 3D runs, Re ∝ ∼ Ra 1/2 , consistent with the supersonic results. The heightened scaling of Re in 2D follows the scaling of velocity (Ma) with Ra, as seen in Fig. 1, and reflects the presence of coherent spinners, which do not exist in 3D.

B. Evolved stratification
In the evolved state, the flows can change the density stratification. In Fig. 4a, we measure the time-and horizontally-averaged density profile in two ways. Empty symbols show the number of density scale heights between the maximum and minimum of the atmospheric density profile. Solid symbols show the number of density scale heights between the top and bottom of the atmosphere. We find that near-sonic and supersonic flows support significant, persistent density inversions in the boundary layers, as was reported previously [6]. This is visible when solid symbols lie below empty symbols. We find this in 2D and 3D, even at very large .
Sample evolved density profiles are displayed in Fig. 4b. The natural log of the temporally and horizontally averaged density profile, ln ρ = ln ρ 0 + ln ρ 1 , is displayed for four cases. At low (dotted green line), the density stratification is, to first order, unchanged from the initial density stratification. At high , in both 2D (solid purple line) and 3D (dashed purple line and dash-dot-dot red line), the evolved stratification differs significantly from the initial state and does not increase monotonically with depth. To measure the number of density scale heights between two points in the atmosphere, z 1 & z 2 , we calculate n ρ (z 1 , z 2 ) = ln ρ(z 2 ) − ln ρ(z 1 ). Thus, the values plotted in Fig. 4a for the cases in Fig. 4b can be directly read off. For example, at = 1 and Ra/Ra crit ∼ 10 4 (dash-dot-dot red line), measuring the stratification between the boundaries retrieves n ρ (L z , 0) ≈ −0.3, but measuring between the maximum and minimum value of the profile retrieves n ρ (min, max) ≈ 1. 6.
Surprisingly, the evolved n ρ is always less than the initial n ρ = 3, and turbulent pressure support plays a larger role than atmospheric slumping. This appears to arise as a result of convection making the interior isentropic in the presence of fixed-temperature boundary conditions; we expect the behavior of the stratification to be dependent on the choice of thermal boundary conditions. The agreement of Nu & Re across (Fig. 3), particularly at low Ra in which all four of our cases collapse onto a single power law, is striking in light of the vastly different evolved stratifications felt by the flows.

IV. DISCUSSION & FUTURE WORK
We have found that the evolved flow properties of stratified, compressible convection scale in a manner reminiscent of incompressible, Boussinesq Rayleigh-Bénard convection. We argue that polytropically stratified atmospheres are the natural extension of the RBC problem with an additional control parameter, , whose primary role is to set the Ma of the flows. We show that other properties of the evolved solutions (Nu, Re) are nearly identical at vastly different values of , except for where there is a transition between the subsonic and supersonic regimes. We also see that Nu scales similarly in 2D and 3D, and that Ma in 3D solutions appears to be a function of alone, allowing for simple specification of the evolved Ma using input parameters. The stratification of these polytropic atmospheres evolves in a complex manner, and future work should aim to understand the importance of stratification on convective heat transport and other flow properties.
Our studies here will serve as a foundation for comparing heat transport in stratified convection to that in RBC [16] and for better quantifying transport in stratified convection. These results can be used to determine if fluid properties scale appropriately in simplified equation sets, such as the anelastic equations. This work will also be useful in coming to understand more realistic systems, such as rapidly rotating atmospheres [22], atmospheres bounded by stable regions [23], and regions with realistic profiles of κ.

ACKNOWLEDGMENTS
EHA acknowledges the support of the University of Colorado's George Ellery Hale Graduate Student Fellowship. This work was additionally supported by NASA LWS grant number NNX16AC92G. Computations were conducted with support by the NASA High End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center on Pleiades with allocations GID s1647 and GID g26133. We thank Jon Aurnou, Axel Brandenburg, Keith Julien, Mark Rast, and Jeff Oishi for many useful discussions. We also thank the three anonymous referees whose careful comments greatly improved the quality of this paper.