Mathematics

The energy budget equations are as follows:

\begin{align*} \frac{\partial A_Z}{\partial t} &= -C_A - C_Z + BA_Z + \Delta G_Z \\ \frac{\partial A_E}{\partial t} &= C_A - C_E + BA_E + \Delta G_E \\ \frac{\partial K_Z}{\partial t} &= C_K + C_Z + BK_Z + \Delta R_Z \\ \frac{\partial K_E}{\partial t} &= C_E - C_K + BK_E + \Delta R_E \end{align*}

In these equations, available potential energy (APE) is divided into zonal (\(A_Z\)) and eddy (\(A_E\)) components, as is kinetic energy (\(K_Z\) and \(K_E\), respectively). The transformations between these forms of energy are denoted by \(C\), with subscripts \(Z\) and \(E\) for conversions between zonal and eddy forms, and \(A\) and \(K\) indicating conversions between APE and kinetic energy, respectively. Thus, \(C_A\) represents the conversion between \(A_Z\) and \(A_E\), \(C_E\) denotes the conversion from \(A_E\) to \(K_E\), \(C_K\) signifies the transformation from \(K_E\) to \(K_Z\), and \(C_Z\) describes the conversion from \(A_Z\) to \(K_Z\). The residual terms are defined as:

\begin{align*} \Delta R_Z &= B \Phi_Z + D_Z + \epsilon_{KZ} \\ \Delta R_E &= B \Phi_E + D_E + \epsilon_{KE} \\ \Delta G_Z &= G_Z + \epsilon_{GZ} \\ \Delta G_E &= G_E + \epsilon_{GE} \end{align*}

Where APE generation and dissipation of kinetic energy are indicated by \(G\) and \(D\), with \(G_Z\) and \(G_E\) marking the generation of \(A_Z\) and \(A_E\), and \(D_Z\) and \(D_E\) representing the dissipation of \(K_Z\) and \(K_E\), respectively.

Note

Sign convention. These budgets follow Brennan and Vincent (1980, Eqs. 5a-6b). All conversion terms are signed rates, and a positive value means transfer in the direction named by the term: positive \(C_A\) transfers \(A_Z \rightarrow A_E\), positive \(C_Z\) transfers \(A_Z \rightarrow K_Z\), positive \(C_E\) transfers \(A_E \rightarrow K_E\), and positive \(C_K\) transfers \(K_E \rightarrow K_Z\). \(D_Z\) and \(D_E\) are likewise signed frictional terms (negative for a loss), not positive-definite loss magnitudes; note that Michaelides (1987) writes them with the opposite sign convention, as \(-D\) in the kinetic budgets.

Residuals. \(\Delta R_Z\) and \(\Delta R_E\) are what the code exports as RKz and RKe. They are composite residuals: each one lumps together the boundary pressure work, the frictional dissipation, any unresolved-scale transfer, and the accumulated numerical error. They must not be interpreted as friction alone. Likewise RGz and RGe (\(\Delta G_Z\), \(\Delta G_E\)) are generation plus numerical error. This matches the residual definitions of Dias Pinto and da Rocha (2011, their Eqs. 8-11).

The code implements exactly these budgets in src/utils/calc_budget_and_residual.py, which rearranges them to

\[\begin{split}\begin{align*} RG_Z &= \frac{\partial A_Z}{\partial t} + C_Z + C_A - BA_Z \\ RK_Z &= \frac{\partial K_Z}{\partial t} - C_Z - C_K - BK_Z \\ RG_E &= \frac{\partial A_E}{\partial t} - C_A + C_E - BA_E \\ RK_E &= \frac{\partial K_E}{\partial t} - C_E + C_K - BK_E \end{align*}\end{split}\]

Firstly, we define the zonal mean of a variable \(X\), between longitudes \(\lambda_{1}\) and \(\lambda_{2}\):

\[\begin{align*} [X]_\lambda &= \frac{1}{\lambda_2 - \lambda_1} \int_{\lambda_2}^{\lambda_1} X d\lambda \end{align*}\]

The eddy component of this variable is its deviation from the zonal mean:

\[\begin{align*} (X)_\lambda &= X - [X]_\lambda \end{align*}\]

The domain mean of the variable \(X\), defined over the computational domain bounded by longitudes \(\lambda_1\) and \(\lambda_2\), and latitudes \(\varphi_1\) and \(\varphi_2\), is given by:

\[\begin{align*} [X]_{\lambda\phi} &= \left(\frac{1}{\lambda_2 - \lambda_1}\right) \left(\frac{1}{\sin\phi_2 - \sin\phi_1}\right) \int_{\lambda_2}^{\lambda_1} X \cos\phi d\lambda d\phi \end{align*}\]

Similarly, we define the deviation of the zonal mean from the domain mean:

\[\begin{align*} ([X]_\lambda)_\phi &= [X]_\lambda - [X]_{\lambda\phi} \end{align*}\]

Note

The averages above are evaluated on the model grid with the trapezoidal rule, and are normalised by that same rule rather than by the analytic \(\lambda_2-\lambda_1\) and \(\sin\phi_2-\sin\phi_1\). The two differ by \(O(\Delta\phi^2)\), about \(1.6\times10^{-4}\) on a 2.5 degree grid. Matching them makes the average of a constant exact, so the departures \(([X]_\lambda)_\phi\) and \((X)_\lambda\) satisfy their defining identity to roundoff. This matters for the geopotential, whose zero point is arbitrary: without it the diagnosed boundary fluxes would depend on the reference level chosen for \(\Phi\).

From the definitions above, the four energy components used in the LEC computation are defined as follows:

\[\begin{split}\begin{align*} A_Z &= \int_{p_t}^{p_b} \frac{([(T)_\lambda ])_\phi^{2}]_{\lambda \phi}} {2[\sigma]_{\lambda \phi}} dp \\ A_E &= \int_{p_t}^{p_b} \frac{[(T)_\lambda^{2}]_{\lambda \phi}]} {2[\sigma]_{\lambda \phi}} dp \\ K_Z &= \int_{p_t}^{p_b} \frac{[[u]_\lambda^2 + [v]_\lambda^2]_{\lambda \phi}}{2g} dp \\ K_E &= \int_{p_t}^{p_b} \frac{[(u)_\lambda^2 + (v)_\lambda^2]_{\lambda \phi}}{2g} dp \end{align*}\end{split}\]

where \(p\) is the atmospheric pressure, with subscripts \(b\) and \(t\) denoting the lower (base) and upper (top) pressure boundaries of the atmosphere, respectively. \(T\) represents temperature, \(g\) is the acceleration due to gravity, and \(u\) and \(v\) are the zonal and meridional wind components, respectively. The static stability parameter \(\sigma\) is defined as:

\[\begin{align*} \sigma &= \left[\frac{gT}{c_p}-\frac{pg}{R}\frac{\partial T}{\partial p}\right]_{\lambda \phi} \end{align*}\]

where \(c_p\) is the specific heat at constant pressure, and \(R\) is the ideal gas constant for dry air.

The four conversion terms are defined as follows, integrating over the atmospheric column from the base (\(p_b\)) to the top (\(p_t\)) pressures:

\[\begin{split}\begin{aligned} &C_Z = \int_{p_t}^{p_b} - [\left([T]_\lambda)_\phi ([\omega]_\lambda\right)_\phi]_{\lambda\phi} \ \frac{R}{gp} \ dp \label{eq:CZ} \\ &C_E = \int_{p_t}^{p_b} - [(T)_\lambda (\omega)_\lambda]_{\lambda\phi} \ \frac{R}{gp} \ dp \label{eq:CE} \\ &C_{\mathrm{overturning}} = -\int_{p_t}^{p_b} [\omega]_{\lambda\phi}\,\alpha\,\frac{dp}{g}, \qquad \alpha=\frac{R[T]_{\lambda\phi}}{p} \label{eq:Coverlap} \\ &C_A = \int_{p_t}^{p_b} - \left( \frac{1}{a\sigma} \left[ (v)_\lambda (T)_\lambda \frac{\partial ([T]_\lambda)_\phi}{\partial \phi} \right]_{\lambda\phi} + \frac{1}{\sigma} \left[ (\omega)_\lambda (T)_\lambda \frac{\partial ([T]_\lambda)_\phi}{\partial p} \right]_{\lambda\phi} \right) dp \label{eq:CA} \\ &C_K = \int_{p_t}^{p_b} \frac{1}{g} \left(\left[ \frac{\cos\phi}{a} (u)_\lambda (v)_\lambda \frac{\partial}{\partial\phi} \left(\frac{[u]_\lambda}{\cos\phi}\right)\right]_{\lambda\phi} + \left[ \frac{(v)_\lambda^2}{a} \frac{\partial [v]_\lambda}{\partial\phi} \right]_{\lambda\phi} \right. \left. + \left[ \frac{\tan\phi}{a} (u)_\lambda^2 [v]_\lambda \right]_{\lambda\phi} + \left[ (\omega)_\lambda (u)_\lambda \frac{\partial [u]_\lambda}{\partial p} \right]_{\lambda\phi} + \left[ (\omega)_\lambda (v)_\lambda \frac{\partial [v]_\lambda}{\partial p} \right]_{\lambda\phi} \right) dp \label{eq:CK} \end{aligned}\end{split}\]

where \(a\) is the Earth’s radius and \(\omega\) is the vertical velocity in isobaric coordinates. The separately exported \(C_{\mathrm{overturning}}\) diagnoses the strength of domain-mean overturning. It is not a missing conversion in the Lorenz-cycle budget: in the exact pressure-work identity it cancels the \(\overline{\Phi}\,\overline{\omega}\) part of the top/bottom geopotential flux. It therefore remains outside \(RG_Z\) and \(RK_Z\). Under the sign convention above, domain-mean ascent (\(\omega<0\)) produces a positive value. It vanishes for the global cycle, where mass continuity requires the horizontal mean \([\omega]_{\lambda\phi}=0\), but need not vanish in a limited area. It belongs neither to \(C_Z\), which uses the area anomaly \(\omega^*\), nor to \(C_E\), which uses the zonal eddy \(\omega'\).

Mass-continuity diagnostic

The toolkit also exports the pressure-coordinate mass-continuity residual

\[M(p)=\left\langle\frac{u|_{\lambda_e}-u|_{\lambda_w}} {a\cos\varphi\,\Delta\lambda}\right\rangle +\left\langle\frac{1}{a\cos\varphi} \frac{\partial([v]\cos\varphi)}{\partial\varphi}\right\rangle +\frac{\partial\overline{\omega}}{\partial p}.\]

The two horizontal terms are evaluated in their telescoped boundary forms and normalised by the same trapezoidal measure \(\int\cos\varphi\,d\varphi\) used by \(\langle\cdot\rangle\) elsewhere, so all three summands of \(M\) share one area-averaging convention; otherwise the mismatch between them would itself be \(O(\Delta\varphi^2)\), the same order as the inconsistency \(M\) exists to measure.

The per-level file M_<pressure-coordinate>.csv contains \(M(p)\) in \(\mathrm{s^{-1}}\); the M column in the main result is \(\int M\,dp/g\) in \(\mathrm{kg\,m^{-2}\,s^{-1}}\). The continuum value is zero. A nonzero value measures the mismatch between the archived horizontal winds and archived \(\omega\) and sets the numerical noise floor for geopotential-flux diagnostics that use full \(\Phi\).

The APE generation and K dissipation terms are defined as:

\[\begin{split}\begin{align*} G_Z &= \int_{p_t}^{p_b} \frac{[([q]_\lambda)_\phi ([T]_\lambda)_\phi]_{\lambda \phi}}{c_p[\sigma]_{\lambda \phi}} dp \\ G_E &= \int_{p_t}^{p_b} \frac{[(q)_\lambda (T)_\lambda]_{\lambda \phi}}{c_p[\sigma]_{\lambda \phi}} dp \\ D_Z &= - \int_{p_t}^{p_b} \frac{1}{g} [[u]_\lambda [F_\lambda]_\lambda + [v]_\lambda [F_\phi]_\lambda]_{\lambda \phi} dp \\ D_E &= - \int_{p_t}^{p_b} \frac{1}{g} [(u)_\lambda (F_\lambda)_\lambda + (v)_\lambda (F_\phi)_\lambda]_{\lambda \phi} dp \end{align*}\end{split}\]

Here, \(F_{\lambda}\) and \(F_{\varphi}\) represent the zonal and meridional frictional components, respectively, and \(q\) is the diabatic heating term, computed as a residual from the thermodynamic equation:

\[\begin{align*} \frac{q}{c_p} &= \frac{\partial T}{\partial t} + \vec{V}_H \cdot \vec{\nabla}_p T - S_p\omega \end{align*}\]

where \(\vec{V}_H \cdot \vec{\nabla}_p T\) represents the horizontal advection of temperature and \(S_p\) approximates the static stability, given by:

\[\begin{align*} S_p &\equiv -\frac{T}{\theta}\frac{\partial \theta}{\partial p} \end{align*}\]

where \(\theta\) is the potential temperature.

The boundary terms are given by:

\[\begin{split}\begin{aligned} & \mathrm{BAZ}=c_1 \int_{p_1}^{p_2} \int_{\varphi_1}^{\varphi_2} \frac{1}{2[\sigma]_{\lambda_{\varphi}}}\left(2\left([T]_\lambda\right)_{\varphi}(T)_\lambda u+\left([T]_{\lambda_{\varphi}}\right)_{\varphi}^2 u\right)_{\lambda_1}^{\lambda_2} \nonumber \\ & \times d \varphi d p+c_2 \int_{p_1}^{p_2} \frac{1}{2[\sigma]_{\lambda \varphi}}\left(2\left[(v)_\lambda(T)_\lambda\right]_\lambda\left([T]_\lambda\right)_{\varphi} \cos \varphi \right. \left.+\left([T]_\lambda\right)_{\varphi}^2[v]_\lambda \cos \varphi\right)_{\varphi_1}^{\varphi_2} d p \nonumber \\ & -\frac{1}{2[\sigma]_{\lambda \varphi}}\left(\left[2(\omega)_\lambda(T)_\lambda\right]_\lambda\left([T]_\lambda\right)_{\varphi}+\left[[\omega]_\lambda\left([T]_\lambda\right)_{\varphi}^2\right]_{\lambda_{\varphi}}\right)_{p_1}^{p_2} \\ & \mathrm{BAE}=c_1 \int_{p_1}^{p_2} \int_{\varphi_1}^{\varphi_2} \frac{1}{2[\sigma]_{\lambda \varphi}}\left[u(T)_\lambda^2\right]_{\lambda_1}^{\lambda_2} d \varphi d p \nonumber \\ & +c_2 \int_{p_1}^{p_2} \frac{1}{2[\sigma]_{\lambda \varphi}}\left(\left[(T)_\lambda^2 v\right]_\lambda \cos \varphi\right)_{\varphi_1}^{^{\varphi_2}} d p \\ & -\left(\frac{\left[\omega(T)_\lambda^2\right]_{\lambda \varphi}}{2[\sigma]_{\lambda \varphi}}\right)_{p_1}^{p_2} \nonumber \\ & \mathrm{BKZ}=c_1 \int_{p_1}^{p_2} \int_{\varphi_1}^{\varphi_2} \frac{1}{2 g}\left(u\left[u^2+v^2-(u)_\lambda^2-(v)_\lambda^2\right]\right)_{\lambda_1}^{\lambda_2} \nonumber \\ & \times d \varphi d p+c_2 \int_{p_1}^{p_2} \frac{1}{2 g}\left(\left[v \cos \varphi \left[u^2+v^2\right.\right.\right. \left.\left.-(u)_\lambda^2-(v)_\lambda^2\right]\right]_{\varphi_1}^{\varphi_2} d p \\ & -\left(\frac{1}{2 g}\left[\omega\left[u^2+v^2-(u)_\lambda^2-(v)_\lambda^2\right]\right]_{\lambda \varphi}\right)_{p_1}^{p_2} \nonumber \\ & \mathrm{BKE}=c_1 \int_{p_1}^{p_2} \int_{\varphi_1}^{\varphi_2} \frac{1}{2 g}\left(u\left[(u)_\lambda^2+(v)_\lambda^2\right]\right)_{\lambda_1}^{\lambda_2} d \varphi d p \nonumber \\ & +c_2 \int_{p_1}^{p_2} \frac{1}{2 g}\left(\left[v \cos \varphi\left[(u)_\lambda^2+(v)_\lambda^2\right]\right]_\lambda\right)_{\varphi_1}^{\varphi_2} d p \\ & -\left(\frac{1}{2 g}\left[\omega\left[(u)_\lambda^2+(v)_\lambda^2\right]\right]_{\lambda \varphi}\right)_{p_1}^{p_2} \nonumber \end{aligned}\end{split}\]

where \(c_1=-\left[a\left(\lambda_2-\lambda_1\right)\left(\sin \varphi_2-\sin \varphi_1\right)\right]^{-1}, c_2=-\left[a\left(\sin \varphi_2-\sin \varphi_1\right)\right]^{-1}\).

Lastly, the terms \(B\Phi_Z\) and \(B\Phi_E\) are given by:

\[\begin{split}\begin{aligned} \mathrm{B}\Phi_Z= &\ c_1 \int_{p_t}^{p_b}\!\!\int_{\varphi_s}^{\varphi_n} \frac{1}{g}\left([u]_\lambda\,\Delta_{EW}\Phi' + \Phi^*\,\Delta_{EW}u'\right) d\varphi\, dp \\ &+c_2 \int_{p_t}^{p_b} \frac{1}{g} \left(\cos\varphi\,[v]_\lambda\,\Phi^*\right)_{\varphi_s}^{\varphi_n} dp -\frac{1}{g}\left(\left[\omega^*\Phi^*\right]_{\lambda\varphi}\right)_{p_t}^{p_b} \\ \mathrm{~B} \Phi \mathrm{E}= & c_1 \int_{p_1}^{p_2} \int_{\varphi_1}^{\varphi_2} \frac{1}{g}\left((u)_\lambda(\Phi)_{\lambda_\lambda}\right)_{\lambda_1}^{\lambda_2} d \varphi d p \nonumber \\ & +c_2 \int_{p_1}^{p_2} \frac{1}{g}\left(\left[(v)_\lambda(\Phi)_{\lambda_\lambda}\right]_\lambda \cos \varphi\right)_{\varphi_1}^{\varphi_2} d p \\ & -\frac{1}{g}\left(\left[(\omega)_\lambda(\Phi)_\lambda\right]_{\lambda_{\varphi}}\right)_{p_1}^{p_2} \nonumber \end{aligned}\end{split}\]

Here \(\Phi^*=[\Phi]_\lambda-[\Phi]_{\lambda\varphi}\), \(\omega^*=[\omega]_\lambda-[\omega]_{\lambda\varphi}\), \(\Phi'=\Phi-[\Phi]_\lambda\), \(u'=u-[u]_\lambda\) and \(\Delta_{EW}X=X|_{\lambda_2}-X|_{\lambda_1}\).

The east and west faces follow the structure of Brennan and Vincent (1980, pp. 964-965), who write that wall as the flux \(u\Phi-u'\Phi'\) evaluated from west to east; expanding the east-minus-west difference gives \([u]_\lambda\Delta_{EW}\Phi'\) together with a second contribution carrying \(\Delta_{EW}u'\). That second contribution is taken here with the area departure \(\Phi^*\) rather than the zonal mean \([\Phi]_\lambda\): a limited domain carries a net mass flux through its walls, so with the full geopotential the term would change if an arbitrary constant were added to \(\Phi\), whereas \(\Phi^*\) leaves it unchanged. Both contributions vanish on a periodic (global) domain, where the two walls coincide.

Michaelides (1987, p. 25) writes this wall as \(c_1\iint([v]_\lambda\Phi^*)|_{\lambda_1}^{\lambda_2}d\varphi\,dp/g\). Both factors are zonal means and therefore independent of longitude, which makes that expression identically zero.

The north/south and top/bottom faces follow Muench (1965) and Michaelides (1987), using the area departures \(\Phi^*\) and \(\omega^*\). Brennan and Vincent (1980) write those two faces with the full geopotential instead; in a limited domain that form additionally carries the domain-mean overturning conversion proportional to \(\int\overline{\omega}\,\alpha\,dp/g\), which is an interior conversion rather than a boundary flux and is diagnosed separately as \(C_{\mathrm{overturning}}\).