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:

  1. Salt budget

  2. Salinity budget

  3. 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 (\theta), the budget is described with salinity (S).

\frac{\partial S}{\partial t} = -\nabla \cdot (S \mathbf{u}) - \nabla\cdot\mathbf{F}_\textrm{diff}^{S} + {F}_\textrm{forc}^{S}

The total tendency (\frac{\partial S}{\partial t}) is equal to advective convergence (-\nabla \cdot (S \mathbf{u})), diffusive flux convergence (-\nabla \cdot \mathbf{F}_\textrm{diff}^{S}) and a forcing term {F}_\textrm{forc}^{S}.

In the case of ECCOv4, salt is strictly a conserved mass and can be described as

\underbrace{\frac{\partial(s^* S)}{\partial t}}_{G^{Slt}_\textrm{total}} = \underbrace{-\nabla_{z^{*}} \cdot(s^* S \, \mathbf{v}_{res}) - \frac{\partial(S\,w_{res})}{\partial z^* }}_{G^{Slt}_\textrm{advection}}\underbrace{- s^* ({\nabla\cdot\mathbf{F}_\textrm{diff}^{S}})}_{G^{Slt}_\textrm{diffusion}} + \underbrace{s^* {F}_\textrm{forc}^{S}}_{G^{Slt}_\textrm{forcing}}

The change in salt content over time (G^{Slt}_\textrm{total}) is equal to the convergence of the advective flux (G^{Slt}_\textrm{advection}) and diffusive flux (G^{Slt}_\textrm{diffusion}) plus a forcing term associated with surface salt exchanges (G^{Slt}_\textrm{forcing}). As with the heat budget, we present both the horizontal (\mathbf{v}_{res}) and vertical (w_{res}) components of the advective term. Again, we have \mathbf{u}_{res} 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 z^* and the scale factor s^* 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, G^{Slt}_\textrm{total} (i.e., the left side of the salt budget equation) can be extended as follows

\frac{\partial(s^* S)}{\partial t} = s^* \frac{\partial S}{\partial t} + S \frac{ \partial s^* }{\partial t}

When substituting G^{Slt}_\textrm{total} with the right hand side of the above equation, we can solve for the salinity tendency (\frac{\partial S}{\partial t}):

\frac{\partial S}{\partial t} = -\frac{1}{s^*} \,\left[S\,\frac{\partial s^* }{\partial t} + \nabla_{z^*}\cdot(s^* S\,\mathbf{v}_{res}) + \frac{\partial(S\,w_{res})}{\partial z^*}\right] - \nabla \cdot \mathbf{F}_\textrm{diff}^{S} + F_\textrm{forc}^{S}

Since s^* = 1 + \frac{\eta}{H} we can define the temporal change in s^* as

\frac{\partial s^*}{\partial t} =  \frac{1}{H}\,\frac{\partial \eta}{\partial t}

This constitutes the conservation of volume in ECCOv4, which can be formulated as

\frac{1}{H}\,\frac{\partial \eta}{\partial t} = -\nabla_{z^* } \cdot (s^*\mathbf{v}) - \frac{\partial w}{\partial z^* } + \mathcal{F}

You can read more about volume conservation and the z^* coordinate system in another tutorial. \mathcal{F} denotes the volumetric surface fluxes and can be decomposed into net atmospheric freshwater fluxes (i.e., precipitation minus evaporation, P - E), continental runoff (R) and exchanges due to sea ice melting/formation (I). Here \mathbf{v} = (u, v) and w are the resolved horizontal and vertical velocities, respectively.

Thus, the conservation of salinity in ECCOv4 can be described as

\underbrace{\frac{\partial S}{\partial t}}_{G^{Sln}_\textrm{total}} = \underbrace{\frac{1}{s^* }\,\left[S\,\left(\nabla_{z^* } \cdot (s^* \mathbf{v}) + \frac{\partial w}{\partial z^* }\right) - \nabla_{z^* } \cdot (s^* S \, \mathbf{v}_{res}) - \frac{\partial(S\,w_{res})}{\partial z^* }\right]}_{G^{Sln}_\textrm{advection}} \underbrace{- \nabla \cdot \mathbf{F}_\textrm{diff}^{S}}_{G^{Sln}_\textrm{diffusion}} + \underbrace{F_\textrm{forc}^{S} - S\,\mathcal{F}}_{G^{Sln}_\textrm{forcing}}

Notice here that, in contrast to the salt budget equation, the salinity equation explicitly includes the surface forcing (S\,\mathcal{F}). \mathcal{F} represents surface freshwater exchanges (P-E+R-I) and F_{\textrm{forc}}^{S} represents surface salt fluxes (i.e., addition/removal of salt). Besides the convergence of the advective flux (\nabla\cdot(S\mathbf{u}_{res})), the salinity equation also includes the convergence of the volume flux multiplied by the salinity (S\,(\nabla\cdot\mathbf{u})), 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 S from a given reference salinity S_{ref}. Thus, within a control volume V the freshwater content is defined as a volume (V_{fw}):

V_{fw} = \iiint_V\frac{S_{ref} - S}{S_{ref}}dV

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:

\underbrace{\frac{\partial V_{fw}}{\partial t}}_{G^{fw}_\textrm{total}} = \underbrace{-\nabla \cdot \mathbf{F}_\textrm{adv}^{fw}}_{G^{fw}_\textrm{advection}} \underbrace{- \nabla \cdot \mathbf{F}_\textrm{diff}^{fw}}_{G^{fw}_\textrm{diffusion}} + \underbrace{\mathcal{F}}_{G^{fw}_\textrm{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
import sys
import glob

from os.path import expanduser,join
user_home_dir = expanduser('~')


# indicate whether you are working in a cloud instance (True if yes, False otherwise)
incloud_access = True
[2]:
# Suppress warning messages for a cleaner presentation
import warnings
warnings.filterwarnings('ignore')
[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.

[ ]:
## if working in the AWS cloud, access datasets needed for this tutorial

ShortNames_list = ["ECCO_L4_GEOMETRY_LLC0090GRID_V4R4",\
                   "ECCO_L4_OCEAN_3D_SALINITY_FLUX_LLC0090GRID_MONTHLY_V4R4",\
                   "ECCO_L4_OCEAN_3D_VOLUME_FLUX_LLC0090GRID_MONTHLY_V4R4",\
                   "ECCO_L4_OCEAN_BOLUS_STREAMFUNCTION_LLC0090GRID_MONTHLY_V4R4",\
                   "ECCO_L4_FRESH_FLUX_LLC0090GRID_MONTHLY_V4R4",\
                   "ECCO_L4_SSH_LLC0090GRID_SNAPSHOT_V4R4",\
                   "ECCO_L4_TEMP_SALINITY_LLC0090GRID_SNAPSHOT_V4R4"]
if incloud_access == True:
    # apply fixes to xgcm (you might have to do this manually if not on a cloud instance,
    # see "Horizontal convergence of advective heat flux" section)
    import xgcm
    xgcm_path = xgcm.__path__[0]
    if xgcm.__version__ == '0.8.1':
        import os
        os.remove(join(xgcm_path,'grid_ufunc.py'))
        os.system('cp ../misc/xgcm_fixes/grid_ufunc.py '+xgcm_path+'/')
        os.remove(join(xgcm_path,'padding.py'))
        os.system('cp ../misc/xgcm_fixes/padding.py '+xgcm_path+'/')

        # re-import xgcm
        import importlib
        importlib.reload(xgcm)

    from ecco_s3_retrieve import ecco_podaac_s3_get_diskaware
    files_dict = ecco_podaac_s3_get_diskaware(ShortNames=ShortNames_list,\
                                              StartDate='1993-01',EndDate='2016-12',\
                                              max_avail_frac=0.5,\
                                              download_root_dir=ECCO_dir)
[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.

sys.path.append(join(user_home_dir,'ECCOv4-py'))   # change to directory hosting ecco_v4_py as needed

import ecco_v4_py as ecco
[7]:
## Load the model grid
if incloud_access == True:
    ecco_grid = xr.open_dataset(files_dict[ShortNames_list[0]])
else:
    ecco_grid = xr.open_dataset(glob.glob(join(ECCO_dir,'*GEOMETRY*','*.nc'))[0])

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
if incloud_access == True:
    ds_etan_snaps = xr.open_mfdataset(files_dict[ShortNames_list[-2]],\
                                      data_vars='minimal',coords='minimal',compat='override')
    ds_salt_snaps = xr.open_mfdataset(files_dict[ShortNames_list[-1]],\
                                      data_vars='minimal',coords='minimal',compat='override')
else:
    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']
if incloud_access == True:
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 \Delta t, 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

G^{Slt}_\textrm{total} = G^{Slt}_\textrm{advection} + G^{Slt}_\textrm{diffusion} + G^{Slt}_\textrm{forcing}

The total tendency of salt (G^{Slt}_\textrm{total}) is the sum of the salt tendencies from advective convergence (G^{Slt}_\textrm{advection}), diffusive heat convergence (G^{Slt}_\textrm{diffusion}) and total forcing (G^{Slt}_\textrm{forcing}).

We present calculations sequentially for each term starting with G^{Slt}_\textrm{total} 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 s^*S.

[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^{Slt}_\textrm{total} (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 calling diff_2d_vector below. If you have this issue, as a temporary fix you can locate the xgcm directory

import xgcm
print(xgcm.__path__)

and then substitute padding.py and grid_ufunc.py with the files linked below:

padding.py

grid_ufunc.py

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) 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.

[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
_images/ECCO_v4_Salt_and_salinity_budget_72_1.png

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
_images/ECCO_v4_Salt_and_salinity_budget_78_1.png

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
_images/ECCO_v4_Salt_and_salinity_budget_80_1.png

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()
_images/ECCO_v4_Salt_and_salinity_budget_85_0.png

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()
_images/ECCO_v4_Salt_and_salinity_budget_90_0.png

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 G^{Slt}_{advection} and G^{Slt}_{diffusion} 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()
_images/ECCO_v4_Salt_and_salinity_budget_93_0.png

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.

G^{Sln}_\textrm{total} = G^{Sln}_\textrm{advection} + G^{Sln}_\textrm{diffusion} + G^{Sln}_\textrm{forcing}

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 (s^*) to be used in our calculations below.

s^* = 1 + \frac{\eta}{H}

This requires following model output: - Depth: Ocean depth, H (m) - ETAN: Surface Height Anomaly, \eta (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.

G^{Sln}_\textrm{total} = \frac{\Delta S}{\Delta t}

[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
_images/ECCO_v4_Salt_and_salinity_budget_121_1.png

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
_images/ECCO_v4_Salt_and_salinity_budget_127_1.png

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
_images/ECCO_v4_Salt_and_salinity_budget_129_1.png

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()
_images/ECCO_v4_Salt_and_salinity_budget_134_0.png

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()
_images/ECCO_v4_Salt_and_salinity_budget_137_0.png

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()
_images/ECCO_v4_Salt_and_salinity_budget_140_0.png

Evaluating the freshwater budget

Fresh water is defined as:

fw = \frac{S_{ref} - S}{S_{ref}},

where S is salinity, S_{ref} is a reference value.

As with the salt and a salinity budget we will evaluate each term in the freshwater budget.

G^{fw}_\textrm{total} = G^{fw}_\textrm{advection} + G^{fw}_\textrm{diffusion} + G^{fw}_\textrm{forcing}

Each term is largely analagous to the salinity budget (including the scale factor s^* to account for volume changes), except that we must calculate G^{fw}_\textrm{diffusion} as the residual of the other terms. We will then compare with G^{Sln}_\textrm{diffusion}.

[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:

\mathbf{\mathcal{F}_{adv}} = \iint_A\mathbf{u}_{res} \cdot \left(\frac{S_{ref} - S}{S_{ref}}\right)dA

u_{res} 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:

\frac{\partial fw}{\partial t} = -\frac{1}{S_{ref}}\frac{\partial S}{\partial t}

so G^{fw}_\textrm{forcing}= -\frac{1}{S_{ref}} G^{Sln}_\textrm{forcing}

[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
_images/ECCO_v4_Salt_and_salinity_budget_165_1.png

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
_images/ECCO_v4_Salt_and_salinity_budget_167_1.png

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 G^{fw}_\textrm{diffusion} is to the scaled Salinity diffusive flux G^{Sln}_\textrm{diffusion}

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()
_images/ECCO_v4_Salt_and_salinity_budget_171_0.png

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()
_images/ECCO_v4_Salt_and_salinity_budget_175_0.png

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()
_images/ECCO_v4_Salt_and_salinity_budget_178_0.png

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>
[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
_images/ECCO_v4_Salt_and_salinity_budget_196_1.png
[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
_images/ECCO_v4_Salt_and_salinity_budget_197_1.png
[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
_images/ECCO_v4_Salt_and_salinity_budget_199_1.png
[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
_images/ECCO_v4_Salt_and_salinity_budget_200_1.png
[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
_images/ECCO_v4_Salt_and_salinity_budget_201_1.png
[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
_images/ECCO_v4_Salt_and_salinity_budget_202_1.png
[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
_images/ECCO_v4_Salt_and_salinity_budget_203_1.png

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>
_images/ECCO_v4_Salt_and_salinity_budget_211_1.png

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>
_images/ECCO_v4_Salt_and_salinity_budget_213_1.png