Energy Budget

Definitions¶

There are four types of energy (per unit mass) that we need to keep track of in the fluid:

$I = c_v T$ is the internal energy (valid for an ideal gas, see e.g. Vallis (2017)) -- this is a measure of the disorganized molecular-scale motion in the fluid
$\Phi = g z$ is the gravitational potential energy
$L = L_v q$ is the latent energy associated with the evaporated water in the air ( $L_v$ is the latent heat of vaporization and $q$ is the specific humidity)
$K = \frac{1}{2} \left(u^2 + v^2 + w^2 \right)$ is the kinetic energy of organized motion of the fluid

We’ll define the total energy per unit mass in the atmosphere as

\begin{align*} E &= I + \Phi + L + K \\ &= c_v T + gz + L_v q + \frac{1}{2} \left(u^2 + v^2 + w^2 \right) \end{align*}

(1)

Budget equations for energy components¶

Potential energy¶

Latent energy¶

Internal energy¶

Kinetic energy¶

A total energy equation and global conservation¶

After some work and cancellation of various conversion terms:

\frac{DE}{Dt} = -\alpha \nabla \cdot \left( \vec{v} p \right) - \alpha \nabla \cdot \vec{F}_{rad} + D_h

(2)

with $D_h = c_p D_T + L_v D_q$ the net surface enthalpy flux and $\vec{F}_{rad}$ the radiative flux.

Alternate form and Moist Static Energy¶

We can rewrite the first law of thermodynamics as

\begin{align*} dQ &= dI + p d\alpha \\ &= dI + d(p\alpha) - \alpha dp \\ &= d(c_v T + RT) - \alpha dp \\ &= d(c_p T ) - \alpha dp \end{align*}

(3)

where $c_p T$ is refered to as the sensible energy $S$

Using this, we can define another form of total energy

\mathcal{E} = S + \Phi + L + K

(4)

and we’ll find it useful to define the moist static energy $m$ as the first three terms of $\mathcal{E}$ :

m = S + \Phi + L = c_p T + g z + L_v q

(5)

We’ll see from the data that $K$ is typically a small fraction of $\mathcal{E}$ so is often ignored in the total energy budget, and we will talk instead about fluxes of moist static energy rather than fluxes of total energy.

Merdional energy transport¶

Write the zonal- and time-averaged energy budget in pressure coordinates and flux form:

\frac{\partial}{\partial t} [\overline{\mathcal{E}}] = -\frac{1}{a \cos\phi} \frac{\partial}{\partial \phi} \left( [\overline{v\mathcal{E}}] \cos\phi \right) -\frac{\partial}{\partial p} \left( [\overline{\omega \mathcal{E}}] \right) + g \left[ \frac{\partial \overline{\vec{F}_{rad}}}{\partial p} \right] + [\overline{D_h}]

(6)

and we can decompose the meridional energy flux into its components:

[\overline{v\mathcal{E}}] = c_p [\overline{vT}] + g [\overline{vz}] + L_v [\overline{vq}] + [\overline{vK}]

(7)

Finally we can decompose each component into mean meridional, stationary eddy, and transient eddy terms, e.g.

[\overline{vT}] = [\overline{v}][\overline{T}] + [\overline{v^*} ~\overline{T^*}] + [\overline{v^\prime T^\prime}]

(8)

and likewise for the other components.

The observed energy budget¶

We are going to make some plots of northward energy fluxes from the ERA5 data.

import numpy as np
import xarray as xr
import dask
import matplotlib.pyplot as plt
from arraylake import Client as alClient
from dask.distributed import Client as dClient

alclient = alClient()
alclient.login()

repo = alclient.get_repo("earthmover-public/era5")
session = repo.readonly_session("main")
pressure_temporal = xr.open_zarr(session.store, group="pressure/temporal", chunks={})
pressure_temporal

  2026-10-28T22:23:53.513790Z  WARN aws_smithy_http_client::client::tls::rustls_provider::build_connector: native root CA certificate loading errors: [Error { context: "failed to read PEM from file", kind: Io { inner: Os { code: 13, kind: PermissionDenied, message: "Permission denied" }, path: "/etc/ssl/certs/sendmail.pem" } }, Error { context: "failed to read PEM from file", kind: Io { inner: Os { code: 13, kind: PermissionDenied, message: "Permission denied" }, path: "/etc/pki/tls/certs/sendmail.pem" } }]
    at /home/conda/feedstock_root/build_artifacts/bld/rattler-build_icechunk_1785432413/build_env/.cargo/registry/src/index.crates.io-1949cf8c6b5b557f/aws-smithy-http-client-1.1.12/src/client/tls/rustls_provider.rs:73

# ## Option 1: run this on a laptop
# ##  Comment out this cell if using option 2
# dclient = dClient()
# dclient

## Option 2: If running this notebook on jupyterlab.its.albany.edu, uncomment and run this cell
##  This sets up a cluster of workers to parallelize the calculations

from dask_jobqueue import SLURMCluster

cluster = SLURMCluster(cores=64, # size of a single job -- typically one node of the HPC cluster
                       processes=8, # 8 workers x 8 threads
                       memory="374GB",  # the max memory on the 64-core burst-daes nodes, I believe
                       walltime="12:00:00",
                       queue="burst-daes",
                      )

cluster.scale(4)
dclient = dClient(cluster)
dclient

2026-10-28 22:54:24,567 - tornado.application - ERROR - Uncaught exception GET /status/ws (10.16.2.49)
HTTPServerRequest(protocol='http', host='169.226.65.130:8787', method='GET', uri='/status/ws', version='HTTP/1.1', remote_ip='10.16.2.49')
Traceback (most recent call last):
  File "/network/rit/lab/snowclus/pixi/aug26/.pixi/envs/default/lib/python3.14/site-packages/tornado/websocket.py", line 965, in _accept_connection
    open_result = handler.open(*handler.open_args, **handler.open_kwargs)
  File "/network/rit/lab/snowclus/pixi/aug26/.pixi/envs/default/lib/python3.14/site-packages/tornado/web.py", line 3409, in wrapper
    return method(self, *args, **kwargs)
  File "/network/rit/lab/snowclus/pixi/aug26/.pixi/envs/default/lib/python3.14/site-packages/bokeh/server/views/ws.py", line 149, in open
    raise ProtocolError("Token is expired.")
bokeh.protocol.exceptions.ProtocolError: Token is expired.

#  We will work with data for one year for demonstration purposes
interval = "2021"

# a subset of the full data
d = pressure_temporal[["u", "v", "t", "q", "z"]].sel(valid_time=interval)

The flux of sensible energy $vT$ ¶

%%time
# Lazy seasonal means over all four seasons
vbar = d.v.groupby("valid_time.season").mean("valid_time")
Tbar = d.t.groupby("valid_time.season").mean("valid_time")
vT_bar = (d.v * d.t).groupby("valid_time.season").mean("valid_time")

CPU times: user 2.22 s, sys: 115 ms, total: 2.33 s
Wall time: 2.38 s

%%time
# Single pass over the data -- the expensive step
vbar, Tbar, vT_bar = dask.compute(vbar, Tbar, vT_bar)

/network/rit/lab/snowclus/pixi/aug26/.pixi/envs/default/lib/python3.14/site-packages/distributed/client.py:3429: UserWarning: Sending large graph of size 62.99 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

CPU times: user 21min 2s, sys: 1min 17s, total: 22min 19s
Wall time: 32min 43s

def compute_flux_components(vbar, Abar, vA_bar, Aref=0.):
    '''For northward velocity vbar and scalar field Abar, 
    compute the components of the northward flux decomposition.'''
    
    # Zonal means and stationary-eddy deviations
    vbar_z = vbar.mean("longitude")
    Abar_z = Abar.mean("longitude")
    vstar = vbar - vbar_z
    Astar = Abar - Abar_z
    
    vA_total   = vA_bar.mean("longitude") - vbar_z * Aref
    vA_mmc   = vbar_z * (Abar_z - Aref)
    vA_stat  = (vstar * Astar).mean("longitude")
    vA_trans = (vA_bar - vbar * Abar).mean("longitude")
    
    # Pack results into a dictionary
    flux = {}
    flux['total'] = vA_total
    flux['transient'] = vA_trans
    flux['stationary'] = vA_stat
    flux['MMC'] = vA_mmc
    return flux

sensible_flux = compute_flux_components(vbar, Tbar, vT_bar, Aref=273.15)

sensible_labels = {'total':      r'$[\overline{vT}]$', 
          'transient':  r'$[\overline{v^\prime T^\prime}]$',
          'stationary': r'$[\overline{v}^* \overline{T}^*]$',
          'MMC':        r'$[\overline{v}][\overline{T}]$',
         }

lat = pressure_temporal.latitude

# Define a mask so our plot excludes the Antarctic continent and high ice sheet
antarctic_mask = (lat < -70.) & (pressure_temporal.pressure_level>700.)

def cross_section_plot(fluxdict, labeldict, title, scaling_factor=1.):
    fig, axes = plt.subplots(4,1, figsize=(8,12), layout="constrained")

    for num, flux_name in enumerate(fluxdict):
        flux = fluxdict[flux_name].mean(dim="season") * scaling_factor
        flux_masked = flux.where(~antarctic_mask)
        CS = flux_masked.plot.contour(ax=axes[num],
                                yincrease=False, 
                              )
        axes[num].clabel(CS)
        axes[num].set_title("{}, {}".format(flux_name, labeldict[flux_name]))
    fig.suptitle(title)
    return fig, axes

fig, axes = cross_section_plot(sensible_flux, sensible_labels, 
                               title="Northward sensible heat flux for {} (ºC m/s)".format(interval)
                              )

Plot the mass-weighted vertical averages by season¶

def column_mean(data):
    return data.integrate(coord="pressure_level")/(-1000.)

def vertical_mean_plot(fluxdict, labeldict, title, scaling_factor=1.):
    fig, axes = plt.subplots(4,1, figsize=(6,8), layout="constrained")

    for num, flux_name in enumerate(fluxdict):
        flux = fluxdict[flux_name] * scaling_factor
        flux_masked = flux.where(~antarctic_mask)
        flux_integrated = column_mean(flux_masked)
        ax = axes[num]
        ax.plot(lat, flux_integrated.mean(dim='season'), linewidth=4, label='Annual')
        for s in ["DJF", "JJA"]: 
            ax.plot(lat, flux_integrated.sel(season=s), linewidth=1, label=s)
        ax.set_title("{}, {}".format(flux_name, labeldict[flux_name]))
        ax.grid()
        ax.legend()
        ax.set_xlim(-70,90)
    fig.suptitle(title)
    return fig, axes

fig, axes = vertical_mean_plot(sensible_flux, sensible_labels,
                   title="Vertical mean northward sensible heat flux for {} (ºC m/s)".format(interval)
                  )

Latent energy fluxes¶

%%time
# Lazy seasonal means over all four seasons
qbar = d.q.groupby("valid_time.season").mean("valid_time")
vq_bar = (d.v * d.q).groupby("valid_time.season").mean("valid_time")

CPU times: user 793 ms, sys: 6.77 ms, total: 800 ms
Wall time: 796 ms

%%time
# Single pass over the data -- the expensive step
qbar, vq_bar = dask.compute(qbar, vq_bar)

/network/rit/lab/snowclus/pixi/aug26/.pixi/envs/default/lib/python3.14/site-packages/distributed/client.py:3429: UserWarning: Sending large graph of size 53.62 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

CPU times: user 16min 30s, sys: 53.9 s, total: 17min 24s
Wall time: 27min 53s

latent_flux = compute_flux_components(vbar, qbar, vq_bar)

latent_labels = {'total':      r'$[\overline{vq}]$', 
          'transient':  r'$[\overline{v^\prime q^\prime}]$',
          'stationary': r'$[\overline{v}^* \overline{q}^*]$',
          'MMC':        r'$[\overline{v}][\overline{q}]$',
         }

fig, axes = cross_section_plot(latent_flux, latent_labels, 
                               scaling_factor = 1000,  # convert to (g/kg)(m/s)
                               title="Northward water vapor flux for {} (g/kg m/s)".format(interval)
                              )

fig, axes = vertical_mean_plot(latent_flux, latent_labels,
                               scaling_factor=1000., # convert to (g/kg)(m/s)
                               title="Vertical mean northward water vapor flux for {} (g/kg m/s)".format(interval)
                              )

Potential energy fluxes¶

%%time
# Lazy seasonal means over all four seasons
Zbar = d.z.groupby("valid_time.season").mean("valid_time")
vZ_bar = (d.v * d.z).groupby("valid_time.season").mean("valid_time")

CPU times: user 822 ms, sys: 16.8 ms, total: 839 ms
Wall time: 833 ms

%%time
# Single pass over the data -- the expensive step
Zbar, vZ_bar, = dask.compute(Zbar, vZ_bar)

/network/rit/lab/snowclus/pixi/aug26/.pixi/envs/default/lib/python3.14/site-packages/distributed/client.py:3429: UserWarning: Sending large graph of size 53.62 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

CPU times: user 16min 28s, sys: 53.9 s, total: 17min 22s
Wall time: 27min 15s

potential_flux = compute_flux_components(vbar, Zbar, vZ_bar)

potential_labels = {'total':      r'$[\overline{vZ}]$', 
          'transient':  r'$[\overline{v^\prime Z^\prime}]$',
          'stationary': r'$[\overline{v}^* \overline{Z}^*]$',
          'MMC':        r'$[\overline{v}][\overline{Z}]$',
         }

g = 9.8  # m/s2

fig, axes = vertical_mean_plot(potential_flux, potential_labels,
                               scaling_factor= 1./g/100., # convert to 10^2 m/s gpm
                               title="Vertical mean northward potential energy flux for {} (100 gpm m/s)".format(interval)
                              )

Moist static energy transport¶

Now compute the total flux

[\overline{vm}] = c_p [\overline{vT}] + g [\overline{vZ}] + L_v [\overline{vq}]

(9)

where we have defined the moist static energy

m \equiv c_p T + g Z + L_v q

(10)

cp = 1004.     # specific heat at constant pressure for dry air (J / kg / K)
Lhvap = 2.5E6    # Latent heat of vaporization (J / kg)

def compute_mse(T,gz,q):
    return cp*T + gz + Lhvap*q

mse_labels = {'total':      r'$[\overline{vm}]$', 
          'transient':  r'$[\overline{v^\prime m^\prime}]$',
          'stationary': r'$[\overline{v}^* \overline{m}^*]$',
          'MMC':        r'$[\overline{v}][\overline{m}]$',
         }

mse_flux = {}
for name in mse_labels:
    mse_flux[name] = compute_mse(sensible_flux[name], 
                                 potential_flux[name],
                                 latent_flux[name])

fig, axes = cross_section_plot(mse_flux, mse_labels, 
                               scaling_factor = 1/cp, # convert to ºC m/s
                               title="Northward moist static energy flux for {} (J/kg m/s)".format(interval)
                              )

fig, axes = vertical_mean_plot(mse_flux, mse_labels,
                               scaling_factor = 1/cp, # convert to ºC m/s
                               title="Vertical mean northward moist static energy flux for {} (ºC m/s)".format(interval)
                              )

Contribution of kinetic energy to total flux¶

Is the term $[\overline{vK}]$ actually small enough to neglect compared to the other terms in the energy flux

[\overline{v\mathcal{E}}] = c_p [\overline{vT}] + g [\overline{vz}] + L_v [\overline{vq}] + [\overline{vK}]

(11)

In other words, how good is the moist static energy flux approximation?

Let’s actually compute and plot $[\overline{vK}]$

Kinetic energy formula¶

We have the definition

K = \frac{1}{2} \left(u^2 + v^2 + w^2 \right)

(12)

Since the vertical component $w$ is typically several orders of magnitude smaller than $u$ and $v$ , we neglect it here. Thus we will compute the flux

v K \approx \frac{v}{2} \left(u^2 + v^2 \right)

(13)

along with its averages and decompositions.

Computing kinetic energy flux from ERA5¶

def kinetic_energy(u,v):
    return (u**2 + v**2)/2

%%time
K = kinetic_energy(d.u, d.v)
Kbar = K.groupby("valid_time.season").mean("valid_time")
vK_bar = (d.v * K).groupby("valid_time.season").mean("valid_time")

CPU times: user 4.19 s, sys: 29.9 ms, total: 4.22 s
Wall time: 4.21 s

%%time
# Single pass over the data -- the expensive step
Kbar, vK_bar, = dask.compute(Kbar, vK_bar)

/network/rit/lab/snowclus/pixi/aug26/.pixi/envs/default/lib/python3.14/site-packages/distributed/client.py:3429: UserWarning: Sending large graph of size 53.62 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

CPU times: user 18min 58s, sys: 1min 3s, total: 20min 1s
Wall time: 29min 41s

kinetic_flux = compute_flux_components(vbar, Kbar, vK_bar)

kinetic_labels = {'total':      r'$[\overline{vK}]$', 
          'transient':  r'$[\overline{v^\prime K^\prime}]$',
          'stationary': r'$[\overline{v}^* \overline{K}^*]$',
          'MMC':        r'$[\overline{v}][\overline{K}]$',
         }

fig, axes = cross_section_plot(kinetic_flux, kinetic_labels, 
                               title="Northward kinetic energy flux for {} (m^3 / s^3)".format(interval)
                              )

fig, axes = vertical_mean_plot(kinetic_flux, kinetic_labels,
                               scaling_factor = 1/cp, # convert to ºC m/s
                               title="Vertical mean northward kinetic energy flux for {} (ºC m/s)".format(interval)
                              )

These are plotted in the same units as the moist static energy fluxes above for comparison. The MSE flux is of order 10 ºC m/s, while the kinetic energy flux is of order 0.1 ºC m/s.

The neglecting of the kinetic energy contribution to the total energy flux seems to be well justified.

Total energy transport¶

Here we will express the northward flux of energy in terms of the total integrated northward energy transport $\mathcal{H}$ as we defined in an earlier chapter.

For completeness, since we have already computed it, we will include the kinetic energy component.

a = 6.373E6      # Radius of Earth (m)
ps = 1E5       # approximate surface pressure (Pa)

scale_factor = 2*np.pi*a*np.cos(np.deg2rad(lat))*ps/g / 1E15 # transform flux in J/kg m/s to integrated transport in PW

total_energy_transport = {}
for name in mse_flux:
    total_energy_transport[name] = (mse_flux[name] + kinetic_flux[name]) * scale_factor
    
total_energy_labels = {'total':      r'$[\overline{v\mathcal{E}}]$', 
          'transient':  r'$[\overline{v^\prime \mathcal{E}^\prime}]$',
          'stationary': r'$[\overline{v}^* \overline{\mathcal{E}}^*]$',
          'MMC':        r'$[\overline{v}][\overline{\mathcal{E}}]$',
         }

fig, axes = vertical_mean_plot(total_energy_transport, total_energy_labels,
                               title="Vertical mean northward total energy transport for {} (PW)".format(interval)
                              )

cluster.close()

References¶

Vallis, G. K. (2017). Atmospheric and Oceanic Fluid Dynamics: fundamentals and large-scale circulation (2nd ed.). Cambridge University Press.

Definitions¶

Budget equations for energy components¶

Potential energy¶

Latent energy¶

Internal energy¶

Kinetic energy¶

A total energy equation and global conservation¶

Alternate form and Moist Static Energy¶

Merdional energy transport¶

The observed energy budget¶

The flux of sensible energy vTvTvT¶

Plot the mass-weighted vertical averages by season¶

Latent energy fluxes¶

Potential energy fluxes¶

Moist static energy transport¶

Contribution of kinetic energy to total flux¶

Kinetic energy formula¶

Computing kinetic energy flux from ERA5¶

Total energy transport¶

The flux of sensible energy $vT$ ¶