A scaling theory for the macroturbulence of weakly supercritical planetary atmospheres

Generated using the official AMS LaTeX template v6.1

Corresponding author: Ryan Eusebi, [email protected]

SIGNIFICANCE STATEMENT: The pole-to-equator temperature gradient is a defining characteristic of a planet’s climate. Its strength, which regulates mid-latitude weather variability, is determined by a balance between radiative heating and poleward turbulent heat transport by large-scale weather systems. In this work, using hundreds of simulations of planetary atmospheres across a wide spectrum of climate states, we show that classical theories cannot explain the heat transport by these weather systems and derive new scalings that capture the behavior. The scalings are implemented as a turbulence closure in a simple one-dimensional energy balance model to explain the pole-to-equator temperature gradients observed across the simulations. The results offer a new framework for understanding the large-scale turbulence of the atmosphere.

Introduction

A complete understanding of the time- and zonal-mean state of a planetary atmosphere, given forcing conditions, relies on understanding the emergent macroturbulence of the system. Such macroturbulence consists of the synoptic scale weather systems, pervasive on Earth, that are responsible for the stormy, variable conditions characteristic of the mid-latitudes. These weather systems, or eddies, are a form of large-scale turbulence responsible for the bulk of atmospheric energy transport poleward of the tropics. Macroturbulence is primarily a product of baroclinic instability, through which eddies convert potential energy from horizontal temperature gradients into kinetic energy [38]. Heat fluxes associated with these eddies control the meridional temperature gradient of the mid-latitudes; eddy momentum fluxes affect the strength and location of the jet streams, surface winds, and the strength of the Hadley circulation [31, 64, 49]. A long-standing goal of the atmospheric community is a description of eddy fluxes based on mean state conditions. While turbulence is a famously difficult problem to understand, the scale separation between the domain size and the eddies in rapidly rotating planets, the relatively weak nonlinearity of the turbulence [53], and the fact that the central goal is understanding spatially inhomogeneous and hence non-trivial mean turbulent fluxes (second moments) rather than an entire hierarchy of higher-order moments, makes understanding their macroturbulence a more tractable problem [21, 62].

The near-surface eddy heat flux, which is the focus of this paper, is a particularly important aspect of atmospheric macroturbulence. Its importance is best appreciated by consideration of the equations of motion in isentropic coordinates, in which the atmosphere may conceptually be decomposed into a surface layer and an interior layer. Various studies have shown that the equatorward mass transport in the surface layer is proportional to the near-surface eddy heat flux [1, 21, 25, 33, 48]. The poleward mass transport in the interior of the atmosphere above the surface layer is proportional to the eddy potential vorticity (PV) flux [25, 21]. Since the mass transports of the two layers must balance, a closed theory for either the near-surface eddy heat flux or the interior eddy PV flux could form the basis of a complete general circulation theory. Held (1999) argues that the near-surface eddy heat flux is a more suitable target because its dynamics more convincingly support a downgradient diffusive representation. Beyond its role in general circulation theory, the near-surface eddy heat flux is widely used in simplified diffusive 1D energy balance models (EBM) to emulate planetary surface climates and climate feedbacks. In these models, the eddy heat flux is crudely represented with a constant diffusivity parameter; they would benefit from a physically motivated parameterization tied to the mean state ([21]; [9]).

Theories for eddy flux parameterizations have primarily been developed in the quasi-geostrophic (QG) framework, often using 2-layer models on an $f$- or $\beta$-plane with homogeneous forcing conditions ([45]; [19]; [24]; [42]; [35]; [23]; [34]; [59, 60]; [2]; [7]; [15, 16]). In these QG simulations, the mean vertical shear (and thus, the meridional temperature gradient) and the stratification are imposed. Thus, the criticality of the flow is also externally imposed ([20]; [47]; [52]). Early theories view the turbulence in terms of energy cascades in baroclinic and barotropic modes, whereby baroclinic potential energy is converted to eddy kinetic energy at the deformation scale from which it cascades to larger scales through the barotropic mode ([45]; [19]; [35]; [23]; [34]). In these theories, the mixing length is typically set by the Rhines scale ([44]) and the eddy velocity scale depends on the criticality. Later theories, such as the vortex gas scaling theory, have offered an alternative view that recognizes the role of bottom drag ([59, 60]; [15, 16]).

While existing theories explain much of the observed behavior in 2-layer QG systems, the application of these theories to more realistic settings, such as those described by the primitive equations on rotating spheres for atmospheres or oceans, is questionable (see section 3). In QG, the stratification and mean shear are externally imposed, whereas in reality, the stratification and shear are emergent features of the system’s equilibration to the external forcing through macroturbulent processes. As suggested by [47] and [52], this limitation of QG ignores the possibility of eddy-mean flow interactions that, for instance, might inhibit the nonlinear eddy-eddy interactions that support cross-scale energy transfer. Similarly, instability in the 2-layer model ([43]) requires that the criticality, a non-dimensional measure of the vertical scale of baroclinic eddies relative to the tropopause depth, be greater than 1, which implies the eddies extend from the bottom layer into the top layer. However, in a continuously stratified system, no such criterion exists as the eddies can set their own dynamical vertical scale ([10]; [20]; [61]). In fact, while there are exceptions ([65]; [66]; [28]; [29, 30]; [6]; [11]), idealized simulations across a wide parameter regime indicate that planetary atmospheres tend to reside in marginally critical states and exhibit no inverse energy cascade ([52]; [53]). Earth’s atmosphere likewise exhibits no inertial range with a substantial inverse energy cascade ([5]; [52]). Additionally, the assumptions of QG might break down in the surface layer of the atmosphere, for instance, at fronts with steep near-surface temperature gradients ([25]; [51]; [48]). Finally, many QG theories have been developed in the homogeneously forced setting, which does not include eddy momentum flux convergences and intense zonal jets that can modulate eddy behaviors([17]; [53]; [14]). Nikurashin 2010), particularly in the nonlinear decay stage of the eddy life cycle (e.g., [55]; [27]; [40]).

Despite these caveats, QG theories have been applied extensively with mixed results in more realistic settings, such as primitive equation models and comprehensive GCMs. [28] and [29] apply the [23] scalings with empirical success for supercritical Boussinesq channel simulations. However, [28] notes a degradation in the scaling for marginally critical simulations. [4] similarly invoke the inverse cascade scalings, and demonstrate a close relation between the Rhines scale and the mixing length in a series of comprehensive GCM simulations with forcing conditions similar to Earth’s. [8] apply the scalings from [34] and [4] to explain domain-averaged eddy diffusivities for vertically integrated eddy heat fluxes in GCM simulations. Other studies have shown that kinematic effects, associated with a suppression in eddy mixing when eddy phase speeds and the mean zonal flow differ, impact eddy mixing rates ([14]; [56]); this process has traditionally been neglected in classical scalings. [53], who deviate from the classic QG theories, introduce bulk scaling laws that accurately describe domain-averaged atmospheric macroturbulence across a spectrum of climates; however, the scalings exhibit an unexplained power law dependence on rotation rate and planetary radius and depend on the baroclinic zone width. The bulk nature of these scalings and their dependence on baroclinic zone width, which itself is defined in terms of the equilibrated eddy flux statistics, complicate their application as a local closure in an EBM. Overall, no study has demonstrated an accurate local closure valid across a wide regime of planetary atmospheres.

This study introduces new, closed scalings for the near-surface eddy heat fluxes based on the analysis of hundreds of dry GCM simulations spanning a broad range of rotation rates, meridional temperature gradients, seasonality, and convective lapse rates. Many of these simulations, which are very similar to those of [52] and [53], reside in the marginally critical regime, like Earth, where classical scalings fail. The new theory amounts to a diffusive closure that depends on an eddy mixing length, an eddy velocity scale, and a correlation coefficient. The theory directly accounts for suppressive mean flow effects ([14]) on the mixing length and correlation coefficient and uses scaling arguments on the non-dimensional zonal momentum balance to determine the characteristic eddy velocity scale. Importantly, the results depend not on the criticality parameter, which QG theories suggest, but instead on a non-dimensional number akin to the thermal Rossby number ([22]).

The rest of the paper is organized as follows: section 2 describes the idealized GCM, along with diagnostic results and calculation conventions. Section 3 develops the scaling theory. In section 4 we demonstrate the skill of the proposed heat flux closure by implementing the new eddy diffusivity in a simple 1D EBM and accurately reproducing global surface temperature distributions across a wide range of climates. Section 5 discusses implications for general circulation theory and section 6 presents the conclusions.

Idealized GCM simulations and diagnostics

Model details

We use the idealized GCM from [47] and [52] with differences in the configuration described below. The GCM, with modifications, has been used in various idealized studies concerning the general circulation (e.g., [64]; [50]; [63]). It is a primitive equation, spectral-transform model built on the Flexible Modeling System (FMS) available from the Geophysical Fluid Dynamics Laboratory (GFDL). The dynamical core is identical to the spectral dynamical core described in [26]. The horizontal resolution is varied between T42, T85, T127, and T213 depending on rotation rate (see Table 1). The vertical discretization uses a centered difference scheme with 30 vertical levels in a $\sigma$-coordinate system (chapter 7.6 of [12]; [47]) with $\sigma= p/p_s$ for the instantaneous surface pressure $p_s$. All simulations are run to steady-state, with a spin-up period of at least 200 days and eddy statistics calculated over a subsequent period of at least 200 days. No significant changes to the eddy statistics are observed with longer averaging or spinup periods.

Parameter	Values
$\Omega_\ast$	0.5, 1, 2, 4
$\Delta_h$ (K)	40, 60, 80, 100, 120, 140, 160, 180, 200
$\phi_0$ (°)	0, 5, 10, 15, 20
$\gamma$	0.7, 0.8, 0.9, 1.0
$\tau_r$ (days)	30, 60

As in [26], diabatic processes are represented by Newtonian relaxation towards a prescribed radiative equilibrium (RE) profile. Unlike [26], a statically unstable RE profile is used ([47]). Such an equilibrium profile is more realistic, and in conjunction with a dry convection scheme, allows the atmosphere to equilibrate with an emergent combination of the stabilizing effects of both convection and dynamical heat transport. We use the RE profile described in [13], which is similar to that used in [50]. The surface RE temperature profile is described by

Here, $T_s^{\mathrm{avg}}$ is the area-weighted global average surface temperature, $\Delta_h$ is the surface RE pole-to-equator surface temperature contrast, and $\phi_0$ is the latitude of maximum $T_s^e$. As in [47], the vertical structure of the RE profile is described by

Here, $T_t^e$ is the skin temperature at the top of the atmosphere (set to 200 K), $p_0$ is the reference surface pressure (set to 1000 hPa), and the exponent $\alpha= 3.5$ controls the lapse rate of the RE state. This vertical profile is a solution of the two-stream equations for a semigray atmosphere ([47]). Temperatures are relaxed toward the RE profile according to a spatially homogeneous radiative timescale $\tau_r$, which is varied across experiments as described in the next section.

The simulations do not explicitly represent moisture. However, a dry convection scheme, mimicking the stabilizing effect of latent heat release from moist convection, relaxes the atmosphere to a rescaled dry adiabat $\gamma\Gamma_d$ (dry adiabatic lapse rate $\Gamma_d$, rescaling parameter $\gamma\le1$) over a timescale of $\tau_c = 1$ day whenever lapse rates exceed the convective lapse rate [47]. The scheme operates globally, as in [64] and [50].

Boundary layer turbulence that mimics the effect of turbulent dissipation is represented as quadratic drag with an inverse length scale drag coefficient $\mu(\sigma)$ that varies linearly from $\mu_s$ at the surface to 0 at $\sigma_b = 0.85$, the top of the boundary layer [47]. We choose the value $\mu_s = 5 \times10^{-6}\,\mathrm{m}^{-1}$ which is common in the literature (e.g., [47]; [52]) and realistic for terrestrial atmospheric boundary layers. The results presented are not sensitive to modest changes of the bottom drag. Biharmonic hyperdiffusion with a coefficient of $6 \times10^{-5}\,\mathrm{m}^{4}\,\mathrm{s}^{-1}$ is included to damp small-scale instabilities.

Simulations across a wide range of parameters

The simulations survey a wide parameter regime of planetary atmospheres and forcing conditions. Rotation rates $\Omega_\ast= \Omega/\Omega_E$ are varied between 0.5 and 4 times Earth’s rotation rate $\Omega_E$; equator-to-pole RE temperature differences $\Delta_h$ vary between 40 and 200 K; $\phi_0$ (a metric of seasonality) varies between 0 and 20°; the convective lapse rate rescaling parameter $\gamma$ varies between 0.7 and 1.0; and the radiative relaxation timescale $\tau_r$ varies between 30 and 60 days. The exact parameter values are summarized in Table 1. Every combination of $\Omega_\ast$, $\Delta_h$, and $\gamma$ was run with $\phi_0 = 0^\circ$ and $\tau_r = 30$ days. For the remaining values of $\phi_0$ and $\tau_r$, each unique value was run with every combination of $\Omega_\ast$ and $\Delta_h$ with the default values of $\phi_0 = 0^\circ$ (if $\phi_0$ is not varied), $\tau_r = 30$ days (if $\tau_r$ is not varied), and $\gamma= 0.7$. The $\Omega_\ast= 4$ and 2 simulations with $\gamma= 1.0$ and $\Delta_h = 40$ K are omitted from the results because vanishingly small static stabilities complicated their analysis. In total, this yields 323 GCM simulations across a broad range of climate states. Such a large ensemble, like in [53] and [36], is necessary to fully test the generality of any theory. We are unaware of a previous study that has included seasonality as a variable parameter in a comprehensive study of planetary atmosphere eddy fluxes. [28] show that the planetary scale $f/\beta= a \tan\phi$ affects criticality and eddy flux behavior in a baroclinic zone; this scale is set by their channel simulation configuration, since $f$ and $\beta$ are prescribed. However, on a rotating sphere, $f/\beta$ is set by the latitudinal band over which the equilibrated baroclinic zone occurs. All else equal, varying $\phi_0$ can change the latitude of the equilibrated baroclinic zone due to variations in the thermal forcing and interactions with the Hadley circulation.

Simulation eddy energies and criticalities

Eddy energy statistics, averaged over the winter hemisphere, summarize the range of atmospheric behaviors represented in the 323 simulations (Figure 1). For almost all simulations, energy is equipartitioned between eddy available potential energy (EAPE) and EKE, except for some of the $\gamma= 1.0$ simulations (Figure 1). Here, EAPE is defined using Lorenz’s approximate definition in terms of potential temperature variations on pressure surfaces [38] vs. the total eddy kinetic energy EKE$_{\mathrm{TOT}}$. The dashed line indicates the relation EAPE $=$ EKE$_{\mathrm{TOT}}$. The EAPE and EKE$_{\mathrm{BC}}$ exhibit a linear relationship (Figure 1) with a proportionality factor of 2.25, which is consistent with [52] and similar to the expected equipartitioning for the linearly most unstable waves. If there was an inverse cascade of energy to large scales in the simulations, the EAPE would scale with the barotropic EKE (EKE$_{\mathrm{BT}}$), which would be much greater than EKE$_{\mathrm{BC}}$ ([23]; [52]). This relationship is consistent with the absence of an inverse energy cascade, which is common in simulations of planetary atmospheres ([47]; [52]). The exceptions include simulations with low $\Delta_h$ and $\gamma= 1.0$, in which the deformation radii become exceedingly small, criticalities become large, and there is spectral room for a modest inverse cascade to operate. Additionally in these simulations, the vanishingly small static stabilities can introduce errors in the calculation of the EAPE. The EKE$_{\mathrm{BT}}$ and EKE$_{\mathrm{BC}}$ maintain a remarkably linear relationship across many simulations, but the linear fit varies across different rotation rates (Figure 1) consistent with [53]. Faster $\Omega_*$ simulations exhibit a greater proportion of EKE in the barotropic mode.

Analysis of the baroclinic zone-averaged criticality $\xi$ indicates that the vast majority of the simulations reside in a marginally critical or subcritical state (Figure 2), consistent with the eddy energies in Figure 1. Only some low-$\Delta_h$ simulations with high $\gamma$ attain modestly larger $\xi$. Here, $\xi$ is approximated in terms of its near-surface quantity (similar to [28]) as:

The planetary radius $a$ is equal to Earth’s radius for all simulations. The Coriolis parameter is indicated by $f$ and its meridional derivative by $\beta$. Across all simulations, we set $H_t = 10$ km. While the tropopause height does vary slightly across simulations, its variations are small compared to all other quantities in our analyses so the approximation only minimally affects the results. The notation $(\overline{\cdot})$ represents a time and zonal mean along a pressure surface, $(\overline{\cdot})^{s}$ indicates an average evaluated over the near-surface pressure range, and the subscript $(\cdot)_s$ indicates a near-surface-averaged quantity.¹ All near-surface quantities are calculated as averages between 700 hPa and 1000 hPa. This pressure range corresponds to the surface layer in isentropic coordinates referred to in [25]. For the remainder of the manuscript, the baroclinic zone refers to all latitudes in the winter hemisphere where the near-surface eddy heat flux is greater than 50% of the simulation maximum.

The eddy heat flux and diffusivity

Figure 3 shows an example of the zonal and temporally averaged eddy heat flux profile for an Earth-like simulation. The upper boundary of our near-surface pressure range is indicated by the thick black dotted line. Note that the majority of the eddy heat fluxes tend to be confined to the near-surface layer. It is useful to understand these eddy heat fluxes in terms of the diffusivity $D$, defined through a downgradient closure as:

In Earth-like simulations, the near-surface diffusivity $D$ varies with latitude and attains its largest value in the vicinity of the strongest heat fluxes (Figure 3b). The remainder of this work will aim to understand what controls $D$.

Following [59], $D$ can be expressed in terms of an eddy mixing length $\ell_m$ and an eddy velocity scale $v_e$:

The eddy velocity scale $v_e$ is defined as the near-surface rms velocity fluctuations:

The mixing length $\ell_m$ is defined from the rms temperature fluctuations and the mean temperature gradient:

The correlation coefficient $c$ is defined as:

The correlation $c$ is assumed constant in many QG studies ([59], [60]). However, in this study $c$ varies appreciably between different climate states and latitudinally within a given climate state (section 3c). Its effects must be accounted for to obtain an accurate scaling. The remainder of this manuscript will develop scalings for $\ell_m$, $v_e$, and $c$ to obtain a fully closed scaling for $D$.

Theory

We seek scalings for the equilibrated near-surface eddy heat flux diffusivity $D$, which unlike 2-layer QG theories, does not require a priori assumptions of the criticality, the existence of an inverse energy cascade, and the partitioning of the eddy kinetic energy between the barotropic and baroclinic modes. The eddy velocity scale $v_e$ is simply its near-surface value. To make contact with previous theories, we will often focus on non-dimensional quantities. The non-dimensional eddy velocity is given by $v_e^* = v_e/U$, where $U$ in primitive equations is defined as the thermal wind shear over the depth of the troposphere:

We calculate $U$ in terms of $\partial_\phi\overline{\theta}_s$ with applications to 1D EBMs in mind. The non-dimensional eddy mixing length is given by $\ell_m^* = \ell_m/\lambda$, with the outer Rossby deformation radius calculated as:

where $N^2 = \frac{g}{\overline{\theta}_s}\partial_z\overline{\theta}^s$ is defined by its near-surface quantity for simplicity. Note that with these definitions, it follows that the criticality in (2) can be expressed as $\xi= U/\beta\lambda^2$.

The mixing length and its suppression by the mean flow

A popular choice for the mixing length is the Rhines scale:

as suggested by [23, 4, 8]. For slower rotation rates, $\ell_{\mathrm{Rh}}$ systematically overestimates $\ell_m$; also, biases persist between different rotation rates (fig:4a). The scaling only provides a good fit if the planetary parameters are not varied sufficiently (particularly the rotation rate, as in [4]). To rectify these errors, we consider the kinematic effects that mean flow advection of potential temperature anomalies will have on $\ell_m$.

Following other studies ([8]), dimensional considerations suggest that the characteristic decorrelation timescale for eddy mixing should depend on the eddy velocity scale $v_e$ and $\beta$ through:

where $\gamma_\beta$ is an inverse decorrelation timescale. The product $v_e\gamma_\beta^{-1}$ yields the length scale:

This is similar to the Rhines scale, except that the length scale is expressed in terms of the near-surface eddy velocity $v_e$, not the barotropic eddy velocity. Unlike other studies, we do not assume that the barotropic component dominates the heat flux (as it does in 2-layer QG, where the baroclinic mode is itself equivalent to the temperature anomalies). In fact, in the simulations the contribution of the baroclinic and barotropic modes to the total eddy heat flux is proportional to the EKE in that mode (not shown). As such, the near-surface velocity $v_e$ is the relevant velocity scale.

Various studies identify a suppression of the meridional eddy diffusivity when the zonal mean flow differs from the zonal eddy phase speed ([14, 32, 56]). Under these conditions, tracer filaments may be peeled out of an eddy by the mean zonal flow before substantial meridional mixing can occur, which reduces the mixing length and suppresses the diffusivity. Accounting for suppressive effects and following [14], the mixing length is expressed in terms of the characteristic eddy decorrelation timescale $\gamma^{-1}$ and velocity scale $v_e$ as:

where $k$ and $l$ are the zonal and meridional wavenumber of the eddy respectively, $\kappa^2 = k^2+l^2$, $c_w$ is the phase speed of the eddy, and $\overline{u}(p)$ is the mean zonal velocity at pressure $p$. First, assume the wavenumber ratio $k^2/\kappa^2$ is an O(1) constant. We substitute $\gamma=\gamma_\beta$ and assume that the characteristic eddy scale is set by similar dimensional reasoning such that $k\sim\ell_\beta^{-1}$. To understand the eddy phase speed $c_w$, consider that a precursor for baroclinic instability is phase-locking of two Dopplershifted counter-propagating Rossby waves in the vertical direction. The phase speed at which this phase-locking occurs will take on a value characteristic of the mid-troposphere mean zonal flow [14]. Consequently, the eddy phase speeds of the baroclinically unstable upper and lower tropospheric Rossby waves will differ from the upper and lower tropospheric mean zonal wind, assuming a non-zero vertical shear $U$. Therefore, at the near-surface pressure level, we can express the difference between the phase speed and mean near-surface zonal flow $\bar{u}_s$ as some fraction $\alpha$ of the mean shear $U$:

Plugging in these quantities into (12) and dividing the numerator and denominator by $\gamma_\beta$ yields our mixing length scaling accounting for mean flow suppression:

which in non-dimensional form is expressed as:

The suppression factor becomes negligible when $v_e \gg\alpha U$. The Charney model of baroclinic instability [10] suggests steering levels (the vertical level at which the eddy phase speed equals the mean zonal flow) below the mid-troposphere level. Indeed, setting $\alpha= 0.1$ in (15), consistent with a steering level around 700 hPa, yields accurate mixing length predictions across the simulations (Figure 4b). The resulting suppression factor is significant and varies in magnitude across all simulations, becoming negligible for high $\Omega_*$ simulations and increasing to as strong as 50% for some low rotation rate simulations. Overall, the scaling explains 98% of the variance across all simulations in log-space (Fig. 4b). Figure 4c shows the same scaling (15), but applied locally at various individual latitudes (indicated by the marker color) throughout all simulations shown in panel b. The scaling is able to capture the local mixing length behavior accurately: it explains 95% of the variability in the non-dimensional mixing length across more than 3,000 points plotted in log-space. The scaling over-predicts the suppression at some low latitude points, where heat fluxes are already weak and induced by surface-trapped waves.

The eddy velocity scale

[52] and [53] demonstrate that across a wide regime of planetary atmospheres, the vertically and domain-integrated (across the baroclinic zone) EKE scales with the EAPE, which in turn scales with mean available potential energy (MAPE) integrated over the baroclinic zone. However, since the baroclinic zone is itself a feature of the equilibrated turbulence, a fully closed scaling should ideally be independent of this feature. Additionally, these scalings were developed to explain the bulk turbulence properties of atmospheres. For a local scaling, a similar equipartitioning argument can be invoked to relate the local EKE to the MAPE integrated over the characteristic mixing length $\ell_m$ of eddies at a given latitude. However, in Appendix A, we show that the functional form of the mixing length expression (15) complicates using this assumption to derive a closed velocity scaling. Additionally, the GCM results indicate quantitative errors when applying this assumption locally.

Instead, we use scaling arguments on the zonal momentum balance to determine an eddy velocity scale. Consider the time- and zonal-average zonal momentum balance in the upper troposphere:

Here we have neglected a frictional term (friction primarily acts in the boundary layer through the quadratic drag formulation) and the vertical eddy momentum flux convergence term. We scale the variables as $\bar{u}=U\hat{u}$, $\bar{v}=V\hat{v}$, $\bar{w}=VH_t/a\hat{w}$, $z=H_t\hat{z}$, and $\overline{u'v'}=V_e^2\hat{u}^*\hat{v}^*$, where terms with $\hat{(\cdot)}$ are non-dimensional. The scaling for $\bar{w}$ derives from scaling the zonally averaged continuity equation. The zonal and meridional eddy velocities are assumed to scale isotropically as $V_e$. The upper troposphere mean zonal velocity scale $U$ may be expressed in terms of the meridional temperature gradient using the thermal wind balance (7). We define the non-dimensional parameter $\chi$:

Dividing through by $fV$ and plugging in the scalings, the zonal momentum balance (16) becomes:

The parameter $\chi$ is akin to a planetary scale Rossby number comparing the thermal wind shear $U$ against the planetary velocity scale $fa$. For rapidly rotating planets like Earth, $\chi\ll1$. This number is similar to the “thermal Rossby number” defined in [22] for describing axisymmetric Hadley circulations. For $\chi\ll1$, the eddy momentum flux divergence (EMFD) term must balance the Coriolis term at leading order. For the EMFD term to be O(1), it must be true that the non-dimensional upper tropospheric eddy velocity scales as:

In QG theories, the ageostrophic meridional velocity $V$ is assumed to scale with $U$ multiplied by the Rossby number Ro. Given $\chi$ is similar to Ro for the planetary scale primitive equations, it is tempting to assume $V\sim U\chi$. Doing so renders the non-dimensional eddy velocity $V_e/U$ a constant, yet the non-dimensional eddy velocity varies by an order of magnitude across all simulations.

Instead, consider that the zonal momentum balance of the upper troposphere (18) can be expressed in terms of just three non-dimensional numbers $\chi$, $V_e/U$, and $V/U$. The latter number $V/U$ does not appear at leading order anywhere else in the equations of motion, continuity, or heat; therefore, we expect it to be a function of the other two non-dimensional parameters in the zonal momentum balance. If a power law functional form is assumed for this relationship, then $V_e/U$ must also be related to $\chi$ through some power law. If it is furthermore assumed that the upper tropospheric eddy velocity $V_e$ is proportional to the near-surface eddy velocity $v_e$ (which is expected, given that the baroclinic kinetic energy primarily resides in the upper and lower troposphere), then the non-dimensional eddy velocity $\hat{v}_e^*$ should also follow a power law relationship with $\chi$. Indeed, an excellent fit to the simulation data is:

The dimensional near-surface eddy velocity may be expressed as:

The dimensional eddy velocity exhibits a slight dependence on $f$ and $a$ and exhibits a slightly stronger than square root dependence on the dimensional meridional temperature gradient. Ultimately, the dependence on $\chi$ is physically motivated, but the power law fit itself is an empirical result. We speculate the inverse dependence on $\chi$ in (20) captures a reduction in EKE due to the effects of EMFD across strong zonal jets, whose strength tend to increase with decreasing $\Omega_*$ through thermal wind balance. The strong zonal shears in these stronger jets promote faster conversion rates through EMFD of EKE to mean flow kinetic energy in the upper troposphere [38]. This can broadly act to reduce eddy lifetime and dissociate lower troposphere waves from their upper troposphere counterparts, possibly explaining the higher ratios of $\mathrm{EKE}_{\mathrm{BC}}$ to $\mathrm{EKE}_{\mathrm{BT}}$ for lower $\Omega^*$ in fig:1c.

The scaling (20) is able to explain variations in $v_e^*$ within a given rotation rate and across different rotation rates, including 96% of the $v_e^*$ variance in log-space across all simulations with $\Delta_h \ge80$ K (Figure 5). It explains 93% of the variability in local $v_e^*$ in log-space across all simulations with $\Delta_h \ge80$ K (fig:5c). The scaling suggests $v_e^*$ scales with latitude like $f^{4/5}$ (accounting for the $f$-dependence in $U$)—the GCM results corroborate this prediction, as $v_e^*$ generally increases with latitude. The scaling degrades for the lowest $\Delta_h$, $\gamma= 0.7$ simulations and locally for some low latitude points. Appendix B explains the cause for this discrepancy: in these weakly forced simulations, the thermal stratification of the troposphere is primarily maintained by convection ([47]), which quickly damps near-surface temperature anomalies before any significant baroclinic instability can ensue. The scalings are only expected to hold in the regime where baroclinic eddy fluxes modulate the thermal structure of the atmosphere. Note that the low scaling coefficient of 0.07 (shown in fig:5b) is likely because the near-surface eddy velocity is plotted, but the scaling is intended for the much higher upper-troposphere eddy velocity.

The correlation coefficient

The correlation coefficient $c$, which appears in the expression (4) for $D$, is considered negligible in many studies or is implicit in the definition of the mixing length. Thompson and Young (2006) diagnose the correlation coefficient and find it is approximately constant across the vast majority of their 2-layer QG simulations. However, $c$ varies appreciably across our simulations as well as spatially in a given climate state. This variation was noted by Schneider and Walker (2008), but a quantitative scaling behavior for $c$ has not been proposed. Swanson and Pierrehumbert (1997) estimate near-surface eddy heat fluxes in Earth’s mid-latitude storm tracks using reanalysis data. and propose a simple stochastic Lagrangian damped-advective model for understanding turbulent heat transport by eddies. In this model, Lagrangian parcels are advected by a stochastically forced velocity field from an initial state with a prescribed mean temperature gradient, velocity variance, and Lagrangian integral decorrelation timescale $\tau_d$. The temperature anomalies of the parcels relative to the environment are damped according to some timescale $\tau_Q$. This model of the turbulent motions of stirred parcels yields the following analytical expression for $c$ consistent with the functional form plotted in their Figure 10:

Swanson et al. (1997) fit their model parameters to the data and no closed formula is presented in terms of mean flow parameters. However, we can estimate the Lagrangian integral decorrelation timescale using our eddy mixing decorrelation timescale and the velocity variance with our velocity scaling in the previous section. The decorrelation timescale follows from the velocity and mixing length scalings as $\tau_d = \ell_m/v_e$, yielding:

Variations in $\tau_Q$ are likely to depend intimately on model-specific details of convective and radiative processes that might not generalize well to real planetary systems like Earth. On the contrary, variations in $\tau_d$ should be related to dynamical mechanisms general to all planetary systems. Therefore, our proposed scaling for $c$ depends only on $\tau_d$, and we set $\tau_Q = 6$ days. Figure 6 plots the GCM-calculated correlation coefficient against $\tau_d$ averaged over the baroclinic zones in each simulation in panel a and applied locally in panel b. The plotted lines are the functional relationship for $c$ in (22) with different values of $\tau_Q$ for each line as indicated in the plot legend. The upper and lower lines represent upper and lower limits based on our convective and radiative relaxation timescales of 1 and 30 days, respectively (note that some simulations have a 60 day relaxation timescale, and consequently tend to have lower $c$). The dashed line indicates the empirically determined best estimate of $\tau_Q = 6$ days.

The scaling (22) captures the general increase in $c$ with rotation rate observed in the GCM simulations. However, the accuracy of the scaling, especially when applied locally, is limited by the complicating effects of convection, whose sporadic and intense nature cannot be accounted for with a simple relaxation timescale as with radiative processes. In lower latitudes and lower $\Delta_h$ simulations (indicated with the lighter shading in panel a), we expect the points to fall closer to the $\tau_Q = 1$ day line due to the increased prevalence of convection, as is indeed the case in 6a and 6b. Also, as expected, simulations with higher $\gamma$ tend to have lower $c$ for a given $\tau_d$, since convection is less prevalent. It is reassuring that maximum and minimum values of $c$ correspond to the theoretical upper and lower bounds (Figure 6b). The complicated interplay between convective and radiative processes will ultimately introduce errors in local predictions of the diffusivity. Fortunately, over-predictions in $v_e$ somewhat compensate under-predictions of $c$, since convection has opposite effects on each quantity.

The diffusivity

Combining the three component scalings of the previous sections for $\ell_m^*$, $v_e^*$, and $c$ in (15), (20), (22) yields an expression for the near-surface eddy diffusivity $D$. The non-dimensional diffusivity $D^* = D/U\lambda= cv_e^*\ell_m^*$ is given by:

The dimensional diffusivity $D$ is given by:

We have expressed (25) to facilitate direct comparisons with the Held and Larichev (1996) scaling that predicts $D \sim(U\xi)^{3/2}\beta^{-1/2}$. The numerator of (25) is analogous to this, except our eddy velocity scales like $v_e \sim(U\chi^{0.4})^{3/2} = U^{0.9}(fa)^{0.6}$. The scaling (25) also features an additional term in the denominator accounting for the direct suppressive effects of kinematic mean flow advection on the near-surface mixing length, as in Ferrari and Nikurashin (2010). Importantly, the dimensional eddy diffusivity depends on $\chi$, not the criticality $\xi$.

The scaling (25) explains 97% of the variability in $D^{*}$ in log-space across the $\Delta_{h} > 60$ K simulations (Figure 7a). As expected from the velocity scaling results, $D^{*}$ is over-estimated for the low-$\Delta_{h}$, low-$\gamma$ simulations where the thermal structure is primarily maintained by convection (see Appendix B). The scaling is also able to capture local variations across three orders of magnitude of $D^{*}$ (Figure 7b). Panel c shows the baroclinic zone-averaged dimensional diffusivity predictions, for which it explains 87% of the variability across all simulations. For both $D^{*}$ and $D$, the primary source of error outside of the low-$\Delta_{h}$ simulations is the scaling for $c$. The scaling overestimates $c$ for faster rotation rate simulations, causing a slight over-prediction of $D$ in the faster rotation rate simulations and an under-prediction in some slower rotation rate simulations.

The dimensional near-surface eddy heat flux is the essential quantity that describes eddy effects on the large-scale climate. Figure 8 shows the predictions from the proposed scaling against GCM-calculated values averaged over baroclinic zones in panel a and applied locally in panel b. The scaling predictions are obtained by multiplying the scaling for the dimensional diffusivity $D$ in (25) by $\partial_{y}\overline{\theta}_{s}\cos\phi$. The scaling is able to explain 91% of the variance in simulation-averaged eddy heat fluxes and 79% of the variance in local eddy heat fluxes across all simulations. Again, the primary source of error in $D$ comes from the scaling for $c$, especially in local predictions. We show in section 5 that the diffusivity scaling (25), when applied as a local parameterization in an EBM, can explain the variability in pole-to-equator temperature contrasts across the GCM simulations.

Applying the scalings in an EBM to explain global climate

The EBM

We approximate the near-surface energy balance using a simplified diffusive EBM in which the near-surface potential temperature in steady-state is set by a balance between diffusive meridional heat transport and radiative heating. Using the near-surface equilibrium temperature profile $\theta_E$, we can express the time-dependent zonal-mean near-surface potential temperature $\bar{\theta}_s$ as:

This EBM can be run to steady state and compared against the near-surface GCM $\bar{\theta}_s$. In the low-latitude tropical region, the Hadley circulation is the dominant mechanism of heat transport and must be parameterized to obtain any accurate global distribution of temperature. While the mean advection of heat by the Hadley circulation is not a diffusive process, the effect of its heat transport can be approximated diffusively as in [39] with a so-called “Hadley adjustment.” With this adjustment, the full diffusivity $D$ in (26) is expressed as:

where $w(\phi) \in[0,1]$ is a weighting function that smoothly transitions the dominant heat transport mechanism from the Hadley cell to mid-latitude transient eddies between the tropics and the mid-latitudes. The weighting $w$ is a function of both local latitude $\phi$ and Hadley cell edge latitude $\phi_h$, which is diagnosed from the GCM output. Such Hadley adjustments are common in the literature ([37]; [54]; [3]) and correct for the distinct dynamical differences between the Hadley circulation and the mid-latitudes with respect to their heat transport behavior. We adopt a similar functional form as in [54] and [3] for the weighting function:

The weighting function has a Gaussian shape with width $\phi_h/2$. Similar to [39], $D_h$ is set to a high value, $8 \times10^6\ \mathrm{m^2\,s^{-1}}$, to emulate the efficient heat transport and flatter temperature gradients observed in the tropics. While this is a crude treatment of the Hadley cell heat transport, the main goal of the EBM here is to assess the ability of the scalings to reproduce the mid-latitude heat transport by transient eddies. This functional form of the Hadley correction emulates sufficiently accurate Hadley cell transport, which allows us to focus on the effect of our scalings on the atmospheric heat transport poleward of the tropics. For the eddy diffusivity $D_e$, we use the scaling (25).

The EBM is run for all equinox ($\phi_0 = 0^\circ$) simulations for every combination of $\Delta_h$ and $\Omega_*$ in Table 1 for the $\gamma$ values of 0.7, 0.8, 0.9, and 1.0. Since (25) exhibits no dependence on stratification, the EBM solutions for different $\gamma$ only differ due to differences in $\phi_h$ between the GCM simulations. We do not experiment with the $\phi_0 > 0^\circ$ simulations because upper tropospheric easterlies develop near the equator in the cross-equatorial winter Hadley cells, inducing up-gradient heat transport by the Hadley circulation in the winter hemisphere near the equator. To model these asymmetrically forced climates, a more sophisticated treatment of the Hadley circulation heat transport is necessary.

EBM Results

The EBM accurately reproduces the modeled global temperature contrasts, and can accurately reproduce the observed behaviors for varying $\Omega_*$ and varying $\Delta_h$, with some slight biases with rotation rate due to biases in the scaling for $c$ (Figure 9). The EBM captures the rate of increase in pole-to-equator temperature contrast with increasing $\Delta_h$ as well as the faster increase of the pole-to-equator temperature gradient with $\Delta_h$ at faster rotation rates than at slower rotation rates. Slight discrepancies with varying $\gamma$ are primarily due to changes in the Hadley circulation heat transport. The Hadley circulation heat transport generally weakens with increasing $\gamma$ (see Figure 4 in [64]). The $D_h$ chosen to emulate the Hadley circulation heat transport was chosen to optimize the $\gamma= 0.7$ simulations. Therefore, using a constant $D_h$ for all simulations leads to an over-estimation of Hadley circulation heat transport for higher $\gamma$ simulations, thus causing an under-estimation of the pole-to-equator temperature gradient in the EBM.

The EBM reproduces well the overall structure of the modeled $\bar{\theta}_s$ profiles. Examples for $\Omega_* = 1$ and $\Delta_h$ varying between 60, 120, and 180 K are shown in fig:10a, and examples for $\Delta_h = 120$ K simulations for four different rotation rates are shown in fig:10b. For the more strongly forced high $\Delta_h$ simulations, advection by the mean circulation, including the Hadley, Ferrel, and polar cells, is stronger, causing stronger variations in the slope of the temperature profile with latitude. A simple diffusive EBM cannot capture such variations, but it still captures the global equator-to-pole temperature contrast.

Discussion

Held (1999) discusses how an understanding of the near-surface eddy heat flux can form the basis of a general circulation theory. The zonally averaged meridional mass transport in isentropic coordinates can be approximated in terms of the near-surface eddy heat flux as (see [1]; [25]):

The diffusive closure for the eddy heat flux developed in section 3 can be used to predict $\overline{v'\theta'^{s}}$, and therefore $V$, as a function of the mean flow. What remains is an a priori prediction of the vertical stratification of the extratropical troposphere, which may consist of either a prediction of the lapse rate or the tropopause height. Schneider (2004) and Schneider and Walker (2008) propose one such explanation, namely that atmospheres tend to organize themselves into an O(1) criticality state across a broad range of atmospheres. However, it remains unclear how general this relation should be ([28]; [30]). An alternative route to a complete understanding of the mid-latitude circulation might be to predict the macroturbulent transport of PV in the upper troposphere directly through some diffusive closure, which is also challenging since eddy PV fluxes more strongly distort the mean PV gradient than near-surface eddy heat fluxes distort the temperature gradient (Held 1999). The results of this paper suggest that kinematic effects may be important in our understanding of these upper tropospheric eddy fluxes. Also based on the scaling analysis in section 3, the parameter $\chi$ should be relevant for the upper troposphere PV fluxes, too. In addition to a complete understanding of the mid-latitude circulation, a closed theory of the tropical Hadley circulation and how it interacts with eddies of mid-latitude origin will be necessary for a closed general circulation theory.

While this paper focuses primarily on dry dynamics, with the stabilizing effect of latent heat release encoded through the dry convection scheme, the scalings should be relevant for describing moist atmospheres, such as Earth. The scaling for the dimensional diffusivity in section 3 exhibits no dependence on the stratification, and the GCM simulations with varied convective lapse rates corroborate this scaling prediction. This suggests that the effects of moisture on tropospheric stability may not be important for near-surface eddy heat fluxes in Earth-like climates. Indeed, O’Gorman and Schneider (2008) find that EKE in the baroclinic zones of idealized moist atmospheres scales with the dry MAPE, suggesting that the dry dynamics might generalize well to moist atmospheres. However, as the effects of convection in the dry simulations indicate, the scalings might become erroneous in regions and climates where convection dominates the maintenance of the vertical stratification. Future work will further examine the relevance of the proposed scalings in moist atmospheres.

Conclusions

We present simulation results of near-surface eddy heat fluxes from 323 idealized GCM simulations across a wide parameter regime of pole-to-equator temperature gradients, rotation rates, seasonality, and vertical stratification. We propose a new scaling theory for the near-surface macroturbulent heat transport that accurately reproduces the simulated heat flux, even in parameter regimes where classical 2-layer QG theories fail. There are two key distinctions between our proposed theory and other classical theories, such as that proposed by Held and Larichev (1996):

1. We incorporate the kinematic effect of mean flow suppression arising from a difference in eddy phase speeds and the near-surface mean zonal flow (Ferrari and Nikurashin 2010). When this mean flow suppression is weak, the mixing length scaling reduces to the Rhines scale, as in Held and Larichev (1996). However, for the majority of simulations, this kinematic reduction of the mixing length cannot be ignored to obtain accurate scalings.

2. Our eddy velocity scales with $U\chi^{0.4}$ instead of $U\xi$. The parameter $\chi$ is a planetary scale Rossby number comparing the strength of the thermal wind shear to the planetary velocity scale $fa$. A power law relationship between the eddy velocity and $\chi$ is suggested from nondimensionalizing the zonal momentum balance for the upper troposphere and examining the leading order balance. This scaling implicitly accounts for non-linear interactions between eddies and the mean zonal jet (through eddy momentum fluxes), a process neglected in homogeneously forced QG theories.

Classical QG theories predict that the near-surface eddy heat flux diffusivity scales cubically with the temperature gradient (Held and Larichev 1996). The scalings developed in this work indicate that the diffusivity depends nearly linearly on the temperature gradient in the limit of weak suppression and exhibits nearly zero dependence on the temperature gradient in the limit of strong suppression. Earth-like climates reside between these two limits. The scalings (24) and (25) also feature a latitudinal dependence which is absent in other theories. The scalings do not exhibit any dependence on stratification, so we expect they should apply to moist atmospheres, such as Earth’s. The complicated dependence of the eddy diffusivity on the thermal wind shear and latitude might also explain phenomena in Earth’s storm tracks, such as the midwinter suppression in EKE noted in the Pacific storm track region (Nakamura 1992; Schemm and Schneider 2018; Hadas and Kaspi 2021).

The proposed scalings can be used as a diffusive closure in a simple EBM to understand the global distribution of surface temperatures across different climate states (Section 5). Integrating a realistic treatment of Hadley circulation heat transport into the EBM would provide a simple, fully closed 1D model that could explain the global distribution of temperature, as well as interactions between the zonally averaged extratropical and the Hadley circulations. Altogether, the proposed scalings can form the basis of a general circulation theory, upon which an understanding of the thermal stratification of the extratropical troposphere and the dynamics of the Hadley circulation may be built.

Acknowledgments. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. 2139433. R.E. acknowledges very helpful discussions with Raffaele Ferrari for the mixing length formulation.

Data availability statement. The idealized GCM used is publicly available and scripts for generating the simulation data will be made publicly available on Github after manuscript acceptance.

Appendix A

The relationship between the eddy velocity and mixing length

[53] approximate the bulk EAPE integrated over a baroclinic zone using a Taylor series expansion over the meridional width of the baroclinic zone temperature profile. They then assume this EAPE is equipartitioned with EKE. However, for a local closure, this assumption can lead to errors (evidenced by local mixing lengths that are substantially smaller than the deformation radii, and therefore the baroclinic zone width, in many of the simulations in fig:4c). Instead, we may use a local definition of MAPE in terms of our locally defined mixing length to estimate the EAPE and assume it is proportional to the local EKE at a given latitude. Doing so will yield a statement of equipartitioning between the local eddy velocity scale $v_e$ and the local eddy mixing length $\ell_m$. To see this, first consider Lorenz’s original approximate expression for EAPE in terms of potential temperature variance on pressure surfaces ([38]). Similar to [53], we estimate the vertically averaged quantity in terms of its near-surface value. The vertically and zonally averaged EAPE per unit mass is then approximated as (with units $\mathrm{m^2\,s^{-2}}$):

The heat capacity at constant pressure is denoted by $c_p$ and $\kappa= R_d/c_p$. The reference pressure $p_0 = 1000\ \mathrm{hPa}$ is approximately equal to the near-surface pressure $p_s$, so their ratio is neglected in what follows. The potential temperature variance may be exactly expressed in terms of the temperature gradient using our definition of the mixing length in (5). Substituting (5) into eq:A1 gives:

Converting the vertical derivative from pressure to height coordinates using the hydrostatic relation, expressing this gradient in terms of the deformation radius (8), and utilizing the ideal gas law, the EAPE can be succinctly expressed in terms of $U$ by:

This relation asserts that the EAPE is proportional to the MAPE over the characteristic mixing length of eddies at a given latitude. Assuming equipartitioning between EKE and EAPE furthermore yields:

where we assumed that the near-surface EKE is proportional to the total EKE (this is an accurate assumption for our simulations). This can be rewritten in non-dimensional terms as $v_e^* \sim\ell_m^*$. This is exactly the relation derived in [23] using inverse energy cascade arguments. If $\ell_m^*$ scales like $\xi$, then we recover the result of Held and Larichev (1996). Unfortunately, the scaling for $\ell_m^*$ (15) itself depends on $v_e^*$ in a complicated way through the mean flow suppression, so it cannot simply be plugged into (A4) to obtain a closed relation for $v_e^*$ as in [23]. Additionally, the eddy velocity scale will depend not only on the mixing length near the surface, but also the vertical structure and extent of the eddy heat flux. While the simulation results indicate that $v_e^*$ and $\ell_m^*$ are qualitatively similar, their proportionality has quantitative errors, particularly between different rotation rates. Instead, in the main text, we use scaling arguments in section 3b to determine the eddy velocity scale $v_e$.

Appendix B

Convective interference with the eddy velocity scaling in weakly forced climates

The inaccuracies apparent for the $\gamma= 0.7$ simulations with $\Delta_h \le60$ K can be explained by the effects of convection. Consider that if the atmosphere were in a state of radiative-convective equilibrium, there would be no meridional circulation and therefore no horizontal EKE to maintain the principal zonal momentum balance assumed in section 3b. Generally, the scaling (20) should only be valid if radiative heating and cooling are predominantly balanced by eddy heat fluxes. These heat fluxes act as a source of wave activity that propagates into the upper troposphere and induces the upper tropospheric EMFD whose zonal momentum tendency must be balanced by the Eulerian-mean Ferrel circulation. For low-$\Delta_h$ simulations, we expect weak temperature gradients, and therefore weaker eddy heat fluxes due to the lower EKE suggested from the scaling in (20). To stabilize the vertically unstable RE profile in the extratropics, the dry convection scheme compensates for the weak eddy activity to maintain a statically stable vertical thermal structure. This transition between eddy-dominated and convection-dominated maintenance of the thermal stratification of the troposphere for weakly forced climates has been identified by [47] and [52].

Analysis of the potential temperature variance budget, which is closely related to the EAPE budget, confirms the interference of convection in the eddy energy budget for weakly forced simulations. The steady-state potential temperature variance budget, averaged over a baroclinic zone such that the transport terms can be neglected, is:

The first term on the left is the rate of production of $\theta$ variance, and the two terms on the right are the conversion from $\theta$ variance to EKE and destruction by diabatic processes, either through radiative ($Q_r$) or convective ($Q_c$) processes. Time- and zonal-anomalous diabatic heating is denoted by $Q'$. The atmospheric thermal structure will be eddy-controlled if generation of temperature variance is largely balanced by baroclinic conversion to EKE; it will be convection-controlled if temperature variance generation is balanced by diabatic destruction. fig:B1a plots the ratio of destruction of temperature variance by convection to the rate of conversion of temperature variance to velocity variance on the y-axis against the near-surface lapse rate $\Gamma_s$ relative to the convectively neutral lapse rate $\gamma\Gamma_d$ for each simulation.

For many low-$\Delta_h$ simulations, convection destroys a large fraction of the produced temperature variance, thus limiting the baroclinic conversion of potential to kinetic energy (fig:B1a). This fractional rate of destruction rapidly increases when the mean state near-surface lapse rate approaches a convectively neutral state ($\Gamma_s \approx\gamma\Gamma_d$). This ratio of destruction to baroclinic conversion can be represented by a simple power law in terms of $\Gamma_s$ and $\gamma\Gamma_d$. In principle, one might be able to correct for this convective interference with a scaling predicting the diabatic destruction of temperature variance given the near-surface lapse rate. However, this calculation is beyond the scope of this work. The scaling in (20) overestimates $v_e$, and therefore $D$, for the simulations with $\Delta_h \le60$ K and $\gamma= 0.7$. There is a strong latitudinal dependence indicating that convective damping is more prevalent at lower latitudes (fig:B1b); this is expected given that near-surface heating rates are stronger at lower latitudes.