Coordinate transformation from z levels to density levels

Rebinning ty_trans to ty_trans_rho density levels in MOM5

This transformation is commonly used for the purpose of decomposing the residual meridional overturning streamfunction into mean and eddy components. The mean is taken to be the time-mean Eulerian transport (time-mean in z coordinates), whilst the residual transport is the time-mean in density coordinates (density surface move in time). The eddy transport is the difference between the residual overturning streamfunction, and the Eulerian transport transformed from depth to density coordinates via the time-mean density.

Binning is the discrete version of this transformation from depth to density coordinates. We define target density bins, and for each bin we add the quantity to be binned from all cells with a density that satisfies that bin. In its simplest form, binning is creating a histogram.

We use three different binning methods:

  1. xhistogram, which performs binning in exactly the way MOM5 does and is thus most appropriate for calculating eddy quantities (define edge of bins)

  2. xgcm conservative binning (define edge of bins, but vertically interpolates so looks smoother than xhistogram)

  3. Using density coordinate binning method of Lee et al. (2007) (define isopycnals to bin onto i.e. centre)

Compute times were calculated using the XXLargeMem (28 cpus, 252 Gb mem) Jupyter Lab on NCI’s ARE, using conda environment analysis3-26.06.

This notebook performs direct coordinate searches using the access-nri intake catalog, and so will require conda environment analysis3-25.02 or later.

A starting point for conversion to MOM6

This notebook uses MOM5 variables. To convert this notebook to MOM6, the following variables are required. However, note that MOM6 pan-Antarctic and ACCESS-OM3 models currently do not output 3D volume transports in both z and density coordinates so this particular example won’t be able to be replicated without additional model runs.

Important note In MOM6, online density-binned output diagnostics may have the same variable name as the z-coordinate binned output. In this case, they should have different filenames to distinguish them.

MOM5 variable

MOM6 variable and file

ty_trans

vmo (no COSIMA MOM6 models output this yet)

ty_trans_rho

vmo, filename = '*.ocean_month_rho2.nc' for panan and filename = 'access-om3.mom6.3d.vmo+rho2.1mon.mean.*.nc' for ACCESS-OM3

pot_rho_2

rhopot2

yu_ocean

yq

xt_ocean

xh

yt_ocean

yh

st_ocean

z_l

st_edges_ocean

z_i

potrho

rho2_l

potrho_edges

rho2_i

 
[1]:

import intake

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import glob
import cmocean.cm as cmocean
import xgcm
from xhistogram.xarray import histogram

from dask.distributed import Client

[2]:

client = Client(threads_per_worker = 1)
client

[2]:

Client

Client-db0872fc-7509-11f1-90fd-000003bdfe80

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

Cluster Info

 
[3]:

experiment = '01deg_jra55v13_ryf9091'

catalog = intake.cat.access_nri
expt_datastore = catalog[experiment]

rho_0 = 1035.0 # kg/m^3 reference density
g = 9.81

# reduce computation by choosing only Southern Ocean latitudes
lat_range = slice(-70, -34.99)

Step 1: Load density and quantity being binned

(you can choose any other quantity instead of ty_trans, e.g. dzt. Interpolate density to whatever grid that variable is on)

For the simplicity of the demonstrations here, we will only load two months of monthly ty_trans, ty_trans_rho and pot_rho_2, then take the time average after we weight with the number of days in each month.

Notes on searching for data

  • We can search for multiple variables at once: variable=['ty_trans','ty_trans_rho','pot_rho_2'] will get all three in the same search, which is faster than searching multiple times

  • The access-nri intake catalog doesn’t support start and end dates yet. To work around this, we use ‘regex’ (regular expression) to search for dates. This regex: start_date='2170-0[1-3].*' says “find me all strings that start with either ‘1970-01’, ‘1970-02’, or ‘1970-03’”. For more info on regex’s, see here.

[4]:

%%time
start_time = '2170-01-01'
end_time = '2170-03-01'
time_slice = slice(start_time, end_time)

ds = expt_datastore.search(
    variable=['ty_trans','ty_trans_rho','pot_rho_2'],
    start_date='2170-0[1-3].*',
).to_dask(xarray_open_kwargs={
    'chunks' : 'auto',
})

ty_trans, ty_trans_rho, pot_rho_2 = ds['ty_trans'], ds['ty_trans_rho'], ds['pot_rho_2']

/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='ty_trans' → variable=['ty_trans','ty_trans']
  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='ty_trans_rho' → variable=['ty_trans_rho','ty_trans_rho']
  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='pot_rho_2' → variable=['pot_rho_2','pot_rho_2']
  norm: dict[str, Any] = self._normalise_kwargs(kwargs)

CPU times: user 3.4 s, sys: 1.67 s, total: 5.07 s
Wall time: 9.28 s

[5]:

%%time
ty_trans = ty_trans.sel(time = time_slice).sel(yu_ocean = lat_range)

# weighted time-mean by month length
days_in_month = ty_trans.time.dt.days_in_month
total_days = days_in_month.sum()
ty_trans = (ty_trans * days_in_month).sum('time') / total_days

CPU times: user 118 ms, sys: 69.8 ms, total: 188 ms
Wall time: 137 ms

[6]:

%%time

pot_rho_2 = pot_rho_2.sel(time = time_slice).sel(yt_ocean = lat_range)
pot_rho_2 = (pot_rho_2 * days_in_month).sum('time') / total_days

CPU times: user 62.3 ms, sys: 36.4 ms, total: 98.7 ms
Wall time: 74.6 ms

[7]:

%%time

ty_trans_rho = ty_trans_rho.sel(time = time_slice).sel(grid_yu_ocean = lat_range)
ty_trans_rho = (ty_trans_rho * days_in_month).sum('time') / total_days

CPU times: user 60.5 ms, sys: 46.1 ms, total: 107 ms
Wall time: 77.8 ms

Create an xgcm grid for interpolation, and then interpolate density onto the meridional transport grid

[8]:

ds = xr.Dataset({'ty_trans': ty_trans, 'pot_rho_2': pot_rho_2})
grid = xgcm.Grid(ds, coords = {'Y': {'center': 'yt_ocean', 'right': 'yu_ocean'}}, periodic = False,
                autoparse_metadata = False)

# interpolate density (t-grid) to ty_trans grid (xt_ocean and yu_ocean)
pot_rho_2 = grid.interp(pot_rho_2, 'Y', boundary = 'extend')

Method 1: xhistogram

We use xhistogram which is an xarray-aware method for computing histograms; see its documentation.

Computation in xhistogram occurs via the same method as in MOM5 online binning. It is thus most appropriate for an comparisons between offline and online binned quantities.

Output will be an array with coordinate density the linear centre of these bins. If we choose potrho_edges, the end result will have coordinates potrho, which is the same as online binned ty_trans_rho.

[9]:

potrho_edges = expt_datastore.search(
    variable='potrho_edges',
    start_date='2170-0[1-3].*',
    frequency='1mon'
).to_dask(
    xarray_open_kwargs= {
        'chunks' : 'auto',
        'decode_timedelta': False,
    }
)['potrho_edges']

targetbins = potrho_edges.values

We include the variable we want to bin in weights. This quantity should be extensive, since grid cells vary in sizes, which is true because ty_trans is multiplied by the cell size in x and z directions.

[10]:

# Make sure the variables have a name, otherwise xhistogram doesn't know what to call the bins
ty_trans = ty_trans.rename('ty_trans')
pot_rho_2 = pot_rho_2.rename('pot_rho_2')

ty_trans_mean = histogram(pot_rho_2,
                          bins = [targetbins],
                          dim = ['st_ocean'],
                          weights = ty_trans).rename({pot_rho_2.name + '_bin': 'potrho',
                                                     'xt_ocean': 'grid_xt_ocean',
                                                     'yu_ocean': 'grid_yu_ocean'})

[11]:

%%time
ty_trans_mean = ty_trans_mean.load()

CPU times: user 21.1 s, sys: 14.2 s, total: 35.3 s
Wall time: 1min 1s

[12]:

def cumsum_from_bottom(residual):
    cumsum = residual.cumsum('potrho') - residual.sum('potrho')
    return cumsum

[13]:

%%time
psi_avg = cumsum_from_bottom(ty_trans_rho.sum('grid_xt_ocean') / (1e6 * rho_0)).load() # .load() since we'll use it a few times in this notebook

psi_avg_mean = cumsum_from_bottom(ty_trans_mean.sum('grid_xt_ocean') / (1e6 * rho_0))

CPU times: user 2.56 s, sys: 1.62 s, total: 4.19 s
Wall time: 3.53 s

[14]:

%%time
fig, axes = plt.subplots(nrows = 1, ncols = 3, figsize = (15, 4))

levels = np.linspace(-25, 25, 26)

psi_avg.plot.contourf(ax = axes[0], x = 'grid_yu_ocean', levels = levels, add_colorbar = False)

psi_avg_mean.plot.contourf(ax = axes[1], x = 'grid_yu_ocean', levels = levels, add_colorbar = False)

p = (psi_avg - psi_avg_mean).plot.contourf(ax = axes[2], x = 'grid_yu_ocean', levels = levels, add_colorbar = False)

cbar_ax = fig.add_axes([0.92, 0.15, 0.01, 0.7])
fig.colorbar(p, cax=cbar_ax, label='Transport (Sv)')

axes[0].set_title('Residual')
axes[1].set_title('Mean')
axes[2].set_title('Eddy')

for ax in axes:
    ax.set_ylim(1037.5, 1032)
    ax.set_xlim(-70, -35);

CPU times: user 305 ms, sys: 150 ms, total: 454 ms
Wall time: 353 ms
../_images/03-Advanced-Recipes_Transformation_from_Depth_to_Potential_Density_20_1.png

Method 2: xgcm

Use xgcm’s conservative binning, described in the tutorial available in the xgcm documentation.

This method results in a smoother vertical distribution than xhistogram, as it is not quite a histogram but does some interpolation to the top and bottom of each vertical cell. We have found that the computation currently has issues with interior land boundaries: xgcm isn’t able to compute partial cells at the bottom of the ocean in MOM5. We thus add a correction after the computation to ensure that vertical integrals are preserved and thus that streamfunctions calculated are closed.

[15]:

%%time

ds = expt_datastore.search(
    variable=['st_ocean','st_edges_ocean'],
    frequency='1mon', # and `frequency` gets us down to a single dataset
    start_date='2170-0[1-3].*',
).to_dask(
    xarray_open_kwargs={
        'chunks' : 'auto', # This is the easiest way to get a dataset to open relatively quickly
    }
)

st_ocean = ds['st_ocean']
st_edges_ocean = ds['st_edges_ocean']

/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='st_ocean' → variable=['st_ocean','st_ocean']
  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='st_edges_ocean' → variable=['st_edges_ocean','st_edges_ocean']
  norm: dict[str, Any] = self._normalise_kwargs(kwargs)

CPU times: user 1.11 s, sys: 584 ms, total: 1.7 s
Wall time: 2.48 s

[16]:

ds = expt_datastore.search(
    variable='potrho_edges',
    frequency='1mon',
    start_date='2170-0[1-3].*',
).to_dask(
    xarray_open_kwargs={
        'chunks' : 'auto', # This is the easiest way to get a dataset to open relatively quickly
    }
)


pot_rho_2_target = ds['potrho_edges'].values

[17]:

%%time

ds = xr.Dataset({'ty_trans': ty_trans, 'pot_rho_2': pot_rho_2})

# Note: Quantity must be extensive, e.g., ty_trans, vhrho_nt, uhrho_et, tx_trans but *not*, e.g., v, salt.
#       If not extensive, then we need to multiply by dzt or dzu.
ds = ds.assign_coords({'st_edges_ocean': st_edges_ocean})
ds = ds.chunk({'st_edges_ocean': 76, 'st_ocean': 75}) # xgcm doesn't like it if there is more than 1 chunk in this axis

grid = xgcm.Grid(ds, coords={'Z': {'center': 'st_ocean', 'outer': 'st_edges_ocean'}}, periodic = False,
                autoparse_metadata = False)
ds['pot_rho_2_outer'] = grid.interp(ds.pot_rho_2, 'Z', boundary='extend')
ds['pot_rho_2_outer'] = ds['pot_rho_2_outer'].chunk({'st_edges_ocean': 76})

ty_trans_transformed_cons = grid.transform(ds.ty_trans,
                                           'Z',
                                           pot_rho_2_target,
                                           method='conservative',
                                           target_data=ds.pot_rho_2_outer)

### change name to fit previous bins because default xgcm naming is different
# pot_rho_2_outer is actually the CENTRE of bins (I guess the convention was defined for linear interpolation not conservative)
ty_trans_transformed_cons = ty_trans_transformed_cons.rename({'pot_rho_2_outer': 'potrho'})

CPU times: user 408 ms, sys: 306 ms, total: 714 ms
Wall time: 505 ms

[18]:

%%time
ty_trans_transformed_cons = ty_trans_transformed_cons.load()

CPU times: user 5.55 s, sys: 3.87 s, total: 9.42 s
Wall time: 12.9 s

We do this by finding what is missing from the vertical integral (residual). This is what was in the partial cell. We then add that residual into the densest cell that currently exists in the transformed data. This is only an approximation, because the density of the partial cell may be denser than that of the cell above it. However, this correction means the vertical integral is preserved under the transformation.

[19]:

%%time
#find residual from vertical integral
ty_trans_residual = ty_trans.sum('st_ocean') - ty_trans_transformed_cons.sum('potrho') #this is positive definite

CPU times: user 1.62 s, sys: 1.15 s, total: 2.77 s
Wall time: 1.79 s

[20]:

%%time
# select out bottom values:
ty_trans_transformed_cons2 = ty_trans_transformed_cons.where(ty_trans_transformed_cons != 0)
dens_array = ty_trans_transformed_cons2 * 0 + ty_trans_transformed_cons2.potrho # array of isopycnal value where it exists and nan elsewhere
max_dens = dens_array.max(dim = 'potrho', skipna = True)

CPU times: user 2.27 s, sys: 1.68 s, total: 3.95 s
Wall time: 2.54 s

[21]:

%%time
ty_trans_residual_array = (dens_array.where(dens_array == max_dens)*0 + 1) * ty_trans_residual
ty_trans_new = ty_trans_residual_array.fillna(0)+ ty_trans_transformed_cons

CPU times: user 3.59 s, sys: 2.56 s, total: 6.15 s
Wall time: 3.99 s

[22]:

# rename coords to match ty_trans_rho
ty_trans_new = ty_trans_new.rename({'yu_ocean': 'grid_yu_ocean', 'xt_ocean': 'grid_xt_ocean'})

[23]:

%%time
# We need to rechunk here to make the chunks a bit smaller, otherwise the dask workers will die.
# I've halved the chunk size on `grid_xt_ocean` and `potrho`. Sometimes, {'chunks' : 'auto'} can be
# a bit too ambitious!

ty_trans_new = ty_trans_new.chunk(
    chunks = {
        'grid_yu_ocean': 337,
        'grid_xt_ocean': 200, # was 400
        'potrho': 40, # was 80
    }
)
ty_trans_new = ty_trans_new.load()

/g/data/xp65/public/apps/med_conda/envs/analysis3-26.06/lib/python3.12/site-packages/distributed/client.py:3387: UserWarning: Sending large graph of size 2.63 GiB.
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 4.3 s, sys: 5.48 s, total: 9.78 s
Wall time: 11.4 s

[24]:

psi_avg_mean_2 = cumsum_from_bottom(ty_trans_new.sum('grid_xt_ocean') / (1e6 * rho_0))

[25]:

%%time

fig, axes = plt.subplots(nrows = 1, ncols = 3, figsize = (15, 4))

levels = np.linspace(-25, 25, 26)

psi_avg.plot.contourf(ax = axes[0], x = 'grid_yu_ocean',levels = levels, add_colorbar = False)

psi_avg_mean_2.plot.contourf(ax = axes[1], x = 'grid_yu_ocean', levels = levels, add_colorbar = False)

p = (psi_avg - psi_avg_mean_2).plot.contourf(ax = axes[2], x = 'grid_yu_ocean', levels = levels, add_colorbar = False)

cbar_ax = fig.add_axes([0.92, 0.15, 0.01, 0.7])
fig.colorbar(p, cax=cbar_ax, label='Transport (Sv)')

axes[0].set_title('Residual')
axes[1].set_title('Mean')
axes[2].set_title('Eddy')

for ax in axes:
    ax.set_ylim(1037.5, 1032)
    ax.set_xlim(-70, -35);

CPU times: user 239 ms, sys: 120 ms, total: 358 ms
Wall time: 275 ms
../_images/03-Advanced-Recipes_Transformation_from_Depth_to_Potential_Density_39_1.png

Method 3: Bin ty_trans (mean Eulerian overturning) into density bins

This uses the method by Lee et al. (2007) and which is is described as follows. Firstly, cells in the model output with a density \(\rho\) between two prescribed densities (\(\rho_\text{heavy} >\rho> \rho_\text{light}\)) are selected. Each cell is assigned a proximity to the lighter density \(\rho_\text{light}\), which is the ‘bin fraction’

\[f_b = \frac{\rho_\text{heavy}-\rho}{\rho_\text{heavy}-\rho_\text{light}}.\]

Here, a bin fraction of 1 means the cell density \(\rho = \rho_\text{light}\) and \(f_b=0\) means \(\rho = \rho_\text{heavy}\). The quantity being binned, such as the meridional transport \(vh\), is then multiplied by \(f_b\) and added to the lighter density \(\rho_\text{light}\) bin’s meridional transport, followed by the \(vh(1-f_b)\) being added to the heavier bin, \(\rho_\text{heavy}\). This process is repeated for all sets of consecutive bins, meaning that density bins have input from model output cells with density slightly lower and higher than it.

Lee, M., Nurser, A., Coward, A., and De Cuevas, B. (2007). Eddy advective and diffusive transports of heat and salt in the Southern Ocean. Journal of Physical Oceanography, 37(5), 1376–1393.

[26]:

ty_trans = ty_trans.fillna(0)
pot_rho_2 = pot_rho_2.fillna(0)

[27]:

%%time
pot_rho_2 = pot_rho_2.load()
ty_trans = ty_trans.load()

CPU times: user 4.8 s, sys: 4.5 s, total: 9.3 s
Wall time: 10.7 s

[28]:

# choose bins (in this case, default ty_trans_rho pot_rho_2 bins
rho2_bins = ty_trans_rho.potrho.values

# set up a zero array to be filled in by the algorithm
ty_trans_binned = np.zeros((len(rho2_bins), len(pot_rho_2.yu_ocean), len(pot_rho_2.xt_ocean)))

Note: the following cell takes a while.

[29]:

%%time

from tqdm.auto import tqdm

# loop over the bins, performing algorithm by Lee et al. (2007)
# note that it takes time; faster if ty_trans and pot_rho_2 already loaded
for i in tqdm(range(len(rho2_bins)-1)):
    bin_mask = pot_rho_2.where(pot_rho_2 <= rho2_bins[i+1]).where(pot_rho_2 > rho2_bins[i])*0 + 1
    bin_fractions = (rho2_bins[i+1] - pot_rho_2 * bin_mask) / (rho2_bins[i+1] - rho2_bins[i])
    ## bin ty_trans:
    ty_trans_in_lower_bin = (ty_trans * bin_mask * bin_fractions).sum(dim = 'st_ocean')
    ty_trans_binned[i, :, :] += ty_trans_in_lower_bin.fillna(0).values
    del ty_trans_in_lower_bin
    ty_trans_in_upper_bin = (ty_trans * bin_mask * (1 - bin_fractions)).sum(dim = 'st_ocean')
    ty_trans_binned[i+1, :, :] += ty_trans_in_upper_bin.fillna(0).values
    del ty_trans_in_upper_bin

CPU times: user 13min 59s, sys: 11min 51s, total: 25min 50s
Wall time: 16min 39s

[30]:

# convert numpy array into xarray dataarray
ty_trans_binned_array = xr.DataArray(ty_trans_binned,
                                     coords = [rho2_bins, pot_rho_2.yu_ocean, pot_rho_2.xt_ocean],
                                     dims = ['potrho', 'grid_yu_ocean', 'grid_xt_ocean'],
                                     name = 'ty_trans_binned')

[31]:

# find streamfunction
psi_avg_mean_3 = cumsum_from_bottom(ty_trans_binned_array.sum('grid_xt_ocean')) / (1e6 * rho_0)

[32]:

fig, axes = plt.subplots(nrows = 1, ncols = 3, figsize = (15, 4))

levels = np.linspace(-25, 25, 26)

psi_avg.plot.contourf(ax = axes[0], x = 'grid_yu_ocean',levels = levels, add_colorbar = False)

psi_avg_mean_3.plot.contourf(ax = axes[1], x = 'grid_yu_ocean', levels = levels, add_colorbar = False)

p = (psi_avg - psi_avg_mean_3).plot.contourf(ax = axes[2], x = 'grid_yu_ocean', levels = levels, add_colorbar = False)

cbar_ax = fig.add_axes([0.92, 0.15, 0.01, 0.7])
fig.colorbar(p, cax = cbar_ax, label = 'Transport (Sv)')

axes[0].set_title('Residual')
axes[1].set_title('Mean')
axes[2].set_title('Eddy')

for ax in axes:
    ax.set_ylim(1037.5, 1032)
    ax.set_xlim(-70, -35);
../_images/03-Advanced-Recipes_Transformation_from_Depth_to_Potential_Density_51_0.png

The exact numbers resulting from the transformations differ, due to the different numerical methods. xhistogram is exact when compared to snapshots of MOM5 model output, and is hence most appropriate when comparing between online and offline binned quantities. xgcm and the Lee method can still be used to compare quantities binned offline. For example, if we were to instead bin the daily transports onto isopycnals and calculate eddy terms directly from the correlation of velocity and thickness fluctuations, xgcm and the Lee method may be better as the weighted binning onto isopycnals leaves less chance of ‘gaps’ between density layers. The Lee method has no issues with partial cells, but it is much slower than xgcm.

The difference in the mean streamfunction between the three methods is presented below.

[33]:

fig, axes = plt.subplots(nrows = 1, ncols = 3, figsize = (15, 4))

levels = np.arange(-10, 10.1, 0.5)

(psi_avg_mean - psi_avg_mean_2).plot.contourf(ax = axes[0], x = 'grid_yu_ocean',levels = levels, add_colorbar = False)

(psi_avg_mean - psi_avg_mean_3).plot.contourf(ax = axes[1], x = 'grid_yu_ocean', levels = levels, add_colorbar = False)

p = (psi_avg_mean_2 - psi_avg_mean_3).plot.contourf(ax = axes[2], x = 'grid_yu_ocean', levels = levels, add_colorbar = False)

cbar_ax = fig.add_axes([0.92, 0.15, 0.01, 0.7])
fig.colorbar(p, cax=cbar_ax, label='Transport (Sv)')

axes[0].set_title('Method 1 - Method 2')
axes[1].set_title('Method 1 - Method 3')
axes[2].set_title('Method 2 - Method 3')

for ax in axes:
    ax.set_ylim(1037.5, 1032)
    ax.set_xlim(-70, -35)

fig.suptitle('Difference in Mean Streamfunction', fontsize = 16);
../_images/03-Advanced-Recipes_Transformation_from_Depth_to_Potential_Density_53_0.png

An alternative method to calculate the decomposition of the residual meridional overturning circulation into mean and eddy components is to bin ty_trans (or vhrho_nt) into density bins using daily data (both density and transport). Then the transport \(\overline{vh}\) can be separated into a mean component \(\overline{v}\overline{h}\) and an eddy component \(\overline{v^\prime h^\prime}\), where the overline is a time average and primed quantities the deviation from the time average. Quantities are calculated within density layers, and \(h\) is density layer thickness, calculated by binning dzt or dzu. This calculation, since it uses daily data, is computationally expensive and is difficult to do efficiently in a Jupyter notebook. The resulting streamfunctions are similar, but not identical, and may be more intuitive for isopycnal flows. An example of this method can be found at https://github.com/claireyung/Topographic_Hotspots_Upwelling-Paper_Code/blob/main/Figure_Code/Fig5-Overturning.ipynb.