Comparative terrestrial atmospheric circulation regimes in simplified global circulation models: II. energy budgets and spectral transfers

TL;DR

This study investigates how planetary rotation rate influences atmospheric energy budgets and spectral transfer processes using a simplified global circulation model, revealing trends in energy cycles and spectral slopes across different rotation regimes.

Contribution

It provides new insights into the spectral characteristics and energy transfer mechanisms of atmospheric circulation patterns as a function of planetary rotation rate.

Findings

Maximum Lorenz energy cycle occurs at half Earth's rotation rate.

KE and APE spectra show distinct slopes depending on rotation rate.

Enstrophy and APE predominantly cascade downscale across all regimes.

Abstract

The energetics of possible global atmospheric circulation patterns in an Earth-like atmosphere are explored using a simplified GCM based on the University of Hamburg's Portable University Model for the Atmosphere. Results from a series of simulations, obtained by varying planetary rotation rate {\Omega} with an imposed equator-to-pole temperature difference, were analysed to determine the heat transport and other contributions to the energy budget for the time-averaged, equilibrated flow. These show clear trends with {\Omega}, with the most intense Lorenz energy cycle for an Earth-sized planet occurring with a rotation rate around half that of the present day Earth. KE and APE spectra, E_K(n) and E_A(n) (where n is total spherical wavenumber), also show clear trends with \Omega, with n^{-3} enstrophy-dominated spectra around \Omega* = \Omega/\Omega_E = 1, where \Omega_E is the rotation rate of the Earth) and steeper (\sim n^{-5}) slopes in the zonal mean flow with little evidence for the n^{-5/3} spectrum anticipated for an inverse KE cascade. Instead, both KE and APE spectra become almost flat at scales larger than the internal Rossby radius, L_d, and exhibit near-equipartition at high wavenumbers. At \Omega* << 1, the spectrum becomes dominated by KE with E_K(n) \sim 2-3 E_A(n) at most wavenumbers and a slope \sim n^{-5/3} across most of the spectrum. Spectral flux calculations show that enstrophy and APE are almost always cascaded downscale, regardless of {\Omega}. KE cascades are more complicated, however, with downscale transfers across almost all wavenumbers, dominated by horizontally divergent modes, for \Omega* \lesssim 1/4. At higher rotation rates, transfers of KE become increasingly dominated by rotational components with strong upscale transfers (dominated by eddy-zonal flow interactions) for scales larger than L_d and weaker downscale transfers for scales smaller than L_d.

Equations

$$E_Z(n)=C_Z{\beta }^2{n}^{-5},C_Z\sim 0.5,$$

$$E_R(n)=C_K{ϵ}^{2/3}{n}^{-5/3},C_K\sim 5,$$

$$E_K(p)$$

$$E_A(p)$$

$$\gamma (p)=R/[\Lambda (p)p{\partial }_p\hat{\theta }],$$

$${\partial }_t{E}_K(p)$$

$${\partial }_t{E}_A(p)$$

$$G(p)$$

$$C(p)$$

$$F_{A\uparrow }(p)$$

$$F_{K\uparrow }(p)$$

$$S(p)$$

$$J(p)$$

$$\frac{1}{4\pi a^{2}g}\int\int\int_{0}^{p_{s}}Q{\rm d}p\,{\rm d}x\,{\rm d}y=\dfrac{1}{4\pi a^{2}}\iiint Q\text{d}m=\Big{<}Q\Big{>}.$$

$$\frac{\text{dAZ}}{\text{d}t}=\text{GZ}-\text{CZ}-\text{CA},$$

$$\frac{\text{dAE}}{\text{d}t}=\text{GE}+\text{CA}-\text{CE},$$

$$\frac{\text{dKZ}}{\text{d}t}=\text{CZ}-\text{CK}-\text{FZ},$$

$$\frac{\text{dKE}}{\text{d}t}=\text{CK}+\text{CE}-\text{FE},$$

$${\theta }^{\prime }({\mathbf{x}}_h,p)=\sum _{n\ge 0}\sum _{-n\le m\le n}{\theta }_{nm}^{\prime }(p)Y_{nm}({\mathbf{x}}_h)$$

$$⟨\omega \Phi ⟩=\sum _{n\ge 0}\sum _{-n\le m\le n}(\omega ,\Phi {)}_{nm}$$

$$(\omega ,\Phi {)}_{nm}=\text{Re}\{{\omega }_{nm}^{†}{\Phi }_{nm}\}$$

$$\mathbf{u}$$

$$\zeta$$

$$\delta$$

$$\widehat{\mathbf{a}\cdot \mathbf{b}}=\sum _{n\ge 0}\sum _{-n\le m\le n}(\mathbf{a},\mathbf{b}{)}_{nm}$$

$$(\mathbf{a},\mathbf{b}{)}_{nm}$$

$${\mathrm{div}}_h(\mathbf{a}{)}_{nm}^{†}{\mathrm{div}}_h(\mathbf{b}{)}_{nm}\},$$

$$E_A^{nm}$$

$$E_K^{nm}$$

$$\displaystyle E_{KZ}^{n}(t)=\dfrac{1}{4}\dfrac{a^{2}}{n(n+1)}\big{(}|\zeta_{n}(m=0,t)|^{2}+|\delta_{n}(m=0,t)|^{2}\big{)},$$

$$\displaystyle E_{KE}^{n}(t)=\dfrac{1}{4}\dfrac{a^{2}}{n(n+1)}\sum_{m=-n}^{n}\big{(}|\zeta_{n}(m,t)|^{2}+|\delta_{n}(m,t)|^{2}\big{)};m\neq 0,$$

$$\frac{\partial \zeta }{\partial t}=-({\mathbf{u}}_r\cdot \nabla )\zeta -\mathcal{D}$$

$$\frac{\partial G}{\partial t}=-\zeta ({\mathbf{u}}_r\cdot \nabla )\zeta -\zeta \mathcal{D}.$$

$$\frac{\partial G^n}{\partial t}=J^n+{\mathcal{D}}^n$$

$$J^{n}=-\frac{1}{4}\sum_{m=-n}^{n}\big{[}\zeta^{nm^{\dagger}}(\mathbf{u}_{r}\cdot\mathbf{\nabla}\zeta)^{nm}+\zeta^{nm}(\mathbf{u}_{r}\cdot\mathbf{\nabla}\zeta)^{nm^{\dagger}}\big{]}$$

$$\sum _{n=0}^N{J}^n=0.$$

$${\mathcal{H}}^n=\sum _{l=N}^n{J}^l,$$

$${\partial }_t{E}_K^{nm}(p)$$

Full text

1 Introduction

Atmospheric circulation can be said to occur because of the propensity of fluid motion to transfer heat energy from regions of net heating to regions of net cooling. Horizontal heat transfer is also part of the overall processing and conversion of energy within the atmospheric “heat engine”, in which local imbalances between incoming radiant energy from the parent star and thermal emission tend to modulate the internal energy and increase the potential energy of the atmosphere. Dynamical processes then act to convert such potential energy into various bulk forms of motion in the form of kinetic energy before dissipative processes ultimately reconvert this back to heat again. Precisely how potential energy is transformed into kinetic energy depends strongly on the dynamical constraints governing atmospheric motion, and this in turn depends on a number of external factors, such as the planetary size, mass and rotation, the overall mass of the atmosphere and its composition. Heat, momentum and other tracers may be transported in latitude either by direct meridional overturning in an axisymmetric Hadley-type circulation, or via non-axisymmetric eddies through systematic covariances between meridional velocity and temperature fluctuations.

Heat transport, of course, represents only part of the overall cycle of energy conversion within a planetary atmosphere. A common approach to the analysis of energy conversions is the one based on the work of Lorenz (1955), in which energy reservoirs and exchanges are partitioned between kinetic and (available) potential energy, and between zonally averaged and eddy components. Energy exchanges within the Earth’s global circulation have been analysed in this way for many years (Peixóto and Oort, 1974; Li et al., 2007; Boer and Lambert, 2008), although very few studies have examined the Lorenz energy cycle for other planets (e.g. Lee and Richardson, 2010; Pascale et al., 2013; Schubert and Mitchell, 2014; Tabataba-Vakili et al., 2015). In the context of an exploration of how the global circulation regime changes within a simplified model atmosphere, it is of significant interest to examine how the cycle of energy conversions changes throughout parameter space. This is investigated in the present study for an Earth-like planet at various rotation rates, based on the set of simulations presented by Wang et al. (2018) using the Hamburg PUMA model, and the results presented in Section 4.

The Lorenz approach provides insight into how the atmospheric heat engine transfers energy from the planetary scale, zonally-symmetric flow into non-axisymmetric ”eddies”. But this is only a crude measure of how energy passes from scales that are directly energised by solar heating and radiative cooling into other scales of motion, that takes little account of the macroturbulent processes that distribute energy from the forcing scales towards those affected by dissipation. In this context, the concept of geostrophic turbulence, first introduced by Charney (1971), has been an important paradigm for theories of the large-scale planetary atmospheric and oceanic circulations.

The flow in geostrophic turbulence tends to be a highly chaotic, quasi-2D (horizontal), quasi-geostrophic flow, typically featuring an inverse energy cascade if small-scale forcing is present. Planetary rotation, large aspect ratio (between horizontal vertical scales) and statically stable stratification all act to bring planetary atmospheric flows into quasi-horizontal (quasi-2D) motion. Small-scale forcing is usually envisaged as being provided either by baroclinic instability occuring at scales comparable to the Rossby deformation radius ($L_{D}$) or by small-scale convection, as is possibly in the case of Jovian planets (see e.g. Ingersoll et al., 2000; Read et al., 2007, 2015)). The energy generated through such processes then becomes a small-scale “agitator” of the inverse energy cascade in the barotropic mode, though the precise mechanism for energising this mode is still not well understood. It is a typical feature of 2-D isotropic (in a 2-D planar sense, or horizontally isotropic) turbulence that energy goes from small scale to large scale through a spectrally-local inverse cascade. The direct consequence of such an inverse energy cascade is the emergence of large circular eddies with no preferred directionality. In the presence of a non-negligible background vorticity gradient (e.g. $\beta$-effect), however, it was shown by Rhines (1975) that such large-scale eddies becomes anisotropic, causing an elongation of structures in the zonal direction and ultimately leading to the formation of zonal jets.

Galperin and Sukoriansky (Sukoriansky et al. (2002), Galperin et al. (2006)) recently proposed the paradigm of zonostrophic turbulence as an attempt to characterise universally the regime of eddy-driven multiple zonal jets on a $\beta$-plane. Under a strong $\beta$-effect, it is proposed that flows can develop into the regime of zonostrophic turbulence which is characterised by a strongly anisotropic KE spectrum with a steep ($-5$) slope for the zonally symmetric flow component and a classic Kolmogorov-Batchelor-Kraichnan (KBK) $-5/3$ slope in the non-axisymmetric eddy/residual modes. The segments of the spectra in this regime take the universal form (when appropriately non-dimensionalised):

[TABLE]

where $\epsilon$ is the energy pumping rate of the small-scale excitation (which, in previous 2-D numerical studies of zonostrophic turbulence, is represented as an artificial energy input at a specific wavenumber $n_{\xi}$, see e.g. Huang et al. (2001) and Galperin et al. (2004). In a real planetary atmosphere, this can be due to barotropic or baroclinic eddies, or in the case of gas giants, possibly from small-scale moist convection as well.). $C_{K}$ is the universal Kolmogorov-Kraichnan constant, while barotropic simulations (e.g. Chekhlov et al., 1996; Huang et al., 2001) suggest that $C_{Z}$ can vary between $0.1$ and $1.0$.

Precisely how these and other regimes emerge and under what conditions has been largely unexplored in general until recently, leaving open many questions as to the nature of the circulation of various planets within and beyond our Solar System. In the present work, therefore, we analyse a set of numerical model simulations, obtained using a simplified global circulation model in which atmospheric flows in an Earth-like planetary atmosphere are driven by simple linear relaxation towards a prescribed (steady) zonally-symmetric temperature field (on a timescale $\tau_{R}$) and dissipated by a linear Rayleigh drag (with prescribed timescale $\tau_{fr}$. We vary various planetary parameters (especially the planetary rotation rate, but also the surface friction timescale) and allow the simulation to equilibrate over a timescale of order 20 Earth years. The basic model and the phenomenology of the circulation regimes were described in a companion paper (Wang et al., 2018), which clearly demonstrated a systematic sequence as $\Omega$ was varied from $\Omega^{\ast}=1/16$ to $\Omega^{\ast}=8$. The regimes obtained ranged from a super-rotating, barotropically unstable cyclostrophic atmosphere at the lowest values of $\Omega^{\ast}$ to a highly geostrophically turbulent circulation with multiple zonal jets at $\Omega^{\ast}>>1$ via more Earth-like, geostrophic states with simpler patterns of jets and baroclinic eddies that were either regular and periodic or chaotic in nature.

In this paper, however, we focus on analysing the budgets of kinetic and potential energy and associated heat transport. We begin with an analysis of the global exchanges of energy within the well known framework of the Lorenz energy cycle, but then extend the analysis to consider the more detailed exchange of energy and enstrophy between different scales via the spectra of kinetic and available potential energy and the principal spectral fluxes as a function of spherical harmonic total wavenumber. The computation of spectral fluxes provides arguably the most detailed and precise means of evaluating the direction and intensity of turbulent cascades by directly computing the exchanges of various forms of energy and enstrophy between different scales, as represented in a decomposition of flow structure projected onto a spectrum of spherical harmonics. This approach has been applied for several years to studies of kinetic energy exchanges within the Earth’s atmospheric circulation in numerical simulations and assimilated analyses (e.g. Boer and Shepherd, 1983; Shepherd, 1987; Koshyk and Hamilton, 2001; Burgess et al., 2013), but relatively few such analyses have been extended to include potential energy exchanges and conversions (Lambert, 1984; Augier and Lindborg, 2013; Malardel and Wedi, 2016). They have proved able to provide important insights, however, into how the atmosphere transfers key properties between different scales through nonlinear interactions, and in particular the potential impacts of various parameterization schemes on the modelled cascades of energy and enstrophy (Malardel and Wedi, 2016). A similar approach has recently been applied to the Earth’s oceans, at least on a local scale in the context of mesoscale eddies (Scott and Wang, 2005; Scott and Arbic, 2007), and even to the kinetic energy budget of Jupiter’s atmosphere (Young and Read, 2017), in which both systems reveal the existence of a double cascade (involving both up- and down-scale segments), energised on scales close to the Rossby deformation radius.

Section 3 presents the framework for analysis of the budgets of kinetic and potential energy and the spectral transfers of energy and enstrophy. Results for the various terms in the Lorenz energy budget as a function of planetary rotation rate are presented and discussed in Section 4 while Section 5 provides an overview of trends in the spectra of kinetic and potential energy. The spectral fluxes of energy and enstrophy as a function of $\Omega^{\ast}$ are presented in Section 6 and the overall results are discussed in Section 7.

2 Model setup and experiment design

The model used is the Portable University Model of the Atmosphere (PUMA; e.g. see Fraedrich et al. (1998); Frisius et al. (1998); von Hardenberg et al. (2000)), consisting of a spectral dynamical core solving the dry primitive equations on a sphere, based on the code developed by Hoskins and Simmons (1975). Temperature, divergence, vorticity, and $\ln p_{s}$ (where $p_{s}$ is the surface pressure) are the prognostic variables. The model domain uses finite difference discretization in the vertical using 10 equally spaced $\sigma$ levels (where $\sigma=p/p_{s}$). The integration in time was carried out with a filtered leap-frog semi-implicit scheme (Robert, 1966).

Thermal forcing was applied via a linear Newtonian relaxation towards a prescribed (axisymmetric) temperature field which was constant in time, with a relaxation timescale $\tau_{R}$. The complete restoration temperature field (with equator-to-pole temperature difference of $60\ \text{K}$) was intended to represent a distribution similar to the Earth and is shown in Fig 1.

The radiative timescale, $\tau_{R}$, was set to 30 Earth days in the free atmosphere, decreasing to $2.5$ Earth days at $\sigma=1.0$.

A smooth, spherical planet was assumed in each case with no surface topography. Dissipation consisted of a combination of linear Rayleigh drag towards rest in the lowest two model levels (with a timescale decreasing from zero to $\tau_{F}=0.6$ Earth day at the surface) and a $\nabla^{8}$ hyperdiffusion acting separately on temperature, vorticity and divergence.

Simulations were run from an isothermal state at rest at a series of planetary rotation rates, from $\Omega^{\ast}=\Omega/\Omega_{E}=1/16$ to $\Omega^{\ast}=8$, for a period equivalent to 10 Earth years. Horizontal resolution was set to T42 for slowly rotating simulations ($\Omega^{\ast}\leq 1$), T127 for faster rotating simulations with $\Omega^{\ast}=1,2,4$, and with T170 reserved for simulations with $\Omega^{\ast}=8$. The computed diagnostics were averaged over the final model year from each run. Further details on the model setup and experiment design are presented by Wang et al. (2018).

3 Analysis of energy budgets

The kinetic energy, $E_{K}$, (KE) and available potential energy, $E_{A}$, (APE) of an atmosphere can be defined (e.g. Augier and Lindborg, 2013) in pressure coordinates by

[TABLE]

where ${\bf u}$ is the horizontal component of the total velocity ${\bf v}=({\bf u},\omega)$, $\omega=Dp/Dt$ is the vertical velocity in pressure coordinates, $\theta$ is potential temperature, $\widehat{(.)}$ denotes an average over an entire pressure level and $(.)^{\prime}$ departures therefrom. $\gamma(p)$ is defined as

[TABLE]

where $R$ is the gas constant, $\Lambda(p)=(p_{R}/p)^{\kappa}$, $p_{R}$ is a reference pressure and $\kappa=R/c_{p}$.

The energy budget of atmospheric KE and APE can be derived from combining Eqns. (2) and (3) with equations of motion and conservation of potential temperature (e.g. Augier and Lindborg, 2013) to obtain

[TABLE]

Here $G(p)$ is an APE generation term e.g. due to differential heating, $C(p)$ is the conversion from APE to KE, $F_{K\uparrow}(p)$ and $F_{A\uparrow}(p)$ are vertical fluxes of KE and APE, respectively, and $D_{K}(p)$, $D_{A}(p)$ are diffusion terms with:

[TABLE]

$p_{s}$ and $\Phi_{s}$ are surface pressure and surface geopotential, respectively, and $\delta_{ps}$ is one when $p=p_{s}$ and zero otherwise. The $S(p)$ and $J(p)$ terms are adiabatic processes but which do not conserve total available energy $E_{K}+E_{A}$. However, these terms have been shown to be negligible in the global mean (Siegmund, 1994; Augier and Lindborg, 2013), and so will not be considered further in this analysis.

3.1 Formulation of the Lorenz energy cycle

The generation and growth of non-axisymmetric waves and other disturbances (designated as “eddies”) within terrestrial planetary atmospheres requires the conversion into eddy kinetic and potential energy from other forms of energy in the background environment, the ultimate source of which is solar or stellar irradiation. This process of energy conversion can be illustrated and quantified most simply with the classical Lorenz energy cycle (Lorenz (1955)) which provides a framework for formulating a global mean atmospheric kinetic energy and available potential energy budget, as well as the conversion rates between the zonal mean (indicated by $[.]$ and “eddy” ($.^{*}$) components of these energy forms. In addition, we designate time-averaged quantities by $\overline{(.)}$ and mass-weighted, vertically-integrated areal averages (e.g. of a quantity $Q$) by

[TABLE]

Following e.g. Peixóto and Oort (1974) and James (1995), the conservation equations for kinetic and potential energy can be written as

[TABLE]

where KZ, KE, AZ, and AE refer to zonal mean kinetic energy, eddy kinetic energy, zonal mean available potential energy, and eddy available potential energy respectively. The conversion rates among these components are as shown in Fig. 2 and are as defined in pressure coordinates e.g. by James (1995). $GZ$ and $GE$ are diabatic generation terms (if positive) for AZ and AE, while $FZ$ and $FE$ are dissipation terms for KZ and KE.

There are at least two principal mechanisms of eddy generation within planetary atmospheres: barotropic instability and baroclinic instability. The physical pictures of these two different eddy-generation mechanisms can be clearly distinguished from the viewpoint of energy conversion. Eddies generated through barotropic instability are fed directly by the zonal mean kinetic energy, implying the conversion route of KZ$\rightarrow$KE to be important. Baroclinic instability, on the other hand, converts available potential energy into eddy kinetic energy through the so-called ‘sloping convection’, which corresponds to the conversion route of AZ$\rightarrow$AE$\rightarrow$KE. Therefore, the dominating mechanism of eddy generation can be appreciated by comparing the relative intensity and direction of CK and CE.

3.2 Spherical harmonic transformation

The horizontal structure of the flow and other quantities (such as energy conversion rates) can be further decomposed into spectra by projection onto suitable sets of eigenfunctions. A scalar function (e.g. $\theta^{\prime}$) on a sphere can be transformed into spherical harmonic spectral space via

[TABLE]

where $n$ is the total and $m$ the zonal wavenumber indices. $Y_{nm}$ are spherical eigenfunctions with $\nabla_{h}^{2}Y_{nm}=-n(n+1)Y_{nm}/a^{2}$ (where $\nabla_{h}$ is the horizontal gradient operator and $\nabla_{h}^{2}$ the corresponding Laplacian). The horizontal mean of the product of two scalar variables is then

[TABLE]

with

[TABLE]

where $\text{Re}\{X\}$ is the real part and $X^{\dagger}$ is the complex conjugate of a complex number $X$ (Augier and Lindborg, 2013).

For the horizontal velocity field, $\bf u$, a decomposition into divergent and rotational (non-divergent) components can be performed via the Helmholtz decomposition

[TABLE]

with reference to the horizontal streamfunction $\psi({\bf x}_{h},p)$ and the horizontal velocity potential $\chi({\bf x}_{h},p)$. Using this decomposition, we can obtain the vorticity $\zeta$ and the horizontal divergence $\delta$:

[TABLE]

This decomposition is then used to calculate the horizontal mean of a scalar product between two horizontal vector fields ${\bf a}$ and ${\bf b}$

[TABLE]

with

[TABLE]

(see e.g. Augier and Lindborg, 2013).

Using Eqn. 16 for scalars and Eqn. 21 for vector fields, the spectral versions of APE and KE can be obtained respectively as:

[TABLE]

The APE or KE spectrum can be further decomposed into a zonal mean spectrum and an eddy (or residual) spectrum as the following (e.g. for KE):

[TABLE]

where $E_{KZ}^{n}$ and $E_{KE}^{n}$ represent the zonal and eddy (residual) part of the spectrum respectively, such that $E_{K}^{n}=E_{KZ}^{n}+E_{KE}^{n}$.

3.3 Calculation of spectral enstrophy fluxes

The nonlinear spectral enstrophy transfer flux (Boer and Shepherd (1983), Shepherd (1987), Burgess et al. (2013))can provide more detailed insights into the enstrophy redistribution among different wavenumbers through nonlinear eddy-eddy interactions. Starting from the vorticity equation

[TABLE]

where $\mathbf{u}_{r}=(u_{r},v_{r})$ is the rotational velocity and $\mathcal{D}$ represents the effects on vorticity evolution due to divergence and other vorticity sources and sinks. Multiply by $\zeta$ to obtain the equation for enstrophy ($G=\frac{1}{2}\zeta^{2}$)

[TABLE]

In spectral space this can be rewritten as

[TABLE]

where the interaction term $J^{n}$ is

[TABLE]

Note that interaction terms only redistribute enstrophy among wavenumbers so

[TABLE]

We can then define the enstrophy spectral flux as

[TABLE]

where the sign is adopted conventionally such that a positive value corresponds to a forward cascade while a negative value corresponds to an inverse cascade.

The interaction terms and spectral fluxes of enstrophy can be further decomposed into contributions from eddy-eddy interactions and eddy-zonal flow interactions (Burgess et al. (2013)). Nonlinear interaction terms of enstrophy due to purely eddy-eddy interactions, $J^{n(e)}$, can be obtained from Eq (30) but carrying out the sum in $m$ for $m\neq 0$ only. The contribution to $J^{n}$ through eddy-zonal mean interactions is then simply $J^{n(z)}=J^{n}-J^{n(e)}$. In this way, the spectral flux of enstrophy can be decomposed as $\mathcal{H}^{n}=\mathcal{H}^{n(z)}+\mathcal{H}^{n(e)}$.

3.4 Spectral energy budget

The spectrally-resolved energy budget can be obtained by inserting Eqn. 23 and Eqn. 24 into Eqns. 5 and 6, resulting in (Augier and Lindborg, 2013):

[TABLE]

where $G^{nm}$ is the spectral APE generation term, $C^{nm}$ is the spectral conversion term, $T_{K}^{nm}$ and $T_{A}^{nm}$ are the spectral transfer terms (of KE and APE, respectively) due to non-linear interactions. $L^{nm}$ is a spectral transfer term due to Coriolis forces and $F_{K\uparrow}^{nm}$ and $F_{A\uparrow}^{nm}$ are vertical fluxes. $D_{K}^{nm}$ and $D_{A}^{nm}$ are diffusion terms. These terms are computed via

[TABLE]

The spectral energy and tendency terms obtained in the previous sections are functions of time, zonal wavenumber $m$, total wavenumber $n$, and pressure $p$. To obtain a one dimensional spectrum or spectral flux from these terms, a dependence upon $n$ alone would be preferable for reasons of simplicity of interpretation. The time dependency is removed by averaging the resulting spectral quantities over multiple time steps. Following a summation over zonal wavenumbers and a vertical integration over a pressure range from $p_{b}$ at the lowest level (usually the surface) to $p_{t}$ at the top level, the vertically integrated KE spectrum is obtained via

[TABLE]

and the vertically integrated KE spectral flux via

[TABLE]

where $\sum_{k\geq n}\sum_{-k\leq m\leq k}$ denotes a cumulative sum (from large to small wavenumbers). Other spectral quantities can be similarly vertically integrated and summed. Note that the cumulative summation is performed from large wavenumbers to small wavenumbers and that all spectral fluxes (barring conversion and vertical fluxes) are conserved, meaning the cumulative sum over all wavenumbers $n$ should add up to zero.

4 Lorenz energy cycles as a function of $\Omega^{\ast}$

In this section we compute the various terms in the Lorenz energy budget for each of the 8 rotation rate simulations spanning $1/16\leq\Omega^{\ast}\leq 8$ and explore the main trends in energies and conversion rates. Fig. 3 and Fig. 4 show how the terms in the globally- and time-averaged Lorenz energy cycles vary with rotation rate, $\Omega^{\ast}$, and thermal Rossby number, $\mathcal{R}o_{T}$, defined as in Wang et al. (2018) by

[TABLE]

where $\Delta\theta_{h}$ is the equator-to-pole potential temperature difference, $a$ the planetary radius and $R$ the specific gas constant. Fig. 3(a) plots the magnitudes of the various energy reservoirs, expressed in units of 100 kJ m*-2*. This shows the zonal mean potential energy, AZ, increasing monotonically with $\Omega^{\ast}$ and decreasing with $\mathcal{R}o_{T}$, though with only a shallow variation over the range computed. KZ, on the other hand, exhibits a maximum around $\Omega^{\ast}=1/8$ but then decreasing rapidly with $\Omega^{\ast}$ for higher rotation rates. The maximum in KZ corresponds to $\mathcal{R}o_{T}\simeq 5$. The eddy kinetic and available potential energies follow similar trends to each other, rising with $\Omega^{\ast}$ to a shallow maximum around $\Omega^{\ast}=1/2$ ($\mathcal{R}o_{T}\simeq 0.3$ and then decreasing with $\Omega^{\ast}$ quite sharply at higher rotation rates. KE is seen to dominate over AE by a factor $\sim 2-3$ until $\Omega^{\ast}=4$, beyond which AE $>$ KE.

Variations in the mean energy conversion rates (in units of W m*-2* are shown in Fig. 3(b). CZ represents the direct conversion of zonally-averaged APE to KE by zonal mean overturning circulations, and essentially reflect the relative strengths of the (thermally-direct) Hadley cells and the (thermally-indirect) Ferrel cells within the global circulation. This behaves more or less as one might expect, with strongly positive values of CZ at low rotation rates ($\Omega^{\ast}\lesssim 0.4$; $\mathcal{R}o_{T}\gtrsim 1$) where the direct Hadley circulation is dominant, but becoming negative at higher rotation rates, where the circulation becomes geostrophic and the Ferrel cells become stronger. This largely reflects the increasing dominance of baroclinic eddies at high rotation rates. The barotropic conversion term, CK, also undergoes a similar reversal of sign around ($\Omega^{\ast}\simeq 0.3$; $\mathcal{R}o_{T}\simeq 1$), indicative of a transition from barotropically energised eddies (with CK $>0$) at low rotation rates to baroclinically energised eddies (with CE $>0$) at higher rotation rates. This interpretation is consistent with $\mathcal{R}o_{T}\simeq 1$ as the criterion since $\mathcal{R}o_{T}\lesssim 1$ is also a criterion for strong baroclinic instability. Both conversion rates (CZ and CK) evidently become vanishingly small as $\Omega^{\ast}$ becomes large. This likely reflects the tendency for geostrophic velocities to decrease with increasing $\Omega$, together with the corresponding vertical velocities and the most energetic length scales for eddies.

As anticipated in the previous subsection, the strength of the baroclinic conversion rates, CE and CA, reach their maximum at $\Omega^{\ast}=1/2$. This also agrees with results from Del Genio and Suozzo (1987) in which baroclinic eddies peak in energy conversion efficiency at $\Omega^{\ast}=1/2$, and is also consistent with the peak of meridional eddy heat flux at $\Omega^{\ast}=1/2$ as suggested by the peak in CA $+$ CK in Fig. 4 (see also Fig. 7 of Wang et al. (2018)). Pascale et al. (2013), however, found that the rotation rate corresponding to this peak in CE was sensitive also to other parameters, notably the strength of the bottom friction, which may account for slight differences in this peak being seen in other studies (e.g. Kaspi and Showman (2015)).

5 Kinetic and available potential energy spectra

In the following section, the KE and APE spectra are employed as a diagnostic to reveal insights into processes such as the jet formation mechanism and the transfer of energy and enstrophy across different scales.

5.1 Slow rotation ($\Omega^{\ast}<1$)

The energy spectra for rotation rates $\Omega^{\ast}<1$ are shown in Fig. 5(a)-(d). These simulations were carried out at a horizontal resolution of T42, so it is likely that there are some artifacts in the spectra due to model diffusion at the highest wavenumbers. This is apparent as a steepening of the spectrum for wavenumbers $k\gtrsim 30$. But for the rest of the spectrum it appears to trend towards a self-similar form with a slope close to $n^{-5/3}$ from a $n^{-3}$ spectrum as $\Omega^{\ast}$ decreases. This trend towards a KBK-like slope would appear to suggest a trend towards a spectrum dominated by an energy-dominated cascade, though it is not immediately apparent whether this would entail upscale or downscale transfers of KE. This will be discussed further in the next section on spectral fluxes, but assuming that kinetic energy is converted from potential energy by baroclinic instabilities at scales close to the Rossby deformation wavenumber, $n_{D}$, and given the trend in $n_{D}$ towards small $n$ as $\Omega^{\ast}$ decreases (e.g. see Wang et al. (2018); Kaspi and Showman (2015)), it would seem likely that the cascade would be predominantly downscale if injection of KE is taking place mainly at near-planetary scales.

Even numbered wavenumbers generally appear stronger in KE than odd numbered modes, with the opposite trend in APE, which reflects the symmetries of the winds and temperature structure about the equator with annual mean forcing that is symmetric about the equator. The spectra are also dominated on average by kinetic energy for these rotation rates with a typical ratio $E_{K}(n)/E_{A}(n)\sim 2-3$.

In contrast to these changes to the eddy part of the spectrum, the zonal components of both of the energies follows an $n^{-5}$ slope for the mid-range in wavenumber space. This seems to be a generic feature of the spectrum at almost all rotation rates and is discussed further below. Most spectra (at all values of $\Omega^{\ast}$) show some evidence of a steep drop-off in energy at the highest wavenumbers, indicative of a region where the model dissipation is active, although this is more apparent in some runs than in others. This may indicate that some runs (notably the case for $\Omega^{\ast}=1$) could be somewhat under-dissipated, which may have possible implications for the shape of the spectrum at other wavenumber ranges. Given limitations on computational resources available to us for this study, however, some compromises were necessary in our attempts to capture both realistic inertial ranges and an adequate range of dissipation (we are very far, of course, from the conditions appropriate for Direct Numerical Simulation). Such problems are not unique to this study (see, e.g. Hamilton et al., 2008, for examples at much higher resolution), though the sensitivity of the KE and APE spectra on dissipation should ideally be investigated further in future work.

5.2 Spectra at $\Omega^{\ast}\geq 1$

Figure 5(e)-(h) shows the KE and APE spectra of simulations with $1\Omega_{E}\leq\Omega\leq 8\Omega_{E}$. These simulations were performed at T170 resolution (except for the $\Omega^{*}=1$ simulation at T127). The Earth-like run at $\Omega^{\ast}=1$ (Figure 5e) exhibits a $n^{-3}$ slope between wavenumbers 20 and 90 as well as a fairly consistent $n^{-5}$ slope in the zonal component. It is interesting to note that both KE and APE behave fairly similarly in this region. At lower wavenumbers ($k\lesssim 10$) the spectrum flattens with a hint of a segment tending towards the KBK $n^{-5/3}$ form between $2-3\lesssim k\lesssim 10$, suggestive of an energy-dominated upscale cascade.

This is broadly consistent with various previous studies based on observational/reanalysis datasets of Earth’s atmosphere (e.g. see Baer (1974), Boer and Shepherd (1983), Koshyk et al. (1999)), indicating the probable existence of a forward enstrophy cascade inertial range. The zonal spectrum in both APE and KE is characterised by a much steeper $-5$ slope, however, which is still not well-understood despite the prediction of a $-5$ slope in the early work e.g. of Rhines (1975). Rhines showed that, near the cross-over scale from Rossby waves to turbulence (i.e. near the Rhines wavenumber $k_{R}\simeq(\beta/U)^{1/2}$, where $k$ represents the total dimensional wavenumber $k=n/a$), the typical wind speed is $U\sim\beta/k^{2}$. Since $E=kE(k)\simeq 1/2U^{2}$, the $-5$ power law can then be revealed as $E(k)\sim\beta^{2}k^{-5}$. But this does not explain the extended $-5$ slope solely in the zonal KE spectrum. A study by Huang et al. (2001) (and further developed e.g. by Sukoriansky et al., 2002; Galperin et al., 2004, 2006, 2010) identified the $-5$ slope associated with the zonal KE spectrum with a so-called zonostrophic regime, prevalent with strong planetary rotation and weak dissipation. They further suggested that the shape of the zonal spectrum can be qualitatively explained by the stabilising effect of $\beta$ on the zonal jets, According to Huang et al. (2001), Sukoriansky et al. (2002), Galperin et al. (2006) and others, the $\beta$-effect modifies the necessary condition for barotropic instability from $\partial_{yy}[u]=0$ to $\beta-\partial_{yy}[u]=0$ somewhere within the flow domain. This means the stability criterion is eased by the $\beta$-effect from $\partial_{yy}[u]\neq 0$ to $\beta-\partial_{yy}[u]\neq 0$, which allows the existence of velocity inflection points ($\partial_{yy}[u]$=0) and enables more energy to reside in the zonal modes without violating the stability criterion. But this argument is somewhat heuristic, leading some (e.g. Danilov and Gurarie, 2004) to question whether a well-defined power-law scaling relationship even exists for the zonal KE spectrum.

Figs. 5f-h show the APE and KE spectra of the regime of multiple zonal jets. The strong zigzag feature of the zonal KE spectrum at small spherical wavenumbers is due to the hemispheric symmetry of the predominantly zonal structures across the globe. The classic $-5/3$ slope of an inverse energy cascade in KE cannot be found, indicating that neither the ‘classical’ picture of 2-D isotropic turbulence, nor the zonostrophic regime of Sukoriansky et al. (2002) and Galperin et al. (2010), are applicable to the multiple jet flows found in this regime (at least under the conditions explored here). As shown by Wang et al. (2018), the energy-containing wavenumber, Rossby deformation wavenumber, and the Rhines wavenumber are not widely separated in this regime, although the KE spectrum is likely to be energised on scales close to $L_{D}$ through conversion of APE by baroclinic instabilities. Such a lack of scale separation might suggest that the inverse energy cascade, initiated around the deformation wavenumber, $n_{D}$ (Wang et al., 2018), through eddy-eddy interactions, has been significantly suppressed because of a lack of room to develop an energy-conserving inertial range, although this will be further investigated below with respect to the computed spectral fluxes. Such closeness of the Rhines and deformation scales also suggests that the simulated flows are far from the conditions necessary to observe fully developed zonostrophic dynamics (e.g. Galperin et al., 2006, 2010).

With increasing rotation rate (see Figs. 5f-h), the maximum of the zonal component moves to higher wavenumbers and the $n^{-3}$ slope that could be so well identified at $\Omega^{*}=1$ becomes inclined towards an even steeper slope at higher wavenumbers. This is likely due to the effect of the model-inherent hyperdiffusion, as a result of which the region over which a $n^{-3}$ slope can be discerned becomes smaller and smaller. The same hyperdiffusion effect can be identified for the $n^{-5}$ slope of the zonal component (see also the discussion of this power law above). Overall, however, the slope of the zonal kinetic energy spectrum in these regimes is not well understood. Nevertheless, our work can report that the zonal component of the APE spectrum does have the same slope.

The ratio of KE to APE also varies significantly with wavenumber in this regime, with APE dominating over KE at $n=2$ and with APE/KE ranging from $\sim 30$ at $\Omega^{\ast}=1$ to more than 108 at $\Omega^{\ast}=8$. At higher wavenumbers, however, within the $n^{-3}$ region then KE is seen to dominate with a KE/APE ratio that ranges from around 3 at $\Omega^{\ast}=1$ down to O(1) at $\Omega^{\ast}=8$.

6 Spectral transfer fluxes of energy and enstrophy for varying rotation rates

In this section, we present enstrophy and energy spectral fluxes of the PUMA-S simulations discussed previously. This is intended to provide a more detailed view of the general spectral transfer pathways within our simulated atmospheric circulations across a range of parameter space, in particular using the spectral energy budget formulation of Augier and Lindborg (2013). From such a spectral energy budget, we can answer the question of how the energy of macroturbulent fluid motion is transported between scales and converted between APE and KE. More specifically, we are interested in seeing at which scale kinetic energy is inserted into the system and where this energy ends up. We distinguish between two modes of transfer between scales, depending upon whether transfer is spectrally local, between nearby scales (representing a conventional energy cascade), or non-local, in which energy is directly transferred from one scale to another across large wavenumber intervals (akin to a “waterfall” - M. E. McIntyre, priv. comm.).The latter can occur, for instance, between disturbances of arbitrary wavenumber and the ($m=0$) zonal flow.

To identify this interaction between eddies and the zonal flow (the eddy-mean flow interaction) we perform an additional decomposition into zonal and eddy components. This decomposition was achieved by taking the eddy component (via $X_{eddy}=X-[X]$ of each input variable (i.e. $u,v,\omega,\Phi,T$) and recalculate all the fluxes from this (for which the $\gamma$ value obtained from the initial flux calculation is used). The zonal component is then obtained as the residual. For spectral fluxes, the “eddy” component encompasses eddy-eddy interactions, while the “zonal” component consists of residual interactions with the zonal mean flow (i.e. combining eddy-zonal and zonal-zonal interactions). In the text below the terms wavenumber, total wavenumber and $k$ are used synonymously.

6.1 Spectral enstrophy fluxes

Figure 6 shows a sequence of profiles of spectral enstrophy fluxes, $\mathcal{H}_{n}$, covering the full range of $\Omega^{*}$. Theoretical discussions of quasi-geostrophic turbulence (e.g. Charney, 1971; Salmon, 1978, 1980) suggest that the flux of enstrophy should be downscale throughout the range of scales, although this depends upon the flow satisfying conditions for quasi-geostrophy. In the cases shown in Fig. 6, this trend is more or less consistent with this expectation, but with some exceptions. At low rotation rates ($\Omega^{*}\leq 1/4$ or $\mathcal{R}o_{T}\gtrsim 1$), where the quasi-geostrophic approximation is not likely to be valid, enstrophy fluxes are generally quite weak at all scales, though with a slight trend to increase towards the smallest resolved scales where dissipation becomes significant. At these small scales, the zonal-eddy interaction appears to dominate.

For $\Omega^{*}\gtrsim 1/2$ (or $\mathcal{R}o_{T}\lesssim 1$), however, enstrophy fluxes become significantly larger and positive at moderate to small scales, as anticipated for quasi-geostrophic turbulence. At larger scales, however, there is a small tendency for $\mathcal{H}_{n}$ to become negative, indicating a weak upscale cascade range and implying a spectrally-local net source of enstrophy at the wavenumber $n_{s}$ at which $\mathcal{H}_{n}$ changes sign. This wavenumber gradually increases with $\Omega^{*}$ from around $n=5-6$ at $\Omega^{*}=1/2$ to $n\simeq 75$ at $\Omega^{*}=8$. The magnitude and distribution of enstrophy flux between eddy-eddy and zonal-eddy terms also changes with $\Omega^{*}$. Enstrophy fluxes appear relatively weak for $\Omega^{*}\lesssim 1$ and dominated by zonal-eddy interactions. At higher rotation rates, the eddy-eddy interactions become more dominant, especially at higher wavenumbers, indicating a conventional, spectrally-local cascade of enstrophy, which becomes much stronger at $\Omega^{*}=2$ and 4 and self-similar in shape as it moves towards higher wavenumbers as $\Omega^{*}$ increases. Fluxes become weaker at the highest $\Omega^{*}$ as $n_{s}$ approaches the resolution limit and dissipation presumably acts to damp the dynamics. This should be investigated further, though will require model resolutions that were beyond the scope of what was feasible in the present study.

6.2 Spectral energy fluxes

Figure 7 shows spectral fluxes for KE ($\Pi_{K}$), APE ($\Pi_{A}$), and the total energy $\Pi=\Pi_{K}+\Pi_{K}$ as well as the cumulative conversion $\mathcal{C}$ from APE to KE for simulations across the full range of $\Omega^{*}$. The fluxes have been decomposed into eddy-eddy (“eddy”) and residual zonal (“zonal”) interaction components (as well as between rotational (non-divergent) and divergent components of the flow, shown in more detail for KE in Figure 8 below). The terms presented here are integrated over the whole pressure range of the simulated atmospheres (see Section 3.4). The figure shows that in all cases the total energy flux $\Pi$ is always positive, signifying a downscale transfer (towards higher wavenumbers) of total energy. Potential energy fluxes $\Pi_{A}$ are also uniformly positive, indicative of downscale transfers, with some indications of inertial ranges (with fluxes independent of wavenumber) in some cases.

6.2.1 Spectral energy fluxes: $\Omega^{*}=1$

For the Earth equivalent simulation at $\Omega^{*}=1$ with T127 resolution and normal friction (Fig. 7(e)), the total energy flux $\Pi$ (black solid line) rises sharply at wavenumbers $n=2$ and 3 to a value of $\sim 1.6$ W m*-2*, stays roughly constant until $n=7$ before falling rapidly between wavenumbers 8 and 12 and then decreasing more slowly towards zero at the highest wavenumbers. $\Pi$ consists of two main components, of which the APE component, $\Pi_{A}$, dominates over the KE component, $\Pi_{K}$, up to a wavenumber of 50. Because of its larger magnitude, the trend of the APE flux $\Pi_{A}$ is similar to that of $\Pi$, except that its slope within the $\sim$constant region between $n=3$ and 8 is less steep. This difference between $\Pi$ and $\Pi_{A}$ is the result of an upscale energy transfer, $\Pi_{K}$, of KE between $n=3$ and $\sim$10-12. At around $n=11$ there is an inflection point in the KE spectrum where $\Pi_{K}$ changes sign. This implies that kinetic energy is being transported towards smaller scales (i.e. larger wavenumbers) for $n\gtrsim 10-12$ and towards larger scales for $n\lesssim 10$. In this region of the spatial spectrum, the baroclinic conversion, $\mathcal{C}$, has a steeply descending slope with wavenumber. This is a cumulative term (c.f. Eqn.  43) and so such a strong negative slope in $\mathcal{C}$ denotes a conversion of APE to KE in this wavenumber range of magnitude C_{n}=\mathcal{C}(n=15)-\mathcal{C}(n=7)\simeq 0.9\>\text{W m{}^{-2}}.

Regarding the partitioning between eddy-eddy and residual zonal components, the zonal components evidently dominate $\Pi_{K}$ and $\Pi_{A}$ at smaller wavenumbers ($1\leq n\lesssim 20$) in Fig. 7(e), while eddy-eddy components gain in relative importance at higher wavenumbers ($n>20$), although the total fluxes there are relatively low. On the other hand, the main component of the conversion term occurs in the eddy-eddy component. Taking all of these points together the Earth-like case is evidently consistent with the defining behaviour of idealized baroclinic turbulence (see e.g. Vallis, 2006). At the injection wavenumber for KE (around the Rossby deformation radius wavenumber $n_{D}\sim 8-12$; see Fig. 8(e)), APE is converted into KE via the eddy-eddy component of $C_{n}$ (which is related to the baroclinic $C_{E}$ component of the Lorenz energy budget; see Fig. 3). The resulting KE is transported mostly upscale into the zonal component by an inverse barotropic conversion (cf $C_{K}$ in Lorenz budget), with a smaller amount of KE being transported downscale where the eddy-eddy interaction component dominates (c.f. Fig 8(e)).

At smaller wavenumbers, we see a zonal component also contributes to $\mathcal{C}$. The entire cumulative sum of $\mathcal{C}$ (i.e. the value depicted at wavenumber $n=1$) is comparable in sign and magnitude to $C_{Z}$ in the Lorenz energy budget (Fig. 3).This conversion shows that $C_{Z}$ is negative for wavenumbers $n=4-7$ and positive at the smallest wavenumbers. It is likely that the zonal-zonal components (associated with the Eulerian mean Hadley and Ferrell circulations) dominate at low wavenumbers, while zonal-eddy components are more significant at higher wavenumbers. The thermally-direct Hadley circulation, which dominates at low latitudes, leads to positive $C_{Z}$, which may account for the behaviour of $C_{n}$ for $n\leq 4$. But the Ferrell cell is generally thermally indirect, so would be expected to make a negative contribution to $C_{Z}$, as apparent for $n=4-7$.

Segments of the spectrum where the spectral fluxes $\Pi_{K}$, $\Pi_{A}$ or $\Pi$ itself are approximately constant are identified as inertial ranges, and two such regions can be discerned in this case. Firstly, the region between wavenumbers $n=3$ and 8, where $\Pi_{A}$ is constant (and $\Pi$ and $\Pi_{K}$ are almost constant) may describe an inertial range characterised by forward baroclinic APE transfers and an inverse (rotational) barotropic KE flux. The second region lies at $n\simeq 30-80$ (see Fig 8e for a close-up of $\Pi_{K}$) with both forward APE and KE cascades. This second wavenumber region can be identified in the energy spectrum with an $n^{-3}$ slope in Fig. 5e. The first region, however, is not so easily identifiable with features in the spectrum. Further comparison with Fig. 5e, however, confirms an association of the $n^{-5}$ slope in zonal energies with a downscale flux of both KE and APE. The narrow inertial range around $3\leq n\leq 8$ in KE occurs mainly in the zonal-eddy component, which means that most of the energy jumps directly between the zonal mean wind and a range of energy-significant non-axisymmetric wavenumbers (more like a “waterfall” than a “cascade”?). For the other inertial range for $30\leq n\leq 80$, however, the transfers are more dominated by eddy-eddy interactions. This means that the latter inertial range involves not only scales at which no dissipation occurs while cascading, but also scales where only weak interactions between the zonally symmetric flow and the respective eddy-scales occur.

6.2.2 Spectral energy fluxes: rapidly rotating cases ($\Omega^{*}>1$)

With increasing rotation rate (Figs. 7f-h) the Rossby deformation radius decreases (and $n_{D}$ increases) so that the wavenumber at which most baroclinic conversion occurs increases commensurately. At $\Omega^{*}\geq 2$, the KE inertial range at large wavenumbers, seen for $\Omega^{*}=1$ (see Fig. 8(e)) can no longer be discerned because it starts to close in on the resolution cut-off in these simulations at wavenumber $n=170$ where it becomes affected by the hyperdiffusion. However, the inertial range in APE flux widens and flattens with increasing rotation rate in a region that corresponds to a positive slope in the energy spectra (cf Figs. 5(f)-(h)).

Figure 8(e)-(h) shows the spectral kinetic energy flux for simulations with $\Omega^{*}=1,2,4,8$ in detail. As mentioned above, for $\Omega^{*}=1$ (Fig. 8(e)) we can see that the energy injected at the Rossby deformation lengthscale through conversion from APE by baroclinic instability is transferred both upscale and downscale (indicated by the red solid line). The upscale component can be identified with an upscale barotropic transfer of KE by the rotational part of the flow, which is dominated by the zonal interaction components. The downscale component at higher wavenumbers, however, is dominated by the divergent eddy-eddy interactions. With increasing rotation rate (Figs. 8(f)-(h)), however, in contrast to the $\Omega^{*}=1$ case, the divergent mode decreases sharply in magnitude so that, at the highest rotation rates, only the rotational part of the flux transfers energy in either direction. In addition, the contribution of the eddy-eddy interaction terms at larger wavenumbers becomes stronger. At the highest values of $\Omega^{*}$, therefore, the macroturbulent interactions are almost entirely dominated by the rotational flow, with the divergent eddies playing little role.

6.2.3 Spectral energy fluxes: slowly rotating cases ($\Omega^{*}<1$)

Figure 7(a)-(d) shows the spectral energy fluxes for slowly-rotating simulations ($\Omega^{*}=\frac{1}{16},\frac{1}{8},\frac{1}{4},\frac{1}{2}$). With decreasing rotation rate, the baroclinically active region (i.e. with downscale $\Pi_{A}$, and negative slope in $\mathcal{C}$), identified in the previous section, moves towards smaller wavenumbers. Between $\Omega^{\ast}=frac{1}{8}$ and $\frac{1}{4}$, however, this baroclinically dominated behaviour is suppressed, giving way to a quite different pattern of fluxes at the lowest values of $\Omega^{*}$. This trend is consistent with that found by Mitchell and Vallis (2010), who observed that their super-rotating simulated circulations, unlike Earth-like cases, were not dominated by baroclinic zonal-eddy interactions, as indicated by a lack of divergence of the vertical component of the EP-fluxes (c.f. Mitchell and Vallis, 2010, their Fig. 7). In addition, the zonal components of $\mathcal{C}$, which were comparatively small at higher values of $\Omega^{*}$, now begin to dominate at all lengthscales. This occurs because, at smaller rotation rates (larger values of $\mathcal{R}o_{T}\gtrsim 10$), the Rossby deformation lengthscale exceeds the planetary radius and APE is then injected directly into the KE reservoir at very low wavenumbers directly via interactions with the zonal mean flow. $\mathcal{C}_{zonal}$ at $n=1$ is again very similar to $C_{Z}$ in the corresponding Lorenz budget (cf Fig. 3), which points towards a strong influence of zonal-zonal interactions in this conversion term.

At $\Omega^{\ast}=\frac{1}{16}$ and $\frac{1}{8}$, the qualitative structure of the fluxes is therefore entirely different to the more quasi-geostrophic cases at higher $\Omega^{\ast}$. Conversion from APE to KE now occurs at the smallest wavenumbers, principally via zonal interactions. In addition both $\Pi_{K}$ and $\Pi_{A}$ now feature a well developed inertial range in the form of a forward transfer with an approximately constant spectral flux between wavenumbers $n=6$ and 30. This is indicative of a forward barotropic “waterfall”. In both cases the zonal-eddy interactions dominate. However, the influence of eddy-eddy interactions is still evident and still increases in magnitude at larger wavenumbers.

Figure 8(a)-(d) again features the kinetic energy flux in detail. In the case of decreasing $\Omega^{\ast}$, it is the rotational component that diminishes and the divergent component of the flux that controls the forward energy cascade. This suggests a much greater role for gravity and equatorial inertia-gravity planetary waves as these do not possess a rotational component. This would not be unduly surprising given that the equatorial waveguide grows in width at low rotation rates to span much of the planet.

The behaviour identified in this section fits well with other results obtained for the large thermal Rossby number regime ($\mathcal{R}o_{T}>>1$). For $\Omega^{\ast}=\frac{1}{8}$ and $\frac{1}{16}$, the baroclinic conversion becomes weak and barotropic effects become stronger (also apparent in the corresponding Lorenz energy budgets; Fig. 3), such that $\Pi_{K}>\Pi_{A}$ in this regime. The flow becomes largely zonal and super-rotating flow emerges. Unfortunately, this analysis does not help directly in identifying the mechanism of formation and maintenance of the equatorial superrotation as this occurs mostly in the zonal component in a specific region of the globe, whereas this analysis computes over a global mean and focusses on the non-zonal spherical wavenumber spectrum. What we can learn, however, is that the kinetic energy in the zonal mode of super-rotating cases dissipates via a downscale cascade that involves both zonal-eddy and eddy-eddy interactions, with the latter dominating at high wavenumbers.

7 Discussion

This study has explored how the dynamical transfers of energy and vorticity between different horizontal scales depend upon the planetary rotation rate, at least as represented in a highly simplified, but nevertheless fully nonlinear and generic, numerical circulation model of a prototypical terrestrial planetary atmosphere. Such explorations are important sources of insight into the factors that determine the form, structure and intensity of atmospheric circulations under various conditions, thereby helping us to understand and quantify the similarities and differences between different planets of our own Solar System (and beyond), as well as indicating how aspects of any atmospheric circulation will scale with key planetary parameters.

7.1 Lorenz energy budgets

Heat transport by both eddies and zonally-symmetric meridional overturning form important contributions to the overall energy budget of an atmosphere. This is commonly analysed using the framework originally developed by Lorenz (1967) and is still used as a source of insight for understanding the atmospheres of Earth and other planets (e.g. Peixóto and Oort, 1974; James, 1995; Schubert and Mitchell, 2014; Tabataba-Vakili et al., 2015).

In the present work, we have computed how the various terms in the Lorenz energy budget for a simple, dry, Earth-like atmosphere vary with $\Omega^{\ast}$. Although the magnitudes of the zonal mean energy reservoirs vary monotonically with $\Omega^{\ast}$, with increasing dominance of APE over KE as $\Omega^{\ast}$ increases, the eddy energies rise to a maximum around the value of $\Omega^{\ast}$ where $\mathcal{R}o_{T}\sim 1$. This is also reflected in most of the conversion rates, which also peak in magnitude around a value of $\mathcal{R}o_{T}$ between 1 and 0.1. The trends in energy conversion rate also demonstrate the change in character of the dominant eddy generation processes from mainly barotropic processes at low rotation rates towards predominantly baroclinic processes at more rapid rotation rates, consistent with the onset of strong and deep baroclinic instabilities when $\mathcal{R}o_{T}\lesssim 1$. At much higher rotation rates, however, even baroclinic instability becomes less effective at energy conversion as the Rossby deformation radius becomes much smaller than the planetary radius, leading to a decrease in the intensity of the whole Lorenz energy cycle as $\Omega^{\ast}\rightarrow\infty$.

This tendency of the Lorenz energy cycle to peak in intensity around conditions not too far from the Earth has been noted before, e.g. by Pascale et al. (2013). In their study, this was associated with a maximum in entropy production rates, although the precise conditions were found to depend not just on rotation rate but also on the strength of dissipation in the system. We have not sought to explore this in detail in the present study, but this would be of interest to investigate further in future work.

7.2 Spectral energy budgets

The results shown in Section 6 present for the first time a reasonably comprehensive overview of how the pattern of spectral fluxes of enstrophy and various forms of energy changes between different planetary circulation regimes. The simulation span a broad range of parameter space, extending from an extreme quasi-geostrophic limit through to a highly ageostrophic, super-rotating regime at very low rotation rates, over which the pattern of enstrophy and energy cascades changes significantly.

Despite the use of a highly simplified global circulation model, the results for Earth-like conditions capture a circulation regime with a pattern of enstrophy and energy cascades that compares reasonably well with results from much more realistic models (e.g. Burgess et al., 2013; Augier and Lindborg, 2013; Malardel and Wedi, 2016), at least qualitatively. The enstrophy fluxes at $\Omega^{\ast}=1$ indicate a predominantly forward cascade over most wavenumbers with a flux that increases towards high wavenumbers in the baroclinically active troposphere. The magnitude of the enstrophy flux in the PUMA simulations is generally smaller than found e.g. by Burgess et al. (2013) in their reanalysis data by a factor of $\sim 5$, but this likely reflects differences in the way the models are energised as well as effects of finite spatial resolution. Energy fluxes at $\Omega^{\ast}=1$ are broadly comparable to those found by Augier and Lindborg (2013) in their analyses of AFES and ECMWF simulations, though again somewhat smaller in magnitude. Vertically integrated, upscale rotational KE fluxes are around half the magnitude of those in both AFES and ECMWF simulations, while the forward cascade for $n>20$, which is dominated in all simulations by divergent components, is weaker in the PUMA simulations by a factor $\sim 5$. Available potential and total energy spectral fluxes, however, were quite comparable in magnitude in the PUMA Earth-like simulation to the NWP models at low-moderate wavenumber, though follows the AFES model more closely at high wavenumbers with positive (downscale) fluxes down to the resolution limit. This is consistent with the results of Malardel and Wedi (2016), who also found spectral fluxes to be significantly weaker in their cases with Held-Suarez (linear relaxation) forcing, suggesting that this approach underestimates the realistic energetic forcing of the simulated circulation.

Resolution is likely to be a limiting factor for various features in the circulation. The KE and APE spectra for $\Omega^{\ast}=1$ exhibit some features in common with the Earth’s spectra (e.g. Nastrom and Gage, 1985; Burgess et al., 2013) in following an $n^{-3}$ trend over most of the spectrum for $n\lesssim 100$ in both KE and APE. With normal levels of surface drag, however, there is little evidence in the vertically-averaged spectrum from the PUMA simulations for the mesoscale break in the KE spectrum towards $n^{-5/3}$ around $n\gtrsim 20$. Dissipation may also play an important role in these simulations in removing energy at scales similar to where baroclinic energy conversion is taking place. These and similar effects from sub-gridscale parameterizations were noted in the study by Malardel and Wedi (2016) and it would be of significant interest to explore this further under conditions that are significantly different from Earth.

As $\Omega^{\ast}$ is increased, the results show a gradual transition from the Earth-like pattern of spectral fluxes, in which both rotation and divergent KE components contribute to the cascades, towards a more rotationally dominated KE cascade at all scales. The magnitude of such fluxes quickly become quite weak as $\Omega^{\ast}$ is increased, probably due to relatively strong bottom friction. Nevertheless, the upscale segment at relatively low wavenumbers is dominated by the rotational flow, but at the highest rotation rates (and smaller values of $Ro_{T}$) the forward KE cascade also becomes dominated by rotational components. Such a pattern resembles more closely the distribution of spectral fluxes found in both the Earth’s oceans (Scott and Wang, 2005; Scott and Arbic, 2007) and in Jupiter’s atmosphere (Young and Read, 2017). Limited spatial resolution probably restricts the ability of the simulated flows to develop fully inertial ranges, so the KE and APE spectra are only marginally consistent with expectations of observing clear enstrophy and KE-dominated cascades. But the results are broadly consistent in this regime with some of the predictions of classical geostrophic turbulence theory, as summarised schematically in Figure 9(a) (though with some modifications discussed further below). These results include uniformly downscale transfers of APE and total energy, excitation of the KE spectrum around the deformation scale $n_{D}$ in association with baroclinic instabilities and barotropization, and near-equipartition between the APE and KE spectra at the highest values of $\Omega^{\ast}$. Under Earth-like conditions, however, some of these classical predictions are not borne out, in particular because of the significant role of divergent components of KE and because the deformation scale is not sufficiently well separated from the planetary scale.

7.3 Zonal jet formation

The results shown in Sections 5 and 6 also examine the jet formation mechanism in terms of KE spectra, with particular attention to the paradigm of “zonostrophic turbulence” recently proposed by Galperin and coworkers as a potential candidate for a universal regime for jet formation in various geophysical fluids, including planetary atmospheres (e.g. Sukoriansky et al., 2002; Galperin et al., 2006, 2010). The experiments presented here demonstrate that, provided surface friction is not too strong, the atmosphere does develop strongly coherent zonal jets with a highly anisotropic KE spectrum that shares at least some features in common with the idealised zonostrophic turbulence regime (e.g. with $-5$ and $-5/3$ slopes respectively in the zonal and eddy KE spectra; Sukoriansky et al. (2002); Galperin et al. (2006, 2010)). However, with relatively stronger surface friction, zonal jets appear in a weaker and more meandering/wavy form. This is consistent with the “barotropic governor” mechanism (e.g. James, 1995) which indicates that weak frictional damping leads to stronger barotropic shear in the atmosphere, thereby suppressing the growth of baroclinic instability, and making the circulation equilibrate into a more zonally symmetric state. This also evidently results in the accumulation of KE in predominantly barotropic zonal flows, mainly through non-local upscale transfers of KE directly from eddies into zonally symmetric flow components. This is an aspect of the transfer of energy between scales that was not considered in the early work of Charney (1971) or Salmon (1980), but is indicated in Fig. 9(a) to make the point that the upscale cascade is more complicated and anisotropic than initially considered.

7.4 An Inertio-stratified regime

At values of $\Omega^{\ast}<1$, both the spectra and pattern of energy transfers undergo a significant change of regime, taking place around $Ro_{T}\simeq 1$. The resulting regime at the highest values of $Ro_{T}$, illustrated schematically in Fig. 9(b), has an entirely different character to any of the quasi-geostrophic sub-regimes. In this regime (which might be termed inertio-stratified turbulence), the system is energized by differential heating at the planetary scale in zonally symmetric modes, which immediately begin cascading energy and enstrophy uniformly towards smaller scales following a direct conversion of zonally symmetric APE into KE. For $n\gtrsim 4$, the overall flow develops a clear inertial range that even our low spatial resolution simulations are able to represent quite well, whereby APE and KE cascade uniformly towards small scales with a near-constant spectral flux. In this regime, divergent KE components dominate the KE spectral flux, which in turn dominates over the APE flux. Both the KE and APE spectra also adopt a “classical” KBK form with a clear $n^{-5/3}$ slope until the resolution limit is approached. KE dominates the total energy spectrum, but the ratio of KE to APE appears to tend towards a fixed value $\simeq 2-3$. This is clearly distinct from any of the quasi-geostrophic regimes, but the precise details of which wave modes govern the properties of the cascade in both the horizontal and vertical directions remain to be explored. This regime exhibits a number of similarities to the mesoscale and sub-mesoscale regimes in the Earth’s atmosphere and oceans, recently identified by Callies et al. (2014) and Callies and Ferrari (2013) respectively, in which relatively fast inertia-gravity waves take over the role of Rossby waves in classical QG turbulence. The KE and APE spectra, however, do not conform very closely to what we find in our model spectra. As discussed also by Wang et al. (2018), these very slowly rotating circulations are dominated by strongly super-rotating zonal flows. Like the rapidly rotating, quasi-geostrophic regimes, therefore, they are likely to be characterized by highly anisotropic spectra and energy transfers which should be explored in more detail in future work.

Finally, together with the results presented on analyses of NWP models by Augier and Lindborg (2013), our results demonstrate that spectral fluxes of energy and enstrophy provide a very clear and insightful approach to diagnosing the performance of numerical models. Even the limited results shown so far indicate some significant differences between different model formulations, which might be expected to lead to some helpful advances in model design, based on sound physical principles.

\ack

PLR, FT-V, AV and RMBY acknowledge support from the UK Science and Technology Facilities Council during the course of this research under grants ST/K502236/1, ST/I001948/1 and ST/K00106X/1. The collaboration with Drs Pierre Augier and Erik Lindborg was made possible through the International Network on Waves and Turbulence, funded by the Leverhulme Trust. We are also grateful to Dhruv Balwada and an anonymous referee for their constructive comments on an earlier version of this paper.