Viscous Dissipation and Dynamics in Simulations of Rotating, Stratified Plane-layer Convection

Convection in stars and planets must be maintained against viscous and Ohmic dissipation. Here, we present the first systematic investigation of viscous dissipation in simulations of rotating, density-stratified plane layers of convection. Our simulations consider an anelastic ideal gas, and employ the open-source code Dedalus. We demonstrate that when the convection is sufficiently vigorous, the integrated dissipative heating tends towards a value that is independent of viscosity or thermal diffusivity, but depends on the imposed luminosity and the stratification. We show that knowledge of the dissipation provides a bound on the magnitude of the kinetic energy flux in the convection zone. In our non-rotating cases with simple flow fields, much of the dissipation occurs near the highest possible temperatures, and the kinetic energy flux approaches this bound. In the rotating cases, although the total integrated dissipation is similar, it is much more uniformly distributed (and locally balanced by work against the stratification), with a consequently smaller kinetic energy flux. The heat transport in our rotating simulations is in good agreement with results previously obtained for 3D Boussinesq convection, and approaches the predictions of diffusion-free theory.


INTRODUCTION
Convection occurs in the interior of every main-sequence star and in many planets, and must be maintained against a finite amount of viscous and Ohmic dissipation.In a steady state, the dynamics and the dissipation are therefore linked; constraints on one yield constraints on the other.
Many authors have explored this link.For example, the widelyemployed theory of Rayleigh-Bénard convection developed by Grossmann & Lohse (2000) relies on the exact relationship between viscous dissipation and heat transport (Shraiman & Siggia 2000); in their model, the heat transport depends crucially on whether the viscous and thermal dissipation occur primarily in the bulk of the convective domain or in the boundary layers that form at its top and base.Jones et al. (2022) have recently explored an extension of this theory to the density-stratified case, with the spatial distribution of the dissipation again playing a vital role.In the stellar context, Anders et al. (2022) have shown that the magnitude of the dissipation within a convection zone strongly influences the amount of convective overshooting into adjacent stable layers.The form and magnitude of the dissipation is likewise crucial in a variety of efforts to go "beyond mixing length theory" (MLT) (e.g., Kupka et al. 2022, Canuto 1997, Viallet et al. 2013, Meakin & Arnett 2010, Arnett et al. 2015).In the Sun, where the form and magnitude of convective flows in ★ E-mail: s.lance@exeter.ac.uk † E-mail: laura.currie@durham.ac.uk ‡ E-mail: m.k.m.browning@exeter.ac.uk the deep convection zone are currently the subject of much debate (e.g., Vasil et al. 2021), the total dissipation may provide important constraints on the flows (Ginet 1994).Ohmic dissipation in particular is thought to limit the depth of zonal winds in Jupiter (Liu et al. 2008, Kaspi et al. 2018, Kaspi et al. 2020), may constrain magnetism in the interiors of low-mass stars (Browning et al. 2016), and could influence the radii of hot Jupiters (Batygin & Stevenson 2010).The purpose of this paper is to provide new constraints on the magnitude and spatial distribution of the viscous dissipation that may be occurring in stellar convection zones.Few prior works have systematically investigated this in the astrophysically-relevant case of a gas with density and temperature stratification, and none have done so when rotation is also present.Here we study this issue within one of the simplest possible systems that captures convection, rotation, and stratification, by conducting a series of hydrodynamic simulations of stratified (anelastic) convection in a rotating Cartesian domain, situated at a fixed latitude.The vast majority of the simulations presented here are 2D, although we compare some of these results to a very small number of 3D simulations.Many elements that are important in real stars -including, crucially, magnetic fieldsare thus absent here.However, this setup has the great advantage that it allows us to sample parameter regimes that would be difficult or impossible to probe in equivalent detail in a full 3D spherical geometry.
In particular, we are able to assess how the dissipation scales with luminosity, rotation rate, and stratification in the limit where the diffusivities are small (i.e., when the convective supercriticality is high).In what follows, we argue that in this regime the dissipation rate (integrated over the convection zone) depends only on the luminosity and the stratification, and is (at fixed supercriticality) independent of rotation.However, the spatial distribution of this dissipation -and with it, many other aspects of the dynamics -does depend on rotation, as detailed below.
In the remainder of this introduction, we summarise prior bounds on the viscous and Ohmic dissipation, and describe how our work extends these.In Section 2 we detail our simulation setup.In Section 3 we provide a brief, qualitative overview of the dynamics in our simulations.In Section 4 we examine the magnitude and spatial distribution of the dissipation in these simulations, and how these scale with the convective driving, the stratification, and the rotational influence.In Section 5 we explore the links between the dissipation, dynamics, and heat transport.We show there that knowledge of the dissipation provides novel constraints on the kinetic energy flux.We close in Section 6 with a summary of our results and their possible astrophysical implications.

Overview of prior work: bounds and constraints on dissipative heating
In the interior of a star, the microphysical diffusion of momentum, heat, or magnetic fields is typically very small compared to other physical processes, so that the relevant non-dimensional numbers (e.g., the Reynolds, Rayleigh, and magnetic Reynolds numbers) are usually very large (e.g., Kulsrud 2005, Brun & Browning 2017, Jermyn et al. 2022).This need not imply, however, that viscous and Ohmic dissipation are negligible.
To place our discussion on a firmer footing, and to highlight some of the aims of our work, we briefly describe the thermodynamic constraints on the dissipation here.More complete discussions can be found in Hewitt et al. (1975) (hereafter HMW75), in Backus (1975), Alboussière & Ricard (2013, 2014), and Alboussière et al. (2022).
Consider a volume of convecting fluid with an associated magnetic field , enclosed by some surface .Assume this surface is impenetrable and either stress-free or no-slip, so that the normal component of the fluid velocity , and either all components of or the tangential stress vanish on .The local rate of change of total energy can be expressed by with the fluid's internal energy, its density, Ψ the gravitational potential satisfying g = ∇Ψ, the pressure, the contribution to the total stress tensor from irreversible processes, the thermal conductivity, the temperature, the rate of internal heat generation (e.g., by nuclear fusion or radioactive decay) or cooling (e.g., by any processes not included in the conductive term), and 0 the permeability of free space.We have assumed the MHD approximation holds, so that E = −u × B + ∇ × B, where = 1/( 0 ) is the magnetic diffusivity and the electrical conductivity (e.g., Priest 2014).Physically, the rate of total energy change at a point is given by the sum of the net inward flux of energy (the divergence terms in eqn. 1) and the rate of internal heat generation.
The first global constraint is that total energy is conserved, but this yields little insight into the magnitude of the dissipative heating.Integrating (1) over gives assuming both a steady state and that the electric current, j = (∇ × B)/ 0 , vanishes everywhere outside .Equation (2) implies that the net flux out of is equal to the total rate of internal heating and cooling.But dissipative terms do not appear in this equation; viscous and ohmic heating do not contribute to the overall heat flux.
To constrain the dissipation, we turn instead to the internal energy equation, which can be written as: Assuming a steady state and integrating over the fluid volume V it can be shown that where the total dissipative heating rate, Φ, is defined as The first and second terms inside the integral represent the contributions due to viscous and Ohmic effects respectively.Equation ( 4) implies that the total dissipative heating, integrated over the volume, is exactly balanced by the work done against the background stratification (Currie & Browning 2017).Equivalently, from the first law of thermodynamics, we have where is the specific entropy, implying that Following HMW75, we can divide this equation by , and integrate over volume to find (8) Physically, this equation expresses the fact that there is a flux of entropy in and out of the domain (first term), and that entropy can be generated in the bulk by conduction (second term) or by heating within the domain (third term).If the inward flux of entropy at the bottom is less than the outward flux of entropy out the top -as occurs if there is a temperature contrast across the domain-then the difference must be made up by entropy generation in the convection zone (either by conduction or dissipation).
In HMW75 this equation is used to derive an upper limit on the total amount of dissipative heating that can occur in a convective layer of depth .This is given by where top and 0 denote the upper and lower boundary values of the temperature respectively and is the luminosity through the layer.This upper limit corresponds to the case in equation ( 8) where there is negligible entropy generation by conduction or heating, and where the dissipation occurs at the highest possible temperature (i.e., at the bottom of the domain).In this case, as discussed in HMW75, the total dissipative heating rate is bounded not by , but by / , with a suitably-defined temperature scale height.In general, however, it cannot be assumed that dissipation occurs at the highest possible temperatures.For example, the dissipation could be distributed more uniformly throughout the layer, or it could be concentrated predominantly in boundary layers.In these situations, could in principle be much smaller than the upper bound of equation ( 9).
Prior simulations have shown that in certain circumstances convection can approach a version of this bound.HMW75 demonstrated that for the specific case of a Boussinesq liquid without magnetism, the integrated dissipation approached a value of order the bound at high enough Rayleigh numbers (measuring the ratio of buoyancy driving to viscous and thermal dissipation).Jarvis & McKenzie (1980) expanded on this by investigating the case of compressible convection in the infinite Prandtl number (defined as the ratio of viscous to thermal diffusivities) regime, appropriate for convection within the Earth's mantle.Currie & Browning (2017) (hereafter CB17) extended these results to a gas at finite , as appropriate for convection in stellar interiors.In a series of 2D hydrodynamic simulations without rotation, they found that the total dissipative heating in their calculations obeyed a tighter, but purely empirical bound, specifically (defining ≡ Φ/ , with Φ the total viscous heating and the luminosity) where is a modified thermal scale height involving the scale height at the top and bottom boundaries ( ,top and ,0 respectively), and some vertical height * , defined such that half of the fluid mass lies above and below * .They showed that for sufficiently high supercriticalities, the value of appeared to approach equation (10) asymptotically.
Yet not all convective systems actually approach these upper bounds.Recently Alboussière et al. (2022), studying 2D convection with an unusual equation of state in which entropy was a function solely of density, found much lower levels of dissipation than suggested by equation (10) in most cases.They attributed the difference in part to the different boundary conditions adopted in their work; in particular, they showed that for their equation of state, high levels of dissipation (approaching the bound in eqn.10) were only realised in cases with rigid walls (as employed in CB17), and not in those with periodic boundary conditions.
Together, these prior results demonstrate that different values of the total dissipative heating are possible in stratified convection.A central aim of this paper is to provide constraints on how much dissipation actually occurs, for the astrophysically-relevant case of an ideal gas with rotation.

Model Setup
We model a layer of fluid contained between impermeable, free-slip boundaries at = 0 and = and assume that the horizontal boundaries are periodic.Our coordinate system is such that the horizontal coordinates, and , correspond to longitudinal and latitudinal directions respectively, and the vertical coordinate, , corresponds to the radial axis.In the majority of our simulations, we retain all three components of velocity but assume all variables are independent of , so that those simulations are 2D (however, for a few cases in §5.5 we relax this constraint and consider the fully 3D case).Gravity acts in the negative direction.To drive convection, we impose a flux at the bottom boundary and fix entropy at the top.Note that for our 2D simulations, has units of energy per time per length (rather than per area, as in the 3D cases) and is related to luminosity, , by = / where, in 2D cases, is the horizontal box length and in 3D cases is the horizontal cross-sectional area.
To investigate the effects of rotation on such a convective layer we consider a tilted f-plane, where the rotation vector takes the form = (0, Ω cos , Ω sin ) where Ω is the rotation rate, and is the latitude.We conducted cases at = 90 • and = 45 • ; for clarity, in almost all of our discussion below we focus on cases at = 90 • , which corresponds to a vertical rotation vector, aligned with gravity and representative of polar latitudes on a spherical body.
To allow for the effects of a density stratification, we use the anelastic equation set under the Lantz-Braginsky-Roberts (LBR) approximation (Lantz 1992, Braginsky & Roberts 1995).This is expected to be valid when the flows are sufficiently subsonic and the stratification is nearly adiabatic.We diffuse entropy instead of temperature (see discussions in, e.g., Lecoanet et al. 2014).We also consider only the hydrodynamical problem, so there is no Lorentz force and all dissipation is viscous.
The governing equations (in dimensionless form) are where = ( , , ) is the fluid velocity, is the specific entropy, is the pressure, ¯ and ¯ are the reference state density and temperature (defined in (17) below) and These quantities are all dimensionless and were obtained from their dimensional counterparts using as the characteristic length scale, 2 / as the characteristic time scale (where is the kinematic viscosity), and / as the characteristic velocity.Specific entropy has characteristic scale 0 0 , where 0 and 0 are, respectively, the values of the reference state density and temperature at the bottom of the domain and is the thermal diffusivity.The characteristic scales for , density, and temperature are 0 2 / 2 , 0 , and 0 respectively.Luminosity has scale , where has either characteristic scale in 2D or 2 in 3D.
Equations ( 12) -( 14) contain several dimensionless parameters defined as follows: where is the specific heat capacity at constant pressure and is the acceleration due to gravity.In this work we take , , and to be constant (i.e., they do not vary with depth).
is the Prandtl number and is taken to be unity throughout this study.
is the usual Taylor number (quantifying Coriolis forces relative to viscous effects) and is a flux-based Rayleigh number.Alongside it will also be useful to consider the traditional Rayleigh number defined as = 3 Δ , where Δ is the entropy difference across the layer.and will be varied to examine solutions at different rotation rates and at different levels of convective driving.
The reference state is taken to be a time-independent, hydrostatic, ideal gas given by where = 3 2 is the polytropic index.We will refer to the number of density scale heights across our layer, , to quantify the degree of stratification.This is defined as In dimensionless terms the boundary conditions amount to enforcing = −1 at = 0 and = 0 at = 1.The impermeable and stress-free boundary conditions become respectively = 0 and = 0 on = 0, 1 (for the 3D simulations in §5.5 , we we also have = 0).We solve this system using the pseudo-spectral code Dedalus (Burns et al. 2020).Our simulations at low and moderate generally use 256 grid points in the horizontal and 128 in the vertical with an aspect ratio of two (i.e., the box is twice as wide as it is tall); at higher , higher resolutions were required (with up to 640 grid points in the vertical direction), and in a few of the highest-resolution cases (at = 10 11 ) we have considered an aspect ratio of 1.075 instead.The 3D cases in §5.5 have an aspect ratio of two.
The simulations were initialised either by imposing very small entropy perturbations on a motionless base state, or (for some cases at higher ) from an evolved state at lower .It is convenient to specify the buoyancy driving in a convective system by reference to the value of at which convection first occurs, the critical Rayleigh number .To find , we constructed an eigenvalue problem (EVP) solver using Dedalus (Burns et al. 2020), from which we obtain a grid of growth rates for a given input range of and values, where is the horizontal wavenumber.We then used the open-source Eigentools package (Oishi et al. 2021) to find (taking into account only those modes that would fit into the finite computational domain).Rotation and stratification both modify the values of (see, e.g., Chandrasekhar 1967, Mizerski & Tobias 2011).For the parameters studied here, then varies from of order 100 (for, e.g., cases at = 4, = 0) to nearly 10 8 (for cases at = 10 11 ).

OVERVIEW OF RESULTING DYNAMICS
The convective flows in this system are influenced by rotation, by stratification, and by the level of buoyancy driving.We have conducted simulations that sample a wide variety of possible states within this multi-dimensional parameter space.We consider cases ranging from the nearly-Boussinesq limit ( = 0.2, with a density contrast from top to bottom of only 1.22) up to stronger stratifications with = 4 (density contrast of 55).The energy passing through the system is quantified by the flux-based Rayleigh number , as defined above; our simulations sample both laminar flows near convective onset (with close to ) and more turbulent states that have ∼ 10 6 .The rotation rate in our simulations is quantified by the Taylor number defined above, which varies between = 10 and = 10 11 .The Ekman number is also commonly used to quantify the influence of rotation relative to viscosity; = −1/2 , so here varies from 3.16 × 10 −1 to 3.16 × 10 −6 .We conducted simulations at latitudes of 90 • and 45 • , but in almost all the figures below have chosen to focus on cases at 90 • for clarity.(None of the key quantities reported in this paper, or their scalings with and , appeared to depend significantly on the choice of  1 lists the input parameters and key derived quantities for a small number of these simulations; the full table is available online.At each , simulations were performed at a range of logarithmically spaced supercriticalities.For cases performed at a fixed supercriticality (e.g., Figures 3 and 6), they were instead logarithmically spaced in .
Increasing the rotation rate stabilises the system against convection, increasing the value of .Thus for simulations at constant , increasing in isolation would eventually result in a system that no longer convects.In much of our discussions below we therefore choose to compare simulations at varying but constant supercriticality, / .We also quantify rotation using Ro = /( ) (as in, e.g., Hindman et al. 2020; see also Gilman 1978), which assesses the buoyancy driving relative to the Coriolis force (see Anders et al. 2019 for a discussion of how this relates to other measures of rotation).We sample both rapidly-rotating cases (with some having volume-averaged values of < 1) and ones in which rotation has little dynamical role ( ≫ 1).Many different types of flow are possible within this parameter space.Three illustrative examples can be seen in Figure 1, which shows the specific entropy for (top row) a non-rotating case at Ra F = 10 4 Ra with a moderate density stratification ( = 1.4), (middle row) a rotating case ( = 10 8 ) with the same stratification but = 56.2≈ 3.8 × 10 7 , and (bottom row) a rotating case ( = 10 8 ) at = 10 4 ≈ 6.8 × 10 9 .In our non-rotating simulations, for which a typical case is shown in the topmost panel of Figure 1, the convection tends to consist of a small number of convective cells, and to be steady in time.These simple flow patterns persist to surprisingly high values of in the setup investigated here (as also seen in Boussinesq simulations with stress-free boundaries in, e.g., Wang et al. 2020); in some other problem formulations (e.g., with fixed entropy or temperature boundary conditions) the flow tends to become visibly turbulent at lower values of (see examples in, e.g., Anders & Brown 2017;Rogers et al. 2003).When rotation is dynamically significant, as shown in the middle panel, the convective patterns tend to align with the axis of rotation in accordance with the Taylor-Proudman theorem.Rotating cases at the same but even higher (as shown in the bottom panel), in which the Coriolis force is small relative to inertia, exhibit time-dependent flow with structure on many spatial scales.
Most of our simulations behave, qualitatively, like one of the three examples in Figure 1.The non-rotating simulations (as sampled in the top panel) represent one extreme; the rotating, very highcases (as in the bottom panel) are another.The single-celled case is presumably not realised in any actual star, but serves as a useful limit, showing what can occur when a convective plume travels almost unimpeded from the top to the bottom of the domain; in this limit (as we demonstrate below) most of the dissipation occurs in the bottom boundary layer.The cases with rotation are more realistic, exhibiting flow and dissipation throughout the domain.Below, we explore (for several different stratifications) how the dissipation and dynamics vary in between these extremes, as a function of rotational influence.

The maximum value of viscous dissipation at high
Here, we examine whether the high levels of dissipation found in CB17 are realised in rotating cases as well.We find that, for the levels of stratification examined here, the total amount of dissipative heating in the rotating simulations appears to approach a similar upper bound to that realised in non-rotating calculations.The models here were conducted with a different aspect ratio than in CB17 (here the horizontal layer size is twice its depth, whereas in CB17 they were equal) and different boundary conditions (here periodic, impermeable in CB17), so we also indirectly show that these results are, for an ideal gas equation of state, not directly dependent on these factors.
Figure 2 shows for a representative selection of cases at different and , for three different stratifications.For our nondimensional setup, is given by Recall that is now non-dimensional, and so for our 2D simulations it is equal to the aspect ratio of the layer, while for the 3D cases is equal to the aspect ratio squared.The horizontal lines in Figure 2 show the value of equation ( 10) at each value of .At high enough supercriticalities, both the rotating and non-rotating cases appear to approach this limiting value, which is dependent on the layer depth and stratification but independent of (and likewise also independent of viscosity or diffusivity).We have found no cases that exceed this value, but (because it is only an empirical bound) cannot rule out the possibility that it would be exceeded at higher or for other parameter regimes.We have displayed example cases at = 10 8 ; our cases at = 10 11 and latitude 45 • exhibit identical behaviour.Both the rotating and non-rotating cases shown here exceed = 1 (i.e., the total integrated dissipative heating exceeds the imposed luminosity) at high enough .However, we cannot rule out different asymptotic values of for the rotating and non-rotating cases, as discussed in more detail below.
It is clear from Figure 2 that the empirical bound on dissipative heating given by equation ( 10) is approached only for sufficiently high / , and that the value of / needed to reach the upper bound is different for each .The largest cases have not quite reached the asymptotic upper limit as they have not been performed at a high enough supercriticality.
A complementary view but instead focused on the influence of rotation is provided by Figure 3, which shows for a selection of cases at fixed supercriticality (here = 10 2 ) and three different but varying (i.e., with varying rotational influence relative to viscous effects).In the regime probed here, it is clear that the presence of rotation does not greatly alter the volume-integrated magnitude of viscous dissipation despite significant changes in the dynamics.
Note that because rotation stabilises the system against convection, the high-cases shown here have appreciably higher than nonrotating cases at the same supercriticality.We have that Ra ∝ Ta 2 3 (Chandrasekhar 1967), so for example cases at = 10 11 require about 10 7 times higher than non-rotating equivalents to reach the same supercriticality.The convective flow fields in the rotating cases are, at equivalent / , more complex than in the nonrotating cases, but they eventually asymptote to similar levels of viscous dissipation.Further, the cases shown here span a range of convective Rossby numbers, from ≈ 22 in the lowest-cases to ≈ 0.3 at = 10 11 , and so sample both cases in which rotation plays little dynamical role (with ≫ 1) and those in which it is more significant ( < 1).

Entropy generation by dissipation and conduction
In the previous section, we suggested that in both rotating and nonrotating cases, the total viscous dissipation at first increases with increased buoyancy driving (higher ) and then plateaus at or below a fixed value (10) that depends on the layer height and stratification but is independent of the rotation rate or diffusivities.Here we begin to explore how this arises.To do so, we consider entropy generation by conduction and dissipation at varying .In a steady state, the energy entering the convection zone at the bottom boundary (by conduction) must equal the energy leaving at the top boundary (also by conduction).The top boundary is at a lower temperature than the bottom one, so the conductive entropy flux out the top is larger than the entropy flux entering the domain; the difference must be made up by entropy generation within the domain, associated with either conduction or viscous dissipation.For our simulations (employing entropy diffusion and without magnetic fields), this implies that and so which follows from equation ( 14) after integration (using the divergence theorem and the constraint of mass conservation).We have also used For reference, we note that the dimensional equivalent of equation ( 21), retaining entropy diffusion, would be and the dimensional equivalent for temperature diffusion, and including magnetism and internal heating, is equation ( 8).
In Figure 4, we examine the terms in ( 21) for a series of calculations at varying , , and .In all cases the sum of the terms on the right hand side of equation ( 21) ( and ) correctly matches − , the mismatch between the entropy flux at the top and bottom boundaries.As increases, the relative contributions of conduction ( ) and dissipation ( ) change: at low both processes contribute to the entropy balance, whereas at high enough there is negligible entropy generation by conduction within the bulk.
In the non-rotating cases the bulk becomes nearly isentropic at high , so that the conductive entropy generation term is then confined mainly to thin thermal boundary layers whose width (discussed in §5.4) decreases with .Thus scales roughly with the width of these boundary layers; as shown in §5.4, the largest (top) boundary layer width scales as −1/4 in our simulations.We have therefore plotted a corresponding −1/4 dependence in Figure 4 to guide the eye; in the non-rotating cases appears to follow this trend reasonably well.(The line is not a fit; it is chosen to pass through the fourth data point for illustrative purposes.)The behaviour in the rotating cases is more complicated, as discussed below, partly because in these cases the entropy gradient (and hence also ) is nonzero in the bulk.
These trends are linked to the values of explored above.We have overplotted the measured values of at each in Figure 4.The simulations with the highest values are those in which entropy generation by conduction is negligible; further, in non-rotating cases the -dependence of is well-matched by the -dependence of (which in turn is linked to the width of the thermal boundary layers as noted above).
These results also help us understand why simulations must be run at much higher to reach the "dissipative asymptote" when the stratification is strong (i.e., at high ).The purely conductive state has which increases with increasing .However, knowledge of alone is not sufficient to determine the actual value of at each .In the limit of high , when is negligible, we know that , so the highest possible value consistent with the known occurs if all the dissipation is at the highest possible temperature.As noted in §1, the firm upper bound of HMW75 corresponds to this limit.More generally, the value of actually attained depends on both the magnitude of the dissipative entropy generation term and on where it occurs.For example if − / ≡ 0 were constant throughout the domain, we would have = − 2 0 , with = − 2 0 ∫ 1 0 ( / ¯ ), which (upon substituting for ¯ and integrating) gives Equating this to − allows us to solve for 0 in this limit.This in turn allows calculation of for this situation, which reduces to ≈ = 1 − − / , if is small.Our rotating cases at very high , which have > 1 and exhibit intricate flow fields, exhibit close to the value predicted by equation ( 26), though they slightly exceed it at the highest we have probed.They always remain below the limit of equation (10).For example, at = 1.4,equation ( 26) yields ≈ 1.00, whereas equation ( 10) gives ≈ 1.21 and equation ( 9) yields ≈ 1.54; our highest-, = 10 8 simulation at that stratification, shown in Figure 2, has ≈ 1.05.In comparison, many of our non-rotating simulations (which have simpler flow fields) exhibit values of that exceed equation ( 26), the predicted value of for uniform dissipation.For example, at = 2, equation ( 26) would yield ≈ 1.54, whereas our highest-non-rotating simulations at that stratification have ≈ 2.08.This is closer to, but does not exceed, the limit described by equation ( 10), which for the same stratification is 2.15.(We have found no cases that exceed the empirical bound in eqn. 10, which is always tighter than the firm bound of eqn. 9.For example, at = 2, the latter is 2.8, which is significantly larger than in our highestcases.) If the dissipation were uniformly distributed throughout the domain, equation ( 26) would provide a useful bound on .Some of our simulations exceed this bound, so evidently in at least these cases the dissipation is not uniform: a disproportionate amount must occur at higher temperatures, allowing to be higher than suggested by equation ( 26) while still satisfying the entropy constraint that = − .In the following section we explore how and when this occurs.

The spatial distribution of dissipation
Here, we determine where the dissipation occurs in our simulations.We show -in particular by examination of a "dissipation halfheight" (defined to be the height by which half of the dissipation occurs) -that the non-rotating cases at high which approach the CB17 upper bound correspond to situations in which much of the dissipation occurs close to the bottom of the domain and there is negligible entropy generation by conduction in the bulk.In rotating cases, the dissipation is more uniformly distributed throughout the interior.
We begin by defining , where is the horizontal average of − 2 / (the local dissipative heating).Here represents the total dissipative heating up to height , so if the heating were uniformly distributed throughout the interior (with a constant) / would be a linear function of height, increasing from zero at the lower boundary to − at the  top.Qualitatively, we find that ( ) is close to linear for our rotating cases; the dissipative heating is nearly uniform throughout the domain.In our non-rotating cases, by contrast, much of the dissipation occurs near the lower boundary.Examples, and a discussion of how these are linked to the buoyancy work and to the dynamics, can be found in §5.1 below.
A simple, quantitative assessment of the sites where dissipation occurs is provided by Figure 5, which shows what we call the "dissipation half-height" in a range of cases.We define this as the location at which reaches half its maximum (absolute) value.That is, If the convection and dissipation were uniform throughout the domain, would be 0.5; meanwhile if the dissipation occurs predominantly at the lower boundary, will tend towards the width of the lower dynamical boundary layer.In both the rotating and non-rotating cases, declines at first with increasing supercriticality: more and more of the dissipation occurs close to the bottom boundary.In the non-rotating cases it appears to level out (i.e., is approximately constant) at high enough .In the non-rotating cases studied here -that is, in 2D and at = 1 specifically -at high enough the flow consists of a steady cell of convection.Thus there is dissipation around the single downflow plume, and in the top and bottom dynamical boundary layers, but very little elsewhere in the bulk.In this limit, is related to the point at which the flow is bent from the vertical towards the horizontal, and this occurs progressively nearer the lower boundary at moderate .In the rotating cases, the dissipation is more uniformly distributed throughout the domain, so is (at any given ) higher than in the non-rotating cases.However, still declines as increases.Note that at fixed , as sampled here, increasing / implies decreasing rotational influence on the dynamics, so that the cases at high and lower have a higher Rossby number.A complementary view is provided by Figure 6, which examines the value of in a series of cases at the same supercriticality = 10 2 but varying .Here we find that increases with increasing (i.e., with increasing rotation rate).That is, stronger rotation leads to more of the dissipation occurring in the bulk of the domain, far from the lower boundary.There is a fairly sharp transition between a "low-" state at low (high Rossby number) to a "high-" state for higher ; beyond this, increases slowly with increasing (i.e., with increasing rotational influence).This transition is connected to the transition from singlecell states (as achieved in non-rotating cases or at very low ) to much more intricate, time-dependent flows realised at higher and .These changes in the spatial distribution of the dissipation must be linked to changes in the flow field.Since the local dissipative heating is related to the stress tensor / , changes in either the magnitude of , or in the characteristic lengthscales present in the flow, will affect .Both these quantities are expected to depend on rotation rate (e.g., Aurnou et al. 2020;Currie et al. 2020;Nicoski et al. 2023), so it is natural that the dissipation exhibits some dependence on this as well.However, must obey the bounds described in §4.2 at all rotation rates.
Flows in 3D, or in real stars, are bound to exhibit more complexity at all rotation rates, and so we do not expect the numerical values of (for example) to be the same in such cases.Equivalently, the relative proportion of bulk and boundary dissipation might well be different.(In the context of stellar convection, with no impermeable boundaries, the equivalent of "boundary" dissipation might be convective plumes that are dissipated only when they reach adjacent stably-stratified layers.)The extremely low-states seen here at low are also unlikely to be realised in any real star, since they occur only for single-cell flows with little bulk dissipation.However, we expect that both the thermodynamic bounds discussed here, and the general trend towards increasing dissipation in the bulk at higher rotation rates, may be robust.

LINKS BETWEEN DYNAMICS, HEAT TRANSPORT, AND DISSIPATION
In a steady state, dissipation and dynamics are linked, so insight into either one yields constraints on the other.Here, we briefly explore how systematic variations in the governing parameters of this problem (namely , , and ) lead to changes in the energy transport and in the flow fields, and we explore how these are related to changes in the magnitude and spatial distribution of the dissipation.Our discussion here is also intended to help place our work in context with a large body of previous research on heat transport in both non-rotating and rotating convection.
In stellar astrophysics, the main purpose of a convective theory is to provide estimates of the entropy gradient needed to carry a certain luminosity outwards (e.g., Gough & Weiss 1976).For example, in standard stellar evolution theory, the radius of a star depends on its specific entropy, and how this varies with depth (see, e.g., discussions in Ireland & Browning 2018).There is also substantial astrophysical interest in properties of the flow itself -e.g., its magnitude at each depth -since these in turn will affect mixing, the transport of heat and angular momentum, and the generation of magnetic fields.Hence, we focus our discussions here on the heat transport, on the related question of how entropy varies with height in our simulations, and on the magnitude of the flows themselves.In Section 5.2, we demonstrate that a quantity of particular interest, the kinetic energy flux, can be estimated given knowledge of the dissipative heating.

Energy balances and transport terms
In this section we begin to quantify the links between energy transport in our simulations, and where the dissipation occurs.One view of this is provided by Figure 7, which assesses the energy transport in two example calculations.The top panels consider a nonrotating case at = 4, = 10 5 ; the bottom ones are for a rotating case with = 1.4,= 17.8 , = 10 8 .We assess the transport in two complementary ways.The right panels show the terms arising in the following equation, which arises after integration and manipulation of the total energy equation (eqn. 3): where the volume integrals are over the volume enclosed between ′ = 0 and ′ = and the surface integrals are over the bounding surfaces of that volume.The convective heat flux ( ) is defined by the first term; the conductive heat flux ( ) by the second; the third and fourth terms define heating and cooling terms ( and ) arising from the viscous dissipation and from work done against the background stratification, respectively (as also discussed in §4.3).As noted in CB17, these latter two terms must balance when integrated over the entire layer, but they do not have to balance at each depth.
The left panels instead arise from considering the total energy balance (i.e., including kinetic energy as well as internal), which in a steady state may be written as (see e.g., Viallet et al. 2013 defining the enthalpy flux ( ), the kinetic energy flux ( ) and the viscous flux ( ).It is common for global-scale simulations of stellar convection to decompose the transport in this way (e.g., Browning et al. 2004, Featherstone & Hindman 2016).In the notation here, and in Figure 7, positive fluxes are defined to be vertically upwards.
Whether considering the total or internal energy equations, in a steady state the sum of the transport terms must equal , the total luminosity, which is constant throughout the layer.The sum of the transport terms is indicated in Figure 7 by a solid line, which is indeed constant with depth in all the sampled cases.In general, we use ( ) as a gauge of whether a given simulation has been evolved for a long enough time, and averaged over long enough intervals, for the results to be time-independent.We evolved all cases in this paper long enough for ( ) to vary by less than one percent across the layer, and for other aspects of the dynamics (e.g., the kinetic energy evolution) to equilibrate as well.This means that the simulations were evolved for typically tens of viscous diffusion times, and averaged over intervals ranging from 0.1 to several diffusion times.
The energy transport differs substantially in our non-rotating and rotating cases.The top row in Figure 7 is an example of what can occur in non-rotating, stratified cases: here (left panel) the enthalpy flux exceeds the total flux in magnitude; this excess is compensated largely by a negative (inward-directed) kinetic energy flux.Broadly similar transport has been observed in simulations of stratified convection for decades (see, e.g.Hurlburt et al. 1984 in 2D;Stein & Nordlund 1989;Featherstone & Hindman 2016).Transport by conduction is small throughout the layer, outside of narrow boundary layers.
By contrast, in the example rotating case (bottom row) the kinetic energy flux is negligible; the enthalpy flux is approximately equal to the total luminosity, with conductive transport small outside the boundary layers.At this particular , the conductive boundary layers are still relatively large, and there is evident asymmetry between the top boundary layer and the bottom one (which have different widths).
The connection between this transport and the viscous dissipation is made clearer by comparison to the right panels of Figure 7, which considers the internal energy decomposition for the same cases.In all cases the buoyancy work is fairly evenly distributed throughout the domain -that is, outside of the bottom boundary layer rises linearly towards the top domain.The rotating case has fairly uniform dissipative heating, with also nearly linear.But in the non-rotating case, the dissipative heating is non-uniform: more of the dissipation is occurring near the bottom boundary.At the upper boundary, both / and / must approach +/− , respectively, and they do so in both cases; what differs in the rotating and non-rotating simulations is the spatial distribution of .The convective luminosity , as defined and plotted here, is often larger than the total luminosity in the non-rotating cases, in accord with the fact that is greater in magnitude than throughout much of the convection zone.In the rotating cases, where there is approximate local balance (as well as an exact global balance) between the dissipative heating and the "cooling" by buoyancy work, the convective luminosity is closer to unity.

Predicting the kinetic energy flux from the viscous dissipation
The transport revealed here differs in some important ways from that envisioned in MLT, and some of these differences are connected to where the viscous dissipation occurs.In this section we consider the kinetic energy flux, which is not explicitly included in classical MLT (e.g., Gough & Weiss 1976) but is a robust feature of stratified convection in stellar environments.Generally we find that the enthalpy flux exceeds the total luminosity by a considerable amount, and is compensated for by the inward-directed KE flux.However, prior work has not clearly established what sets the amplitude of these fluxes.Might it be possible, for example, for a star like the Sun to have a thousand solar luminosities moving outwards in the enthalpy flux, and 999 moving inwards via the KE flux?In this section we show that knowledge of the viscous dissipation can answer this question, and more generally provide constraints on the magnitude of the kinetic energy flux.Following CB17, we define , where = 1 ∫ .Hence ℎ is equivalent to the steady-state transport associated with processes other than the convective flux as defined above.Outside of the boundary layers, prior work has found that is generally larger in magnitude than or (e.g., Viallet et al. 2013) , so that in the bulk ℎ ≈ .If the local dissipation and buoyancy work terms do balance at each depth, then ℎ is zero.This in turn implies negligible kinetic energy flux.This is approximately the state attained in some of our rotating cases: there, both and are linear in , and of similar magnitude, so that ℎ ≪ .In our non-rotating 2D cases, by contrast, the concentration of much of the dissipative heating near the bottom boundary implies a substantial mismatch between dissipative heating and buoyancy work throughout much of the bulk, so ℎ is no longer negligible.This in turn requires a substantial kinetic energy flux.
We can use these ideas to place simple bounds on the magnitude of the kinetic energy flux.Consider the extreme case in which none of the viscous dissipation occurs in the bulk (i.e., it is all in a comparatively narrow bottom boundary layer).At some depth just above this boundary layer, nearly all the integrated viscous heating will have occurred, but very little of the integrated work will have; in the notation employed here, / will be close to zero, while / will be nearly equal to its value at the top of the domain.The latter is bounded by = / ˆ , as discussed above, so we have ℎ / ≈ just above the bottom boundary.Hence, if the "other" transport is dominated by the KE flux (rather than or ) we expect the maximum absolute value of / to be bounded by the value of at each stratification.
We examine this prediction in Figure 8, which shows the maximum absolute value of the kinetic energy flux in a series of calculations at varying at high supercriticality ( ≥ 10 4 ), along with the limiting value of given by equation ( 10) and the actual value of attained in the simulation.For cases at of two or less, the measured values adhere closely to the limiting value (indicated by the black line); at the highest , the measured values are lower than the theory.The KE flux closely tracks the measured value of at each stratification (lying slightly below it), in keeping with the simple model described above.
These results suggest that at high enough the maximum possible amplitude of the kinetic energy flux may be estimated simply by calculating = / ˆ (eqn.10).Our non-rotating cases, which consist of simple uni-cellular flows, actually approach this limit.However (as noted above) in rotating cases the dissipative heating more nearly balances the buoyancy work at each height, leading to significantly smaller kinetic energy fluxes.Likewise, more complex flows (as likely realised at higher in real stars) likely lead to smaller KE fluxes as well.We speculate, though, that the limits on the KE flux developed here are unlikely to be exceeded by real convective flows.For this to occur, the dissipation and buoyancy work would have to be even more imbalanced than in our single-plume non-rotating cases; for uniform buoyancy work, this would require the dissipation to be concentrated to an even smaller part of the domain than in these simulations.

Entropy profiles and Nusselt number scalings
In previous sections, we saw that the energy transport in our simulations -and in particular the relative contributions of conduction and convection -varied in response to changes in the key controlling parameters , , and .Here we explore how these variations arise, and in particular how they are linked to the entropy gradients established by the convection.
In our work entropy is fixed at the upper boundary; at the lower boundary its gradient is fixed.Recall also that we have assumed that conductive transport is proportional to entropy gradients (rather than temperature gradients).Together, these imply that in the absence of convection, we would expect a linear specific entropy profile, extending from = 0 at the top to Δ at the bottom.For the models considered here, Δ is If convection occurs, a smaller total entropy contrast between the top and bottom boundaries, Δ , is required to carry the same imposed .In our simulations the total Δ , and its variation with height, are functions of , , and .We assess these for two illustrative cases in Figure 9, which shows ( ) (where • denotes a horizontal average) for both a rotating case (at = 10 8 , = 17.8 ) and a rotating one (with = 10 4 ) at = 1.4.In both cases, conduction must carry all the energy within some distance of the boundaries, so there are steep entropy gradients at the top and bottom of the domain; the entropy gradient is smaller in the bulk.The total Δ across the layer is similar in the two cases, but its spatial variation is different: in the non-rotating the interior is nearly isentropic ( is close to a vertical line), whereas in the rotating case there is a visible, finite slope to throughout the bulk.To analyse these trends quantitatively, we turn to an average measure of the heat transport over the domain (rather than to its spatial variation).In studies of Rayleigh-Bénard convection, it is customary to encapsulate this via the Nusselt number , a dimensionless measure of the heat transport relative to that provided by conduction.There is no universally-accepted definition of that makes sense for all boundary conditions, stratifications, and with/without rotation, but sensible definitions have the property that they are large when convection is efficient, and tend to one as the convection vanishes (i.e., as all transport becomes conductive).For the mixed fixed-flux, fixed-entropy boundary conditions here, we choose to adopt as our definition of .This is a global measure of the efficiency of the convective flow -more efficient convection should have a smaller Δ across the layer, and so a higher -but normalised to the Δ that would be required to carry the flux in the absence of convection.This is akin to the definition for Boussinesq convection adopted by, e.g., Kazemi et al. (2022).
The resulting measures of are plotted for a sample of cases with varying , , and in Figure 10.We have also overplotted several previously-proposed scaling relations, as discussed below.We have chosen here not to normalise each case by , primarily because varies so much across the simulations sampled here; in general, each "track" of simulations shown begins with of order ten times critical at that and .First, consider the non-rotating, weakly-stratified cases at = 0.5.These are well-matched by the power law ∝ 1/4 , which is obtained if transport within the bulk is entirely by convection, transport within narrow thermal boundary layers is by conduction, and the width of the boundary layers is set by the requirement that they be marginally stable against convection (Malkus 1954).The scaling at higher appears to be slightly less steep than this.For comparison, we have overplotted ∝ 2/9 .This scaling would arise if ), as has often been reported in non-rotating experiments and simulations (see, e.g., discussions in Grossmann & Lohse 2000;Siggia 1994).None of our data are consistent with the so-called "ultimate regime" scaling ∝ ( ) 1/2 (e.g.Chavanne et al. 1997), which has been conjectured to apply at very high .The rotating cases exhibit steeper ( ) scalings.Figure 10 shows a series of cases at fixed = 10 8 and a smaller number of cases at = 10 11 .Note that the cases at = 10 11 are situated at 45 • rather than 90 • .At moderate (where is small) the dependence of on appear to be reasonably well described by the overplotted scaling ∝ 3/5 .At higher , when rotation is unimportant dynamically (large ), the = 10 8 cases latch on to the non-rotating scalings.The slope of the ( ) relation is similar for the cases at = 10 8 and = 10 11 .Note that for fixed , increasing is equivalent to increasing / , so a plot of as a function of / would exhibit the same slope.This behaviour is consistent with prior results in different settings.For example, the transition from a steep "rotating" ( ) relation to a shallower "non-rotating" one was reported by (King et al. 2009) to occur in plane-layer experiments and accompanying Boussinesq simulations (with fixed temperature and no-slip boundaries); see also discussions in King et al. (2013) and Aurnou et al. (2020).Simulations in spherical shells (see discussions in Gastine et al. 2016, Long et al. 2020b, discussing fixed-temperature and fixed-flux simulations respectively) likewise exhibit similar trends.
Interestingly, the ∝ 3/5 scaling seen in our rotating cases is in accord with the expectations of rotating mixing-length theory (Currie et al. 2020, Barker et al. 2014, Stevenson 1979).The same scaling law also arises in the classical "CIA balance," which supposes a dynamical balance between Coriolis, inertial, and buoyancy terms in the momentum equation (Aurnou et al. 2020, Vasil et al. 2021) and in asymptotic theories of convection at low Rossby number (Julien et al. 2012).In these theories, the transport is predicted to follow ∝ ( / ) 3/5 ∝ ( −2/3 ) 3/5 , which (upon substituting in the definitions for , , and ) is diffusion-free.This suggests that in the rapidly-rotating limit the diffusive boundary layers are playing a less significant role in the heat transport.In dimensional terms, this scaling implies that the entropy gradient in the bulk of the convection zone becomes steeper when rotation is more rapid, scaling as / ∝ Ω 4/5 , with Ω the angular velocity.

Boundary layers and the link to dissipation
The trends explored above arise partly from the varying influence of viscous and thermal boundary layers in our simulations.In this subsection, we explore how the widths and other properties of these boundary layers vary as the supercriticality of the convection, the density stratification, and the rotational influence are changed.We also discuss the manner in which the boundary layers, heat transport, flow amplitudes, and dissipation are linked -and demonstrate explicitly that knowledge of some of these aspects constrains the others.
Many different definitions of the boundary layers have been employed in the literature on Boussinesq convection, but our inclusion of rotation, our use of a fixed-flux thermal boundary condition at the bottom boundary, and our adoption of stress-free velocity boundary conditions together implies that some of these definitions are not relevant (see discussion in Long et al. 2020a).We choose here to adopt the simple method suggested by Long et al. (2020a), defining the width of these layers (near the top and bottom of the domain) to be the points at which the advective and conductive contributions to the heat transport are equal.Inside the boundary layer conduction dominates the heat transport; outside it convection does.Long et al. (2020a) demonstrate that this method gives sensible results in a variety of Boussinesq settings (with and without rotation), though to our knowledge it has not been previously employed to study anelastic convection simulations.
At the top and bottom boundaries conduction must carry all the energy (because in our simulations the vertical convective velocity goes to zero there, and because near-surface radiative cooling has been ignored).The value of the horizontally averaged entropy gradient / at the bottom boundary is therefore determined by the energy flux entering the domain.At the top boundary, the entropy is fixed (rather than its gradient), but in a steady state the simulation must still develop a sufficiently large entropy gradient to carry the same energy flux out the top boundary.Specifically, because we have assumed that conduction diffuses entropy, we must have at both the top and bottom boundaries.(As elsewhere, all variables here are dimensionless; the dimensional version would have a factor of on the right hand side of the equation.)Because we are considering stratified convection, the density at the top of the domain is smaller than at the bottom, so we expect the entropy gradients / that develop will be greater at the outer boundary than at the inner one.
We now suppose that within the conductive boundary layers / is approximately uniform, and equal to where Δ is the entropy jump across the boundary layer and is its width, ¯ and ¯ are evaluated at the top or bottom of the domain (for the top and bottom boundary layers respectively), and where we have assumed conduction carries all the flux within the boundary layer.
In Figure 11 (top) we compare the resulting predictions for Δ / to measurements in example simulations.We show the ratio of Δ / in the bottom boundary layer to that in the top; the horizontal lines denote for each stratification, which this ratio should approach (per eqn.33, and given that is the same at both boundaries).The agreement between the measured values and the estimated value is good at high for all stratifications shown, in both rotating and non-rotating cases.At low the agreement is less good.With increasing the boundary layers become increasingly asymmetric, so that for example in cases with = 2, Δ / is more than 25 times larger in the top boundary layer than in the bottom.This, again, is a consequence of the much smaller densities and temperatures at the top of these stratified domains, which then require a much larger entropy gradient to carry the imposed flux out the top boundary.
In the bottom panel of Figure 11 we show how the top and bottom boundary layer widths vary with for example cases at = 1.4.As expected, the boundary layers grow thinner at higher .The top and bottom boundary layers appear to follow slightly different trends: we have overplotted ∝ −1/4 , which appears to match the top boundary layer data well, and −1/3 , which matches the bottom boundary layer well.These findings are consistent with, and aid in understanding, our findings for the dynamics and heat transport (i.e., ( ) scalings) in previous sections.In the non-rotating cases at high nearly the entire Δ across the whole domain occurs in the top and bottom boundary layers; hence their width determines the overall Nusselt number for the entire domain.(In rotating cases, the entropy contrast across the bulk can approach or exceed that in the boundary layers; see discussions in, e.g., Barker et al. 2014, Currie et al. 2020.)The top boundary layer is, in our stratified calculations, probably the more restrictive of these because it is thicker; we expect thus expect that in these calculations the Nusselt number should scale approximately as the height of the layer divided by the width of this boundary layer.Here, this implies that should scale as 1/4 , which is in agreement with many of our findings for the non-rotating cases in Section 5.3 above.The steeper heat-transport scalings exhibited by rotating cases are linked to where dissipation occurs.This is because the dissipation and the work done against the background stratification must balance; this balance gives rise to the exact relationship in Boussinesq convec-tion between the ( ) relationship and the viscous dissipation (Shraiman & Siggia 2000), and to a more complex analogue of this in the anelastic case (as shown recently by Jones et al. 2022).Changes in where the dissipation occurs -which in turn arise because the convective velocities and lengthscales change in the presence of rotation, as discussed above -thus also give rise to changes in the heat transport.
This link between dissipation and transport has led some prior authors to separate the viscous dissipation into "bulk" and "boundary" contributions -see, e.g., Grossmann & Lohse (2000), Jones et al. (2022), Gastine et al. (2016) -and to posit that transitions in the heat transport correspond to changes between bulk-dominated or boundary-dominated dissipation (Grossmann & Lohse 2000).In our simulations, we cannot usefully divide the dissipation in this way; because the boundary layers as defined here become very thin at high / , the dissipation is almost always "bulk-dominated."We have argued above that provides, for our setup, a more meaningful distinction between cases where the dissipation is concentrated near the boundaries and those where it is distributed throughout the domain.We find that cases that fall on the rotating scaling relation ∝ ( / ) 3/5 systematically have larger than those which follow the non-rotating heat transport scalings ( ∝ 1/4 ).

Comparison to fully three-dimensional flows
The simulations presented in the preceding figures were restricted to two spatial dimensions (i.e., assuming axisymmetry in one dimension).In this section we provide a preliminary view of whether the key quantities assessed in this paper -namely the overall levels of dissipation (as measured by ) and its spatial distribution (as encapsulated by and by the spatial distribution of ) -are different in 2D and 3D cases.
In Figure 12, we examine the energy transport in two example 3D cases and their closest 2D equivalents.We display both non-rotating cases (bottom panel, at ≈ 2.3 × 10 4 ≈ 100 ) and rotating ones (top panel, at = 10 6 , = 3.3 × 10 5 ≈ 10 , at 45 degrees latitude).The 2D/2.5D cases have the same as the 3D ones.All cases have = 1.4,and an aspect ratio of 2:1 (i.e., the horizontal dimensions extend twice as far as the vertical one).The 3D cases were run at resolutions of 256 × 256 × 128, and evolved for more than a diffusion time.In the figure, data from the 2D cases are over-plotted as symbols, with data from the 3D cases as lines.
It is clear that the energy transport in the 3D cases is very similar to that realised in the 2D/2.5Dcases.In the non-rotating case, for example, has a steep slope near the lower boundary, and a smaller one in the bulk; this reflects the fact that much of the dissipation is occurring near the lower boundary.In the rotating case, is less sharply peaked, reflecting the more even distribution of dissipation throughout the bulk.In both cases, the other transport terms ( , , , and the total transport ) are also quite similar in 2D and 3D.
The other bulk quantities of interest to us -and -are also similar in 2D and 3D.The rotating 3D case shown here has (time-averaged) = 0.69 and = 0.61; the corresponding 2.5D case had = 0.66 and = 0.60.The non-rotating 3D case has = 0.80 and = 0.21; the corresponding 2D case had = 0.82 and = 0.17.Crucially, this suggests that in 3D cases exhibits a similar dependence on rotation as it did in 2.5D: namely, it is significantly larger when rotation is present (because more of the dissipation is distributed in the bulk in that case) than in the non- rotating case (where much of the dissipation is concentrated near the boundaries).
These results suggest that the key quantities examined in this paper may not be too sensitive to the assumption of axisymmetry.We defer a more detailed study of the 3D cases to later work.

DISCUSSION AND CONCLUSIONS
We have presented the first systematic investigation of viscous dissipation in a rotating, stratified plane layer of convection, studied here within the anelastic approximation for an ideal gas.We have shown that, for fixed convective supercriticality / and moderate stratification, the total dissipative heating does not depend appreciably on rotation rate.However, the spatial distribution of the dissipation does vary with rotation, and this has a number of important consequences for the dynamics and heat transport.
The total dissipation is thermodynamically bounded by eqn.( 9), which corresponds to the case in which there is negligible entropy generation by conduction in the bulk of the domain and all the dissipation occurs at the highest possible temperature.In practice we have not found any cases in which the dissipation exceeded the tighter, empirical bound of eqn.(10), which does not depend directly on the diffusivities or on rotation rate.Our non-rotating cases, which exhibit simple mono-cellular flows, approach the latter bound at high enough / .These represent an extreme case in which very little dissipation occurs in the bulk of the convection zone, so we regard them as a limit on that is unlikely to be exceeded by more realistic flows.
In rotating cases the viscous dissipation is more uniformly distributed throughout the layer than in corresponding non-rotating cases.In the non-rotating simulations, much of the dissipation occurs near the bottom of the computational domain, so that although there is a global balance between dissipation and work done against the background stratification, these quantities do not balance at each depth.We defined a new quantity , the height at which half the total viscous dissipation has occurred, which encapsulates the spatial distribution of the dissipation in a simple way, and used it to characterise our simulations in different regimes.
We have shown that the heat transport scalings ( ( )) in our rotating cases appear to be consistent with theoretical diffusionfree predictions arising from either "rotating mixing length theory" (Stevenson 1979, Barker et al. 2014, Currie et al. 2020) or, equivalently, from a conjectured balance between Coriolis, inertial, and buoyancy forces (e.g., Aurnou et al. 2020, Gastine et al. 2016, Vasil et al. 2021).Prior work has shown this in other settings (mainly within the Boussinesq approximation).
We have shown that these changes in heat transport are linked to where in the domain the dissipation occurs.This is similar to the case in Boussinesq convection, where prior work (e.g., Grossmann & Lohse 2000) has established that the ( ) heat transport relation varies depending on whether the dissipation is "bulk" or "boundary" dominated, and broadly in line with very recent theoretical predictions for the anelastic case (Jones et al. 2022).
Here, the situation is more complex than in the Boussinesq case because of the background stratification, but the same basic trends appear to hold.In particular, we find that cases which follow the "rotating" heat transport relation are those for which is especially high.
We also established that, for the setup examined here, the thermal boundary layers in our simulations are asymmetric -the top one being considerably larger than the bottom -and that the thicknesses of the top and bottom boundary layers scale differently with .Finally, we have explored the link between dissipation and the kinetic energy flux.We developed a simple model of the kinetic energy flux in our non-rotating cases -based on the idea that dissipation approaches the upper bound at high enough , and that much of the dissipation occurs near the lower boundary -and showed that it provided a reasonably accurate prediction of the maximum (negative) kinetic energy flux attained in our simulations for each stratification at high enough .We have argued that this provides an upper bound on the kinetic energy flux achievable by real convection.
If our results are applicable to real stars, one conclusion is that rapidly-rotating stars should exhibit less convective overshooting into adjacent stably-stratified regions than slowly-rotating ones.From the perspective adopted here, this is because in rotating convection zones the buoyancy work at each depth is more nearly balanced locally (rather than just globally) by dissipation.Equivalently, the more rapidly rotating cases have more dissipation per unit volume in the bulk (all else being equal).For a fixed stratification, this will imply (following the discussion in, e.g., Anders et al. 2022) that convective motions in rotating stars will have less kinetic energy when they reach the boundary of the Schwarzschild-unstable region, so they will penetrate less deeply into the adjacent layers.
Similarly, we expect a smaller kinetic energy flux in rotating stars and planets than in non-rotating ones.The dynamical consequences of this are not yet clear, but we intend to explore them in the future.

Figure 2 .
Figure 2. Calculated values of at a range of stratifications and supercriticalities for sample non-rotating (black) and rotating cases (blue).Values of = 0.5 (circles), 1.4 (stars, and blue squares), and 2.0 (pluses) are used.The horizontal black lines represent the value of equation (10) for each value of .The rotating cases are at latitude 90 • with = 10 8 .Note that due to the effect of rotation on the critical Rayleigh number, the rotating cases have considerably larger values of for a given stratification.

Figure 3 .
Figure 3.The values of E for a range of rotation rates.All cases have a fixed supercriticality of / = 10 2 ; we consider stratifications of = 0.5, 1, and 1.4.The labels show the range of the convective Rossby number for increasing values of Ta at each stratification.All simulations shown are at latitude 90 • , except for the point at = 10 11 which is at 45 • .

Figure 4 .
Figure 4. Entropy generation terms as given in equation (21) for example non-rotating (top panel) and rotating (bottom panel) simulations at = 1.4 and varying.Shown for each case are the entropy production due to conduction () and viscous dissipation ( ), together with their sum ().Also shown (as horizontal green line) is the entropy mismatch − for each case; in all cases the sum of and matches − .For comparison, we also plot the value of for each case, and indicative −1/4 scaling over the values; see text for discussion of this scaling.

Figure 5 .
Figure 5. Values of for a range of non-rotating and rotating simulations at a range of stratifications, .All rotating cases shown are at latitude 90 • .

Figure 6 .
Figure 6.Values of for a range of simulations at fixed supercriticality, = 10 2 for stratifications of = 0.5, 1.0, and 1.4 and a range of Ta.All cases shown are at latitude 90 • , except for the point at = 10 11 which is at 45 • .

Figure 7 .
Figure 7. Fluxes of energy provided by different transport terms in the total energy equation (left panels) and internal energy equation (right panels), for example 2D cases at (top) = 4 without rotation and (bottom) = 1.4 with rotation ( = 10 8 ).

Figure 8 .
Figure 8.Comparison between KE flux amplitudes in theory and simulation.Shown as a solid line is the estimate of eqn.(10) for in non-rotating convection at very high , which we argue in the text provides a bound on the maximum KE flux amplitude.Also shown are the values of in sample non-rotating simulations (at high but finite supercriticality) and the maximum absolute values of the KE flux in the same simulations.

Figure 9 .
Figure 9. Horizontally-averaged specific entropy as a function of depth in two example cases, both at = 1.4.The solid (orange) line is for a rotating case at = 10 8 and the dashed (blue) line is for a non-rotating case ( = 0).

Figure 10 .
Figure 10.Dimensionless heat transport (Nusselt number) as a function of for a sample of rotating and non-rotating cases at different .

Figure 11 .
Figure 11.Properties of top and bottom boundary layers at varying , , and .Top panel shows the ratios of the dimensionless entropy gradients across the bottom and top boundary layers (as defined in the text); bottom panel shows the width of the top and bottom boundary layers in example nonrotating cases at = 1.4.The horizontal lines show the result of evaluating equation (34) at each stratification.

Figure 12 .
Figure 12.Internal energy transport in rotating (top) and non-rotating (bottom) cases in both 2D and 3D.The rotating cases are at latitude 45 • .

Table 1 .
Input parameters and selected output quantities for example simulations presented in this paper.Indicated are the number of density scale heights across the layer ( ), the supercriticality of the simulation / , the critical Rayleigh number , the Taylor number , the convective Rossby number , the latitude (for rotating cases only), and the output quantities , and.Full machine-readable table available online.