Ctrl+K

Decompositions of the circulation

4. Decompositions of the circulation#

Introduction#

In the previous chapter, we looked at time averaged quantities.

It’s often going to be useful to consider this total flow as consisting of a mean and a fluctuating part.

Time mean and zonal mean#

We consider two means that are often used: time mean and zonal mean, along with deviations from these means.

We define the time mean as

\[ \overline{A} = \frac{1}{T} \int_0^T A dt \]

and the zonal mean as

\[ \left[ A \right] = \frac{1}{2\pi} \int_0^{2\pi} A d\lambda \]

where \(t\) is time, \(T\) is the total time interval over which we are averaging, \(\lambda\) is longitude, and \(A\) is whatever variable we are averaging.

This introduces a standardized notation, where the overbar will always indicate time averages and the square brackets will be reserved for zonal averages.

Example#

Prove that the zonal and time mean mass flux across any latitude circle must vanish, for time scales greater than a year.

Let’s start with the continuity (mass) equation in pressure coordinates:

\[ \nabla \cdot \vec{u} = 0 \]

which we can write out in spherical coordinates like so:

\[ \frac{1}{a \cos\phi}\frac{\partial u }{\partial \lambda} + \frac{1}{a \cos\phi} \frac{\partial \left( v \cos\phi \right)}{\partial \phi} + \frac{\partial \omega}{\partial p} \]

So if we take the time and zonal average of the continuity equation we get

\[ \left[ \frac{1}{a \cos\phi}\frac{\partial \overline{u} }{\partial \lambda} \right] + \left[ \frac{1}{a \cos\phi} \frac{\partial \left( \overline{v} \cos\phi \right)}{\partial \phi} \right] + \left[ \frac{\partial \overline{\omega}}{\partial p} \right] \]

which we can also write as

\[ \frac{1}{2 \pi a \cos\phi}\int_0^{2\pi} \frac{\partial \overline{u} }{\partial \lambda} d\lambda + \frac{1}{a \cos\phi} \frac{\partial }{\partial \phi} \left( \left[ \overline{v} \right] \cos\phi \right) + \frac{\partial \left[ \overline{\omega} \right]}{\partial p} \]

The first term is zero. Why?

Let’s now integrate over the depth of the atmosphere, from surface pressure \(p_s\) to zero:

\[ \frac{1}{a \cos\phi} \frac{\partial }{\partial \phi} \left( \cos\phi \int_0^{p_s} \left[\overline{v}\right] dp \right) + \int_0^{p_s} \frac{\partial \left[\overline{\omega}\right]}{\partial p} dp \]

The second term is the integral of a derivative so resolves to only boundary contributions. At the TOA \(\omega(0) = 0\), and so the term simplifies to just

\[\left[\overline{\omega}\right] |_{p_s}\]

We cannot have flow in or out of the surface. The boundary condition is complicated by the fact that the surface pressure is not constant in time, but over timescales of a year or more we do not expect any time-averaged changes in surface pressure. Therefore

\[\left[\overline{\omega}\right] |_{p_s} \approx 0\]

and the average mass budget reduces to

\[ \frac{\partial}{\partial \phi} \left( \cos\phi \int_0^{p_s} \left[\overline{v}\right]\right) = 0 \]

so

\[ \cos\phi \int_0^{p_s} \left[\overline{v}\right] = C \]

for some constant \(C\)

The integral above is proportional to the next mass flux across latitude circles. So what is the constant \(C\)?

One way to rationalize this is to consider the polar limit where \(\phi \rightarrow \frac{\pi}{2}\) and thus \(\cos\phi \rightarrow 0\). This would then require that the integral becomes unboundedly large if \(C\ne 0\).

We conclude that \(C=0\) for all \(\phi\), and thus the zonal and time mean mass flxu across any latitude circle must vanish.

Departures from the mean#

Let’s introduce some more standardized notation for departures from our two means

Departure from time mean#

\[ A^\prime \equiv A - \overline{A} \]

Departure from zonal mean#

\[ A^* \equiv A - \left[ A \right] \]

Some identities#

Note that

\[ \overline{A^\prime} = 0 \]

and

\[ \left[ A^\prime \right] = 0 \]

as long as the averaging is the same as that used to find the perturbations.

Also the averages permute:

\[\begin{align*} \left[ \overline{A} \right] &= \overline{ \left[ A \right] } & \overline{ \overline{A} } &= \overline{A} & \left[ \left[A \right] \right] \end{align*}\]

Combining zonal and time expansions#

We can decompose the full field into four terms:

\[ A = \overline{A} + A^\prime = \left[\overline{A}\right] + \overline{A}^* + \left[ A^\prime \right] + A^{\prime *} \]

where our four terms are

  1. Zonally symmetric part of the steady time-averaged quantities, e.g. trade winds and midlatitude westerlies \(\left[ \overline{u} \right]\)

  2. Asymmetric part of the time-averaged quantities, e.g. monsoon circulations, PNA, NAO (patterns which vary on only low frequencies)

  3. Instantaneous fluctutations of symmetric part of quantities, e.g. time fluctuations in Hadley cells, mean jet locations

  4. Instantaneous, zonally asymmetric part of quantities, e.g. transient cyclones and anticyclones

Examples using some real data#

We are going to load in some zonal wind data \(u\) from the CFSR dataset.

For simplicity we’ll just work with a few months of data here from a single year (only need to open one data file), and focus on a single pressure level (500 hPa, mid-troposphere).

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs

cfsr_path = '/cfsr/data/'
year = '2021'

ds = xr.open_dataset(cfsr_path + year + '/u.' + year + '.0p5.anl.nc')
ds

<xarray.Dataset>
Dimensions:  (time: 1460, lat: 361, lon: 720, lev: 32)
Coordinates:
  * time     (time) datetime64[ns] 2021-01-01 ... 2021-12-31T18:00:00
  * lat      (lat) float32 -90.0 -89.5 -89.0 -88.5 -88.0 ... 88.5 89.0 89.5 90.0
  * lon      (lon) float32 -180.0 -179.5 -179.0 -178.5 ... 178.5 179.0 179.5
  * lev      (lev) float32 1e+03 975.0 950.0 925.0 900.0 ... 50.0 30.0 20.0 10.0
Data variables:
    u        (time, lev, lat, lon) float32 ...
Attributes:
    description:    u 1000-10 hPa
    year:           2021
    source:         http://nomads.ncdc.noaa.gov/data.php?name=access#CFSR-data
    references:     Saha, et. al., (2010)
    created_by:     User: ab473731
    creation_date:  Sat Jan  2 06:00:32 UTC 2021
# instantaneous u on one day
u500 = ds.u.sel(lev=500.)
u500.sel(time='2021-01-01T00')

<xarray.DataArray 'u' (lat: 361, lon: 720)>
[259920 values with dtype=float32]
Coordinates:
    time     datetime64[ns] 2021-01-01
  * lat      (lat) float32 -90.0 -89.5 -89.0 -88.5 -88.0 ... 88.5 89.0 89.5 90.0
  * lon      (lon) float32 -180.0 -179.5 -179.0 -178.5 ... 178.5 179.0 179.5
    lev      float32 500.0
Attributes:
    level_type:  isobaric_level (hPa)
    units:       m s^-1
    long_name:   zonal wind

Let’s make a helper function to draw maps of lat-lon data:

def make_map(field, levels=np.linspace(-75,75,31), title=None):
    p = field.plot.contourf(levels=levels,
                            subplot_kws={'projection': ccrs.Robinson()},
                                    transform=ccrs.PlateCarree())
    p.axes.coastlines()
    if title is not None:
        plt.title(title)
    return p

And make a map of the instantaneous \(u\) on a particular day and time:

make_map(u500.sel(time='2021-01-01T00'));
../_images/f07c6e0fba1f1fb5eecdce87bb9e4b0d33a6690b296449b7d89906bd057e2979.png

To do some seasonal averaging, we need to pick a time interval. Let’s choose all dates in January, February and March.

This is easy to do by slicing on the time axis:

# A time slice containing all of January, February and March
u500_JFM = u500.sel(time=slice('2021-01-01','2021-03-31'))
u500_JFM

<xarray.DataArray 'u' (time: 360, lat: 361, lon: 720)>
[93571200 values with dtype=float32]
Coordinates:
  * time     (time) datetime64[ns] 2021-01-01 ... 2021-03-31T18:00:00
  * lat      (lat) float32 -90.0 -89.5 -89.0 -88.5 -88.0 ... 88.5 89.0 89.5 90.0
  * lon      (lon) float32 -180.0 -179.5 -179.0 -178.5 ... 178.5 179.0 179.5
    lev      float32 500.0
Attributes:
    level_type:  isobaric_level (hPa)
    units:       m s^-1
    long_name:   zonal wind

This array has 360 points on the time axis, corresponding to 31+28+31 = 90 days @ 4 times daily.

Let’s make a map of \(\overline{u}\) by taking a time average:

ubar = u500_JFM.mean(dim='time')
m = make_map(ubar, title='$\overline{u}$ at 500 hPa for January - March 2021')
../_images/01b837c3f9ddb9c9fa8160e2d2698b5efc074bc5bb65d95b5f1ce6733823e54d.png

If we further take the zonal average, we get \([\overline{u}]\)

Since this quantity is one dimensional (latitude only), we’ll make a line plot instead of a map:

ubar.mean(dim=('lon')).plot(y='lat')
plt.title('$[\overline{u}]$ at 500 hPa for January - March 2021')
plt.grid();
../_images/5e9aa6d35507cc3bf229227686f82c2d7aed3eadcca2d0b0c8370a836b55c5ca.png

Now here’s \(\overline{u^*}\), the time averaged departures from zonal mean:

ustar = u500_JFM - u500_JFM.mean(dim='lon')
ustarbar = ustar.mean(dim='time')
m = make_map(ustarbar, title='$\overline{u^*}$ at 500 hPa for January - March 2021')
../_images/8819b281d6c6514f427548eaacc99fb778860a48eefd556783a4f2ee36681ead.png

And here’s \(u^\prime\) (departures from time average) at one instant:

uprime = u500_JFM - ubar
make_map(uprime.sel(time='2021-01-01T00'));
../_images/46131ebe88f8e742004e3b1c30b11a9234c56be832e033abe83be44c682d60d2.png

and its zonal average \([u^\prime]\) as a line plot:

uprime.sel(time='2021-01-01T00').mean(dim='lon').plot(y='lat')

[<matplotlib.lines.Line2D at 0x7f9d79249360>]
../_images/3b6f8f9861d2e125fbbece447bf1d4c4cb3272262abe29c139b67540aac5e66b.png

And finally the deviations away from that zonal average \(u^{\prime*}\)

make_map((uprime - uprime.mean(dim='lon')).sel(time='2021-01-01T00'))

<cartopy.mpl.contour.GeoContourSet at 0x7f9d7909f610>
../_images/d013b58958744988b17e419a0b6acaff41f319139aafd6c150355f3c6291ee7c.png

Covariances between quantities#

We are often interested in covariances between two quantities, usually a component of the flow and another scalar (e.g. meridional wind and temperature)

\[\begin{align*} \overline{AB} &= \overline{\left( \overline{A} + A^\prime \right) \left(\overline{B} + B^\prime \right)} \\ &= \overline{\overline{A}~\overline{B}} + \overline{\overline{A} B^\prime} + \overline{A^\prime \overline{B}} + \overline{A^\prime B^\prime} \\ &= \overline{A}~\overline{B} + \overline{A} \overline{B^\prime} + \overline{B} \overline{A^\prime} + \overline{A^\prime B^\prime} \\ &= \overline{A}~\overline{B} + \overline{A^\prime B^\prime} \end{align*}\]

The last term is the covariance of \(A\) and \(B\) in time.

Similarly for the zonal average

\[\left[ AB \right] = \left[A\right] \left[B\right] + \left[ A^* B^* \right] \]

where the last term is the covariance of \(A\) and \(B\) is longitude.

We can also combine the time and zonal means

\[\begin{align*} \left[\overline{AB}\right] &= \left[ \overline{A}~\overline{B} + \overline{A^\prime B^\prime} \right] \\ &= \left[\overline{A}\right]\left[\overline{B}\right] + \left[ \overline{A}^* \overline{B}^* \right] + \left[\overline{A^\prime B^\prime} \right] \end{align*}\]

Three terms:

  1. Mean meridional component

  2. Stationary eddy component (zonal variations in time mean

  3. Transient eddy component (time variations in time mean)

Example: the meridional heat flux \(vT\) at 500 hPa#

Let \(A = v\) and \(B = T\)

and let’s consider the decomposition of the heat flux \(vT\)

Actually we are going to subtract the global average \(\left< T \right>\) from the temperature data first so that the numbers are easier to interpret!

v_ds = xr.open_dataset(cfsr_path + year + '/v.' + year + '.0p5.anl.nc')
T_ds = xr.open_dataset(cfsr_path + year + '/t.' + year + '.0p5.anl.nc')

v500_JFM = v_ds.v.sel(lev=500, time=slice('2021-01-01','2021-03-31'))
T500_JFM = T_ds.t.sel(lev=500, time=slice('2021-01-01','2021-03-31'))

T500_JFM

<xarray.DataArray 't' (time: 360, lat: 361, lon: 720)>
[93571200 values with dtype=float32]
Coordinates:
  * time     (time) datetime64[ns] 2021-01-01 ... 2021-03-31T18:00:00
  * lat      (lat) float32 -90.0 -89.5 -89.0 -88.5 -88.0 ... 88.5 89.0 89.5 90.0
  * lon      (lon) float32 -180.0 -179.5 -179.0 -178.5 ... 178.5 179.0 179.5
    lev      float32 500.0
Attributes:
    level_type:  isobaric_level (hPa)
    units:       K
    long_name:   temperature
# Compute the area-weighted global mean on lat-lon grid
T500_JFM_global = T500_JFM.mean(dim='time').weighted(np.cos(np.deg2rad(T500_JFM.lat))).mean(dim=('lat','lon'))
T500_JFM_global

<xarray.DataArray 't' ()>
array(258.15087999)
Coordinates:
    lev      float32 500.0

First we can compute the full field \([\overline{vT}]\):

vTbar = (v500_JFM * (T500_JFM - T500_JFM_global)).mean(dim='time')
levs = np.linspace(-180,180, 19)
make_map(vTbar, levels=levs, title='$\overline{vT}$ at 500 hPa for January - March 2021')

<cartopy.mpl.contour.GeoContourSet at 0x7f9d78e89420>
../_images/ec8f6e8aca0dcdd762a08554d63e4df296929a212530a34731ee83f16d6fff4a.png

How much of the full field is captured by the time mean advection of the time mean temperature?

(Note, the time means here include any eddies that are stationary on our averaging timescale, which is seasonal here)

vbar = v500_JFM.mean(dim='time')
Tbar = T500_JFM.mean(dim='time') - T500_JFM_global

make_map(vbar * Tbar, levels=levs, title='$\overline{v}~\overline{T}$ at 500 hPa for January - March 2021')

<cartopy.mpl.contour.GeoContourSet at 0x7f9d7986bd60>
../_images/0716c7cf4e7e1d26e28578b8804755349216df961ae1a4442e9177a275ab85ff.png

And now let’s look at the transient eddy contribution

vprime = v500_JFM - vbar
Tprime = (T500_JFM - T500_JFM_global) - Tbar
vprime_Tprime_bar = (vprime*Tprime).mean(dim='time')
make_map(vprime_Tprime_bar, title='$\overline{v^\prime T^\prime}$ at 500 hPa for January - March 2021')

<cartopy.mpl.contour.GeoContourSet at 0x7f9d79567bb0>
../_images/65fa33cbcacbcbe7b61b6fe46687cfc52d3994f24b87bfb2c78bca479a824f28.png

Now look at contributions to the zonal average#

Consider the decomposition

\[ [\overline{vT}] = [\overline{v}][\overline{T}] + [\overline{v}^* \overline{T}^*] + [\overline{v^\prime T^\prime}] \]
vTbar_zon = vTbar.mean(dim='lon')
vbar_zon = vbar.mean(dim='lon')
Tbar_zon = Tbar.mean(dim='lon')
vbarzon_Tbarzon = vbar_zon * Tbar_zon
vbar_star = vbar - vbar_zon
Tbar_star = Tbar - Tbar_zon
vbarstar_Tbarstar_zon = (vbar_star * Tbar_star).mean(dim='lon')
vprime_Tprime_bar_zon = vprime_Tprime_bar.mean(dim='lon')

vTbar_zon.plot(y='lat', label='Total, $[\overline{vT}]$')
vbarzon_Tbarzon.plot(y='lat', label='Mean meridional, $[\overline{v}][\overline{T}]$')
vbarstar_Tbarstar_zon.plot(y='lat', label='Stationary eddy, $[\overline{v}^* \overline{T}^*]$')
vprime_Tprime_bar_zon.plot(y='lat', label='Transient eddy, $[\overline{v^\prime T^\prime}]$')
plt.legend()
plt.grid()
plt.title('Zonal mean heat flux decomposition at 500 hPa for January - March 2021');
../_images/3dd5db04ece12fdfe4a78629efbd3ff60e2c34dd0bbb56e6d38e9fda2b39b0cb.png

What does it mean to have a positive “stationary eddy” component?#

The plot above shows that there are maximum in poleward heat flux near 50º latitude in both hemispheres. While the transient eddy term \([\overline{v^\prime T^\prime}]\) (red) is large in both hemispheres, there is a substantial positive (northward) stationary eddy \([\overline{v}^* \overline{T}^*]\) contribution in the Northern Hemisphere.

What is this, and what does it mean?

Map of zonal temperature anomalies \(\overline{T}^*\)#

Let’s take a quick look at the map of \(\overline{T}^*\) to get a sense of its contribution to \([\overline{v}^* \overline{T}^*]\)

make_map(Tbar_star, levels=np.linspace(-10,10,21), title='$\overline{T}^*$ at 500 hPa for January - March 2021');
../_images/ca4b901af67f85edc9a8e28f1f1f014c044a1d2a8f70e8653a0bf5bcbff40245.png

Focussing on the Northern mid-latitudes, there is a prominent wavenumber-2 pattern with alternating warm and cold anomalies. The warm anomalies are centered over the eastern parts of the ocean basins, while the cold anomalies are centered over the eastern continents.

This looks like a “smoothed out” version of the zonal temperature anomalies at 1000 hPa that we’ve already seen.

Map of zonal wind anomalies \(\overline{v}^*\)#

Similarly we can plot \(\overline{v}^*\) and look for spatial covariance between the anomalous meridional flow and the anomalous temperature fields:

make_map(vbar_star, levels=np.linspace(-12,12,25), title='$\overline{v}^*$ at 500 hPa for January - March 2021');
../_images/711d7674ddd000b931d716057678e6a40d45dc64baa6d1524aede3772f763c1d.png

Here, focussing again on the Northern mid-latitudes, we can see that there are regions where the seasonally-averaged meridional flow is anomalously northward (red) and southward (blue).

The key point is that there is covariance between these two fields: the regions of northward flow (over the ocean basins) tend to also be regions of anomalously warm temperatures, while the regions of southward flow (over the continents) tend to also be regions of anomalously cold temperatures.

Thus the product \([\overline{v}^* \overline{T}^*]\) is positive in the zonal average, even though the individual zonal averages \([\overline{v}^*]\), \([\overline{T}^*]\) are both zero by definition!

Zonal mean heat fluxes in the opposite season#

Let’s just repeat the whole calculation for June through August 2021.

First though, we’ll wrap our calculation and plotting code in functions for easier repetition with different time periods.

def calc_zonal_heat_flux(dateslice, lev=500):
    v500 = v_ds.v.sel(lev=lev, time=dateslice)
    T500 = T_ds.t.sel(lev=lev, time=dateslice)
    T500_global = T500.mean(dim='time').weighted(np.cos(np.deg2rad(T500.lat))).mean(dim=('lat','lon'))

    vTbar = (v500 * (T500 - T500_global)).mean(dim='time')
    vbar = v500.mean(dim='time')
    Tbar = T500.mean(dim='time') - T500_global

    vTbar_zon = vTbar.mean(dim='lon')
    vbar_zon = vbar.mean(dim='lon')
    Tbar_zon = Tbar.mean(dim='lon')
    vprime = v500 - vbar
    Tprime = (T500 - T500_global) - Tbar
    vprime_Tprime_bar = (vprime*Tprime).mean(dim='time')

    vbarzon_Tbarzon = vbar_zon * Tbar_zon
    vbar_star = vbar - vbar_zon
    Tbar_star = Tbar - Tbar_zon
    vbarstar_Tbarstar_zon = (vbar_star * Tbar_star).mean(dim='lon')
    vprime_Tprime_bar_zon = vprime_Tprime_bar.mean(dim='lon')
    
    decomp = {'vTbar': vTbar_zon,
              'vbar_Tbar': vbarzon_Tbarzon,
              'vbarstar_Tbarstar': vbarstar_Tbarstar_zon,
              'vprime_Tprime_bar': vprime_Tprime_bar_zon}
    return decomp

def plot_zonal_heat_flux(decomp):
    decomp['vTbar'].plot(y='lat', label='Total, $[\overline{vT}]$')
    decomp['vbar_Tbar'].plot(y='lat', label='Mean meridional, $[\overline{v}][\overline{T}]$')
    decomp['vbarstar_Tbarstar'].plot(y='lat', label='Stationary eddy, $[\overline{v}^* \overline{T}^*]$')
    decomp['vprime_Tprime_bar'].plot(y='lat', label='Transient eddy, $[\overline{v^\prime T^\prime}]$')
    plt.legend()
    plt.grid()
    return 

plot_zonal_heat_flux(calc_zonal_heat_flux(dateslice=slice('2021-06-01','2021-08-31')))
plt.title('Zonal mean heat flux decomposition at 500 hPa for June - August 2021');
../_images/abde03d181370edf3aed5a9cb12220d0ac85c73c29f93192c7da1054aac13072.png

Here we can see a few differences:

  • The NH peak is much smaller in summer than winter (by a factor of ~4)

  • The SH peak is stronger in its summer, but the winter-summer contrast in the SH is smaller than in the NH

  • Unlike in NH winter, there is very little stationary eddy contribution in the SH winter

Dependence on averaging period#

The decomposition between transient and stationary fluxes depend on the chosen averaging period.

For very short time periods, \(\left[\overline{A^\prime B^\prime}\right] \rightarrow 0\). The transient eddies disappear and all fluctuations appear as stationary eddies.

For timescales of ~3 months, the stationary eddies are seasonal features such as the Aleutian Low, and transient eddies are anything with residence times < 3 months (e.g. traveling synoptic systems)

At 12 month timescale, we have mostly transient eddies but still some stationary features due to orographic forcing.

Zonal mean flux over 1 week#

Let’s repeat a calculation but change the averaging period from 3 months down to 1 week:

plot_zonal_heat_flux(calc_zonal_heat_flux(dateslice=slice('2021-01-01','2021-01-07')))
plt.title('Zonal mean heat flux decomposition at 500 hPa for Jan 1 to Jan 7 2021');
../_images/a56fcf8cd06082fa19d58ec54626cb85ba3eeffe86cb252dd21c4be42d8f77c0.png

Zonal mean heat flux over 1 year#

And do the whole year:

plot_zonal_heat_flux(calc_zonal_heat_flux(dateslice=slice('2021-01-01','2021-12-31')))
plt.title('Zonal mean heat flux decomposition at 500 hPa for all of 2021');
../_images/e1aa3be710c48736b649f9c15fab5f03cabb54cefcff43f8294ebadc0d08500e.png

We can see that the partition between stationary eddy and transient eddy is particularly sensitive to our choice of averaging interval.

Contents