Calculating the mean Eulerian salinity budget in ECCOv4r4#
Carenza Williams & Wenrui Jiang, Feb 2026
This notebook demonstrates how to calculate the time-mean Eulerian salinity budget in ECCOv4r4, with a focus on the advection term. The method illustrated in this notebook can be found in our paper titled “Tracer Budgets on Lagrangian Trajectories”.
This is No.1 of a serie of 2 tutorial notebooks. The second one can be found here.
0.Set up#
Let’s start by loading in our dataset and some key packages.
Show code cell source
# -----------------
# Import packages
# -----------------
import warnings
import numpy as np
import xarray as xr
import seaduck as sd
import seaduck.eulerian_budget as sdeb
# Suppress futurewarning
warnings.filterwarnings("ignore", category=FutureWarning,
message="elementwise comparison failed")
# Suppress xgcm padding warning
warnings.filterwarnings("ignore", category=UserWarning,
message="rename 'Z' to 'Z' does not create an index anymore")
# -----------------------------
# Load and inspect the dataset
# -----------------------------
ds_bar = sd.utils.get_dataset('eul_bud_mean')
# Create cell volume
vol = (ds_bar.rA*ds_bar.drF*ds_bar.hFacC).transpose('face','Z','Y','X')
ds_bar['Vol'] = vol
ds_bar
<xarray.Dataset> Size: 368MB
Dimensions: (Z: 50, face: 13, Y: 90, Xp1: 90, Yp1: 90, X: 90, Zl: 50)
Coordinates: (12/20)
* Z (Z) float32 200B -5.0 -15.0 -25.0 ... -5.461e+03 -5.906e+03
PHrefC (Z) float32 200B dask.array<chunksize=(10,), meta=np.ndarray>
drF (Z) float32 200B dask.array<chunksize=(10,), meta=np.ndarray>
k (Z) int64 400B dask.array<chunksize=(10,), meta=np.ndarray>
* face (face) int64 104B 0 1 2 3 4 5 6 7 8 9 10 11 12
* Y (Y) int64 720B 0 1 2 3 4 5 6 7 8 ... 82 83 84 85 86 87 88 89
... ...
XC (face, Y, X) float32 421kB dask.array<chunksize=(13, 45, 45), meta=np.ndarray>
YC (face, Y, X) float32 421kB dask.array<chunksize=(13, 45, 45), meta=np.ndarray>
hFacC (Z, face, Y, X) float32 21MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
maskC (Z, face, Y, X) int8 5MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
rA (face, Y, X) float32 421kB dask.array<chunksize=(13, 45, 45), meta=np.ndarray>
* Zl (Zl) float32 200B 0.0 -10.0 -20.0 ... -5.244e+03 -5.678e+03
Data variables: (12/20)
ADVx_SLT_bar (Z, face, Y, Xp1) float32 21MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
ADVy_SLT_bar (Z, face, Yp1, X) float32 21MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
ADVz_SLT_bar (Zl, face, Y, X) float32 21MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
DFx_SLT_bar (Z, face, Y, Xp1) float32 21MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
DFy_SLT_bar (Z, face, Yp1, X) float32 21MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
DFz_SLT_bar (Zl, face, Y, X) float32 21MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
... ...
spdivup_bar (Z, face, Y, X) float32 21MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
tendSln_bar (Z, face, Y, X) float32 21MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
u_bar (Z, face, Y, Xp1) float32 21MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
v_bar (Z, face, Yp1, X) float32 21MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
w_bar (Zl, face, Y, X) float32 21MB dask.array<chunksize=(10, 13, 45, 45), meta=np.ndarray>
Vol (face, Z, Y, X) float32 21MB dask.array<chunksize=(13, 10, 45, 45), meta=np.ndarray>This dataset has the required variables to calculate the mean salinity budget, as well as some ECCOv4r4 grid variables. These are:
S_bar –> time-mean salinity [PSU]
T_bar –> time-mean potential temperature [°C]
u_bar, v_bar, w_bar –> time-mean total velocity m s [\(^{-1}\)]
ADVx_SLT_bar (etc.) –> time-mean advective flux of salinity [PSU m\(^{3}\) s\(^{-1}\)]
DFx_SLT_bar (etc.) –> time-mean diffusive flux of salinity [PSU m\(^{3}\) s\(^{-1}\)]
forcS_bar –> time-mean salt forcing [PSU s\(^{-1}\)]
forcFW_bar –> time-mean freshwater forcing [PSU s\(^{-1}\)]
tendSln_bar –> time-mean salinity tendency [PSU s\(^{-1}\)]
spdivup_bar –> \(\overline{S'\nabla u'}\) [PSU s\(^{-1}\)]
The variables in the first five lines can be directly downloaded from the ECCOv4r4 website. Since this notebook focuses on the advective process, the Eulerian tendency, the eddy transport, and the forcing terms are precalculated for brevity. Interested readers are directed to Piecuch 2017 for an instructive guide on the calculation of those terms.
# Grid object
grid = sdeb.create_ecco_grid(ds_bar,for_outer = True)
1. Calculate the wall salinity#
A key step to bridge the difference between Eulerian and Lagrangian frameworks is to calculate the “wall salinity”. This involves interpolating the model salinity onto the grid walls so that salinity and velocity can be stored at the same point on the model grid.
First, we will create a Topology object using a seaduck built-in function, which helps us navigate the complex tile connections of ECCO.
# Topology object
tp = sd.Topology(ds_bar)
The ECCO use third order upwind and third order direct space-time scheme for salinity advection in the vertical and horizontal, respectively. third_order_upwind_z and third_order_DST_x/y are built-in seaduck functions that mimic the behavior of those advection schemes.
# ----------------------------------------------
# Create wall salinity (uses seaduck functions)
# ----------------------------------------------
lm = 2 # The extra margin on the left side
rm = 1 # and that on the right side
dx = np.array(ds_bar.dxG)
dy = np.array(ds_bar.dyC)
ds = ds_bar
w = np.array(ds.w_bar)
s = np.array(ds.S_bar.where(ds.maskC!=0))
sz = sdeb.third_order_upwind_z(s,w)# The vertical advection scheme
u = np.array(ds.u_bar)
sx = np.zeros_like(s)
for face in range(13):
xbuffer = sdeb.buffer_x_withface(s,face,lm,rm,tp)
u_cfl = np.array(u[...,face,:,:]/dx[face])
sx[:,face,:,:] = sdeb.third_order_DST_x(xbuffer,u_cfl)
v = np.array(ds_bar.v_bar)
sy = np.zeros_like(s)
for face in range(13):
u_cfl = np.array(v[...,face,:,:]/dy[face])
ybuffer = sdeb.buffer_y_withface(s,face,lm,rm,tp)
sy[:,face,:,:] = sdeb.third_order_DST_y(ybuffer,u_cfl)
ds_bar['sx_bar'] = xr.DataArray(sx.reshape(50,13,90,90),dims = ('Z','face','Y','Xp1'))
ds_bar['sy_bar'] = xr.DataArray(sy.reshape(50,13,90,90),dims = ('Z','face','Yp1','X'))
ds_bar['sz_bar'] = xr.DataArray(sz.reshape(50,13,90,90),dims = ('Zl','face','Y','X'))
Now our dataset also contains the time-mean salinity stored on the x, y and z faces of the model grid. This will be required to calculate the advection term of the salinity budget.
2. Calculating the Salinity Budget#
Now, onto the main event: finding the salinity budget.
The time-mean salinity budget can be expressed as:
where \(\overline{u'\cdot\nabla s'}\) is calculated by:
Terms are defined as:
\(\frac{\partial{\overline{S}}}{\partial{t}}\) –> Eulerian tendency
\(\left( -\overline{u}\cdot\nabla{\overline{S}} \right)\) –> (Mean) advection
\(- \overline{u'\cdot\nabla s'}\) –> Eddy transport
\(-\nabla{\overline{DIF}_{x,y,z}}\) –> Diffusion of salinity
\(\overline{F}_{FW}\) –> Freshwater forcing
\(\overline{F}_{S}\) –> Salt forcing (from run off, ice plume, etc.)
All terms but the advection of salinity are stored in the dataset ‘ds_bar’. We must calculate the advection of salinity term, which we will do one term at a time below.
First, we calculate volume transport from velocity.
Show code cell source
# ---------------------------
# Define transport function
# ---------------------------
def vel_to_trans(u, v, w):
return (
u * ds_bar.drF * ds_bar.dyG,
v * ds_bar.drF * ds_bar.dxG,
w * ds_bar.rA
)
# Calculate volume transport
ut, vt, wt = vel_to_trans(ds_bar.u_bar, ds_bar.v_bar, ds_bar.w_bar)
ds_bar['utrans'] = ut.compute()
ds_bar['vtrans'] = vt.compute()
ds_bar['wtrans'] = wt.compute()
ds_bar['wtrans'][0,:] = 0
We can now calculate the mean component of the salinity fluxes
ds_bar['ubarsbar_x'] = (ut * ds_bar.sx_bar).compute()
ds_bar['ubarsbar_y'] = (vt * ds_bar.sy_bar).compute()
ds_bar['ubarsbar_z'] = (wt * ds_bar.sz_bar).compute()
ds_bar['ADV_x'] = ds_bar.ADVx_SLT_bar.compute()
ds_bar['ADV_y'] = ds_bar.ADVy_SLT_bar.compute()
ds_bar['ADV_z'] = ds_bar.ADVz_SLT_bar.compute()
ds_bar['DIF_x'] = (ds_bar.DFx_SLT_bar).compute()
ds_bar['DIF_y'] = (ds_bar.DFy_SLT_bar).compute()
ds_bar['DIF_z'] = (ds_bar.DFz_SLT_bar).compute()
Pin xgcm=0.6.1
The reason we do .compute() on those fluxes is because of certain breaking changes in xgcm as documented here. We find that xgcm=0.6.1 is a good version to do most ECCO related analysis.
To calculate the divergence of the fluxes, we first create a OceData object
tub = sd.OceData(ds_bar)
Calculate \(\nabla\cdot{(\overline{u}\overline{S})}\)#
# -------------------------------
# Calculate div (uS)
# -------------------------------
divus = sdeb.total_div(tub, grid, 'ubarsbar_x', 'ubarsbar_y', 'ubarsbar_z')
Calculate \(\overline{S}\cdot\nabla{\overline{u}}\)#
# -------------------------------
# Calculate div (uS)
# -------------------------------
divu = sdeb.total_div(tub, grid, 'utrans', 'vtrans', 'wtrans')
# Multiply by mean salinity
# -------------------------------
# Calculate the term -- S div (u)
# -------------------------------
sdivu = ds_bar.S_bar * divu
Now the above terms can be combined to give the required \(\overline{u}\cdot\nabla{\overline{S}}\) term:
Calculate \(\large\overline{u}\cdot\nabla{\overline{S}} = \nabla\cdot{(\overline{u}\overline{S})} - \overline{S}\cdot\nabla{\overline{u}}\)#
# -------------------------------
# Calculate u dot grad(S)
# -------------------------------
ugrads = divus - sdivu
Calculate the divergence of ADV and DIF terms#
# -------------------------------
# Calculate div(ADV)
# -------------------------------
divADV = sdeb.total_div(tub, grid, 'ADV_x', 'ADV_y', 'ADV_z')
# -------------------------------
# Calculate -div(DIF)
# -------------------------------
dif_h = -sdeb.hor_div(tub, grid, 'DIF_x', 'DIF_y')
dif_v = -sdeb.ver_div(tub, grid, 'DIF_z')
3. Check the budget closure#
# ------------------------------
# Calculate the advection term
# ------------------------------
adv = ugrads
# ------------------------------
# Calculate the eddy transport term
# ------------------------------
neg_upgradsp_bar = -(divADV - divus - ds_bar.spdivup_bar)
# ------------------------------
# Calculate the total forcing term
# ------------------------------
forc = ((-ds_bar.forcFW_bar) + ds_bar.forcS_bar)
Once again, we are trying to close this budget:
If we have done this correctly, then the LHS - RHS of the budget equation should equal 0 (to machine precision). Let’s check!
# ------------------------------
# Check that the budget closes
# ------------------------------
res = (ds_bar.tendSln_bar + adv) - (neg_upgradsp_bar + dif_h + dif_v + forc)
# ----------------------------------------------
# Plot the residuals (around the UK and Europe)
# ----------------------------------------------
res[0,2].plot()
<matplotlib.collections.QuadMesh at 0x7f36c0d52ad0>
Success! The fact that the residual is very small and look random suggests that it arise solely from round-off error.
The only thing left to do now is to save our work…
Show code cell source
out = ds[['sx_bar','sy_bar','sz_bar','u_bar','v_bar','w_bar','S_bar','T_bar','dxG','dyG']]
out['neg_upgradsp_bar'] = neg_upgradsp_bar
out['forc_s'] = forc
out['dif_h'] = dif_h
out['dif_v'] = dif_v
out['conv_us'] = -divus
out['tend_s'] = ds_bar.tendSln_bar
out.to_zarr('lag_budg.zarr',mode = 'w')
<xarray.backends.zarr.ZarrStore at 0x7f36d0a37740>
… and now we can use these eulerian budget terms to calculate the Lagrangian budget in the next notebook.