Closed water mass transformation budgets in MOM6

Import packages

%load_ext autoreload
%autoreload 2

import warnings
import matplotlib.pyplot as plt

import numpy as np
import xarray as xr
import xbudget
import regionate
import xwmt
import xwmb
import xgcm

print('xgcm version', xgcm.__version__, '\nxbudget version', xbudget.__version__, '\nregionate version', regionate.__version__, '\nxwmt version', xwmt.__version__, '\nxwmb version', xwmb.__version__)

xgcm version 0.9.0
xbudget version 0.5.1
regionate version 0.5.1
xwmt version 0.1.0
xwmb version 0.5.5

Load example data and model grid

To keep this example dataset reasonably small while still providing accurate results, we have post-processed raw monthly-mean model outputs in the following ways:

• conservatively regridded all diagnostics from their original depth bins to density coordinates

• conservatively regridded from the native eddy-permitting grid (nominally 0.25º) to by a factor of 6x6 grid columns (nominally 1.5º)

• averaged over a year to further reduce the size of the dataset and suppress seasonal variations

from load_example_model_grid import load_MOM6_coarsened_diagnostics
grid = load_MOM6_coarsened_diagnostics()
display(grid)

File 'MOM6_global_example_sigma2_budgets_v0_0_6.nc' already exists at ../data/MOM6_global_example_sigma2_budgets_v0_0_6.nc. Skipping download.

<xgcm.Grid>
X Axis (periodic, boundary='periodic'):
  * center   xh --> outer
  * outer    xq --> center
Y Axis (not periodic, boundary='extend'):
  * center   yh --> outer
  * outer    yq --> center
Z Axis (not periodic, boundary='extend'):
  * center   sigma2_l --> outer
  * outer    sigma2_i --> center

1. Global water mass budget

Collect high-level budget terms

We load in a preset `xbudget <hdrake/xbudget>`__ dictionary that describes the format of the MOM6 mass, heat, and salt budgets and use the collect_budgets utility function to fill out the empty slots in the dictionary with standardized names of terms in the full budget.

import xbudget
xbudget_dict = xbudget.load_preset_budget(model="MOM6")
xbudget.collect_budgets(grid, xbudget_dict)

Compute water mass budget

Terms in the budget are computed according to Drake et al. (2025) for finely spaced density surfaces.

# specific tracer that defines the water mass
lam = "sigma2"

import warnings
with warnings.catch_warnings():
    warnings.simplefilter(action='ignore', category=FutureWarning)

    # Instantiate the water mass budget class
    wmb = xwmb.WaterMassBudget(
        grid,
        xbudget_dict
    )

    wmb.mass_budget(lam, greater_than=True)
    wmt = wmb.wmt.squeeze()
    wmt.load()

Plot global water mass budget

The global water mass budget in density coordinates is given by: \begin{equation} \partial_{t} \mathcal{M}_{\geq} = \mathcal{G}^{(BF)}_{\geq} + \mathcal{G}^{(Mix)}_{\geq} + \mathcal{G}^{(Spu)}_{\geq} + \mathcal{S}_{\geq}, \end{equation} where:

• \(\partial_{t} \mathcal{M}_{\geq}(\sigma_{2}, t)\) quantifies the rate of change of global waters denser than \(\sigma_{2}\)

• \(\mathcal{G}^{(BF)}_{\geq}\) quantifies the rate at which waters lighter than \(\sigma_{2}\) transform into waters denser than \(\sigma_{2}\) due to surface buoyancy fluxes (via heat, salt, or freshwater fluxes)

• \(\mathcal{G}^{(Mix)}_{\geq}\) quantifies the rate at which waters lighter than \(\sigma_{2}\) transform into waters denser than \(\sigma_{2}\) due to surface buoyancy fluxes (via heat, salt, or freshwater fluxes)

• \(\mathcal{G}^{(Spu)}_{\geq}\) quantifies the rate at which waters lighter than \(\sigma_{2}\) transform into waters denser than \(\sigma_{2}\) due to spurious numerical mixing. Note: This term cannot be explicitly computed and is thus indirectly inferred from the residual of all the other terms in the budget.

• \(\mathcal{S}_{\geq}\) quantifies the rate at which surface freshwater fluxes add mass directly to a water mass. Note: This term is generally small and distinct from the effect that freshwater fluxes have on salinity (and hence density), which is instead included in \(\mathcal{G}^{(BF)}_{\geq}\).

plt.figure(figsize=(6,4.5))

kwargs = {"alpha":0.7, "lw":1.5}

plt.plot(-wmt.mass_tendency,                wmt.sigma2_l_target, label=r"$-\partial_{t} \mathcal{M}_{_{\geq}}$",**kwargs)
plt.plot( 0.*wmt.sigma2_l_target,           wmt.sigma2_l_target, label=r"$\Psi_{_{\geq}}$",                   **kwargs)
plt.plot( wmt.mass_source,                  wmt.sigma2_l_target, label=r"$\mathcal{S}_{_{\geq}}$",            **kwargs)
plt.plot( wmt.boundary_fluxes,              wmt.sigma2_l_target, label=r"$\mathcal{G}^{(BF)}_{\geq}$",      **kwargs)
plt.plot( wmt.diffusion,                    wmt.sigma2_l_target, label=r"$\mathcal{G}^{(Mix)}_{\geq}$",     **kwargs)
plt.plot( wmt.spurious_numerical_mixing,    wmt.sigma2_l_target, label=r"$\mathcal{G}^{(Spu)}_{\geq}$", **kwargs)
plt.axhline(34, color="grey", linestyle="dashed")

plt.legend(loc="upper right")
plt.grid(True, alpha=0.25)
plt.ylabel(r"potential density $\sigma_{2}$ [kg/m$^{3}$]");
plt.xlabel("Water mass transformations [kg/s]")
plt.ylim(38, 20);

Scientific Interpretation

Consider the mass budget for \(\sigma_{2} = 34\) \(\text{kg/m}^{3}\), highlighted by the dashed grey line in the figure above:

• Surface buoyancy loss (presumably at high latitudes) form dense waters at a rate of \(\mathcal{G}^{(BF)}_{\geq} \simeq 2.5 \times 10^{10}\) kg/s.

• Interior mixing by various parameterized physical processes destroys dense waters at a rate of \(\mathcal{G}^{(Mix)}_{\geq} \simeq 2.25 \times 10^{10}\) kg/s (by mixing them with lighter waters), nearly balancing the formation rate at the surface.

• A small amount of freshwater is directly added to dense water masses via surface mass fluxes, with a globally-integrated rate of \(\mathcal{S}_{\geq} \simeq 10^{9}\) kg/s.

• Spurious numerical mixing associated with errors in the horizontal advection scheme or vertical advection scheme (Lagrangian remapping-regridding scheme) induce water mass transformations that are generally weaker than those due to parameterized mixing, but are not always of the same sign! We interpret the positive \(\mathcal{G}^{(Spu)}_{\geq}\) at \(\sigma_{2} = 34\) \(\text{kg/m}^{3}\) to reflect spurious numerical entrainment in the overflows associated with North Atlantic Deep Water formation.

• The sum of dense water formation by surface buoyancy fluxes and spurious numerical entrainment is much larger than the rate of destruction by parameterized interior mixing, resulting in an imbalance that causes the mass of dense water to increase at a rate of \(\partial_{t} \mathcal{M}_{\geq} \simeq 10^{10}\) kg/s.

On the attribution of the residual to spurious numerical mixing

While we have here interpreted the residual water mass transformations implied by the non-closure of the water mass budget as spurious numerical mixing, it is possible. However, in Drake et al. (2025) we have performed sensitivity experiments to test several alternative error sources and found them to be either small or random, such that they can not explain the leading-order structure shown in the budget above.

2. Regional water mass budgets

Example 1: specifying a region by its boundary

# Note: the properties of this region are quite different from the rest of the Baltic!
name = "North Atlantic"
lons = np.array([-70, -40,  30,   5])
lats = np.array([ 50,  75,  60,  44])
region = regionate.GriddedRegion(name, lons, lats, grid)

plt.figure(figsize=(7,5.5))
plt.subplot(facecolor=(0.8, 0.8, 0.8))
pc = plt.pcolor(
    grid._ds['geolon_c'],
    grid._ds['geolat_c'],
    grid._ds['deptho'].where(grid._ds['deptho']!=0),
    cmap="Blues",
)
plt.colorbar(pc, label="depth [m]")
plt.plot(region.lons_c, region.lats_c, color="k");
plt.xlabel("longitude");
plt.ylabel("latitude");

import warnings

lam = "sigma2"
with warnings.catch_warnings():
    warnings.simplefilter(action='ignore', category=FutureWarning)
    wmb = xwmb.WaterMassBudget(
        grid,
        xbudget_dict,
        region # specify the region
    )
    wmb.mass_budget(lam, greater_than=True)
    wmt = wmb.wmt.squeeze()
    wmt.load()

plt.figure(figsize=(6,4.5))

kwargs = {"alpha":0.7, "lw":1.5}

plt.plot(-wmt.mass_tendency,                wmt.sigma2_l_target, label=r"$-\partial_{t} \mathcal{M}_{_{\geq}}$",**kwargs)
plt.plot( wmt.convergent_mass_transport,    wmt.sigma2_l_target, label=r"$\Psi_{_{\geq}}$",                   **kwargs)
plt.plot( wmt.mass_source,                  wmt.sigma2_l_target, label=r"$\mathcal{S}_{_{\geq}}$",            **kwargs)
plt.plot( wmt.boundary_fluxes,              wmt.sigma2_l_target, label=r"$\mathcal{G}^{(BF)}_{\geq}$",      **kwargs)
plt.plot( wmt.diffusion,                    wmt.sigma2_l_target, label=r"$\mathcal{G}^{(Mix)}_{\geq}$",     **kwargs)
plt.plot( wmt.spurious_numerical_mixing,    wmt.sigma2_l_target, label=r"$\mathcal{G}^{(Spu)}_{\geq}$", **kwargs)
plt.axhline(wmt.convergent_mass_transport.idxmin(), color="grey", linestyle="dashed", label="$\sigma_{2}$ of $\min(\Psi_{\geq})$ (AMOC)")

plt.legend(loc="upper right")
plt.grid(True, alpha=0.25)
plt.ylabel(r"potential density $\sigma_{2}$ [kg/m$^{3}$]");
plt.xlabel("Water mass transformations [kg/s]")
plt.ylim(37.2, 34);

Example 2: specifying a region by its interior mask

# Note: the properties of this region are quite different from the rest of the Baltic!
name = "Southern Ocean"
mask = grid._ds.geolat < -30.
regions = regionate.MaskRegions(mask, grid, name) # returns a dictionary of all distinct region objects
region = regions.region_dict[0] # in this case, there is only one contiguous region

plt.figure(figsize=(7,5.5))
plt.subplot(facecolor=(0.8, 0.8, 0.8))
pc = plt.pcolor(
    grid._ds['geolon_c'],
    grid._ds['geolat_c'],
    grid._ds['deptho'].where(grid._ds['deptho']!=0),
    cmap="Blues",
)
plt.colorbar(pc, label="depth [m]")
plt.plot(region.lons_c, region.lats_c, color="k");
plt.xlabel("longitude");
plt.ylabel("latitude");

import warnings

lam = "sigma2"
with warnings.catch_warnings():
    warnings.simplefilter(action='ignore', category=FutureWarning)
    wmb = xwmb.WaterMassBudget(
        grid,
        xbudget_dict,
        region # specify the region
    )
    wmb.mass_budget(lam, greater_than=True)
    wmt = wmb.wmt.squeeze()
    wmt.load()

plt.figure(figsize=(6,4.5))

kwargs = {"alpha":0.7, "lw":1.5}

plt.plot(-wmt.mass_tendency,                wmt.sigma2_l_target, label=r"$-\partial_{t} \mathcal{M}_{_{\geq}}$",**kwargs)
plt.plot( wmt.convergent_mass_transport,    wmt.sigma2_l_target, label=r"$\Psi_{_{\geq}}$",                   **kwargs)
plt.plot( wmt.mass_source,                  wmt.sigma2_l_target, label=r"$\mathcal{S}_{_{\geq}}$",            **kwargs)
plt.plot( wmt.boundary_fluxes,              wmt.sigma2_l_target, label=r"$\mathcal{G}^{(BF)}_{\geq}$",      **kwargs)
plt.plot( wmt.diffusion,                    wmt.sigma2_l_target, label=r"$\mathcal{G}^{(Mix)}_{\geq}$",     **kwargs)
plt.plot( wmt.spurious_numerical_mixing,    wmt.sigma2_l_target, label=r"$\mathcal{G}^{(Spu)}_{\geq}$", **kwargs)
sigma2_amoc = wmt.convergent_mass_transport.idxmax().compute()
sigma2_smoc = wmt.convergent_mass_transport.sel(sigma2_l_target=slice(sigma2_amoc, None)).idxmin().compute()
plt.axhline(sigma2_amoc, color="grey", linestyle="dashed", label="$\sigma_{2}$ of $\max(\Psi_{\geq})$ (AMOC)")
plt.axhline(sigma2_smoc, color="grey", linestyle="dotted", label="$\sigma_{2}$ of $\max(\Psi_{\geq})$ (SMOC)")


plt.legend(loc="upper right")
plt.grid(True, alpha=0.25)
plt.ylabel(r"potential density $\sigma_{2}$ [kg/m$^{3}$]");
plt.xlabel("Water mass transformations [kg/s]")
plt.ylim(37.5, 32);

3. More granular control over transformation budgets

We now repeat the same calculation, with the only differences being that we use the following two keyword arguments to further break down the budget.

Example 1. Decomposing surface heat flux contributions

To decompose the nonadvective surface transformations by term, we pass decompose = ["surface_exchange_flux", "nonadvective"] to xwmb.WaterMassBudget

import warnings

lam = "sigma2"
with warnings.catch_warnings():
    warnings.simplefilter(action='ignore', category=FutureWarning)
    wmb = xwmb.WaterMassBudget(
        grid,
        xbudget_dict,
        region,
        decompose = ["surface_exchange_flux", "nonadvective"]
    )
    wmb.mass_budget(lam, greater_than=True)
    wmt = wmb.wmt.squeeze()
    wmt.load()

plt.figure(figsize=(6,4.5))

kwargs = {"alpha":0.7, "lw":1.5}

plt.plot( wmt.boundary_fluxes,              wmt.sigma2_l_target, label=r"$\mathcal{G}^{(BF)}_{\geq}$", color="k", lw=2)

# Plot transformations due to each contributing heat flux term
heat_fluxes = ["latent", "sensible", "longwave", "shortwave"]
heat_flux_dict = {
    heat_flux: f"surface_exchange_flux_nonadvective_{heat_flux}_heat"
    for heat_flux in heat_fluxes
}
plt.plot( wmt[heat_flux_dict["latent"]],              wmt.sigma2_l_target, label=r"$\mathcal{G}^{(LATE)}_{\geq}$",      **kwargs)
plt.plot( wmt[heat_flux_dict["sensible"]],              wmt.sigma2_l_target, label=r"$\mathcal{G}^{(SENS)}_{\geq}$",      **kwargs)
plt.plot( wmt[heat_flux_dict["longwave"]],              wmt.sigma2_l_target, label=r"$\mathcal{G}^{(LW)}_{\geq}$",      **kwargs)
plt.plot( wmt[heat_flux_dict["shortwave"]],              wmt.sigma2_l_target, label=r"$\mathcal{G}^{(SW)}_{\geq}$",      **kwargs)

sigma2_smoc = wmt.convergent_mass_transport.sel(sigma2_l_target=slice(sigma2_amoc, None)).idxmin().compute()
plt.axhline(sigma2_smoc, color="grey", linestyle="dotted", label="$\sigma_{2}$ of $\max(\Psi_{\geq})$ (SMOC)")


plt.legend(loc="upper left")
plt.grid(True, alpha=0.25)
plt.ylabel(r"potential density $\sigma_{2}$ [kg/m$^{3}$]");
plt.xlabel("Water mass transformations [kg/s]")
plt.ylim(37.5, 32);

Example 2. Spatial patterns

To compute spatial patterns of transformation rates, rather than the default of area-integrating, we pass integrate=False to xwmb.WaterMassBudgetmass_budget.

import warnings

lam = "sigma2"
with warnings.catch_warnings():
    warnings.simplefilter(action='ignore', category=FutureWarning)
    wmb = xwmb.WaterMassBudget(
        grid,
        xbudget_dict,
        region,
        decompose = ["surface_exchange_flux", "nonadvective"]
    )
    wmb.mass_budget(lam, greater_than=True, integrate=False)
    wmt = wmb.wmt.squeeze()
    wmt.load()

Supposing we are interested in how these fluxes contribute to the formation of dense waters, we select the isopycnal of maximum SMOC transport.

def plot_smoc_transformation_map(term, ax):
    pc = (wmt[heat_flux_dict[term]]*1e-6).sel(sigma2_l_target=sigma2_smoc).plot(
        ax=ax, x="geolon", y="geolat",
        cmap="RdBu_r", vmin=-10, vmax=10
    )
    (grid._ds.wet == 0).plot.contour(ax=ax, x="geolon", y="geolat", colors="k", levels=[0])
    pc.colorbar.set_label(r"WMT [10$^{6}$ kg/s/grid column]")
    ax.set_title(f"Dense water formation due to {term} heat flux")
    ax.set_ylim(-80, -50)

fig, axes = plt.subplots(2,2, figsize=(12, 6))
for ax, term in zip(axes.flatten(), heat_flux_dict.keys()):
    plot_smoc_transformation_map(term, ax)
plt.tight_layout()