Barotropic Streamfunction

This recipe demonstrates how to compute and plot the barotropic streamfunction (\(\psi\)).

Background

The barotropic streamfunction (\(\psi\)) is obtained from the integration of the depth-integrated transports starting from a physical boundary at which we know the transport is zero - therefore there are different ways to calculate it depending on your choice of boundary for the integration.

This recipe calculates \(\psi\) integrating the depth-integrated zonal transport, \(U\) in the meridional direction, starting from the Antarctic continent:

\[\psi = \int_{y_{\rm Antarctica}}^{y} U \, \mathrm{d}y ,\]

You can see from this expression that the direction of the transport is then parallel to streamlines, and the intensity of that transport is given by the difference between streamlines. For this statement to be valid, the flow must be incompressible (or approximately so).

Requirements

This recipe uses MOM6 output.

The workflow is,

  1. Load the MOM6 depth-integrated zonal mass transport (umo_2d).

  2. Convert it to volume transport in Sverdrups (Sv).

  3. Integrate meridionally (cumulative sum along latitude) to obtain the barotropic streamfunction (\(\psi\)).

  4. Plot a circumpolar map with a circular boundary and land mask.

MOM5 adaptation

To adapt this recipe for MOM5 output, the following diagnostic equivalences may be useful:

MOM6 diagnostic (x-coord,y-coord)

MOM5 diagnostic (x-coord,y-coord)

umo_2d(xq, yh)

tx_trans_int_z(xu_ocean, xt_ocean)

deptho(xh, yh)

ht(xt_ocean, xt_ocean)

 
[1]:

import cartopy.crs as ccrs
import cmocean
import intake
import matplotlib.path as mpath
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import dask.distributed as dask

[2]:

client = dask.Client(threads_per_worker=1)
client

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/distributed/node.py:188: UserWarning: Port 8787 is already in use.
Perhaps you already have a cluster running?
Hosting the HTTP server on port 33711 instead
  warnings.warn(

[2]:

Client

Client-bb4ded8a-7675-11f1-88a2-0000008afe80

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: /proxy/33711/status

Cluster Info

 
[3]:

catalog = intake.cat.access_nri
experiment = "panant-01-zstar-v13"

Load transports - use a preprocessing function to subselect a domain.

[4]:

def subset_domain(ds):
    ds = ds.sel(yh=slice(None, -50))
    return ds

umo = catalog[experiment].search(variable='umo_2d',
                                start_date='200[0,1].*',
                                frequency='1mon'
                                ).to_dask(preprocess=subset_domain)
deptho = catalog[experiment].search(variable='deptho',
                                    frequency='fx'
                                    ).to_dask(preprocess=subset_domain)

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/intake_esm/source.py:314: ConcatenationWarning: Attempting to concatenate datasets without valid dimension coordinates: retaining only first dataset. Request valid dimension coordinate to silence this warning.
  warnings.warn(

Take the time mean of the transport:

[5]:

mass_transport = umo['umo_2d'].mean('time')
mass_transport

[5]:

<xarray.DataArray 'umo_2d' (yh: 666, xq: 3601)> Size: 10MB
dask.array<mean_agg-aggregate, shape=(666, 3601), dtype=float32, chunksize=(666, 3601), chunktype=numpy.ndarray>
Coordinates:
  * yh       (yh) float64 5kB -81.11 -81.07 -81.02 ... -50.16 -50.09 -50.03
  * xq       (xq) float64 29kB -280.0 -279.9 -279.8 -279.7 ... 79.8 79.9 80.0
Attributes:
    units:          kg s-1
    long_name:      Ocean Mass X Transport Vertical Sum
    cell_methods:   yh:sum xq:point time: mean
    time_avg_info:  average_T1,average_T2,average_DT
    standard_name:  ocean_mass_x_transport_vertical_sum
    interp_method:  none

Looking at the mass_transport attributes, we can confirm that it is a mass transport from its units. To convert to volume transoprt we need to divide by the reference density.

[6]:

ρ0 = 1030  # reference seawater density, kg m⁻³
volume_transport = mass_transport / ρ0   # convert kg s⁻¹ -> m³ s⁻¹

To integrate the transport in the meridional direction, we just need to sum cumulatively northwards - this is because the transport is the total amount of water through a grid face.

We also divide by 1e6 to convert from m³ s⁻¹ to Sv (1 Sv = 10⁶ m³ s⁻¹)

[7]:

psi = volume_transport.cumsum('yh') / 1e6 # convert m³ s⁻¹ -> Sv

Now lets create a land mask for plotting:

[8]:

land_mask = xr.where(np.isnan(deptho['deptho']), 1, np.nan)

[9]:

fig = plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.SouthPolarStereo())
ax.set_extent([-180, 180, -80, -50], crs=ccrs.PlateCarree())
ax.set_facecolor('lightgrey')

# Map the plot boundaries to a circle
theta = np.linspace(0, 2 * np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)

psi.plot.contourf(ax=ax,
                  levels=np.arange(-50, 160, 10),
                  extend='both',
                  cmap=cmocean.cm.speed,
                  cbar_kwargs={'shrink':.7, 'label':'Sv'},
                  transform=ccrs.PlateCarree())

land_mask.plot.contourf(ax=ax,
                        colors=['lightgrey'],
                        add_colorbar=False,
                        transform=ccrs.PlateCarree())


plt.title("Barotropic streamfunction (MOM6)")
plt.tight_layout()
../_images/02-Easy-Recipes_Barotropic_Streamfunction_14_0.png