Salt, Salinity and Freshwater Budgets¶
Contributors: Jan-Erik Tesdal, Ryan Abernathey, Ian Fenty, Emma Boland, and Andrew Delman.
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://dspace.mit.edu/handle/1721.1/111094?show=full). 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
Salinity 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:
Datasets to download¶
If you don’t have any of the following datasets already, you will need to download them to complete the tutorial. Aside from the grid geometry file (which has no time dimension), you will need 2 monthly datasets for the full time span of ECCOv4r4 output (1992 through 2017). The ShortNames of the datasets are:
ECCO_L4_GEOMETRY_LLC0090GRID_V4R4
ECCO_L4_OCEAN_3D_SALINITY_FLUX_LLC0090GRID_MONTHLY_V4R4 (1993-2016)
ECCO_L4_OCEAN_3D_VOLUME_FLUX_LLC0090GRID_MONTHLY_V4R4 (1993-2016)
ECCO_L4_OCEAN_BOLUS_STREAMFUNCTION_LLC0090GRID_MONTHLY_V4R4 (1993-2016)
ECCO_L4_FRESH_FLUX_LLC0090GRID_MONTHLY_V4R4 (1993-2016)
ECCO_L4_SSH_LLC0090GRID_SNAPSHOT_V4R4 (1993/1/1-2017/1/1, 1st of each month)
ECCO_L4_TEMP_SALINITY_LLC0090GRID_SNAPSHOT_V4R4 (1993/1/1-2017/1/1, 1st of each month)
If you haven’t yet been through the download tutorial or used the ecco_download module, it may help you to review that information before downloading the datasets.
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.
from os.path import expanduser,join
import sys
user_home_dir = expanduser('~')
sys.path.append(join(user_home_dir,'ECCOv4-py')) # change to directory hosting ecco_v4_py as needed
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 a high-level directory for ECCO fields
ECCO_dir = join(user_home_dir,'Downloads','ECCO_V4r4_PODAAC')
Note: Change
ECCO_dir
to the directory path where ECCOv4 output is stored.
[7]:
# Load the model grid
grid_shortname = "ECCO_L4_GEOMETRY_LLC0090GRID_V4R4"
grid_filename = "GRID_GEOMETRY_ECCO_V4r4_native_llc0090.nc"
ecco_grid = xr.open_dataset(join(ECCO_dir,grid_shortname,grid_filename))
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)
Load monthly snapshots¶
If you don’t already have the relevant files needed locally you will need to download them. Here we will only need snapshots at monthly intervals, but snapshots are available from PO.DAAC (via NASA Earthdata Cloud) at daily intervals. Using daily snapshots spanning 24 years would require downloading >8000 files from each dataset! To limit the number of files we download, we’ll create a text file with only the file names that we need and then use wget
to download them.
[9]:
# define first and last dates of snapshots
snap_first_date = '1993-01-01'
snap_last_date = '2018-01-01'
# snapshot dataset shortnames
etan_snaps_shortname = "ECCO_L4_SSH_LLC0090GRID_SNAPSHOT_V4R4"
salt_snaps_shortname = "ECCO_L4_TEMP_SALINITY_LLC0090GRID_SNAPSHOT_V4R4"
# form of snapshot filenames
etan_snaps_fileform = "SEA_SURFACE_HEIGHT_snap_"
salt_snaps_fileform = "OCEAN_TEMPERATURE_SALINITY_snap_"
# define function for generating list of monthly snapshots to download
def snaps_monthly_textlist(out_filename,ShortName,fileform,snap_first_date,snap_last_date):
"Creates text list of snapshots to download at monthly intervals"
podaac_cloud_path = "https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/"
filename_list = []
snap_firstdate = np.datetime64(snap_first_date,'D')
snap_lastdate = np.datetime64(snap_last_date,'D')
snap_currdate = snap_firstdate
while snap_currdate - snap_lastdate < np.timedelta64(0,'D'):
curr_filename = podaac_cloud_path + ShortName + '/' + fileform \
+ str(snap_currdate) + 'T000000_ECCO_V4r4_native_llc0090.nc\n'
filename_list.append(curr_filename)
snap_nextdate = snap_currdate + 1
while str(snap_nextdate)[8:10] != '01':
snap_nextdate += 1
snap_currdate = snap_nextdate
curr_filename = podaac_cloud_path + ShortName + '/' + fileform \
+ str(snap_lastdate) + 'T000000_ECCO_V4r4_native_llc0090.nc\n'
filename_list.append(curr_filename)
with open(out_filename,'w') as f:
f.writelines(filename_list)
# generate list of monthly snapshots
etan_snaps_list = 'etan_snaps_filenames.txt'
snaps_monthly_textlist(etan_snaps_list,etan_snaps_shortname,etan_snaps_fileform,\
snap_first_date,snap_last_date)
salt_snaps_list = 'salt_snaps_filenames.txt'
snaps_monthly_textlist(salt_snaps_list,salt_snaps_shortname,salt_snaps_fileform,\
snap_first_date,snap_last_date)
Now we can use wget to download the list of files in etan_snaps_filenames.txt and salt_snaps_filenames.txt, e.g., wget –no-verbose –no-clobber –continue -i etan_snaps_filenames.txt -P ~/Downloads/ECCO_V4r4_PODAAC/ECCO_L4_SSH_LLC0090GRID_SNAPSHOT_V4R4/
[10]:
# define year bounds [year_start,year_end)...year_end is excluded
year_start = 1993
year_end = 2017
# directories where snapshots are located
etan_snaps_dir = join(ECCO_dir,etan_snaps_shortname)
salt_snaps_dir = join(ECCO_dir,salt_snaps_shortname)
# define function to get list of files in year range
def files_in_year_range(file_dir,year_start,year_end):
"Creates text list of files in the year range [year_start,year_end)"
import os
import re
files_in_range = []
for file in os.listdir(file_dir):
# use regex search to find year associated with file
year_match = re.search('[0-9]{4}(?=-[0-9]{2})',file)
curr_yr = int(year_match.group(0))
if (curr_yr >= year_start) and (curr_yr < year_end):
files_in_range.append(join(file_dir,file))
return files_in_range
# get list of files in year range (including year_end)
etan_snaps_in_range = files_in_year_range(etan_snaps_dir,year_start,year_end+1)
salt_snaps_in_range = files_in_year_range(salt_snaps_dir,year_start,year_end+1)
# open files as two xarray datasets, then merge to create one dataset with the variables we need
ds_etan_snaps = xr.open_mfdataset(etan_snaps_in_range,\
data_vars='minimal',coords='minimal',compat='override')
ds_salt_snaps = xr.open_mfdataset(salt_snaps_in_range,\
data_vars='minimal',coords='minimal',compat='override')
ecco_monthly_snaps = xr.merge([ds_etan_snaps['ETAN'],ds_salt_snaps['SALT']])
# Exclude snapshots after Jan 1 of year_end
years_spanned = year_end - year_start
ecco_monthly_snaps = ecco_monthly_snaps.isel(time=np.arange(0, (12*years_spanned) + 1))
[11]:
# Drop superfluous coordinates (We already have them in ecco_grid)
ecco_monthly_snaps = ecco_monthly_snaps.reset_coords(drop=True)
Load monthly mean data¶
To download these files you can use the ecco_download module or wget. The ShortNames of the needed datasets are below.
[12]:
# ShortNames of monthly mean datasets needed
FW_surf_flux_shortname = "ECCO_L4_FRESH_FLUX_LLC0090GRID_MONTHLY_V4R4"
salt_flux_shortname = "ECCO_L4_OCEAN_3D_SALINITY_FLUX_LLC0090GRID_MONTHLY_V4R4"
vol_flux_shortname = "ECCO_L4_OCEAN_3D_VOLUME_FLUX_LLC0090GRID_MONTHLY_V4R4"
bolus_strmfcn_shortname = "ECCO_L4_OCEAN_BOLUS_STREAMFUNCTION_LLC0090GRID_MONTHLY_V4R4"
etan_monthly_shortname = "ECCO_L4_SSH_LLC0090GRID_MONTHLY_V4R4"
salt_monthly_shortname = "ECCO_L4_TEMP_SALINITY_LLC0090GRID_MONTHLY_V4R4"
# open datasets where monthly mean files are located
datasets_to_open = ['FW_surf_flux','salt_flux','vol_flux','bolus_strmfcn',\
'etan_monthly','salt_monthly']
for dataset in datasets_to_open:
curr_shortname = eval(dataset + '_shortname')
curr_dir = join(ECCO_dir,curr_shortname)
curr_files_in_range = files_in_year_range(curr_dir,year_start,year_end)
exec('ds_' + dataset + ' = xr.open_mfdataset(curr_files_in_range,parallel=True,data_vars=\'minimal\','\
+ 'coords=\'minimal\', compat=\'override\')')
[13]:
# create one dataset with all the monthly mean variables we need
ecco_monthly_mean = xr.merge([ds_FW_surf_flux[['SFLUX','oceFWflx']],\
ds_salt_flux[['oceSPtnd','ADVx_SLT','ADVy_SLT','ADVr_SLT',\
'DFxE_SLT','DFyE_SLT','DFrE_SLT','DFrI_SLT']],\
ds_vol_flux[['UVELMASS','VVELMASS','WVELMASS']],\
ds_bolus_strmfcn[['GM_PsiX','GM_PsiY']],\
ds_etan_monthly['ETAN'],\
ds_salt_monthly['SALT']])
[14]:
# 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, divided into chunks (dask arrays) along the time axes.
[15]:
ds = xr.merge([ecco_monthly_mean,
ecco_monthly_snaps.rename({'time':'time_snp','ETAN':'ETAN_snp', 'SALT':'SALT_snp'})])\
.chunk({'time':1,'time_snp':1})
Predefine coordinates for global regridding of the ECCO output (used in resample_to_latlon
)¶
[16]:
new_grid_delta_lat = 1
new_grid_delta_lon = 1
new_grid_min_lat = -90
new_grid_max_lat = 90
new_grid_min_lon = -180
new_grid_max_lon = 180
Create the xgcm ‘grid’ object¶
[17]:
# Change time axis of the snapshot variables
ds.time_snp.attrs['c_grid_axis_shift'] = 0.5
[18]:
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).
[19]:
delta_t = grid.diff(ds.time_snp, 'T', boundary='fill', fill_value=np.nan)
[20]:
# 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 .
[21]:
# Calculate the s*S term
sSALT = ds.SALT_snp*(1+ds.ETAN_snp/ecco_grid.Depth)
[22]:
# Total tendency (psu/s)
sSALT_diff = sSALT.diff('time_snp')
sSALT_diff = sSALT_diff.rename({'time_snp':'time'})
del sSALT_diff.time.attrs['c_grid_axis_shift'] # remove attribute from DataArray
sSALT_diff = sSALT_diff.assign_coords(time=ds.time) # correct time coordinate values
G_total_Slt = sSALT_diff/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.
Note: when using at least one recent version of
xgcm
(v0.8.1), errors were triggered when callingdiff_2d_vector
below. If you have this issue, as a temporary fix you can locate thexgcm
directoryimport xgcm print(xgcm.__path__)and then substitute
padding.py
andgrid_ufunc.py
with the files linked below:Alternatively, you can revert
xgcm
to an earlier version (e.g., v0.6.0) that should not have this issue.
[23]:
# Set fluxes on land to zero (instead of NaN)
ds['ADVx_SLT'] = ds.ADVx_SLT.where(ecco_grid.hFacW.values > 0,0)
ds['ADVy_SLT'] = ds.ADVy_SLT.where(ecco_grid.hFacS.values > 0,0)
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)
[24]:
# Set fluxes on land to zero (instead of NaN)
ds['ADVr_SLT'] = ds.ADVr_SLT.where(ecco_grid.hFacC.values > 0,0)
# Convergence of vertical advection (psu/s)
adv_vConvS = grid.diff(ds.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
) inWVELMASS
. Thus, to keep the surface forcing term explicitly represented, one needs to zero out the values ofWVELMASS
at the surface so as to avoid double counting (seeECCO_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.
[25]:
# 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.
[26]:
# Set fluxes on land to zero (instead of NaN)
ds['DFxE_SLT'] = ds.DFxE_SLT.where(ecco_grid.hFacW.values > 0,0)
ds['DFyE_SLT'] = ds.DFyE_SLT.where(ecco_grid.hFacS.values > 0,0)
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.
[27]:
# Set fluxes on land to zero (instead of NaN)
ds['DFrE_SLT'] = ds.DFrE_SLT.where(ecco_grid.hFacC.values > 0,0)
ds['DFrI_SLT'] = ds.DFrI_SLT.where(ecco_grid.hFacC.values > 0,0)
# Convergence of vertical diffusion (psu m^3/s)
dif_vConvS = grid.diff(ds.DFrE_SLT, 'Z', boundary='fill') + grid.diff(ds.DFrI_SLT, 'Z', boundary='fill')
Total diffusive salt convergence¶
[28]:
# 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.
[29]:
# Load SFLUX and add vertical coordinate
SFLUX = ds.SFLUX.assign_coords(k=0).expand_dims(dim='k',axis=1)
# Calculate forcing term by adding SFLUX and oceSPtnd (g/m^2/s)
forcS = xr.concat([SFLUX+ds.oceSPtnd,ds.oceSPtnd.isel(k=slice(1,None))], 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.
[30]:
# Forcing (psu/s)
G_forcing_Slt = forcS/rhoconst/(ecco_grid.hFacC*ecco_grid.drF)
Salt budget: Map of residual¶
[31]:
# 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
[32]:
# individual month weights
time_month_weights = delta_t/delta_t.sum()
# Accumulated residual
resSlt = ((((rhs-G_total_Slt)*vol).sum(dim='k')/vol.sum(dim='k'))*time_month_weights)\
.sum(dim='time').compute()
[33]:
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()
-179.9 179.9
-180.0 180.0
-89.9 89.9
-90.0 90.0
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¶
[34]:
# 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}$']
[35]:
# Set an index for the time (t) and depth (k) axis
t, k = 100, 0
Example maps at a particular time and depth level¶
[36]:
axes_class = (GeoAxes,dict(map_projection=cartopy.crs.Robinson(central_longitude=-159)))
fig = plt.figure(figsize=(14,8))
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')
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')
plt.show()
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
Time-mean distribution¶
[37]:
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].isel(k=k).mean('time'),
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()
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
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¶
[38]:
# 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()
[39]:
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¶
[40]:
# 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
[41]:
# compute vertical profiles at selected point and load into memory
G_total_Slt_atpt = G_total_Slt.isel(tile=t,j=j,i=i).compute()
G_advection_Slt_atpt = G_advection_Slt.isel(tile=t,j=j,i=i).compute()
G_diffusion_Slt_atpt = G_diffusion_Slt.isel(tile=t,j=j,i=i).compute()
G_forcing_Slt_atpt = G_forcing_Slt.isel(tile=t,j=j,i=i).compute()
rhs_atpt = G_advection_Slt_atpt + G_diffusion_Slt_atpt + G_forcing_Slt_atpt
[42]:
fig, axs = plt.subplots(2, 2, figsize=(14,8))
plt.sca(axs[0,0])
G_total_Slt_atpt.isel(k=k).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_atpt.isel(k=k).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_atpt.isel(k=k).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_atpt.isel(k=k).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¶
[43]:
fig = plt.subplots(1, 2, sharey=True, figsize=(12,7))
plt.subplot(1, 2, 1)
plt.plot(G_total_Slt_atpt.mean('time'), ecco_grid.Z,
lw=4, color='black', marker='.', label=r'$G^{Slt}_{total}$ (LHS)')
plt.plot(G_advection_Slt_atpt.mean('time'), ecco_grid.Z,
lw=2, color='red', marker='.', label=r'$G^{Slt}_{advection}$')
plt.plot(G_diffusion_Slt_atpt.mean('time'), ecco_grid.Z,
lw=2, color='orange', marker='.', label=r'$G^{Slt}_{diffusion}$')
plt.plot(G_forcing_Slt_atpt.mean('time'), ecco_grid.Z,
lw=2, color='blue', marker='.', label=r'$G^{Slt}_{forcing}$')
plt.plot(rhs_atpt.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_atpt.mean('time'), ecco_grid.Z,
lw=4, color='black', marker='.', label=r'$G^{Slt}_{total}$ (LHS)')
plt.plot(rhs_atpt.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)
[44]:
# 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.
[46]:
# Total tendency (psu/s)
SALT_diff = ds.SALT_snp.diff('time_snp')
SALT_diff = SALT_diff.rename({'time_snp':'time'})
del SALT_diff.time.attrs['c_grid_axis_shift'] # remove attribute from DataArray
SALT_diff = SALT_diff.assign_coords(time=ds.time) # correct time coordinate values
G_total_Sln = SALT_diff/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)
[47]:
# 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
# Set fluxes on land to zero (instead of NaN)
u_transport = u_transport.where(ecco_grid.hFacW.values > 0,0)
v_transport = v_transport.where(ecco_grid.hFacS.values > 0,0)
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
).
[48]:
# 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.
[49]:
# Vertical volume transport (m^3/s)
w_transport = ds.WVELMASS.where(ds.k_l>0,0.) * ecco_grid.rA
[50]:
# Set land values of flux to zero (instead of NaN)
w_transport = w_transport.where(ecco_grid.hFacC.values > 0,0)
# 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
).
[51]:
# Vertical convergence of salinity (psu m^3/s)
adv_vConvSln = (-ds.SALT*vConvV + adv_vConvS)/rstarfac
Total advective salinity convergence¶
[52]:
# 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.
[53]:
# 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)
[54]:
# 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).isel(k=slice(1,None)))],
dim='k')
[55]:
# Sea surface forcing for salinity (psu/s)
G_forcing_Sln = (-ds.SALT*forcV + G_forcing_Slt)/rstarfac
Salinity budget: Map of residual¶
[56]:
# 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
[57]:
# individual month weights
time_month_weights = delta_t/delta_t.sum()
# Accumulated residual
resSln = (((((rhs_Sln-G_total_Sln)*vol).sum(dim='k'))/vol.sum(dim='k'))\
*time_month_weights).sum(dim='time').compute()
[58]:
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()
-179.9 179.9
-180.0 180.0
-89.9 89.9
-90.0 90.0
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¶
[59]:
# 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}$']
[60]:
# Set an index for the time (t) and depth (k) axis
t, k = 100, 0
Example maps at a particular time and depth level¶
[61]:
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()
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
Time-mean distribution¶
[62]:
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].isel(k=k).mean('time'),
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()
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
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¶
[63]:
# 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()
[64]:
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 salinity budget closure¶
[65]:
# 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)
# compute vertical profiles at selected point and load into memory
G_total_Sln_atpt = G_total_Sln.isel(tile=t,j=j,i=i).compute()
G_advection_Sln_atpt = G_advection_Sln.isel(tile=t,j=j,i=i).compute()
G_diffusion_Sln_atpt = G_diffusion_Sln.isel(tile=t,j=j,i=i).compute()
G_forcing_Sln_atpt = G_forcing_Sln.isel(tile=t,j=j,i=i).compute()
rhs_Sln_atpt = G_advection_Sln_atpt + G_diffusion_Sln_atpt + G_forcing_Sln_atpt
12 0 87 16
[69]:
fig, axs = plt.subplots(2, 2, figsize=(14,8))
plt.sca(axs[0,0])
G_total_Sln_atpt.isel(k=k).plot(color='k',lw=2)
rhs_Sln_atpt.isel(k=k).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_atpt.isel(k=k).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_atpt.isel(k=k).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_atpt.isel(k=k).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.
[71]:
fig = plt.subplots(1, 2, sharey=True, figsize=(12,7))
plt.subplot(1, 2, 1)
plt.plot(G_total_Sln_atpt.mean('time'), ecco_grid.Z,
lw=4, color='black', marker='.', label=r'$G^{Sln}_{total}$ (LHS)')
plt.plot(G_advection_Sln_atpt.mean('time'), ecco_grid.Z,
lw=2, color='red', marker='.', label=r'$G^{Sln}_{advection}$')
plt.plot(G_diffusion_Sln_atpt.mean('time'), ecco_grid.Z,
lw=2, color='orange', marker='.', label=r'$G^{Sln}_{diffusion}$')
plt.plot(G_forcing_Sln_atpt.mean('time'), ecco_grid.Z,
lw=2, color='blue', marker='.', label=r'$G^{Sln}_{forcing}$')
plt.plot(rhs_Sln_atpt.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_atpt.mean('time'), ecco_grid.Z,
lw=4, color='black', marker='.', label=r'$G^{Sln}_{total}$ (LHS)')
plt.plot(rhs_Sln_atpt.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¶
Fresh water is defined as:
where is salinity, is a reference value.
As with the salt and a salinity budget we will evaluate each term in the freshwater budget.
Each term is largely analagous to the salinity budget (including the scale factor to account for volume changes), except that we must calculate as the residual of the other terms. We will then compare with .
[72]:
# Reference salinity
Sref = 35.0
Total freshwater tendency¶
As with salinity, we calculated the monthly-averaged time tendancy of freshwater by differencing monthly snapshots.
[76]:
f = (Sref - ds.SALT_snp)/Sref
# Total freshwater tendency (m^3/s)
f_diff = f.diff('time_snp')
f_diff = f_diff.rename({'time_snp':'time'})
del f_diff.time.attrs['c_grid_axis_shift'] # remove attribute from DataArray
f_diff = f_diff.assign_coords(time=ds.time) # correct time coordinate values
G_total_Fw = f_diff*vol/delta_t
Advective flux of freshwater¶
Advective fluxes of freshwater are calculated offline using salinity and velocity fields to calculate the volume convergence of freshwater:
is the residual mean velocity field, which contains both the resolved (Eulerian), as well as the Gent-McWilliams bolus velocity (i.e., the parameterization of unresolved eddy effects).
GM Bolus Velocity¶
[81]:
# Set values on land to zero (instead of NaN)
ds['GM_PsiX'] = ds.GM_PsiX.where(ecco_grid.hFacW.values > 0,0)
ds['GM_PsiY'] = ds.GM_PsiY.where(ecco_grid.hFacS.values > 0,0)
UVELSTAR = grid.diff(ds.GM_PsiX, 'Z', boundary='fill')/ecco_grid.drF
VVELSTAR = grid.diff(ds.GM_PsiY, 'Z', boundary='fill')/ecco_grid.drF
[82]:
GM_PsiXY_diff = grid.diff_2d_vector({'X' : ds.GM_PsiX*ecco_grid.dyG,
'Y' : ds.GM_PsiY*ecco_grid.dxG}, boundary = 'fill')
WVELSTAR = (GM_PsiXY_diff['X'] + GM_PsiXY_diff['Y'])/ecco_grid.rA
Calculate advective freshwater flux¶
[83]:
SALT_at_u = grid.interp(ds.SALT, 'X', boundary='extend')
SALT_at_v = grid.interp(ds.SALT, 'Y', boundary='extend')
SALT_at_w = grid.interp(ds.SALT, 'Z', boundary='extend')
As with the salinity budget, we zero out the top cell of WVELMASS
to avoid double counting the surface flux
[84]:
# Freshwater advective (Eulerian+Bolus) fluxes (m^3/s)
ADVx_FW = (ds.UVELMASS+UVELSTAR)*ecco_grid.dyG*ecco_grid.drF*(Sref-SALT_at_u)/Sref
ADVy_FW = (ds.VVELMASS+VVELSTAR)*ecco_grid.dxG*ecco_grid.drF*(Sref-SALT_at_v)/Sref
ADVr_FW = (ds.WVELMASS.where(ds.k_l>0).fillna(0.)+WVELSTAR)*ecco_grid.rA*(Sref-SALT_at_w)/Sref
[85]:
# set fluxes on land to zero (instead of NaN)
ADVx_FW = ADVx_FW.where(ecco_grid.hFacW.values > 0,0)
ADVy_FW = ADVy_FW.where(ecco_grid.hFacS.values > 0,0)
ADVr_FW = ADVr_FW.where(ecco_grid.hFacC.values > 0,0)
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']))
[86]:
# Convergence of vertical advection (m^3/s)
adv_vConvFw = grid.diff(ADVr_FW, 'Z', boundary='fill')
[87]:
# Sum horizontal and vertical convergences (m^3/s)
G_advection_Fw = (adv_hConvFw + adv_vConvFw)/rstarfac
Freshwater forcing¶
Include salinity forcing as follows:
so
[88]:
# Freshwater forcing (m^3/s)
forcFw = ds.oceFWflx/rhoconst*ecco_grid.rA
# Expand to fully 3d (using G_advection_Fw as template)
forcing_Fw = xr.concat([forcFw.reset_coords(drop=True).assign_coords(k=0).expand_dims('k'),
xr.zeros_like(G_advection_Fw).isel(k=slice(1,None))],
dim='k').where(ecco_grid.hFacC==1)
# Sum FW and Salinity forcing, changing G_forcing_Slt from [m psu/s] to [m^3/s]
G_forcing_Fw = (forcing_Fw-G_forcing_Slt*ecco_grid.rA/rhoconst/Sref)/rstarfac
Diffusive freshwater flux¶
We calculate the diffusive freshwater flux as the residuals of the remaining budget terms as this term is not output by ECCO
[89]:
# Convergence of freshwater diffusion (m^3/s)
G_diffusion_Fw = G_total_Fw - G_forcing_Fw - G_advection_Fw
Freshwater budget: Spatial distributions¶
[90]:
# In order to plot the budget terms in one figure, let's add them in a list
var = [G_total_Fw,G_advection_Fw,G_diffusion_Fw,G_forcing_Fw]
varstrngs = [r'$G^{Fw}_{total}$',r'$G^{Fw}_{advection}$',r'$G^{Fw}_{diffusion}$',r'$G^{Fw}_{forcing}$']
[91]:
# Set an index for the time (t) and depth (k) axis
t, k = 100, 0
Example maps at a particular time and depth level¶
[92]:
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 freshwater 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/1e3, transform=cartopy.crs.PlateCarree(),
vmin=-2, vmax=2, 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$^{3}$ m$^3$ s$^{-1}$', fontsize=14, fontweight='bold')
plt.show()
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
Time-mean distribution¶
[94]:
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 freshwater 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].isel(k=k).mean('time'),
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/1e3, 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$^{3}$ m$^3$ s$^{-1}$', fontsize=14, fontweight='bold')
plt.show()
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0
Freshwater budget through time¶
We cannot evaluate the closure of the fw budget as we have used the residual of available terms to determine the diffusive flux. We can compare how close the derived FW diffusive flux is to the scaled Salinity diffusive flux
Global average budget¶
[99]:
# Take volume-weighted mean of these terms
tmp_a = (G_total_Fw*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_b = (G_advection_Fw*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_d = (G_forcing_Fw*vol).sum(dim=('k','i','j','tile'))/vol.sum()
# tmp_c1 = (G_diffusion_Fw*vol).sum(dim=('k','i','j','tile'))/vol.sum()
tmp_c1 = tmp_a - tmp_d - tmp_b
tmp_c2 = (-G_diffusion_Sln/Sref*vol*vol).sum(dim=('k','i','j','tile'))/vol.sum()
[100]:
fig, axs = plt.subplots(2, 2, figsize=(14,8))
plt.sca(axs[0,0])
tmp_a.plot(color='k')
axs[0,0].set_title(r'a. $G^{Fw}_{total}$ [m$^3$ s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[0,1])
tmp_b.plot(color='r')
axs[0,1].set_title(r'b. $G^{Fw}_{advection}$', fontsize=12)
plt.grid()
plt.sca(axs[1,0])
tmp_c1.plot(color='orange',lw=2)
tmp_c2.plot(color='m')
axs[1,0].set_title(r'c. $G^{Fw}_{diffusion}$ (orange) & -$G^{Sln}_{diffusion}$/$S_{\mathrm{ref}}$ (pink) [m$^3$ s$^{-1}$] [m$^3$ s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[1,1])
tmp_d.plot(color='b')
axs[1,1].set_title(r'd. $G^{Fw}_{forcing}$ [m$^3$ s$^{-1}$]', fontsize=12)
plt.grid()
plt.subplots_adjust(hspace = .5, wspace=.2)
plt.suptitle('Global Freshwater Budget', fontsize=16)
plt.show()
Local freshwater budget closure¶
[101]:
# 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
[102]:
G_total_Fw_atpt = G_total_Fw.isel(tile=t,j=j,i=i).compute()
G_advection_Fw_atpt = G_advection_Fw.isel(tile=t,j=j,i=i).compute()
G_forcing_Fw_atpt = G_forcing_Fw.isel(tile=t,j=j,i=i).compute()
# G_diffusion_Fw_atpt = G_diffusion_Fw.isel(tile=t,j=j,i=i).compute()
negG_diffusion_Sln_vol_divSref_atpt = (-G_diffusion_Sln*vol/Sref).isel(tile=t,j=j,i=i).compute()
G_diffusion_Fw_atpt = G_total_Fw_atpt - G_forcing_Fw_atpt - G_advection_Fw_atpt
[103]:
fig, axs = plt.subplots(2, 2, figsize=(14,8))
plt.sca(axs[0,0])
G_total_Fw_atpt.isel(k=k).plot(color='k')
axs[0,0].set_title(r'a. $G^{Fw}_{total}$ [m$^3$ s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[0,1])
G_advection_Fw_atpt.isel(k=k).plot(color='r')
axs[0,1].set_title(r'b. $G^{Fw}_{advection}$ [m$^3$ s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[1,0])
G_diffusion_Fw_atpt.isel(k=k).plot(color='orange')
negG_diffusion_Sln_vol_divSref_atpt.isel(k=k).plot(color='m')
axs[1,0].set_title(r'c. $G^{Fw}_{diffusion}$ (orange) & -$G^{Sln}_{diffusion}$/$S_{\mathrm{ref}}$ (pink) [m$^3$ s$^{-1}$]', fontsize=12)
plt.grid()
plt.sca(axs[1,1])
G_forcing_Fw_atpt.isel(k=k).plot(color='b')
axs[1,1].set_title(r'd. $G^{Fw}_{forcing}$ [m$^3$ s$^{-1}$]', fontsize=12)
plt.grid()
plt.subplots_adjust(hspace = .5, wspace=.2)
plt.suptitle('Freshwater 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 freshwater budget terms¶
This section illustrates the balance in the freshwater budget along the depth axis.
[104]:
fig = plt.subplots(1, 2, sharey=True, figsize=(12,7))
plt.subplot(1, 2, 1)
plt.plot(G_total_Fw_atpt.mean('time'), ecco_grid.Z,
lw=4, color='black', marker='.', label=r'$G^{Fw}_{total}$ (LHS)')
plt.plot(G_advection_Fw_atpt.mean('time'), ecco_grid.Z,
lw=2, color='red', marker='.', label=r'$G^{Fw}_{advection}$')
plt.plot(G_diffusion_Fw_atpt.mean('time'), ecco_grid.Z,
lw=2, color='orange', marker='.', label=r'$G^{Fw}_{diffusion}$')
plt.plot(G_forcing_Fw_atpt.mean('time'), ecco_grid.Z,
lw=2, color='blue', marker='.', label=r'$G^{Fw}_{forcing}$')
plt.xlabel(r'Tendency [m$^3$ 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_Fw_atpt.mean('time'), ecco_grid.Z,
lw=4, color='black', marker='.', label=r'$G^{Fw}_{total}$ (LHS)')
plt.setp(plt.gca(), 'yticklabels',[])
plt.xlabel(r'Tendency [m$^3$ s$^{-1}$]', fontsize=14)
plt.ylim([-200,0])
plt.suptitle('Vertical profile of freshwater budget at a specific grid point (tile = %i, k = %i, j = %i, i = %i)'%(t,k,j,i), fontsize=16)
plt.show()
Save budget terms¶
Now that we have all the terms evaluated, let’s save them to a dataset. Here are two examples:
Add budget variables to a new dataset¶
[105]:
#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})
[106]:
ds.time.encoding = {}
ds = ds.reset_coords(drop=True)
Save to zarr dataset¶
[107]:
from dask.diagnostics import ProgressBar
[108]:
ds
[108]:
<xarray.Dataset> Dimensions: (i: 90, j: 90, tile: 13, k: 50, time: 288) Coordinates: * i (i) int32 0 1 2 3 4 5 6 7 8 ... 81 82 83 84 85 86 87 88 89 * j (j) int32 0 1 2 3 4 5 6 7 8 ... 81 82 83 84 85 86 87 88 89 * 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 * time (time) datetime64[ns] 1993-01-16T12:00:00 ... 2016-12-16... Data variables: G_total_Sln (time, k, tile, j, i) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray> G_advection_Sln (time, k, tile, j, i) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray> G_diffusion_Sln (time, k, tile, j, i) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray> G_forcing_Sln (time, k, tile, j, i) float32 dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
- i: 90
- j: 90
- tile: 13
- k: 50
- time: 288
- i(i)int320 1 2 3 4 5 6 ... 84 85 86 87 88 89
- axis :
- X
- long_name :
- grid index in x for variables at tracer and 'v' locations
- swap_dim :
- XC
- comment :
- In the Arakawa C-grid system, tracer (e.g., THETA) and 'v' variables (e.g., VVEL) have the same x coordinate on the model grid.
- coverage_content_type :
- coordinate
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])
- j(j)int320 1 2 3 4 5 6 ... 84 85 86 87 88 89
- axis :
- Y
- long_name :
- grid index in y for variables at tracer and 'u' locations
- swap_dim :
- YC
- comment :
- In the Arakawa C-grid system, tracer (e.g., THETA) and 'u' variables (e.g., UVEL) have the same y coordinate on the model grid.
- coverage_content_type :
- coordinate
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])
- tile(tile)int320 1 2 3 4 5 6 7 8 9 10 11 12
- long_name :
- lat-lon-cap tile index
- comment :
- The ECCO V4 horizontal model grid is divided into 13 tiles of 90x90 cells for convenience.
- coverage_content_type :
- coordinate
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12])
- k(k)int320 1 2 3 4 5 6 ... 44 45 46 47 48 49
- axis :
- Z
- long_name :
- grid index in z for tracer variables
- swap_dim :
- Z
- coverage_content_type :
- coordinate
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])
- time(time)datetime64[ns]1993-01-16T12:00:00 ... 2016-12-...
- long_name :
- center time of averaging period
- axis :
- T
- bounds :
- time_bnds
- coverage_content_type :
- coordinate
- standard_name :
- time
array(['1993-01-16T12:00:00.000000000', '1993-02-15T00:00:00.000000000', '1993-03-16T12:00:00.000000000', ..., '2016-10-16T12:00:00.000000000', '2016-11-16T00:00:00.000000000', '2016-12-16T12:00:00.000000000'], dtype='datetime64[ns]')
- G_total_Sln(time, k, tile, j, i)float32dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
Array Chunk Bytes 5.65 GiB 20.08 MiB Shape (288, 50, 13, 90, 90) (1, 50, 13, 90, 90) Count 2918 Tasks 288 Chunks Type float32 numpy.ndarray - G_advection_Sln(time, k, tile, j, i)float32dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
Array Chunk Bytes 5.65 GiB 20.08 MiB Shape (288, 50, 13, 90, 90) (1, 50, 13, 90, 90) Count 5822245 Tasks 288 Chunks Type float32 numpy.ndarray - G_diffusion_Sln(time, k, tile, j, i)float32dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
Array Chunk Bytes 5.65 GiB 20.08 MiB Shape (288, 50, 13, 90, 90) (1, 50, 13, 90, 90) Count 4627904 Tasks 288 Chunks Type float32 numpy.ndarray - G_forcing_Sln(time, k, tile, j, i)float32dask.array<chunksize=(1, 50, 13, 90, 90), meta=np.ndarray>
Array Chunk Bytes 5.65 GiB 20.08 MiB Shape (288, 50, 13, 90, 90) (1, 50, 13, 90, 90) Count 14983 Tasks 288 Chunks Type float32 numpy.ndarray
[110]:
# specify directory to save budget terms in: currently set to ~/Downloads
save_dir = join(user_home_dir,'Downloads')
with ProgressBar():
ds.to_zarr(join(save_dir,'eccov4r4_budg_Sln')) # change name of file as you wish
[########################################] | 100% Completed | 14m 31s
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.
[111]:
# Suppress warning messages for a cleaner presentation
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import xarray as xr
from os.path import expanduser,join
import sys
user_home_dir = expanduser('~')
sys.path.append(join(user_home_dir,'ECCOv4-py')) # change to directory hosting ecco_v4_py as needed
import ecco_v4_py as ecco
import matplotlib.pyplot as plt
%matplotlib inline
# Define a high-level directory for ECCO fields
ECCO_dir = join(user_home_dir,'Downloads','ECCO_V4r4_PODAAC')
# Load the model grid
grid_shortname = "ECCO_L4_GEOMETRY_LLC0090GRID_V4R4"
grid_filename = "GRID_GEOMETRY_ECCO_V4r4_native_llc0090.nc"
ecco_grid = xr.open_dataset(join(ECCO_dir,grid_shortname,grid_filename))
# Volume (m^3)
vol = (ecco_grid.rA*ecco_grid.drF*ecco_grid.hFacC).transpose('tile','k','j','i')
# year range
year_start = 1993
year_end = 2017
# ShortNames of monthly mean datasets needed
FW_surf_flux_shortname = "ECCO_L4_FRESH_FLUX_LLC0090GRID_MONTHLY_V4R4"
vol_flux_shortname = "ECCO_L4_OCEAN_3D_VOLUME_FLUX_LLC0090GRID_MONTHLY_V4R4"
bolus_strmfcn_shortname = "ECCO_L4_OCEAN_BOLUS_STREAMFUNCTION_LLC0090GRID_MONTHLY_V4R4"
salt_monthly_shortname = "ECCO_L4_TEMP_SALINITY_LLC0090GRID_MONTHLY_V4R4"
# open datasets where monthly mean files are located
datasets_to_open = ['FW_surf_flux','vol_flux','bolus_strmfcn',\
'salt_monthly']
for dataset in datasets_to_open:
curr_shortname = eval(dataset + '_shortname')
curr_dir = join(ECCO_dir,curr_shortname)
curr_files_in_range = files_in_year_range(curr_dir,year_start,year_end)
exec('ds_' + dataset + ' = xr.open_mfdataset(curr_files_in_range,parallel=True,data_vars=\'minimal\','\
+ 'coords=\'minimal\', compat=\'override\')')
# create one dataset with all the monthly mean variables we need
ecco_monthly_mean = xr.merge([ds_FW_surf_flux['oceFWflx'],\
ds_vol_flux[['UVELMASS','VVELMASS','WVELMASS']],\
ds_bolus_strmfcn[['GM_PsiX','GM_PsiY']],\
ds_salt_monthly['SALT']])
[112]:
# specify directory budget terms are saved in: currently set to ~/Downloads
save_dir = join(user_home_dir,'Downloads')
# Load terms from zarr dataset
file_to_open = join(save_dir,'eccov4r4_budg_Sln')
time_month_weights = xr.open_zarr(file_to_open).time_month_weights
G_total_Sln = xr.open_zarr(file_to_open).G_total_Sln
G_advection_Sln = xr.open_zarr(file_to_open).G_advection_Sln
G_diffusion_Sln = xr.open_zarr(file_to_open).G_diffusion_Sln
G_forcing_Sln = xr.open_zarr(file_to_open).G_forcing_Sln
[113]:
# set values on land to zero (instead of NaN)
ecco_monthly_mean['GM_PsiX'] = ecco_monthly_mean.GM_PsiX.where(ecco_grid.hFacW.values > 0,0)
ecco_monthly_mean['GM_PsiY'] = ecco_monthly_mean.GM_PsiY.where(ecco_grid.hFacS.values > 0,0)
grid = ecco.get_llc_grid(ecco_monthly_mean)
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 = ADVr_SLT.chunk({'time':1,'tile':13,'k_l':50,'j':90,'i':90})
#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
Comparison between LHS and RHS of the budget equation¶
[114]:
# Total convergence
ConvSln = G_advection_Sln + G_diffusion_Sln
ConvSln2 = G_advection + G_diffusion_Sln
[115]:
# Sum of terms in RHS of equation
rhsSln = ConvSln + G_forcing_Sln
rhsSln2 = ConvSln2 + G_forcing_Sln
Map of residuals¶
[116]:
resSln = (((((G_advection_Sln + G_diffusion_Sln + G_forcing_Sln - G_total_Sln)*vol)\
.sum(dim='k'))/vol.sum(dim='k'))*time_month_weights).sum(dim='time').compute()
resSln2 = (((((rhsSln2 - G_total_Sln)*vol)\
.sum(dim='k'))/vol.sum(dim='k'))*time_month_weights).sum(dim='time').compute()
[117]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, resSln,
cmin=-1e-9, cmax=1e-9, 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()
-179.9 179.9
-180.0 180.0
-89.9 89.9
-90.0 90.0
[118]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, resSln2,
cmin=-1e-9, cmax=1e-9, 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()
-179.9 179.9
-180.0 180.0
-89.9 89.9
-90.0 90.0
[119]:
# depth integral, time mean of DataArray
def depthint_tmean(da):
da = da.chunk({'time':1,'k':50,'tile':13,'j':90,'i':90})
da_depthint_tmean = (((da*vol).sum('k')/vol.sum('k'))*time_month_weights).sum('time')
return da_depthint_tmean
G_advection_depthint_tmean = depthint_tmean(G_advection).compute()
G_advection_Sln_depthint_tmean = depthint_tmean(G_advection_Sln).compute()
G_total_Sln_depthint_tmean = depthint_tmean(G_total_Sln).compute()
G_diffusion_Sln_depthint_tmean = depthint_tmean(G_diffusion_Sln).compute()
G_forcing_Sln_depthint_tmean = depthint_tmean(G_forcing_Sln).compute()
[120]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_advection_depthint_tmean,\
cmin=-1e-9, cmax=1e-9, show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Advection [psu s$^{-1}$]', fontsize=16)
plt.show()
-179.9 179.9
-180.0 180.0
-89.9 89.9
-90.0 90.0
[121]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_advection_Sln_depthint_tmean,\
cmin=-1e-9, cmax=1e-9, show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Advection [psu s$^{-1}$]', fontsize=16)
plt.show()
-179.9 179.9
-180.0 180.0
-89.9 89.9
-90.0 90.0
[122]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_total_Sln_depthint_tmean,
cmin=-1e-9, cmax=1e-9, show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Total tendency [psu s$^{-1}$]', fontsize=16)
plt.show()
-179.9 179.9
-180.0 180.0
-89.9 89.9
-90.0 90.0
[123]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_diffusion_Sln_depthint_tmean,
cmin=-1e-9, cmax=1e-9,show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Diffusion [psu s$^{-1}$]', fontsize=16)
plt.show()
-179.9 179.9
-180.0 180.0
-89.9 89.9
-90.0 90.0
[124]:
plt.figure(figsize=(15,5))
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, G_forcing_Sln_depthint_tmean,
cmin=-1e-9, cmax=1e-9,show_colorbar=True, cmap='RdBu_r',dx=0.2, dy=0.2)
plt.title(r'Forcing [psu s$^{-1}$]', fontsize=16)
plt.show()
-179.9 179.9
-180.0 180.0
-89.9 89.9
-90.0 90.0
Histogram of residuals¶
We can look at the distribution of residuals to get a little more confidence.
[125]:
from xhistogram.xarray import histogram
[126]:
res_closed = np.abs(G_advection_Sln + G_diffusion_Sln + G_forcing_Sln - G_total_Sln)
res_closed.name = 'Residual_closed'
[127]:
res_offline = np.abs(G_advection + G_diffusion_Sln + G_forcing_Sln - G_total_Sln)
res_offline.name = 'Residual_offline'
[128]:
res_adv = np.abs(G_diffusion_Sln + G_forcing_Sln - G_total_Sln)
res_adv.name = 'Residual_advection'
[129]:
res_dif = np.abs(G_advection_Sln + G_forcing_Sln - G_total_Sln)
res_dif.name = 'Residual_diffusion'
[130]:
res_frc = np.abs(G_advection_Sln + G_diffusion_Sln - G_total_Sln)
res_frc.name = 'Residual_forcing'
[131]:
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()
plt.legend(labels=['Closed residual','Res w/o adv','Res w/o diff','Res w/o forcing'])
[131]:
<matplotlib.legend.Legend at 0x18442501eb0>
Note that the fact that the closed residual has the highest number in the low-value bins is a good thing; since all of these arrays have the same number of elements, this means that far fewer elements of the res_closed
array have higher values:
[132]:
fig,ax = plt.subplots(1,1,figsize=(7,7))
histogram(res_closed, bins = [np.linspace(0, 1.e-7,1000)]).plot()
histogram(res_dif, bins = [np.linspace(0, 1.e-7,1000)]).plot()
histogram(res_frc, bins = [np.linspace(0, 1.e-7,1000)]).plot()
histogram(res_adv, bins = [np.linspace(0, 1.e-7,1000)]).plot()
ax.set_yscale('log')
plt.legend(labels=['Closed residual','Res w/o adv','Res w/o diff','Res w/o forcing'])
[132]:
<matplotlib.legend.Legend at 0x18443fe3760>