Relative Vorticity

This notebook shows a simple example of calculation the vertical component of relative vorticity, but is applicable to any curl calculation as well.

Background

Relative vorticity is defined as the curl of the velocity field:

\[\zeta = \partial_x v - \partial_y u.\]

For demonstration purposes we will compute the vorticity at some arbitrary depth for a given month, but the method can be applied to any time-varying, 3D field.

We will demonstrate two methods for computing relative vorticity:

  1. The recommended method, leveraging the functionality of xgcm to handle grid operations.

  2. A naïve method of simply converting degrees of longitude/latitude into metres and differentiating via xarray.

We compare these two methods to the model’s vorticity diagnostic. We show that Method 1 is accurate, and should be used over Method 2.

This recipe is intended to work with output from ACCESS-OM2, where MOM5 uses an Arakawa B-grid and velocities are given at cell corners. For adaptation to MOM6, which uses an Arakawa C-grid, you should take into account that velocities are in cell faces. This means that:

MOM5 diagnostic (x-coord, y-coord)

MOM6 diagnostic (x-coord, y-coord)

u (xu_ocean, yu_ocean)

uo (xq, yh)

v (xu_ocean, yu_ocean)

vo (xh, yq)

vorticity_z (xt_ocean, yt_ocean)

RV (xq, yq)

Most experiments will also have four different options of \(\mathrm{d}x\) and \(\mathrm{d}y\) to account for distances in different grids/faces. For MOM6, we recommend that you approach the calculation using Method 1 by adapting the interpolation steps to suit MOM6 diagnostics and where each flow field lives.

[1]:

import intake
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from dask.distributed import Client

Load a dask client.

[2]:

client = Client(threads_per_worker=1, memory_limit=0)
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 45463 instead
  warnings.warn(

[2]:

Client

Client-5fe66124-7a80-11f1-a984-000003fefe80

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

Cluster Info

Open the ACCESS-NRI default catalog and the first cycle of the 0.1\(^{\circ}\) IAF experiment. This recipe uses this experiment because it has velocity and vorticity diagnostics available at the same time step. This will allow us to compare methods.

[3]:

catalog = intake.cat.access_nri

[4]:

experiment = "01deg_jra55v140_iaf"

Load the diagnostics needed for this recipe - note that intake is not able, currently, to load to_dask datasets in different grids, so we will load them separately.

We will load grid diagnostics (distances in meridional and zonal directions, as well as latitude, longitude), one day of velocities at a certain depth level, and MOM5’s relative vorticity diagnostic, which we will use to compare the three methods.

[5]:

geolon_c = catalog[experiment].search(variable=['geolon_c'],
                                      frequency="fx").to_dask()['geolon_c']
geolon_c = geolon_c.sel(xu_ocean=slice(-90, -40), yu_ocean=slice(20, 60))

geolat_c = catalog[experiment].search(variable=['geolat_c'],
                                      frequency="fx").to_dask()['geolat_c']
geolat_c = geolat_c.sel(xu_ocean=slice(-90, -40), yu_ocean=slice(20, 60))

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/access_nri_intake/aliases.py:192: UserWarning: Value aliasing: variable='geolon_c' → variable=['geolon_c','geolon_c']
  norm: dict[str, Any] = self._normalise_kwargs(kwargs)
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/access_nri_intake/aliases.py:192: UserWarning: Value aliasing: variable='geolat_c' → variable=['geolat_c','geolat_c']
  norm: dict[str, Any] = self._normalise_kwargs(kwargs)

[6]:

# u-grid diagnostics
dxu = catalog["01deg_jra55v140_iaf"].search(variable=["dxu"],
                                            frequency="fx").to_dask().sel(xu_ocean=slice(-90, -40), yu_ocean=slice(20, 60))
dyu = catalog["01deg_jra55v140_iaf"].search(variable=["dyu"],
                                            frequency="fx").to_dask().sel(xu_ocean=slice(-90, -40), yu_ocean=slice(20, 60))

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/access_nri_intake/aliases.py:192: UserWarning: Value aliasing: variable='dxu' → variable=['dxu','dxu']
  norm: dict[str, Any] = self._normalise_kwargs(kwargs)
/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(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/access_nri_intake/aliases.py:192: UserWarning: Value aliasing: variable='dyu' → variable=['dyu','dyu']
  norm: dict[str, Any] = self._normalise_kwargs(kwargs)
/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(

[7]:

# t-grid diagnostics
dxt = catalog["01deg_jra55v140_iaf"].search(variable=["dxt"],
                                            frequency="fx").to_dask().sel(xt_ocean=slice(-90, -40), yt_ocean=slice(20, 60))
dyt = catalog["01deg_jra55v140_iaf"].search(variable=["dyt"],
                                            frequency="fx").to_dask().sel(xt_ocean=slice(-90, -40), yt_ocean=slice(20, 60))

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/access_nri_intake/aliases.py:192: UserWarning: Value aliasing: variable='dxt' → variable=['dxt','dxt']
  norm: dict[str, Any] = self._normalise_kwargs(kwargs)
/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(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/access_nri_intake/aliases.py:192: UserWarning: Value aliasing: variable='dyt' → variable=['dyt','dyt']
  norm: dict[str, Any] = self._normalise_kwargs(kwargs)
/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(

Load velocity as monthly snapshots (indicated by the variable cell methods argument)

[8]:

velocity_diagnostics = ["u", "v"]
depth_level = 30
def subset_ugrid(ds):
    ds = ds.sel(xu_ocean=slice(-90, -40), yu_ocean=slice(20, 60)).sel(st_ocean=depth_level, method='nearest')
    return ds

ds_velocity = catalog[experiment].search(variable=velocity_diagnostics,
                                         frequency="1mon",
                                         variable_cell_methods='time: point'
                                         ).to_dask(preprocess=subset_ugrid)

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/access_nri_intake/aliases.py:192: UserWarning: Value aliasing: variable='u' → variable=['u','u']
  norm: dict[str, Any] = self._normalise_kwargs(kwargs)
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/access_nri_intake/aliases.py:192: UserWarning: Value aliasing: variable='v' → variable=['v','v']
  norm: dict[str, Any] = self._normalise_kwargs(kwargs)
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/intake_esm/source.py:305: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly.
  self._ds = xr.combine_by_coords(

Load vorticity as monthly snapshots (indicated by the variable cell methods argument)

[9]:

def subset_tgrid(ds):
    ds = ds.sel(xt_ocean=slice(-90, -40), yt_ocean=slice(20, 60)).sel(st_ocean=depth_level, method='nearest')
    return ds

ds_vorticity = catalog[experiment].unwrap().search(variable="vorticity_z",
                                                   frequency="1mon",
                                                   variable_cell_methods='time: point'
                                                   ).to_dask(preprocess=subset_tgrid)

Combine all we need in the same dataset and select a random time step

[10]:

ds = xr.merge([ds_velocity, ds_vorticity, dxu, dyu, dxt, dyt], compat='override')
ds = ds.isel(time=-1)

Method 1: Using xgcm

`xgcm <https://xgcm.readthedocs.io/en/stable/>`__ is a package that deals with staggered grids that are typically used in ocean models. An excerpt from xgcm’s docs mentions:

“(in model output datasets), different variables are located at different positions with respect to a volume or area element (e.g. cell center, cell face, etc.) xgcm solves the problem of how to interpolate and difference these variables from one position to another.”

[11]:

import xgcm

We first need to create a grid object that has all the information regarding our staggered grid. For our case, grid needs to know the location of the xt_ocean, xu_ocean points (and same for \(y\)) and their relative orientation to one another, i.e., that xu_ocean is shifted to the right of xt_ocean by \(\frac1{2}\) grid-cell.

xgcm also expects you to pass on grid information in which xt_ocean, xu_ocean are of the same length and staggered in the correct direction (u to the right of t) - and same for the y-direction. Lets check that.

[12]:

ds

[12]:

<xarray.Dataset> Size: 8MB
Dimensions:      (yu_ocean: 550, xu_ocean: 500, yt_ocean: 551, xt_ocean: 500)
Coordinates:
  * yu_ocean     (yu_ocean) float64 4kB 20.08 20.17 20.26 ... 59.87 59.92 59.97
  * xu_ocean     (xu_ocean) float64 4kB -89.9 -89.8 -89.7 ... -40.2 -40.1 -40.0
  * yt_ocean     (yt_ocean) float64 4kB 20.03 20.12 20.22 ... 59.9 59.95 60.0
  * xt_ocean     (xt_ocean) float64 4kB -89.95 -89.85 -89.75 ... -40.15 -40.05
    st_ocean     float64 8B 29.45
    time         datetime64[ns] 8B 2019-01-01
Data variables:
    u            (yu_ocean, xu_ocean) float32 1MB dask.array<chunksize=(171, 260), meta=np.ndarray>
    v            (yu_ocean, xu_ocean) float32 1MB dask.array<chunksize=(171, 260), meta=np.ndarray>
    vorticity_z  (yt_ocean, xt_ocean) float32 1MB dask.array<chunksize=(171, 260), meta=np.ndarray>
    dxu          (yu_ocean, xu_ocean) float32 1MB dask.array<chunksize=(550, 500), meta=np.ndarray>
    dyu          (yu_ocean, xu_ocean) float32 1MB dask.array<chunksize=(550, 500), meta=np.ndarray>
    dxt          (yt_ocean, xt_ocean) float32 1MB dask.array<chunksize=(551, 500), meta=np.ndarray>
    dyt          (yt_ocean, xt_ocean) float32 1MB dask.array<chunksize=(551, 500), meta=np.ndarray>
Attributes:
    title:                                   ACCESS-OM2-01
    grid_type:                               mosaic
    grid_tile:                               1
    intake_esm_attrs:file_id:                ocean.1mon.st_edges_ocean:76.st_...
    intake_esm_attrs:frequency:              1mon
    intake_esm_attrs:variable_cell_methods:  ,,,time: point,,
    intake_esm_attrs:variable_units:         meters,meters,days since 1900-01...
    intake_esm_attrs:realm:                  ocean
    intake_esm_attrs:temporal_label:         point
    intake_esm_attrs:_data_format_:          netcdf
    intake_esm_dataset_key:                  ocean.1mon.st_edges_ocean:76.st_...

The x-direction is correct, but in the y-direction we need to remove the last t-point:

[13]:

ds = ds.isel(yt_ocean=slice(None,-1))

Create a grid object that stores grid information (coordinate names, properties and dimensions):

[16]:

metrics = {
    ("X",): ["dxt", "dxu"],  # X distances
    ("Y",): ["dyt", "dyu"],  # Y distances
}

grid = xgcm.Grid(ds,
                 metrics=metrics,
                 coords={'X': {'center': 'xt_ocean', 'right': 'xu_ocean'},
                         'Y': {'center': 'yt_ocean', 'right': 'yu_ocean'}},
                 autoparse_metadata=False)

We can use this grid object to interpolate diagnostics across grids - for example, grid.interp(v, 'Y') will interpolate v in the y-direction, bringing v from yu_ocean to yt_ocean.

We can also compute derivatives using the .derivative function. For example, \(\partial_x v\) is obtained via grid.derivative(v, 'X') - note that the result will also shift v in the x-direction from corners (xu_ocean) to centres (xt_ocean).

We want to calculate vorticity at the centres of tracer cells in order to compare to the vorticity_z diagnostic - so we will interpolate accordingly.

[18]:

vx = grid.derivative(ds.v, "X", boundary="extend")
vx = grid.interp(vx, "Y", boundary="extend")
uy = grid.derivative(ds.u, "Y", boundary="extend")
uy = grid.interp(uy, "X", boundary="extend")

ζ_xgcm = vx - uy

ζ_xgcm = ζ_xgcm.rename("Relative Vorticity")
ζ_xgcm.attrs["long_name"] = "Relative Vorticity, ∂v/∂x-∂u/∂y"
ζ_xgcm.attrs["units"] = "s-1"

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/xgcm/grid.py:493: UserWarning: Metric at ('yu_ocean', 'xt_ocean') being interpolated from metrics at dimensions ('yu_ocean', 'xu_ocean'). Boundary value set to 'extend'.
  warnings.warn(
/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/xgcm/grid.py:493: UserWarning: Metric at ('yt_ocean', 'xu_ocean') being interpolated from metrics at dimensions ('yu_ocean', 'xu_ocean'). Boundary value set to 'extend'.
  warnings.warn(

Now, let’s plot ζ_xgcm, as well as the difference with the vorticity diagnostic:

[19]:

maxvalue = 5e-5
levels = np.linspace(-maxvalue, maxvalue, 24)

[21]:

fig, axs = plt.subplots(1,2,figsize=(12, 4))

c=axs[0].contourf(geolon_c,
                  geolat_c,
                  ζ_xgcm,
                  levels=levels,
                  cmap='RdBu_r',
                  extend='both')
plt.colorbar(c)
axs[1].set_title("$\zeta_{xgcm}$")

c=axs[1].contourf(geolon_c,
                  geolat_c,
                  np.abs(ζ_xgcm - ds['vorticity_z']),
                  levels=np.linspace(0, 1e-10, 22),
                  cmap='Oranges',
                  extend='max')
plt.colorbar(c)
axs[1].set_title("$\zeta_{xgcm}$ - vorticity_z")

for ax in axs:
    ax.set_xlabel("longitude")
    ax.set_ylabel("latitude")
    ax.set_xlim(-90, -40)
    ax.set_ylim(20, 60);

plt.tight_layout()
../_images/03-Advanced-Recipes_Relative_Vorticity_30_0.png

Note that xgcm doesn’t handle the edges of the domain quite well, and a similar issue could happen around coastlines - you should be extra careful in those areas!

Another thing to consider is that the model uses double precision arithmetic - but these diagnostics are saved as float32. Therefore, the differences between our offline calculation and the diagnostics are not exactly zero - but they are around 5 orders of magnitude smaller which is pretty good!

Method 2 (naïve computation using xarray)

To compute relative vorticity \(\zeta = \partial_x v - \partial_y u\) we simply differentiate the velocity components with respect of lon (here xu_ocean in degrees) and lat (here yu_ocean in degrees). We then convert the derivatives from units of degrees\(^{-1}\) to m\(^{-1}\). To do so, we use the value of the radius of the Earth Rearth and also take into account that as we go towards the poles the lon-grid spacing is scaled by \(\cos(\) lat \()\).

(Note the unicode characters like ζ can be used in python.)

[22]:

# values used by MOM5
Ω = 7.292e-5  # Earth's rotation rate; in radians/s
Rearth = 6371e3  # Earth's radius; in m

Calculate the Coriolis parameter \(f = 2\Omega \sin(\texttt{lat})\).

[23]:

f = 2 * Ω * np.sin(np.deg2rad(geolat_c))  # convert lat in radians
f = f.rename("Coriolis")
f.attrs["long_name"] = "Coriolis parameter"
f.attrs["units"] = "s-1"
f.attrs["coordinates"] = "geolon_c geolat_c"

Now we can calculate relative vorticity:

[24]:

scale_factor_x = np.deg2rad(Rearth * np.cos(np.deg2rad(geolat_c)))
scale_factor_y = (np.pi / 180 * Rearth)
ζ_naive = ds["v"].differentiate("xu_ocean") / scale_factor_x - ds["u"].differentiate("yu_ocean") / scale_factor_y
ζ_naive = ζ_naive.rename("Relative Vorticity").drop_vars('geolat_c')
ζ_naive.attrs["long_name"] = "Relative Vorticity, ∂v/∂x-∂u/∂y"
ζ_naive.attrs["units"] = "s-1"

We now plot ζ_naive and the difference from the vorticity diagnostic - note that for this we actually need to bring ζ_naive to the t-cells.

[25]:

ζ_naive_t = grid.interp(ζ_naive, ["X", "Y"], boundary="extend")

[27]:

fig, axs = plt.subplots(1, 2, figsize=(12, 4))

c=axs[0].contourf(geolon_c,
                 geolat_c,
                 ζ_naive,
                 levels=levels,
                 cmap='RdBu_r',
                 extend='both')
plt.colorbar(c)
axs[1].set_title("$\zeta_{xgcm}$")

c=axs[1].contourf(geolon_c,
                  geolat_c,
                  np.abs(ζ_naive_t - ds['vorticity_z']),
                  levels=np.linspace(0, 1e-5, 22),
                  cmap='Oranges',
                  extend='max')
plt.colorbar(c)
axs[1].set_title("$\zeta_{xgcm}$ - vorticity_z")

for ax in axs:
    ax.set_xlabel("longitude")
    ax.set_ylabel("latitude")
    ax.set_xlim(-90, -40)
    ax.set_ylim(20, 60)

plt.tight_layout()
../_images/03-Advanced-Recipes_Relative_Vorticity_41_0.png

As you can see this method is very different from the vorticity_z diagnostic - the residual is of the same order of magnitude as the diagnostic.

Commentary on the two methods

Below are some brief discussions/conclusions on the two methods showcased in this recipe.

Note 1

We have demonstrated that using xgcm is the most accurate way of calculating relative vorticity - and the simplest one. Using offline naive computations might lead to significantly different results - so you should be careful if you are interested in budget calculations.

Note 2

There are other ways of calculating derivatives with xgcm documented online, namely using .diff() to differentiate and then dividing by the appropriate grid metric. You should be aware that how you choose to do this will impact the end result - you can read a more extensive discussion in this issue. As this recipe demonstrates, the recommendation is to use .derivative().

Note 3

It is always useful to understand how the model actually calculates a diagnostic. In the case of vorticity_z, you can see here that the vorticity of a tracer cell is given by the average of \(\partial_x v\) at its northern and southern faces, and equivalently \(\partial_y u\) is the average of the derivatives at the western and eastern faces. Taking this into account, it is possible to do a naive computation that matches the diagnostic, but you must have understood the model’s code before doing it! Jumping in head first without knowledge of how the model computed the diagnostic did not yield correct results, as exemplified by Method 2.

Note 4

Remember to be careful around coastlines, the edges of your domain, and north of 65 North in the tripolar region!

Example application - the Rossby number

To conclude, let’s visualize the Rossby number

\[\mathrm{Ro} = \frac{\zeta }{f},\]

where \(f=2 \Omega \sin \theta\) is the Coriolis parameter.

[23]:

f = grid.interp(grid.interp(f, "X"), "Y", boundary="extend")

Ro = ζ_xgcm / f
Ro = Ro.rename("Rossby number ζ/f")

[24]:

plt.figure(figsize=(8, 5))

Ro.plot.contourf(levels=np.linspace(-0.15, 0.15, 21),
                 cmap='RdBu_r',
                 extend='both')

plt.title(
    "Rossby number $\zeta/f$"
)
plt.xlabel("longitude")
plt.ylabel("latitude")
plt.xlim(-90, -40)
plt.ylim(20, 60);
../_images/03-Advanced-Recipes_Relative_Vorticity_46_0.png