Salt, Salinity and Freshwater Budgets¶
Contributors: Jan-Erik Tesdal, Ryan Abernathey and Ian Fenty
A major part of this tutorial is based on “A Note on Practical Evaluation of Budgets in ECCO Version 4 Release 3” by Christopher G. Piecuch (https://ecco.jpl.nasa.gov/drive/files/Version4/Release3/doc/v4r3_budgets_howto.pdf). Calculation steps and Python code presented here are converted from the MATLAB code presented in the above reference.
Objectives¶
This tutorial will go over three main budgets which are all related:
- Salt budget
- Salinty budget
- Freshwater budget
We will describe the governing equations for the conservation for both salt, salinity and freshwater content and discuss the subtle differences one needs to be aware of when closing budgets of salt and freshwater content (extensive quantities) versus the budget of salinity (an intensive quantity) in ECCOv4.
Introduction¶
The general form for the salt/salinity budget can be formulated in the same way as with the heat budget, where instead of potential temperature (), the budget is described with salinity (
).
The total tendency () is equal to advective convergence (
), diffusive flux convergence (
) and a forcing term
.
In the case of ECCOv4, salt is strictly a conserved mass and can be described as
The change in salt content over time () is equal to the convergence of the advective flux (
) and diffusive flux (
) plus a forcing term associated with surface salt exchanges (
). As with the heat budget, we present both the horizontal (
) and
vertical (
) components of the advective term. Again, we have
as the “residual mean” velocities, which contain both the resolved (Eulerian) and parameterizing “GM bolus” velocities. Also note the use of the rescaled height coordinate
and the scale factor
which have been described in the volume and
heat budget tutorials.
The salt budget in ECCOv4 only considers the mass of salt in the ocean. Thus, the convergence of freshwater and surface freshwater exchanges are not formulated specifically. An important point here is that, given the nonlinear free surface condition in ECCOv4, budgets for salt content (an extensive quantity) are not the same as budgets for salinity (an intensive quantity). In order to accurately describe variation in salinity, we need to take into account the variation of both salt and volume.
Using the product rule, (i.e., the left side of the salt budget equation) can be extended as follows
When substituting with the right hand side of the above equation, we can solve for the salinity tendency (
):
Since we can define the temporal change in
as
This constitutes the conservation of volume in ECCOv4, which can be formulated as
You can read more about volume conservation and the coordinate system in another tutorial.
denotes the volumetric surface fluxes and can be decomposed into net atmospheric freshwater fluxes (i.e., precipitation minus evaporation,
), continental runoff (
) and exchanges due to sea ice melting/formation (
). Here
and
are the resolved horizontal and vertical velocities, respectively.
Thus, the conservation of salinity in ECCOv4 can be described as
Notice here that, in contrast to the salt budget equation, the salinity equation explicitly includes the surface forcing ().
represents surface freshwater exchanges (
) and
represents surface salt fluxes (i.e., addition/removal of salt). Besides the convergence of the advective flux (
), the salinity equation also includes the convergence of the volume flux multiplied by the
salinity (
), which accounts for the concentration/dilution effect of convergent/divergent volume flux.
The (liquid) freshwater content is defined here as the volume of freshwater (i.e., zero-salinity water) that needs to be added (or subtracted) to account for the deviation between salinity from a given reference salinity
. Thus, within a control volume
the freshwater content is defined as a volume (
):
Similar to the salt and salinity budgets, the total tendency (i.e., change in freshwater content over time) can be expressed as the sum of the tendencies due to advective convergence, diffusive convergence, and forcing:
Prepare environment and load ECCOv4 diagnostic output¶
Import relevant Python modules¶
[1]:
import numpy as np
import xarray as xr
[2]:
# Suppress warning messages for a cleaner presentation
import warnings
warnings.filterwarnings('ignore')
[3]:
## Import the ecco_v4_py library into Python
## =========================================
## -- If ecco_v4_py is not installed in your local Python library,
## tell Python where to find it.
#import sys
#sys.path.append('/Users/jet/git/ECCOv4-py')
import ecco_v4_py as ecco
[4]:
# Plotting
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid
from cartopy.mpl.geoaxes import GeoAxes
import cartopy
%matplotlib inline
Add relevant constants¶
[5]:
# Seawater density (kg/m^3)
rhoconst = 1029
## needed to convert surface mass fluxes to volume fluxes
Load ecco_grid¶
[6]:
## Set top-level file directory for the ECCO NetCDF files
## =================================================================
# Define main directory
base_dir = '/work/noaa/gfdlscr/jtesdal/ECCOv4-release'
# Define a high-level directory for ECCO fields
ECCO_dir = base_dir + '/Release3_alt'
**Note**: Change ``base_dir`` to your own directory path.
[7]:
# Load the model grid
grid_dir = ECCO_dir + '/nctiles_grid/'
ecco_grid = ecco.load_ecco_grid_nc(grid_dir, 'ECCOv4r3_grid.nc')
Volume¶
Calculate the volume of each grid cell. This is used when converting advective and diffusive flux convergences and calculating volume-weighted averages.
[8]:
# Volume (m^3)
vol = (ecco_grid.rA*ecco_grid.drF*ecco_grid.hFacC).transpose('tile','k','j','i')
Load monthly snapshots¶
[9]:
data_dir= ECCO_dir + '/nctiles_monthly_snapshots'
year_start = 1993
year_end = 2017
# Load one extra year worth of snapshots
ecco_monthly_snaps = ecco.recursive_load_ecco_var_from_years_nc(data_dir, \
vars_to_load=['ETAN','SALT'],\
years_to_load=range(year_start, year_end+1))
num_months = len(ecco_monthly_snaps.time.values)
# Drop the last 11 months so that we have one snapshot at the beginning and end of each month within the
# range 1993/1/1 to 2015/1/1
ecco_monthly_snaps = ecco_monthly_snaps.isel(time=np.arange(0, num_months-11))
loading files of ETAN
loading files of SALT
[10]:
# Drop superfluous coordinates (We already have them in ecco_grid)
ecco_monthly_snaps = ecco_monthly_snaps.reset_coords(drop=True)
Load monthly mean data¶
[11]:
data_dir= ECCO_dir + '/nctiles_monthly'
# Find the record of the last snapshot
## This is used to defined the exact period for monthly mean data
year_end = ecco_monthly_snaps.time.dt.year.values[-1]
ecco_monthly_mean = ecco.recursive_load_ecco_var_from_years_nc(data_dir, \
vars_to_load=['SFLUX', 'oceSPtnd', 'ADVx_SLT', 'ADVy_SLT', 'ADVr_SLT',
'DFxE_SLT', 'DFyE_SLT', 'DFrE_SLT', 'DFrI_SLT', 'oceFWflx',
'UVELMASS', 'VVELMASS', 'WVELMASS', 'GM_PsiX', 'GM_PsiY',
'ETAN', 'SALT'], years_to_load=range(year_start, year_end))
loading files of ADVr_SLT
loading files of ADVx_SLT
loading files of ADVy_SLT
loading files of DFrE_SLT
loading files of DFrI_SLT
loading files of DFxE_SLT
loading files of DFyE_SLT
loading files of ETAN
loading files of GM_PsiX
loading files of GM_PsiY
loading files of SALT
loading files of SFLUX
loading files of UVELMASS
loading files of VVELMASS
loading files of WVELMASS
loading files of oceFWflx
loading files of oceSPtnd
[12]:
# Drop superfluous coordinates (We already have them in ecco_grid)
ecco_monthly_mean = ecco_monthly_mean.reset_coords(drop=True)
Merge dataset of monthly mean and snapshots data¶
Merge the two datasets to put everything into one single dataset
[13]:
ds = xr.merge([ecco_monthly_mean,
ecco_monthly_snaps.rename({'time':'time_snp','ETAN':'ETAN_snp', 'SALT':'SALT_snp'})])
Predefine coordinates for global regridding of the ECCO output (used in resample_to_latlon
)¶
[14]:
new_grid_delta_lat = 1
new_grid_delta_lon = 1
new_grid_min_lat = -90+new_grid_delta_lat/2
new_grid_max_lat = 90-new_grid_delta_lat/2
new_grid_min_lon = -180+new_grid_delta_lon/2
new_grid_max_lon = 180-new_grid_delta_lon/2
Create the xgcm ‘grid’ object¶
[15]:
# Change time axis of the snapshot variables
ds.time_snp.attrs['c_grid_axis_shift'] = 0.5
[16]:
grid = ecco.get_llc_grid(ds)
Number of seconds in each month¶
The xgcm grid
object includes information on the time axis, such that we can use it to get , which is the time span between the beginning and end of each month (in seconds).
[17]:
delta_t = grid.diff(ds.time_snp, 'T', boundary='fill', fill_value=np.nan)
[18]:
# convert to seconds
delta_t = delta_t.astype('f4') / 1e9
Evaluating the salt budget¶
We will evalute each term in the above salt budget
The total tendency of salt () is the sum of the salt tendencies from advective convergence (
), diffusive heat convergence (
) and total forcing (
).
We present calculations sequentially for each term starting with which will be derived by differencing instantaneous monthly snapshots of
SALT
. The terms on the right hand side of the heat budget are derived from monthly-averaged fields.
Total salt tendency¶
We calculate the monthly-averaged time tendency of SALT
by differencing monthly SALT
snapshots. Remember that we need to include a scale factor due to the nonlinear free surface formulation. Thus, we need to use snapshots of both ETAN
and SALT
to evaluate .
[19]:
# Calculate the s*S term
sSALT = ds.SALT_snp*(1+ds.ETAN_snp/ecco_grid.Depth)
[20]:
# Total tendency (psu/s)
G_total_Slt = grid.diff(sSALT, 'T', boundary='fill', fill_value=0.0)/delta_t
The nice thing is that now the time values of (
G_total_Slt
) line up with the time values of the time-mean fields (middle of the month)
Advective salt convergence¶
Horizontal advective salt convergence¶
The relevant fields from the diagnostic output here are - ADVx_SLT
: U Component Advective Flux of Salinity (psu m^3/s) - ADVy_SLT
: V Component Advective Flux of Salinity (psu m^3/s)
The xgcm grid
object is then used to take the convergence of the horizontal heat advection.
[21]:
ADVxy_diff = grid.diff_2d_vector({'X' : ds.ADVx_SLT, 'Y' : ds.ADVy_SLT}, boundary = 'fill')
# Convergence of horizontal advection (psu m^3/s)
adv_hConvS = (-(ADVxy_diff['X'] + ADVxy_diff['Y']))
Vertical advective salt convergence¶
The relevant field from the diagnostic output is - ADVr_SLT
: Vertical Advective Flux of Salinity (psu m^3/s)
[22]:
# Load monthly averages of vertical advective flux
ADVr_SLT = ds.ADVr_SLT.transpose('time','tile','k_l','j','i')
**Note**: For ``ADVr_SLT``, ``DFrE_SLT`` and ``DFrI_SLT``, we need to make sure that sequence of dimensions are consistent. When loading the fields use ``.transpose('time','tile','k_l','j','i')``. Otherwise, the divergences will be not correct (at least for ``tile = 12``).
[23]:
# Convergence of vertical advection (psu/s)
adv_vConvS = grid.diff(ADVr_SLT, 'Z', boundary='fill')
**Note**: In case of the volume budget (and salinity conservation), the surface forcing (``oceFWflx``) is already included at the top level (``k_l = 0``) in ``WVELMASS``. Thus, to keep the surface forcing term explicitly represented, one needs to zero out the values of ``WVELMASS`` at the surface so as to avoid double counting (see ``ECCO_v4_Volume_budget_closure.ipynb``). The salt budget only balances when the sea surface forcing is not added to the vertical salt flux (at the air-sea
interface).
Total advective salt convergence¶
We can get the total convergence by simply adding the horizontal and vertical component.
[26]:
# Sum horizontal and vertical convergences and divide by volume (psu/s)
G_advection_Slt = (adv_hConvS + adv_vConvS)/vol
Diffusive salt convergence¶
Horizontal diffusive salt convergence¶
The relevant fields from the diagnostic output here are - DFxE_SLT
: U Component Diffusive Flux of Salinity (psu m^3/s) - DFyE_SLT
: V Component Diffusive Flux of Salinity (psu m^3/s)
As with advective fluxes, we use the xgcm grid
object to calculate the convergence of horizontal salt diffusion.
[27]:
DFxyE_diff = grid.diff_2d_vector({'X' : ds.DFxE_SLT, 'Y' : ds.DFyE_SLT}, boundary = 'fill')
# Convergence of horizontal diffusion (psu m^3/s)
dif_hConvS = (-(DFxyE_diff['X'] + DFxyE_diff['Y']))
Vertical diffusive salt convergence¶
The relevant fields from the diagnostic output are - DFrE_SLT
: Vertical Diffusive Flux of Salinity (Explicit part) (psu m^3/s) - DFrI_SLT
: Vertical Diffusive Flux of Salinity (Implicit part) (psu m^3/s) > Note: Vertical diffusion has both an explicit (DFrE_SLT
) and an implicit (DFrI_SLT
) part.
[28]:
# Load monthly averages of vertical diffusive fluxes
DFrE_SLT = ds.DFrE_SLT.transpose('time','tile','k_l','j','i')
DFrI_SLT = ds.DFrI_SLT.transpose('time','tile','k_l','j','i')
# Convergence of vertical diffusion (psu m^3/s)
dif_vConvS = grid.diff(DFrE_SLT, 'Z', boundary='fill') + grid.diff(DFrI_SLT, 'Z', boundary='fill')
Total diffusive salt convergence¶
[29]:
# Sum horizontal and vertical convergences and divide by volume (psu/s)
G_diffusion_Slt = (dif_hConvS + dif_vConvS)/vol
Salt forcing¶
There are two relevant model diagnostics: - SFLUX
: total salt flux (match salt-content variations) (g/m^2/s) - oceSPtnd
: salt tendency due to salt plume flux (g/m^2/s)
The local forcing term reflects surface salt exchanges. There are two relevant model diagnostics here, namely the total salt exchange at the surface (SFLUX
), which is nonzero only when sea ice melts or freezes, and the salt plume tendency (oceSPtnd
), which vertically redistributes surface salt input by sea ice formation. We will merge SFLUX
and oceSPtnd
into a single data array (forcS
) and convert it to units of psu per second.
[30]:
# Load SFLUX and add vertical coordinate
SFLUX = ds.SFLUX.assign_coords(k=0).expand_dims('k').transpose('time','tile','k','j','i')
# Calcualte forcing term by adding SFLUX and oceSPtnd (g/m^2/s)
forcS = xr.concat([SFLUX+ds.oceSPtnd,ds.oceSPtnd[:,:,1:]], dim='k')
SFLUX
and oceSPtnd
is given in g/m^2/s. Dividing by density and corresponding vertical length scale (drF
) results in g/kg/s, which is the same as psu/s.
[31]:
# Forcing (psu/s)
G_forcing_Slt = forcS/rhoconst/(ecco_grid.hFacC*ecco_grid.drF)
Salt budget: Map of residual¶
[40]:
# Total convergence (psu/s)
ConvSlt = G_advection_Slt + G_diffusion_Slt
# Sum of terms in RHS of equation (psu/s)
rhs = ConvSlt + G_forcing_Slt
[41]:
# Accumulated residual
resSlt = (rhs-G_total_Slt).sum(dim='k').sum(dim='time').compute()
[42]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, resSlt,
cmin=-1e-9, cmax=1e-9, show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Accumulated residual (RHS - LHS) [psu s$^{-1}$]', fontsize=16)
plt.show()

The above map confirms that the residual (summed over depth and time) is essentially zero everywhere, and the ECCOv4 salt budget can be closed to machine precision.
Salt budget: Spatial distributions¶
[43]:
# In order to plot the budget terms in one figure, let's add them in a list
var = [G_total_Slt,G_advection_Slt,G_diffusion_Slt,G_forcing_Slt]
varstrngs = [r'$G^{Slt}_{total}$',r'$G^{Slt}_{advection}$',r'$G^{Slt}_{diffusion}$',r'$G^{Slt}_{forcing}$']
[44]:
# Set an index for the time (t) and depth (k) axis
t, k = 100, 0
Example maps at a particular time and depth level¶
[48]:
axes_class = (GeoAxes,dict(map_projection=cartopy.crs.Robinson(central_longitude=-159)))
fig = plt.figure(figsize=(14,8))
fig.suptitle('Spatial distribution at z = %i m of salt budget components in '\
%np.round(-ecco_grid.Z[k].values)+str(ds.time[t].dt.strftime("%b %Y").values),
fontsize=16, fontweight='bold')
axgr = AxesGrid(fig, 111, axes_class=axes_class, nrows_ncols=(2, 2), axes_pad=(0.1 ,0.5),
share_all=True, label_mode='')
for i, ax in enumerate(axgr):
new_grid_lon, new_grid_lat, field_nearest_1deg =\
ecco.resample_to_latlon(ecco_grid.XC, ecco_grid.YC,
var[i].isel(time=t,k=k),
new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,
new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,
fill_value = np.NaN, mapping_method = 'nearest_neighbor',
radius_of_influence = 120000)
ax.coastlines(linewidth=1.0)
ax.add_feature(cartopy.feature.LAND,color='lightgrey')
ax.set_title(varstrngs[i],fontsize=16)
p = ax.contourf(new_grid_lon, new_grid_lat, field_nearest_1deg*1e7, transform=cartopy.crs.PlateCarree(),
vmin=-5, vmax=5, cmap='RdBu_r', levels=np.linspace(-5, 5, 51), extend='both')
cax = fig.add_axes([0.92, 0.2, 0.015, 0.6])
cb = fig.colorbar(p, cax=cax, orientation='vertical',ticks=np.linspace(-5, 5, 11))
cb.ax.tick_params(labelsize=12)
cb.set_label(r'10$^{-7}$ psu s$^{-1}$', fontsize=14, fontweight='bold')
plt.show()

Time-mean distribution¶
[49]:
axes_class = (GeoAxes,dict(map_projection=cartopy.crs.Robinson(central_longitude=-159)))
fig = plt.figure(figsize=(14,8))
fig.suptitle('Spatial distribution at z = %i m of salt budget components averaged over period %i-%i,'\
%(np.round(-ecco_grid.Z[k].values),year_start,year_end),
fontsize=16, fontweight='bold')
axgr = AxesGrid(fig, 111, axes_class=axes_class, nrows_ncols=(2, 2), axes_pad=(0.1 ,0.5),
share_all=True, label_mode='')
for i, ax in enumerate(axgr):
new_grid_lon, new_grid_lat, field_nearest_1deg =\
ecco.resample_to_latlon(ecco_grid.XC, ecco_grid.YC,
var[i].mean('time').isel(k=k),
new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,
new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,
fill_value = np.NaN, mapping_method = 'nearest_neighbor',
radius_of_influence = 120000)
ax.coastlines(linewidth=1.0)
ax.add_feature(cartopy.feature.LAND,color='lightgrey')
ax.set_title(varstrngs[i],fontsize=16)
p = ax.contourf(new_grid_lon, new_grid_lat, field_nearest_1deg*1e7, transform=cartopy.crs.PlateCarree(),
vmin=-5, vmax=5, cmap='RdBu_r', levels=np.linspace(-5, 5, 51), extend='both')
cax = fig.add_axes([0.92, 0.2, 0.015, 0.6])
cb = fig.colorbar(p, cax=cax, orientation='vertical',ticks=np.linspace(-5, 5, 11))
cb.ax.tick_params(labelsize=12)
cb.set_label(r'10$^{-7}$ psu s$^{-1}$', fontsize=14, fontweight='bold')
plt.show()

From the maps above, we can see that the balance in the salt budget is mostly between the advective and diffusive convergence, and the forcing term is only relevant close to the sea ice edge.
Salt budget closure through time¶
Global average budget closure¶
[50]:
# Take volume-weighted mean of these terms
tmp_a1 = (G_total_Slt*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_a2 = (rhs*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_b = (G_advection_Slt*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_c = (G_diffusion_Slt*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_d = (G_forcing_Slt*vol).sum(dim=('k','i','j','tile'))/vol.sum()
[51]:
fig, axs = plt.subplots(2, 2, figsize=(14,8))
plt.sca(axs[0,0])
tmp_a1.plot(color='k',lw=2)
tmp_a2.plot(color='grey')
axs[0,0].set_title(r'a. $G^{Slt}_{total}$ (black) / RHS (grey) [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[0,1])
tmp_b.plot(color='r')
axs[0,1].set_title(r'b. $G^{Slt}_{advection}$ [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[1,0])
tmp_c.plot(color='orange')
axs[1,0].set_title(r'c. $G^{Slt}_{diffusion}$ [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[1,1])
tmp_d.plot(color='b')
axs[1,1].set_title(r'd. $G^{Slt}_{forcing}$ [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.subplots_adjust(hspace = .5, wspace=.2)
plt.suptitle('Global Salt Budget', fontsize=16)
plt.show()

The globally-averaged salt budget is driven by the forcing term, which mostly represents the input/output of salt from sea ice melting/freezing.
Local salt budget closure¶
[52]:
# Pick any set of indices (tile, k, j, i) corresponding to an ocean grid point
t,k,j,i = (12,0,87,16)
print(t,k,j,i)
12 0 87 16
[53]:
fig, axs = plt.subplots(2, 2, figsize=(14,8))
plt.sca(axs[0,0])
G_total_Slt.isel(tile=t,k=k,j=j,i=i).plot(color='k',lw=2)
rhs.isel(tile=t,k=k,j=j,i=i).plot(color='grey')
axs[0,0].set_title(r'a. $G^{Slt}_{total}$ (black) / RHS (grey) [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[0,1])
G_advection_Slt.isel(tile=t,k=k,j=j,i=i).plot(color='r')
axs[0,1].set_title(r'b. $G^{Slt}_{advection}$ [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[1,0])
G_diffusion_Slt.isel(tile=t,k=k,j=j,i=i).plot(color='orange')
axs[1,0].set_title(r'c. $G^{Slt}_{diffusion}$ [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[1,1])
G_forcing_Slt.isel(tile=t,k=k,j=j,i=i).plot(color='b')
axs[1,1].set_title(r'd. $G^{Slt}_{forcing}$ [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.subplots_adjust(hspace = .5, wspace=.2)
plt.suptitle('Salt Budget at a specific grid point (tile = %i, k = %i, j = %i, i = %i)'%(t,k,j,i), fontsize=16)
plt.show()

The balance looks very different for the local salt budget of a specific grid point. We see much greater magnitudes, mostly in the advective and diffusive part. The forcing component is an order of magnitude smaller than and
and only relevant when sea ice is melting/freezing.
Vertical profiles of the salt budget terms¶
[45]:
fig = plt.subplots(1, 2, sharey=True, figsize=(12,7))
plt.subplot(1, 2, 1)
plt.plot(G_total_Slt.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
lw=4, color='black', marker='.', label=r'$G^{Slt}_{total}$ (LHS)')
plt.plot(G_advection_Slt.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
lw=2, color='red', marker='.', label=r'$G^{Slt}_{advection}$')
plt.plot(G_diffusion_Slt.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
lw=2, color='orange', marker='.', label=r'$G^{Slt}_{diffusion}$')
plt.plot(G_forcing_Slt.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
lw=2, color='blue', marker='.', label=r'$G^{Slt}_{forcing}$')
plt.plot(rhs.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z, lw=1, color='grey', marker='.', label='RHS')
plt.xlabel(r'Tendency [psu s$^{-1}$]', fontsize=14)
plt.ylim([-200,0])
plt.ylabel('Depth (m)', fontsize=14)
plt.gca().tick_params(axis='both', which='major', labelsize=12)
plt.legend(loc='lower left', frameon=False, fontsize=12)
plt.subplot(1, 2, 2)
plt.plot(G_total_Slt.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
lw=4, color='black', marker='.', label=r'$G^{Slt}_{total}$ (LHS)')
plt.plot(rhs.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z, lw=1, color='grey', marker='.', label='RHS')
plt.setp(plt.gca(), 'yticklabels',[])
plt.xlabel(r'Tendency [psu s$^{-1}$]', fontsize=14)
plt.ylim([-200,0])
plt.show()

The above examples illustrate that we can close the salt budget globally/spatially averaged, locally (for each grid point) at a specific time or averaged over time.
Given the nonlinear free surface condition, budgets for salt content (an extensive quantity) are not the same as budgets for salinity (an intensive quantity). The surface freshwater exchanges do not enter into the salt budget, since such fluxes do not affect the overall salt content, but rather make it more or less concentrated. However, a budget for salinity can be derived based on the conservation equations for salt and volume, and estimated using diagnostic model output. Such details are given in the section below.
Evaluating the salinity budget¶
In this section, we demonstrate how to estimate the salinity budget using output from the ECCOv4 solution. Each term in the following salinity budget equation will be evaluated.
Scale factor¶
Closing the salinity budget requires accurate estimates of volume changes for each grid cell. Thus, we need to explicitly calculate the scale factor () to be used in our calculations below.
This requires following model output: - Depth
: Ocean depth, (m) -
ETAN
: Surface Height Anomaly, (m)
[52]:
# Scale factor
rstarfac = ((ecco_grid.Depth + ecco_monthly_mean.ETAN)/ecco_grid.Depth)
Total salinity tendency¶
We calculate the monthly-averaged time tendency of salinity by differencing monthly SALT
snapshots. This operation includes dividing by the number of seconds between each snapshot.
[57]:
# Total tendency (psu/s)
G_total_Sln = grid.diff(ds.SALT_snp, 'T', boundary='fill', fill_value=0.0)/delta_t
Advective salinity convergence¶
Based on the derivation in the Introduction section, the salinity budget requires terms from both the volume and salt budgets. For the advective convergence of salinity, we first need to derive the convergence of volume.
Horizontal convergence¶
Relevant model output: - UVELMASS
: U Mass-Weighted Component of Velocity (m/s) - VVELMASS
: V Mass-Weighted Component of Velocity (m/s)
[58]:
# Horizontal volume transports (m^3/s)
u_transport = ds.UVELMASS * ecco_grid.dyG * ecco_grid.drF
v_transport = ds.VVELMASS * ecco_grid.dxG * ecco_grid.drF
uv_diff = grid.diff_2d_vector({'X' : u_transport, 'Y' : v_transport}, boundary = 'fill')
# Convergence of the horizontal flow (m^3/s)
hConvV = -(uv_diff['X'] + uv_diff['Y'])
Advective convergence of salinity has two parts: the advective salt flux (adv_hConvS
), and the tendency due to volume convergence (hConvV
).
[63]:
# Horizontal convergence of salinity (m^3/s)
adv_hConvSln = (-ds.SALT*hConvV + adv_hConvS)/rstarfac
Vertical convergence¶
Relevant model output: - WVELMASS
: Vertical Mass-Weighted Component of Velocity (m/s) > Note: WVELMASS[k=0] == -oceFWflx/rho0
. If we don’t zero out the top cell, we end up double counting the surface flux.
[61]:
# Vertical volume transport (m^3/s)
w_transport = ds.WVELMASS.where(ds.k_l>0).fillna(0.) * ecco_grid.rA
[62]:
# Convergence of the vertical flow (m^3/s)
vConvV = grid.diff(w_transport, 'Z', boundary='fill')
Again, to get the vertical convergence of salinity we need both the vertical salt flux (adv_vConvS
) and convergece of vertical flow (vConvV
).
[64]:
# Vertical convergence of salinity (psu m^3/s)
adv_vConvSln = (-ds.SALT*vConvV + adv_vConvS)/rstarfac
Total advective salinity convergence¶
[65]:
# Total convergence of advective salinity flux (psu/s)
G_advection_Sln = (adv_hConvSln + adv_vConvSln)/vol
Diffusive salinity convergence¶
The diffusive flux of salinity is pretty much the same as for salt. The only step is dividing the convergence of salt diffusion by the scale factor.
[66]:
# Horizontal convergence
dif_hConvSln = dif_hConvS/rstarfac
# Vertical convergence
dif_vConvSln = dif_vConvS/rstarfac
# Sum horizontal and vertical convergences and divide by volume (psu/s)
G_diffusion_Sln = (dif_hConvSln + dif_vConvSln)/vol
Salinity forcing¶
The forcing term is comprised of both salt flux (forcS
) and volume (i.e., surface freshwater) fluxes (forcV
). We now require monthly mean salinity SALT
to convert forcV
to appropriate units.
Volume forcing¶
oceFWflx
: net surface Fresh-Water flux into the ocean (kg/m^2/s)
[67]:
# Load monthly averaged freshwater flux and add vertical coordinate
oceFWflx = ds.oceFWflx.assign_coords(k=0).expand_dims('k')
# Sea surface forcing on volume (1/s)
forcV = xr.concat([(oceFWflx/rhoconst)/(ecco_grid.hFacC*ecco_grid.drF),
xr.zeros_like(((oceFWflx[0]/rhoconst)/\
(ecco_grid.hFacC*ecco_grid.drF)).transpose('time','tile','k','j','i'))[:,:,1:]],
dim='k').transpose('time','tile','k','j','i')
[68]:
# Sea surface forcing for salinity (psu/s)
G_forcing_Sln = (-ds.SALT*forcV + G_forcing_Slt)/rstarfac
Salinity budget: Map of residual¶
[69]:
# Total convergence (psu/s)
ConvSln = G_advection_Sln + G_diffusion_Sln
# Sum of terms in RHS of equation (psu/s)
rhs_Sln = ConvSln + G_forcing_Sln
[70]:
# Accumulated residual
resSln = (rhs_Sln-G_total_Sln).sum(dim='k').sum(dim='time').compute()
[78]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, resSln,
cmin=-2e-8, cmax=2e-8, show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Accumulated residual (RHS - LHS) [psu s$^{-1}$]', fontsize=16)
plt.show()

The residual in the salinity budget are more extensive compared to the salt budget. Here errors occur that are mostly found in the continental shelves and high latitudes. However, given that the above map shows the accumulated residual, the errors are very small compared to the salinity tendencies’ overall range of values.
Salinity budget: Spatial distributions¶
[89]:
# In order to plot the budget terms in one figure, let's add them in a list
var = [G_total_Sln,G_advection_Sln,G_diffusion_Sln,G_forcing_Sln]
varstrngs = [r'$G^{Sln}_{total}$',r'$G^{Sln}_{advection}$',r'$G^{Sln}_{diffusion}$',r'$G^{Sln}_{forcing}$']
[90]:
# Set an index for the time (t) and depth (k) axis
t, k = 100, 0
Example maps at a particular time and depth level¶
[93]:
axes_class = (GeoAxes,dict(map_projection=cartopy.crs.Robinson(central_longitude=-159)))
fig = plt.figure(figsize=(14,8))
fig.suptitle('Spatial distribution at z = %i m of salinity budget components in '\
%np.round(-ecco_grid.Z[k].values)+str(ds.time[t].dt.strftime("%b %Y").values),
fontsize=16, fontweight='bold')
axgr = AxesGrid(fig, 111, axes_class=axes_class, nrows_ncols=(2, 2), axes_pad=(0.1 ,0.5),
share_all=True, label_mode='')
for i, ax in enumerate(axgr):
new_grid_lon, new_grid_lat, field_nearest_1deg =\
ecco.resample_to_latlon(ecco_grid.XC, ecco_grid.YC,
var[i].isel(time=t,k=k),
new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,
new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,
fill_value = np.NaN, mapping_method = 'nearest_neighbor',
radius_of_influence = 120000)
ax.coastlines(linewidth=1.0,zorder=2)
ax.add_feature(cartopy.feature.LAND,color='lightgrey',zorder=1)
ax.set_title(varstrngs[i],fontsize=16)
p = ax.contourf(new_grid_lon, new_grid_lat, field_nearest_1deg*1e6, transform=cartopy.crs.PlateCarree(),
vmin=-1, vmax=1, cmap='RdBu_r', levels=np.linspace(-1, 1, 51), extend='both',zorder=0)
cax = fig.add_axes([0.92, 0.2, 0.015, 0.6])
cb = fig.colorbar(p, cax=cax, orientation='vertical',ticks=np.linspace(-1, 1, 11))
cb.ax.tick_params(labelsize=12)
cb.set_label(r'10$^{-6}$ psu s$^{-1}$', fontsize=14, fontweight='bold')
plt.show()

Time-mean distribution¶
[99]:
axes_class = (GeoAxes,dict(map_projection=cartopy.crs.Robinson(central_longitude=-159)))
fig = plt.figure(figsize=(14,8))
fig.suptitle('Spatial distribution at z = %i m of salinity budget components averaged over period %i-%i,'\
%(np.round(-ecco_grid.Z[k].values),year_start,year_end),
fontsize=16, fontweight='bold')
axgr = AxesGrid(fig, 111, axes_class=axes_class, nrows_ncols=(2, 2), axes_pad=(0.1 ,0.5),
share_all=True, label_mode='')
for i, ax in enumerate(axgr):
new_grid_lon, new_grid_lat, field_nearest_1deg =\
ecco.resample_to_latlon(ecco_grid.XC, ecco_grid.YC,
var[i].mean('time').isel(k=k),
new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,
new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,
fill_value = np.NaN, mapping_method = 'nearest_neighbor',
radius_of_influence = 120000)
ax.coastlines(linewidth=1.0,zorder=2)
ax.add_feature(cartopy.feature.LAND,color='lightgrey',zorder=1)
ax.set_title(varstrngs[i],fontsize=16)
p = ax.contourf(new_grid_lon, new_grid_lat, field_nearest_1deg*1e7, transform=cartopy.crs.PlateCarree(),
vmin=-5, vmax=5, cmap='RdBu_r', levels=np.linspace(-5, 5, 51), extend='both',zorder=0)
cax = fig.add_axes([0.92, 0.2, 0.015, 0.6])
cb = fig.colorbar(p, cax=cax, orientation='vertical',ticks=np.linspace(-5, 5, 11))
cb.ax.tick_params(labelsize=12)
cb.set_label(r'10$^{-6}$ psu s$^{-1}$', fontsize=14, fontweight='bold')
plt.show()

Unlike with the salt budget, we now see a clear spatial pattern in the forcing term, which resembles surface freshwater flux.
Salinity budget closure through time¶
This section illustrates that we can close the salinity budget globally and locally (i.e., at any given grid point).
Global average budget closure¶
[94]:
# Take volume-weighted mean of these terms
tmp_a1 = (G_total_Sln*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_a2 = (rhs_Sln*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_b = (G_advection_Sln*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_c = (G_diffusion_Sln*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_d = (G_forcing_Sln*vol).sum(dim=('k','i','j','tile'))/vol.sum()
[95]:
fig, axs = plt.subplots(2, 2, figsize=(14,8))
plt.sca(axs[0,0])
tmp_a1.plot(color='k',lw=2)
tmp_a2.plot(color='grey')
axs[0,0].set_title(r'a. $G^{Sln}_{total}$ (black) / RHS (grey) [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[0,1])
tmp_b.plot(color='r')
axs[0,1].set_title(r'b. $G^{Sln}_{advection}$ [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[1,0])
tmp_c.plot(color='orange')
axs[1,0].set_title(r'c. $G^{Sln}_{diffusion}$ [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[1,1])
tmp_d.plot(color='b')
axs[1,1].set_title(r'd. $G^{Sln}_{forcing}$ [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.subplots_adjust(hspace = .5, wspace=.2)
plt.suptitle('Global Salinity Budget', fontsize=16)
plt.show()

Local salt budget closure¶
[96]:
# Pick any set of indices (tile, k, j, i) corresponding to an ocean grid point
t,k,j,i = (12,0,87,16)
print(t,k,j,i)
12 0 87 16
[97]:
fig, axs = plt.subplots(2, 2, figsize=(14,8))
plt.sca(axs[0,0])
G_total_Sln.isel(tile=t,k=k,j=j,i=i).plot(color='k',lw=2)
rhs_Sln.isel(tile=t,k=k,j=j,i=i).plot(color='grey')
axs[0,0].set_title(r'a. $G^{Sln}_{total}$ (black) / RHS (grey) [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[0,1])
G_advection_Sln.isel(tile=t,k=k,j=j,i=i).plot(color='r')
axs[0,1].set_title(r'b. $G^{Sln}_{advection}$ [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[1,0])
G_diffusion_Sln.isel(tile=t,k=k,j=j,i=i).plot(color='orange')
axs[1,0].set_title(r'c. $G^{Sln}_{diffusion}$ [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[1,1])
G_forcing_Sln.isel(tile=t,k=k,j=j,i=i).plot(color='b')
axs[1,1].set_title(r'd. $G^{Sln}_{forcing}$ [psu s$^{-1}$]', fontsize=12)
plt.grid()
plt.subplots_adjust(hspace = .5, wspace=.2)
plt.suptitle('Salinity budget at a specific grid point (tile = %i, k = %i, j = %i, i = %i)'%(t,k,j,i), fontsize=16)
plt.show()

Vertical profiles of the salinity budget terms¶
This section illustrates the balance in the salinity budget along the depth axis.
[98]:
fig = plt.subplots(1, 2, sharey=True, figsize=(12,7))
plt.subplot(1, 2, 1)
plt.plot(G_total_Sln.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
lw=4, color='black', marker='.', label=r'$G^{Slt}_{total}$ (LHS)')
plt.plot(G_advection_Sln.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
lw=2, color='red', marker='.', label=r'$G^{Slt}_{advection}$')
plt.plot(G_diffusion_Sln.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
lw=2, color='orange', marker='.', label=r'$G^{Slt}_{diffusion}$')
plt.plot(G_forcing_Sln.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
lw=2, color='blue', marker='.', label=r'$G^{Slt}_{forcing}$')
plt.plot(rhs_Sln.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z, lw=1, color='grey', marker='.', label='RHS')
plt.xlabel(r'Tendency [psu s$^{-1}$]', fontsize=14)
plt.ylim([-200,0])
plt.ylabel('Depth (m)', fontsize=14)
plt.gca().tick_params(axis='both', which='major', labelsize=12)
plt.legend(loc='lower left', frameon=False, fontsize=12)
plt.subplot(1, 2, 2)
plt.plot(G_total_Sln.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z,
lw=4, color='black', marker='.', label=r'$G^{Slt}_{total}$ (LHS)')
plt.plot(rhs_Sln.isel(tile=t,j=j,i=i).mean('time'), ecco_grid.Z, lw=1, color='grey', marker='.', label='RHS')
plt.setp(plt.gca(), 'yticklabels',[])
plt.xlabel(r'Tendency [psu s$^{-1}$]', fontsize=14)
plt.ylim([-200,0])
plt.show()

Evaluating the freshwater budget - Work in progress¶
As with the salt and a salinity budget we will evaluate each term in the freshwater budget.
Suggestion by Matt Mazloff:
The reference is just a multiplying factor
So the budget is
And the user can specify the later or easily convert. I would just use 35 for the tutorial. And you use the exact salinity budget terms. The diffusive flux of freshwater is minus the diffusive flux of salt.
[46]:
# Reference salinity
Sref = 35.0
Total freshwater tendency¶
[47]:
f = (Sref - ecco_monthly_snaps.SALT)/Sref
G_total_Fw_month = f.isel(time=range(1,num_months)).values - f.isel(time=range(0,num_months-1)).values
[48]:
# Convert numpy array to an xarray DataArray with matching dimensions as the monthly mean fields
G_total_Fw_month = xr.DataArray(G_total_Fw_month, coords={'time': ecco_monthly_mean.time.values,
'tile': ecco_monthly_mean.tile.values,
'k': ecco_monthly_mean.k.values,
'j': ecco_monthly_mean.j.values,
'i': ecco_monthly_mean.i.values},
dims=['time','tile','k','j','i'])
[49]:
# Freshwater tendency (m^3/s)
G_total_Fw = G_total_Fw_month*vol / secs_per_month
Advective flux of freshwater¶
Advective fluxes of freshwater are calculated offline using salinity and velocity fields:
GM Bolus Velocity¶
[50]:
UVELSTAR = grid.diff(ecco_monthly_mean.GM_PsiX, 'Z', boundary='fill')/ecco_grid.drF
VVELSTAR = grid.diff(ecco_monthly_mean.GM_PsiY, 'Z', boundary='fill')/ecco_grid.drF
[51]:
GM_PsiXY_diff = grid.diff_2d_vector({'X' : ecco_monthly_mean.GM_PsiX*ecco_grid.dyG,
'Y' : ecco_monthly_mean.GM_PsiY*ecco_grid.dxG}, boundary = 'fill')
WVELSTAR = (GM_PsiXY_diff['X'] + GM_PsiXY_diff['Y'])/ecco_grid.rA
Calculate advective freshwater flux¶
[52]:
SALT_at_u = grid.interp(ecco_monthly_mean.SALT, 'X', boundary='extend')
SALT_at_v = grid.interp(ecco_monthly_mean.SALT, 'Y', boundary='extend')
SALT_at_w = grid.interp(ecco_monthly_mean.SALT, 'Z', boundary='extend')
[53]:
# Freshwater advective (Eulerian+Bolus) fluxes (m^3/s)
ADVx_FW = (ecco_monthly_mean.UVELMASS+UVELSTAR)*ecco_grid.dyG*ecco_grid.drF*(Sref-SALT_at_u)/Sref
ADVy_FW = (ecco_monthly_mean.VVELMASS+VVELSTAR)*ecco_grid.dxG*ecco_grid.drF*(Sref-SALT_at_v)/Sref
ADVr_FW = WVELMASS*ecco_grid.rA*(Sref-SALT_at_w)/Sref
[54]:
ADVxy_diff = grid.diff_2d_vector({'X' : ADVx_FW, 'Y' : ADVy_FW}, boundary = 'fill')
# Convergence of horizontal advection (m^3/s)
adv_hConvFw = (-(ADVxy_diff['X'] + ADVxy_diff['Y']))
[55]:
# Convergence of vertical advection (m^3/s)
adv_vConvFw = grid.diff(ADVr_FW, 'Z', boundary='fill')
[56]:
# Sum horizontal and vertical convergences (m^3/s)
G_advection_Fw = adv_hConvFw + adv_vConvFw
Freshwater forcing¶
[57]:
# Freshwater forcing (m^3/s)
forcFw = ecco_monthly_mean.oceFWflx/rhoconst*ecco_grid.rA
# Expand to fully 3d (using G_advection_Fw as template)
G_forcing_Fw = xr.concat([forcFw.reset_coords(drop=True).assign_coords(k=0).expand_dims('k')\
.transpose('time','tile','k','j','i'),
xr.zeros_like(G_advection_Fw.transpose('time','tile','k','j','i'))[:,:,1:]],
dim='k').transpose('time','tile','k','j','i').where(ecco_grid.hFacC==1)
Diffusive freshwater flux¶
[58]:
# Convergence of freshwater diffusion (m^3/s)
G_diffusion_Fw = G_total_Fw - G_forcing_Fw - G_advection_Fw
Save budget terms¶
Now that we have all the terms evaluated, let’s save them to a dataset. Here are two examples:
Add all variables to a new dataset¶
[59]:
#varnames = ['G_total_Slt','G_advection_Slt','G_diffusion_Slt','G_forcing_Slt',
# 'G_total_Sln','G_advection_Sln','G_diffusion_Sln','G_forcing_Sln',
# 'G_total_Fw', 'G_advection_Fw', 'G_diffusion_Fw', 'G_forcing_Fw']
varnames = ['G_total_Sln','G_advection_Sln','G_diffusion_Sln','G_forcing_Sln']
ds = xr.Dataset(data_vars={})
for varname in varnames:
ds[varname] = globals()[varname].chunk(chunks={'time':1,'tile':13,'k':50,'j':90,'i':90})
[60]:
ds.time.encoding = {}
ds = ds.reset_coords(drop=True)
Save to zarr dataset¶
[61]:
from dask.diagnostics import ProgressBar
[62]:
ds
[62]:
<xarray.Dataset> Dimensions: (i: 90, j: 90, k: 50, tile: 13, time: 264) Coordinates: * time (time) datetime64[ns] 1993-01-16T12:00:00 ... 2014-12-16... * tile (tile) int32 0 1 2 3 4 5 6 7 8 9 10 11 12 * k (k) int32 0 1 2 3 4 5 6 7 8 ... 41 42 43 44 45 46 47 48 49 * j (j) int32 0 1 2 3 4 5 6 7 8 ... 81 82 83 84 85 86 87 88 89 * i (i) int32 0 1 2 3 4 5 6 7 8 ... 81 82 83 84 85 86 87 88 89 Data variables: G_total_Sln (time, tile, k, j, i) float64 dask.array<chunksize=(1, 13, 50, 90, 90), meta=np.ndarray> G_advection_Sln (time, tile, k, j, i) float32 dask.array<chunksize=(1, 13, 50, 90, 90), meta=np.ndarray> G_diffusion_Sln (time, tile, k, j, i) float32 dask.array<chunksize=(1, 13, 50, 90, 90), meta=np.ndarray> G_forcing_Sln (time, tile, k, j, i) float32 dask.array<chunksize=(1, 13, 50, 90, 90), meta=np.ndarray>
- i: 90
- j: 90
- k: 50
- tile: 13
- time: 264
- time(time)datetime64[ns]1993-01-16T12:00:00 ... 2014-12-...
- long_name :
- center time of averaging period
- standard_name :
- time
- bounds :
- time_bnds
- axis :
- T
array(['1993-01-16T12:00:00.000000000', '1993-02-15T12:00:00.000000000', '1993-03-16T12:00:00.000000000', ..., '2014-10-16T12:00:00.000000000', '2014-11-16T12:00:00.000000000', '2014-12-16T12:00:00.000000000'], dtype='datetime64[ns]')
- tile(tile)int320 1 2 3 4 5 6 7 8 9 10 11 12
- standard_name :
- tile_index
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], dtype=int32)
- k(k)int320 1 2 3 4 5 6 ... 44 45 46 47 48 49
- long_name :
- z-dimension of the t grid
- standard_name :
- z_grid_index
- swap_dim :
- Z
- axis :
- Z
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49], dtype=int32)
- j(j)int320 1 2 3 4 5 6 ... 84 85 86 87 88 89
- long_name :
- y-dimension of the t grid
- standard_name :
- y_grid_index
- swap_dim :
- YC
- axis :
- Y
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89], dtype=int32)
- i(i)int320 1 2 3 4 5 6 ... 84 85 86 87 88 89
- long_name :
- x-dimension of the t grid
- standard_name :
- x_grid_index
- swap_dim :
- XC
- axis :
- X
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89], dtype=int32)
- G_total_Sln(time, tile, k, j, i)float64dask.array<chunksize=(1, 13, 50, 90, 90), meta=np.ndarray>
Array Chunk Bytes 11.12 GB 42.12 MB Shape (264, 13, 50, 90, 90) (1, 13, 50, 90, 90) Count 265 Tasks 264 Chunks Type float64 numpy.ndarray - G_advection_Sln(time, tile, k, j, i)float32dask.array<chunksize=(1, 13, 50, 90, 90), meta=np.ndarray>
Array Chunk Bytes 5.56 GB 21.06 MB Shape (264, 13, 50, 90, 90) (1, 13, 50, 90, 90) Count 64384 Tasks 264 Chunks Type float32 numpy.ndarray - G_diffusion_Sln(time, tile, k, j, i)float32dask.array<chunksize=(1, 13, 50, 90, 90), meta=np.ndarray>
Array Chunk Bytes 5.56 GB 21.06 MB Shape (264, 13, 50, 90, 90) (1, 13, 50, 90, 90) Count 35658 Tasks 264 Chunks Type float32 numpy.ndarray - G_forcing_Sln(time, tile, k, j, i)float32dask.array<chunksize=(1, 13, 50, 90, 90), meta=np.ndarray>
Array Chunk Bytes 5.56 GB 21.06 MB Shape (264, 13, 50, 90, 90) (1, 13, 50, 90, 90) Count 3093 Tasks 264 Chunks Type float32 numpy.ndarray
[63]:
with ProgressBar():
ds.to_zarr(base_dir + '/eccov4r3_budg_Slt_Sln_Fw')
[########################################] | 100% Completed | 12min 42.1s
Load budget variables from file¶
After having saved the budget terms to file, let’s restart the kernel and load only the relevant data and Python modules.
[1]:
# Suppress warning messages for a cleaner presentation
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import xarray as xr
import ecco_v4_py as ecco
import matplotlib.pyplot as plt
%matplotlib inline
base_dir = '/mnt/efs/data/ECCOv4-release'
ECCO_dir = base_dir + '/Release3_alt'
grid_dir= ECCO_dir + '/nctiles_grid/'
ecco_grid = ecco.load_ecco_grid_nc(grid_dir, 'ECCOv4r3_grid.nc')
grid = ecco.get_llc_grid(ecco_grid)
# Volume (m^3)
vol = (ecco_grid.rA*ecco_grid.drF*ecco_grid.hFacC).transpose('tile','k','j','i')
data_dir= ECCO_dir + '/nctiles_monthly'
year_start = 1993
year_end = 2015
ecco_monthly_mean = ecco.recursive_load_ecco_var_from_years_nc(data_dir, \
vars_to_load=['UVELMASS', 'VVELMASS', 'WVELMASS', 'GM_PsiX', 'GM_PsiY', 'SALT', 'oceFWflx'],
years_to_load=range(year_start, year_end))
loading files of GM_PsiX
loading files of GM_PsiY
loading files of SALT
loading files of UVELMASS
loading files of VVELMASS
loading files of WVELMASS
loading files of oceFWflx
[2]:
# Load terms from zarr dataset
G_total_Sln = xr.open_zarr(base_dir + '/eccov4r3_budg_Slt_Sln_Fw').G_total_Sln
G_advection_Sln = xr.open_zarr(base_dir + '/eccov4r3_budg_Slt_Sln_Fw').G_advection_Sln
G_diffusion_Sln = xr.open_zarr(base_dir + '/eccov4r3_budg_Slt_Sln_Fw').G_diffusion_Sln
G_forcing_Sln = xr.open_zarr(base_dir + '/eccov4r3_budg_Slt_Sln_Fw').G_forcing_Sln
[4]:
UVELSTAR = grid.diff(ecco_monthly_mean.GM_PsiX, 'Z', boundary='fill')/ecco_grid.drF
VVELSTAR = grid.diff(ecco_monthly_mean.GM_PsiY, 'Z', boundary='fill')/ecco_grid.drF
GM_PsiXY_diff = grid.diff_2d_vector({'X' : ecco_monthly_mean.GM_PsiX*ecco_grid.dyG,
'Y' : ecco_monthly_mean.GM_PsiY*ecco_grid.dxG}, boundary = 'fill')
WVELSTAR = (GM_PsiXY_diff['X'] + GM_PsiXY_diff['Y'])/ecco_grid.rA
SALT_at_u = grid.interp(ecco_monthly_mean.SALT, 'X', boundary='extend')
SALT_at_v = grid.interp(ecco_monthly_mean.SALT, 'Y', boundary='extend')
SALT_at_w = grid.interp(ecco_monthly_mean.SALT, 'Z', boundary='extend')
# Remove oceFWflx from WVELMASS
WVELMASS = ecco_monthly_mean.WVELMASS.transpose('time','tile','k_l','j','i')
oceFWflx = ecco_monthly_mean.oceFWflx.assign_coords(k_l=0).expand_dims('k_l').transpose('time','tile','k_l','j','i')
# Seawater density (kg/m^3)
rhoconst = 1029
oceFWflx = (oceFWflx/rhoconst)
WVELMASS = xr.concat([WVELMASS.sel(k_l=0) + oceFWflx, WVELMASS[:,:,1:]],
dim='k_l').transpose('time','tile','k_l','j','i')
# Salinity advective (Eulerian+Bolus) fluxes (psu m^3/s)
ADVx_SLT = (ecco_monthly_mean.UVELMASS+UVELSTAR)*ecco_grid.dyG*ecco_grid.drF*SALT_at_u
#ADVx_SLT = ecco_monthly_mean.UVELMASS*ecco_grid.dyG*ecco_grid.drF*SALT_at_u
ADVy_SLT = (ecco_monthly_mean.VVELMASS+VVELSTAR)*ecco_grid.dxG*ecco_grid.drF*SALT_at_v
#ADVy_SLT = ecco_monthly_mean.VVELMASS*ecco_grid.dxG*ecco_grid.drF*SALT_at_v
ADVr_SLT = (WVELMASS+WVELSTAR)*ecco_grid.rA*SALT_at_w
#ADVr_SLT = WVELMASS*ecco_grid.rA*SALT_at_w
ADVxy_diff = grid.diff_2d_vector({'X' : ADVx_SLT, 'Y' : ADVy_SLT}, boundary = 'fill')
adv_hConvS = (-(ADVxy_diff['X'] + ADVxy_diff['Y']))
adv_vConvS = grid.diff(ADVr_SLT, 'Z', boundary='fill')
G_advection = (adv_hConvS + adv_vConvS)/vol
[98]:

Comparison between LHS and RHS of the budget equation¶
[5]:
# Total convergence
ConvSln = G_advection_Sln + G_diffusion_Sln
ConvSln2 = G_advection + G_diffusion_Sln
[6]:
# Sum of terms in RHS of equation
rhsSln = ConvSln + G_forcing_Sln
rhsSln2 = ConvSln2 + G_forcing_Sln
Map of residuals¶
[7]:
resSln = (G_advection_Sln + G_diffusion_Sln + G_forcing_Sln - G_total_Sln).sum(dim='k').sum(dim='time').compute()
resSln2 = (rhsSln2 - G_total_Sln).sum(dim='k').sum(dim='time').compute()
[8]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, resSln,
cmin=-1e-7, cmax=1e-7, show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Residual [psu s$^{-1}$] (RHS - LHS)', fontsize=16)
plt.show()

[20]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, resSln2,
cmin=-1e-4, cmax=1e-4, show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Residual [psu s$^{-1}$] (RHS - LHS)', fontsize=16)
plt.show()

[21]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_advection,\
cmin=-1e-4, cmax=1e-4, show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Advection [psu s$^{-1}$]', fontsize=16)
plt.show()

[22]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_total,
cmin=-1e-4, cmax=1e-4, show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Total tendency [psu s$^{-1}$]', fontsize=16)
plt.show()

[23]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_diffusion,
cmin=-1e-4, cmax=1e-4,show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Diffusion [psu s$^{-1}$]', fontsize=16)
plt.show()

[24]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_forcing,
cmin=-1e-4, cmax=1e-4,show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Forcing [psu s$^{-1}$]', fontsize=16)
plt.show()

[88]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_forcing_Sln[-1,:,0],
show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Forcing [psu s$^{-1}$]', fontsize=16)
plt.show()

[87]:
[87]:
<xarray.DataArray 'G_forcing_Sln' (tile: 13, j: 90, i: 90)> dask.array<getitem, shape=(13, 90, 90), dtype=float32, chunksize=(13, 90, 90), chunktype=numpy.ndarray> Coordinates: * i (i) int32 0 1 2 3 4 5 6 7 8 9 10 ... 80 81 82 83 84 85 86 87 88 89 * j (j) int32 0 1 2 3 4 5 6 7 8 9 10 ... 80 81 82 83 84 85 86 87 88 89 k int32 0 * tile (tile) int32 0 1 2 3 4 5 6 7 8 9 10 11 12 time datetime64[ns] 2014-12-16T12:00:00
- tile: 13
- j: 90
- i: 90
- dask.array<chunksize=(13, 90, 90), meta=np.ndarray>
Array Chunk Bytes 421.20 kB 421.20 kB Shape (13, 90, 90) (13, 90, 90) Count 266 Tasks 1 Chunks Type float32 numpy.ndarray - i(i)int320 1 2 3 4 5 6 ... 84 85 86 87 88 89
- axis :
- X
- long_name :
- x-dimension of the t grid
- standard_name :
- x_grid_index
- swap_dim :
- XC
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89], dtype=int32)
- j(j)int320 1 2 3 4 5 6 ... 84 85 86 87 88 89
- axis :
- Y
- long_name :
- y-dimension of the t grid
- standard_name :
- y_grid_index
- swap_dim :
- YC
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89], dtype=int32)
- k()int320
- axis :
- Z
- long_name :
- z-dimension of the t grid
- standard_name :
- z_grid_index
- swap_dim :
- Z
array(0, dtype=int32)
- tile(tile)int320 1 2 3 4 5 6 7 8 9 10 11 12
- standard_name :
- tile_index
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], dtype=int32)
- time()datetime64[ns]2014-12-16T12:00:00
- axis :
- T
- bounds :
- time_bnds
- long_name :
- center time of averaging period
- standard_name :
- time
array('2014-12-16T12:00:00.000000000', dtype='datetime64[ns]')
[35]:
tmp = np.abs(G_advection_Sln + G_diffusion_Sln + G_forcing_Sln - G_total_Sln).values.ravel()
[52]:
tmp = res_offline.values.ravel()
[53]:
plt.figure(figsize=(10,3));
plt.hist(tmp[np.nonzero(tmp > 0)],np.linspace(0, 1.2e-9,1000));
plt.grid()

[36]:
plt.figure(figsize=(10,3));
plt.hist(tmp[np.nonzero(tmp > 0)],np.linspace(0, 1.2e-11,1000));
plt.grid()

[12]:

Histogram of residuals¶
We can look at the distribution of residuals to get a little more confidence.
[25]:
from xhistogram.xarray import histogram
[48]:
res_closed = np.abs(G_advection_Sln + G_diffusion_Sln + G_forcing_Sln - G_total_Sln)
res_closed.name = 'Residual_closed'
[50]:
res_offline = np.abs(G_advection + G_diffusion_Sln + G_forcing_Sln - G_total_Sln)
res_offline.name = 'Residual_offline'
[71]:
res_adv = np.abs(G_diffusion_Sln + G_forcing_Sln - G_total_Sln)
res_adv.name = 'Residual_advection'
[72]:
res_dif = np.abs(G_advection_Sln + G_forcing_Sln - G_total_Sln)
res_dif.name = 'Residual_diffusion'
[73]:
res_frc = np.abs(G_advection_Sln + G_diffusion_Sln - G_total_Sln)
res_frc.name = 'Residual_forcing'
[74]:
histogram(res_closed, bins = [np.linspace(0, 1.2e-11,1000)]).plot()
histogram(res_dif, bins = [np.linspace(0, 1.2e-11,1000)]).plot()
histogram(res_frc, bins = [np.linspace(0, 1.2e-11,1000)]).plot()
histogram(res_adv, bins = [np.linspace(0, 1.2e-11,1000)]).plot()
[74]:
[<matplotlib.lines.Line2D at 0x7f502dc650f0>]

[85]:
histogram(G_forcing_Sln, bins = [np.linspace(-2e-10, 2e-10,1000)]).plot()
[85]:
[<matplotlib.lines.Line2D at 0x7f502e71fd30>]

[78]:
G_forcing_Sln.max().values
[78]:
array(3.7624181e-06, dtype=float32)
[49]:
[49]:
[<matplotlib.lines.Line2D at 0x7f507b70ca20>]

[67]:
histogram(res_adv, bins = [np.linspace(0, 5e-11,1000)]).plot()
[67]:
[<matplotlib.lines.Line2D at 0x7f50e440d978>]

[59]:
histogram(res_offline, bins = [np.linspace(0, 1.2e-4,1000)]).plot()
[59]:
[<matplotlib.lines.Line2D at 0x7f507ba18d68>]

[ ]: