ECCOv4 Global Volume Budget Closure

Here we demonstrate the closure of volume budgets in ECCOv4 configurations. This notebook is draws heavily from evaluating_budgets_in_eccov4r3.pdf which explains the procedure with Matlab code examples by Christopher G. Piecuch. See ECCO Version 4 release documents: /doc/evaluating_budgets_in_eccov4r3.pdf).

Objectives

Illustrate how volume budgets are closed globally.

Introduction

ECCOv4 uses the z^* coordinate system in which the depth of the vertical coordinate, z^* varies with time as:

z^* = \frac{z - \eta(x,y,t)}{H(x,y) + \eta(x,y,t)} H(x,y)

With H being the model depth, \eta being the model sea level anomaly, and z being depth.

If the vertical coordinate didn’t change through time then volume fluxes across the ‘u’ and ‘v’ grid cell faces of a tracer cell could be calculated by multiplying the velocities at the face with the face area:

volume flux across ‘u’ face in the +x direction = \mathit{UVEL}(x,y,k) \times \mathit{drF}(k) \times \mathit{dyG}(x,y) \times \mathit{hFacW}(x,y,k)

volume flux across ‘v’ face in the +y direction = \mathit{VVEL}(x,y,k) \times \mathit{drF}(k) \times \mathit{dxG}(x,y) \times \mathit{hFacS}(x,y,k)

With dyG and dxG being the lengths of the ‘u’ and ‘v’ faces, drF being the grid cell height and hFacW and hFacS being the vertical fractions of the ‘u’ and ‘v’ grid cell faces that are open water (ECCOv4 uses partial cells to better represent bathymetry which can allows 0 < hfac \le 1).

However, because the vertical coordiate varies with time in the z^* system, the grid cell height drF varies with time as drF\times s^*(t), with

s^*(x,y,k,t) = 1 + \frac{\eta(x,y,t)}{H}

with s^* > 1 when \eta > 0

Thus, to calculate the volume fluxes grid cell through horizontal faces we must account for the time-varying grid cell face areas:

volume flux across ‘u’ in the +x direction face with z^* coordinates = \mathit{UVEL}(x,y,k) \times \mathit{drF}(k) \times \mathit{dyG}(x,y) \times \mathit{hFacW}(x,y,k) \times s^*(x,y,k,t)

volume flux across ‘v’ in the +y direction face with z^* coordinates = \mathit{VVEL}(x,y,k) \times \mathit{drF}(k) \times \mathit{dxG}(x,y) \times \mathit{hFacS}(x,y,k) \times s^*(x,y,k,t)

To make budget calculations easier we provide the scaled velocities quantities UVELMASS and VVELMASS,

\mathit{UVELMASS}(x,y,k) = \mathit{UVEL}(x,y,k) \times \mathit{hFacW}(x,y,k) \times s^{*}(x,y,k,t)

and

\mathit{VVELMASS}(x,y,k) = \mathit{VVEL}(x,y,k) \times \mathit{hFacS}(x,y,k) \times s^{*}(x,y,k,t)

It is worth noting that the word mass in UVELMASS and VVELMASS is confusing since there is no mass involved here. Think of these terms as simply being UVEL and VVEL multiplied by the fraction of the grid cell height that is open grid cell face across which the volume transport occurs. Partial cell bathymetry can make this fraction (hFacW, hFacS) less than one, and the s^* scaling factor further adjusts this fraction higher or lower through time.

Fully closing the budget requires the vertical volume fluxes across the top and bottom ‘w’ faces of the grid cell and surface freshwater fluxes. Regarding vertical volume fluxes, there are no s^* or hFac equivalent scaling factors that modify our top and bottom grid cell areas. Therefore, vertical volume fluxes through ‘w’ faces are simply:

volume flux across ‘w’ face in the +z direction = \mathit{WVEL}(x,y,k) \times \mathit{rA}(x,y)

Note: Inexplicably, the term WVEL is provided with the silly name WVELMASS. Sometimes it’s difficult to ignore other people’s poor life choices, but please try to do so here. Ignore the confusing name, WVELMASS is identical to WVEL.

In the z^* coordinate system the depth of the surface grid cell is always z^* = 0. In the MITgcm, WELMASS at the top of the surface grid cell is the liquid volume flux out of the ocean surface and is proportional to the vertical ocean mass flux, oceFWflx

ETAN in a Boussinesq Model

ECCOv4 uses a volume-conserving Boussinesq formulation of the MITgcm. Because volume is conserved in Boussinesq formulations, seawater density changes do not change model sea level anomaly, ETAN. The following demonstration ofETAN budget closure considers volumetric fluxes but does not take into consideration expansion/contraction due to changes in density. Furthermore, model ETAN changes with the exchange of water between ocean and sea-ice while in reality the Archimedes principle holds that sea level should not change following the growth or melting of sea-ice because floating sea-ice displaces a volume of seawater equal to its weight. Thus, ETAN is not comparable to observed sea level.

We correct ETAN to make a sea surface height field that is comaprable with observations by making three corrections: with a) the “Greatbatch correction”, a time varying, globally-uniform correction to ocean volume due to changes in global mean density, b) the inverted barometer (IB) correction (see SSHIBC) and c) the ‘sea ice load’ correction to account for the displacement of seawater due to submerged sea-ice and snow (see sIceLoad). A demonstration of these corrections is outside the scope of this tutorial. Here we focus on closing the model budget keeping in mind that we are neglecting sea level changes from changes in global mean density and the fact that ETAN does not account for volume displacement due to submerged sea-ice.

Greatbatch, 1994. J. of Geophys. Res. Oceans, https://doi.org/10.1029/94JC00847

Evaluating the model sea level anomaly ETAN volume budget

We will evalute

\underbrace{\frac{\partial \eta}{\partial t}}_{G_\text{total tendency}} = \underbrace{\int_{-H}^0 \left( -\nabla_{z^*}(s^*\,{\bf v}- \frac{\partial w}{\partial z^*} \right) dz^*}_{G_{\text{volumetric divergence}}} + \underbrace{F}_{{G_{surface fluxes}}}

The total tendency of \eta, G_{\text{total tendency}} is the sum of the \eta tendencies from volumetric divergence, G_{\text{volumetric divergence}}, and volumetric surface fluxes, G_{\text{surface fluxes}}.

In discrete form, using indexes that start from k=0 (surface tracer cell) and running to k=nk-1 (bottom tracer cell)

\frac{\eta(i,j)}{\partial t} = \sum_{k=nk-1}^{0} \underbrace{\left[\mathit{UVELMASS}(i_g,j,k)-\mathit{UVELMASS}(i_g+1,j,k) \right]\, \mathit{dyG}(i_g,j) \, \mathit{drF}(k) }_{\text{volumetric flux in minus out in x direction}}  +
\\
\sum_{k=nk-1}^{0} \underbrace{\left[\mathit{VVELMASS}(i,j_g,k)-\mathit{VVELMASS}(i,j_g+1,k) \right] \, \mathit{dxG}(i,j_g) \, \mathit{drF}(k)}_{\text{volumetric flux in minus out in y direction}} +
\\
\sum_{k_l=nk}^{1} \underbrace{\mathit{WVELMASS}(i,j,k_l)\, \mathit{drA}(i,j)}_{\text{volumetric flux through grid cell bottom surface}} +
\\
\underbrace{\mathit{oceFWflx}(i,j)/ \mathit{rhoConst}}_{\text{volumetric flux through  the top surface of the uppermost tracer cell}}

In the above we intentionally sum WVELMASS fluxes from the BOTTOM surface of the lowermost grid cell (at k_l = 50) to the BOTTOM face of the uppermost grid cell (k_l = 1) so that we can explicitly include the surface volume flux (forcing) term, \mathit{oceFWflx}(i,j)/\mathit{rhoConst}.

We will calculate \partial \eta / \partial t by differencing instantaneous monthly snapshots of \eta as

\frac{\partial \eta}{\partial t} = \frac{\eta(i,j,t+1) - \eta(i,j,t)}{\Delta t}

The UVELMASS, VVELMASS, WVELMASS and oceFWflx terms must be time-average quantities between the monthly \eta snapshots.

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 mean datasets for the full time span of ECCOv4r4 output (1992 through 2017), and 2 sets of snapshots at monthly intervals. The ShortNames of the datasets are:

  • ECCO_L4_GEOMETRY_LLC0090GRID_V4R4

  • ECCO_L4_OCEAN_3D_VOLUME_FLUX_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_OBP_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. If you are downloading the snapshots, you may also find the snaps_monthly_textlist function helpful; it generates a file list of only the snapshot files for the 1st of each month, which you can then download using wget.

Prepare environment and loading the relevant model variables

[1]:
import numpy as np
import sys
import xarray as xr
import glob
from copy import deepcopy
import matplotlib.pyplot as plt
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')

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

# indicate whether you are working in a cloud instance (True if yes, False otherwise)
incloud_access = False
[2]:
# setting up a dask LocalCluster
from dask.distributed import Client
from dask.distributed import LocalCluster
cluster = LocalCluster()
client = Client(cluster)
client
[2]:

Client

Client-fc1ec589-0ca3-11ef-8879-0659db0d1d75

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:8787/status

Cluster Info

[3]:
# Density kg/m^3
rhoconst = 1029

## needed to convert surface mass fluxes to volume fluxes

# lat/lon resolution in degrees to interpolate the model
# fields for the purposes of plotting
map_dx = .2
map_dy = .2
[4]:
## Set top-level file directory for the ECCO NetCDF files
## =================================================================

## currently set to ~/Downloads/ECCO_V4r4_PODAAC,
## the default if ecco_podaac_download was used to download dataset granules
ECCO_dir = join(user_home_dir,'Downloads','ECCO_V4r4_PODAAC')
[5]:
## if working in the AWS cloud, access datasets needed for this tutorial

ShortNames_list = ["ECCO_L4_GEOMETRY_LLC0090GRID_V4R4",\
                   "ECCO_L4_OCEAN_3D_VOLUME_FLUX_LLC0090GRID_MONTHLY_V4R4",\
                   "ECCO_L4_FRESH_FLUX_LLC0090GRID_MONTHLY_V4R4",\
                   "ECCO_L4_SSH_LLC0090GRID_SNAPSHOT_V4R4",\
                   "ECCO_L4_OBP_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 Volume Flux Divergence 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)
[6]:
## 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.  The example below adds
##    ecco_v4_py to the user's path if it is stored in the folder
##    ECCOv4-py under the user's home directory

sys.path.append(join(user_home_dir,'ECCOv4-py'))

import ecco_v4_py as ecco

Load ecco_grid

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

Open 2D MONTHLY \eta snapshots

[8]:
year_start = 1993
year_end = 2017


# open ETAN snapshots (beginning of each month)
if incloud_access == True:
    ecco_monthly_snaps_etan = (xr.open_mfdataset(files_dict[ShortNames_list[-2]],\
                                           data_vars='minimal',coords='minimal',compat='override')\
                                           ['ETAN']).to_dataset()
    ecco_monthly_snaps_obp = xr.open_mfdataset(files_dict[ShortNames_list[-1]],
                                           data_vars='minimal',coords='minimal',compat='override')
else:
    ecco_monthly_snaps_etan = (xr.open_mfdataset(join(ECCO_dir,'*SSH*SNAPSHOT*','*-01T*.nc'),\
                                           data_vars='minimal',coords='minimal',compat='override')\
                                           ['ETAN']).to_dataset()
    ecco_monthly_snaps_obp = xr.open_mfdataset(join(ECCO_dir,'*OBP*SNAPSHOT*','*-01T*.nc'),
                                           data_vars='minimal',coords='minimal',compat='override')
ecco_monthly_snaps = xr.merge((ecco_monthly_snaps_etan,ecco_monthly_snaps_obp))

# time mask for snapshots
time_snap_mask = np.logical_and(ecco_monthly_snaps.time.values >= np.datetime64(str(year_start)+'-01-01','ns'),\
                                ecco_monthly_snaps.time.values < np.datetime64(str(year_end)+'-01-02','ns'))

ecco_monthly_snaps = ecco_monthly_snaps.isel(time=time_snap_mask)
[9]:
# print time range of monthly snapshots
print(ecco_monthly_snaps.time.isel(time=[0, -1]).values)
['1993-01-01T00:00:00.000000000' '2017-01-01T00:00:00.000000000']
[10]:
# find the record of the last ETAN snapshot
last_record_date = ecco_monthly_snaps.time[-1].values
print(last_record_date)
last_record_year = str(last_record_date)[:4]
2017-01-01T00:00:00.000000000

Open MONTHLY mean data

[11]:
## Load ECCO variables
if incloud_access == True:
    ecco_vars_uvw = xr.open_mfdataset(files_dict[ShortNames_list[1]],\
                                      data_vars='minimal',coords='minimal',compat='override')
    ecco_vars_fflx = xr.open_mfdataset(files_dict[ShortNames_list[2]],\
                                      data_vars='minimal',coords='minimal',compat='override')
else:
    ecco_vars_uvw = xr.open_mfdataset(join(ECCO_dir,'*_OCEAN_3D_VOLUME_FLUX*MONTHLY*','*.nc'),\
                                      data_vars='minimal',coords='minimal',compat='override')
    ecco_vars_fflx = xr.open_mfdataset(join(ECCO_dir,'*_FRESH_FLUX*MONTHLY*','*.nc'),\
                                      data_vars='minimal',coords='minimal',compat='override')
ecco_monthly_mean = xr.merge((ecco_vars_uvw,ecco_vars_fflx['oceFWflx']))

# time mask for monthly means
time_mean_mask = np.logical_and(ecco_monthly_mean.time.values >= np.datetime64(str(year_start)+'-01-01','ns'),\
                                ecco_monthly_mean.time.values < np.datetime64(str(year_end)+'-01-01','ns'))

ecco_monthly_mean = ecco_monthly_mean.isel(time=time_mean_mask)
[12]:
# first and last monthly-mean records
print(ecco_monthly_mean.time.isel(time=[0, -1]).values)
['1993-01-16T12:00:00.000000000' '2016-12-16T12:00:00.000000000']
[13]:
# each monthly mean record is bookended by a snapshot.
# we should have one more snapshot than monthly mean record
print('number of monthly mean records: ', len(ecco_monthly_mean.time))
print('number of monthly snapshot records: ', len(ecco_monthly_snaps.time))
number of monthly mean records:  288
number of monthly snapshot records:  289

Create the xgcm ‘grid’ object

the xgcm grid object makes it easy to make flux divergence calculations across different tiles of the lat-lon-cap grid.

[14]:
ecco_xgcm_grid = ecco.get_llc_grid(ecco_grid)
ecco_xgcm_grid
[14]:
<xgcm.Grid>
Z Axis (not periodic, boundary=None):
  * center   k --> left
  * right    k_u --> center
  * left     k_l --> center
  * outer    k_p1 --> center
Y Axis (not periodic, boundary=None):
  * center   j --> left
  * left     j_g --> center
X Axis (not periodic, boundary=None):
  * center   i --> left
  * left     i_g --> center

Calculate LHS: \eta time tendency: G_{\text{total tendency}}

We calculate the monthly-averaged time tendency of ETAN by differencing monthly ETAN snapshots. Subtract the numpy arrays \eta(t+1) - \eta(t). This operation gives us \Delta ETAN / \Delta t (month) records.

[15]:
num_months = len(ecco_monthly_snaps.time)
G_total_tendency_month = \
    ecco_monthly_snaps.ETAN.isel(time=range(1,num_months)).values - \
    ecco_monthly_snaps.ETAN.isel(time=range(0,num_months-1)).values

# The result is a numpy array of 264 months
print('shape of G_total_tendency_month: ', G_total_tendency_month.shape)
shape of G_total_tendency_month:  (288, 13, 90, 90)
[16]:
ecco_monthly_mean.oceFWflx.shape
[16]:
(288, 13, 90, 90)
[17]:
# Convert this numpy array to an xarray Dataarray to take advantage
# of xarray time indexing.
# The easiest way is to copy an existing DataArray that has the
# dimensions and time indexes that we want,
# replace its values, and change its name.

tmp = ecco_monthly_mean.oceFWflx.copy(deep=True)
tmp.values = G_total_tendency_month
tmp.name = 'G_total_tendency_month'
G_total_tendency_month = tmp

# the nice thing is that now the time values of G_total_tendency_month now line
# up with the time values of the time-mean fields (middle of the month)
print('\ntime of first array of G_total_tendency_month');
print(G_total_tendency_month.time[0].values)

time of first array of G_total_tendency_month
1993-01-16T12:00:00.000000000

Now convert \Delta ETAN / \Delta t (month) to \Delta ETAN / \Delta t (seconds) by dividing by the number of seconds in each month.

[18]:
secs_per_month = (ecco_monthly_snaps.time[1:].values - \
ecco_monthly_snaps.time[0:-1].values).astype('float')/(1.e9)   #divide by 1.e9 to convert ns to seconds


# Make a DataArray with the number of seconds in each month,
#time indexed to the times in dETAN_dT_perMonth (middle of each month)
secs_per_month = xr.DataArray(secs_per_month, \
                              coords={'time': G_total_tendency_month.time.values}, \
                              dims='time')

# show number of seconds in the first two months:
print('# of seconds in Jan and Feb 1993 ', secs_per_month[0:2].values)

# sanity check: show number of days in the first two months:
print('# of days in Jan and Feb 1993 ', secs_per_month[0:2].values/3600/24)
# of seconds in Jan and Feb 1993  [2678400. 2419200.]
# of days in Jan and Feb 1993  [31. 28.]
[19]:
# convert dETAN_dT_perMonth to perSeconds
G_total_tendency = G_total_tendency_month / secs_per_month

Plot the time-mean \partial \eta / \partial t, total \Delta \eta, and one example \partial \eta / \partial t field

Time-mean \partial \eta / \partial t

To calculate the time averaged G_{total_tendency} one might be tempted to do the following

> G_total_tendency_mean = G_total_tendency.mean('time')

But that would be folly because the ‘mean’ function does not know that the number of days in each month is different! The result would downweight Februarys and upweight Julys. We have to weight the tendency records by the length of each month. A clever way of doing that is provided in the xarray documents: https://xarray-test.readthedocs.io/en/latest/examples/monthly-means.html

In our case we know the length of each month, we just calculated it above in secs_per_month. We will weight each month by the relative # of seconds in each month and sum to get a weighted average.

[20]:
# the weights are just the # of seconds per month divided by total seconds
month_length_weights = secs_per_month / secs_per_month.sum()

The time mean of the ETAN tendency, \overline{G_{\text{total tendency}}}, is given by

\overline{G_{\text{total tendency}}} = \sum_{i=1}^{nm} w_i G_{\text{total tendency}}

with \sum_{i=1}^{nm} w_i = 1 and nm=number of months

[21]:
# the weights sum to 1
print(month_length_weights.sum())
<xarray.DataArray ()> Size: 8B
array(1.)
[22]:
# the weighted mean weights by the length of each month (in seconds)
G_total_tendency_mean = (G_total_tendency*month_length_weights).sum('time')
[23]:
plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC,\
                              G_total_tendency_mean,show_colorbar=True,\
                              cmin=-1e-9, cmax=1e-9, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);
plt.title('Average $\partial \eta / \partial t$ [m/s]', fontsize=20);
_images/ECCO_v4_Volume_budget_closure_33_0.png

Total \Delta \eta

The time average eta tendency is small, about 1 billionth of a meter per second. The ECCO period is coming up to a billion seconds though… How much did ETAN change over the analysis period?

[24]:
# the number of seconds in the entire period
seconds_in_entire_period = \
    float(ecco_monthly_snaps.time[-1] - ecco_monthly_snaps.time[0])/1e9
print ('seconds in analysis period: ', seconds_in_entire_period)

# which is also the sum of the number of seconds in each month
print('sum of seconds in each month ', secs_per_month.sum().values)
seconds in analysis period:  757382400.0
sum of seconds in each month  757382400.0
[25]:
ETAN_delta = G_total_tendency_mean*seconds_in_entire_period
[26]:
plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              ETAN_delta,show_colorbar=True,\
                              cmin=-1, cmax=1, \
                              cmap='RdBu_r', user_lon_0=-67, \
                              dx=map_dx,dy=map_dy);
plt.title('Predicted $\Delta \eta$ [m] from $\partial \eta / \partial t$', \
          fontsize=20);
_images/ECCO_v4_Volume_budget_closure_37_0.png

We can sanity check the total ETAN change that we found by multipling the time-mean ETAN tendency with the number of seconds in the simulation by comparing that with the difference in ETAN between the end of the last month and start of the first month.

[27]:
ETAN_delta_method_2 = ecco_monthly_snaps.ETAN.isel(time=-1).values - \
    ecco_monthly_snaps.ETAN.isel(time=0).values
[28]:
plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              ETAN_delta_method_2,
                              show_colorbar=True,\
                              cmin=-1, cmax=1, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);

plt.title('Actual $\Delta \eta$ [m]', fontsize=20);
_images/ECCO_v4_Volume_budget_closure_40_0.png
[29]:
plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              ETAN_delta_method_2-ETAN_delta,\
                              show_colorbar=True,\
                              cmin=-1e-6, cmax=1e-6, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);
plt.title('Difference between actual and predicted $\Delta \eta$ [m]', \
          fontsize=20);
_images/ECCO_v4_Volume_budget_closure_41_0.png

That’s a big woo, these are the same to within 10^-6 meters!

Example \partial \eta / \partial t field

[30]:
plt.figure(figsize=(20,8));

curr_t_ind = 100
curr_time = G_total_tendency.time[curr_t_ind].values
print(curr_time)
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              G_total_tendency.isel(time=curr_t_ind),\
                              show_colorbar=True,\
                              cmin=-1e-7, cmax=1e-7,\
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);

plt.title('$\partial \eta / \partial t$ [m/s] during ' +
          str(curr_time)[0:4] +'/' + str(curr_time)[5:7], fontsize=20);
2001-05-16T12:00:00.000000000
_images/ECCO_v4_Volume_budget_closure_44_1.png

For any given month the time rate of change of ETAN is two orders of magnitude smaller than the 1993-2015 mean. In the above we are looking at May 2001. We see positive ETAN tendency due sea ice melting in the northern hemisphere (e.g., Baffin Bay, Greenland Sea, and Chukchi Sea).

Calculate RHS: \eta tendency due to surface fluxes, G_\text{surface fluxes}

Surface mass fluxes are given in oceFWflx. Convert surface mass flux to a vertical velocity by dividing by the reference density rhoConst= 1029 kg m-3

[31]:
# tendency of eta implied by surface volume fluxes (m/s)
G_surf_fluxes = ecco_monthly_mean.oceFWflx/rhoconst

Plot the time-mean, total, and one month average of G_{\text{surface fluxes}}

Time-mean G_{\text{surface fluxes}}

We calculate the time-mean surface flux \eta tendency using the same weights as the total \eta tendency.

[32]:
G_surf_fluxes_mean= (G_surf_fluxes*month_length_weights).sum('time')
G_surf_fluxes_mean.shape
[32]:
(13, 90, 90)
[33]:
plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              G_surf_fluxes_mean,
                              show_colorbar=True,\
                              cmin=-1e-7, cmax=1e-7, \
                              cmap='RdBu_r', user_lon_0=-67,
                              dx=map_dx,dy=map_dy)
plt.title('Average $\partial \eta / \partial t$ [m/s] implied by oceFWflx surface mass fluxes\n Negative = Water out',
          fontsize=20);
_images/ECCO_v4_Volume_budget_closure_50_0.png

Total \Delta \eta due to surface fluxes

If there were no other terms on the RHS to balance surface fluxes, the total change in ETAN between 1993 and 2017 would be order of 10s of meters almost everwhere.

[34]:
ETAN_delta_surf_fluxes = G_surf_fluxes_mean*seconds_in_entire_period
[35]:
plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              ETAN_delta_surf_fluxes,show_colorbar=True,\
                              cmin=-100, cmax=100, \
                              cmap='RdBu_r', user_lon_0=-67, \
                              dx=map_dx,dy=map_dy);
plt.title('$\Delta \eta$ [m] implied by oceFWflx surface mass fluxes',
          fontsize=20);
_images/ECCO_v4_Volume_budget_closure_53_0.png

Example \partial \eta / \partial t impliedy by surface fluxes

[36]:
plt.figure(figsize=(20,8));

curr_t_ind = 100
curr_time = G_surf_fluxes.time[curr_t_ind].values
print(curr_time)
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              G_surf_fluxes.isel(time=curr_t_ind),show_colorbar=True,\
                             cmin=-1e-7, cmax=1e-7, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);
plt.title('$\partial \eta / \partial t$ [m/s] implied by oceFWflx surface mass fluxes ' +
          str(curr_time)[0:4] +'-' + str(curr_time)[5:7], fontsize=20);
2001-05-16T12:00:00.000000000
_images/ECCO_v4_Volume_budget_closure_55_1.png

For any given month the time rate of change of ETAN is almost the same as its 22 year mean. Differences are largest in the high latitudes where sea-ice melt and growth during any particular month induce large changes in ETAN.

Calculate RHS: \eta tendency due to volumetric flux divergence, G_\text{volumetric fluxes}

First we will look at vertical volumetric flux divergence, then horizontal volumetric flux divergence.

Vertical volumetric flux divergence

It turns out that WVELMASS at k_l=0 (the top face of the top tracer cell) is proportional to the ocean surface mass flux oceFWflx. The vertical velocity of the ocean surface is equal to the rate at which water is being added or removed across the top surface of the uppermost grid cell. This is demonstrated by differencing the velocity at the top ‘w’ face of the uppermost tracr cell WVELMASS (k_l = 0) and the velocity equivalent of transporting the surface mass flux term oceFWFlx through this same face.

First, find the time-mean vertical velocity at the liquid ocean surface

[37]:
WVELMASS_surf_mean = \
    (ecco_monthly_mean.WVELMASS.isel(k_l=0)*month_length_weights).sum('time')

Next, find the time-mean vertical velocity implied by the oceFWflx at k_l=0:

[38]:
WVEL_from_oceFWflx_mean = \
    -(ecco_monthly_mean.oceFWflx*month_length_weights).sum('time')/rhoconst
[39]:
plt.figure(figsize=(15,15))
#plt.sca(axs[0,0])
F=ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                                WVELMASS_surf_mean,\
                                show_colorbar=True,\
                                cmin=-1e-7, cmax=1e-7, \
                                cmap='RdBu_r', user_lon_0=-67,\
                                dx=2,dy=2, subplot_grid=[3,1,1]);
F[1].set_title('A: surface velocity from WVEL [m/s]')

F=ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                                WVEL_from_oceFWflx_mean,\
                                show_colorbar=True,\
                                cmin=-1e-7, cmax=1e-7, \
                                cmap='RdBu_r', user_lon_0=-67,\
                                dx=2,dy=2, subplot_grid=[3,1,2])
F[1].set_title("B: surface velocity implied by oceFWflx [m/s]")

F=ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                                WVELMASS_surf_mean-WVEL_from_oceFWflx_mean,\
                                show_colorbar=True,\
                                cmin=-1e-7, cmax=1e-7, \
                                cmap='RdBu_r', user_lon_0=-67,\
                                dx=2,dy=2, subplot_grid=[3,1,3])
F[1].set_title("difference between A and B");

_images/ECCO_v4_Volume_budget_closure_62_0.png

WVELMASS at the surface evidently is the same as the surface velocity implied by the surface mass flux oceFWflx. Thus, we do not actually need oceFWflx to close the volume budget. However, to keep the surface forcing term explicitly represented, we will keep oceFWflx and instead zero out the values of WVELMASS at the surface so as to avoid double counting.

Calculate the vertical volume fluxes at each level: (velocity x grid cell area) [m3 s-1]

[40]:
vol_transport_z = ecco_monthly_mean.WVELMASS * ecco_grid.rA

Set the volume transport at the surface level to be zero because we already counted the fluxes out of the domain with oceFWflx. We will also set the volume transport in land areas (currently NaN values) to zero so the divergence adjacent to the bottom is computed correctly.

[41]:
# vol_transport_z.isel(k_l=0).values[:] = 0
# vol_transport_z.values[ecco_grid.hFacC.values == 0] = 0

# ## For some reason the above commands don't seem to work with dask arrays,
# ## probably due to some complication resetting the values of an existing dask array.
## However, using where works with dask arrays
vol_transport_z = vol_transport_z.where(np.expand_dims(ecco_grid.k_l.values,axis=(0,-1,-2,-3)) != 0,\
                                        0)
vol_transport_z = vol_transport_z.where(ecco_grid.hFacC.values > 0,0)

Each grid cell has a top and bottom surface and therefore WVELMASSshould have 51 vertical levels (one more than the number of tracer cells). For some reason we only have 50 vertical levels, with the bottom of the 50th tracer cell missing. To calculate vertical flux divergence we need to add this 51st WVELMASS which is everywhere zero (no volume flux from the seafloor). The xgcm library handles this in its diff routine by specifying the boundary=’fill’ and fill_value = 0.

[42]:
# re-chunk xarray DataArray for better performance
vol_transport_z = vol_transport_z.chunk({'time':12,'k_l':vol_transport_z.sizes['k_l'],'tile':1})
[43]:
# this computation can take a little while, so we'll loop through 12 time indices (1 year)
# at a time and show when each computation has been completed

import time

G_vertical_flux_divergence_depth_integrated_time_mean = xr.DataArray(\
    data=np.zeros(vol_transport_z.shape[-3:]),dims=['tile','j','i'])
size_superchunk = 12
for time_superchunk in range(int(np.ceil(vol_transport_z.sizes['time']/size_superchunk))):
    time_log = [time.time()]

    curr_time_ind = np.arange(size_superchunk*time_superchunk,\
                              np.fmin(size_superchunk*(time_superchunk+1),\
                                      vol_transport_z.sizes['time']))


    # volume flux divergence into each grid cell, m^3 / s
    vol_vert_divergence = ecco_xgcm_grid.diff(vol_transport_z.isel(time=curr_time_ind),axis='Z',\
                                              boundary='fill',fill_value=0).squeeze()

    # change in eta per unit time due to volumetric vertical convergence
    # at each depth level: m/s
    G_vertical_flux_divergence = vol_vert_divergence / ecco_grid.rA

    # change in eta per unit time due to vertical integral of
    # volumetric horizonal convergence: m/s
    G_vertical_flux_divergence_depth_integrated = G_vertical_flux_divergence.sum('k')

    # Calculate the time-mean surface flux $\eta$ tendency using
    # the same weights as the total $\eta$ tendency.
    G_vertical_flux_divergence_depth_integrated_time_mean.values += \
        (G_vertical_flux_divergence_depth_integrated * month_length_weights[curr_time_ind]).sum('time').values

    time_log = np.append(time_log,time.time())

    print(f'Completed calculation for time ind = {curr_time_ind[[0,-1]]}, '\
            +f'time elapsed = {np.diff(time_log)[0]} s')
Completed calculation for time ind = [ 0 11], time elapsed = 9.466320753097534 s
Completed calculation for time ind = [12 23], time elapsed = 8.571366786956787 s
Completed calculation for time ind = [24 35], time elapsed = 8.247650384902954 s
Completed calculation for time ind = [36 47], time elapsed = 8.209043025970459 s
Completed calculation for time ind = [48 59], time elapsed = 8.605789422988892 s
Completed calculation for time ind = [60 71], time elapsed = 8.250980615615845 s
Completed calculation for time ind = [72 83], time elapsed = 8.65609860420227 s
Completed calculation for time ind = [84 95], time elapsed = 8.249817848205566 s
Completed calculation for time ind = [ 96 107], time elapsed = 8.15756893157959 s
Completed calculation for time ind = [108 119], time elapsed = 8.413262367248535 s
Completed calculation for time ind = [120 131], time elapsed = 8.668389558792114 s
Completed calculation for time ind = [132 143], time elapsed = 7.98816442489624 s
Completed calculation for time ind = [144 155], time elapsed = 8.698549747467041 s
Completed calculation for time ind = [156 167], time elapsed = 8.800422191619873 s
Completed calculation for time ind = [168 179], time elapsed = 8.268079996109009 s
Completed calculation for time ind = [180 191], time elapsed = 9.052908420562744 s
Completed calculation for time ind = [192 203], time elapsed = 7.919287443161011 s
Completed calculation for time ind = [204 215], time elapsed = 9.170875072479248 s
Completed calculation for time ind = [216 227], time elapsed = 8.453101873397827 s
Completed calculation for time ind = [228 239], time elapsed = 9.024392127990723 s
Completed calculation for time ind = [240 251], time elapsed = 8.407300472259521 s
Completed calculation for time ind = [252 263], time elapsed = 9.076092720031738 s
Completed calculation for time ind = [264 275], time elapsed = 8.12934160232544 s
Completed calculation for time ind = [276 287], time elapsed = 8.642892122268677 s

Plot the time-mean G_{\text{vertical flux divergence}}

[44]:
plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC,\
                              G_vertical_flux_divergence_depth_integrated_time_mean,\
                              show_colorbar=True,\
                              cmap='RdBu_r', \
                              cmin=-1e-9, cmax=1e-9,
                              dx=map_dx,dy=map_dy);
plt.title('Average $\partial \eta / \partial t$ [m/s] due to vertical flux divergence\n Negative = Water out',
          fontsize=20);
_images/ECCO_v4_Volume_budget_closure_71_0.png

These values are everywhere essentially zero (numerical noise). On average, the only vertical flux divergence in the column is across the ocean surface. Below the surface, the sum of vertical flux divergence in all tracer cells in the column must be zero because any divergence in any one particular cell is exactly offset by convergence in another cell. Net convergence into the column manifests as a positive vertical velocity at the surface which is equal to oceFWflux in the time-mean. Thanks to Hong Zhang for comments that improved this explanation.

Horizontal Volume Flux Divergence

[45]:
# Volumetric transports in x and y(m^3/s)
vol_transport_x = ecco_monthly_mean.UVELMASS * ecco_grid.dyG * ecco_grid.drF
vol_transport_y = ecco_monthly_mean.VVELMASS * ecco_grid.dxG * ecco_grid.drF
[46]:
# Set fluxes on land to zero (instead of NaN)
vol_transport_x = vol_transport_x.where(ecco_grid.hFacW.values > 0,0)
vol_transport_y = vol_transport_y.where(ecco_grid.hFacS.values > 0,0)

# re-chunk arrays for better performance
vol_transport_x = vol_transport_x.chunk({'time':1,'k':-1,'tile':-1,'j':-1,'i_g':-1})
vol_transport_y = vol_transport_y.chunk({'time':1,'k':-1,'tile':-1,'j_g':-1,'i':-1})
[47]:
def diff_2d_flux_llc90(flux_vector_dict):
    """
    A function (similar to xgcm.Grid.diff_2d_vector) that differences flux variables on the llc90 grid.
    Can be used in case of errors using xgcm's diff_2d_vector.
    """

    u_flux = flux_vector_dict['X']
    v_flux = flux_vector_dict['Y']

    u_flux_padded = u_flux.pad(pad_width={'i_g':(0,1)},mode='constant',constant_values=np.nan)\
                            .chunk({'i_g':u_flux.sizes['i_g']+1})
    v_flux_padded = v_flux.pad(pad_width={'j_g':(0,1)},mode='constant',constant_values=np.nan)\
                            .chunk({'j_g':v_flux.sizes['j_g']+1})

    def da_replace_at_indices(da,indexing_dict,replace_values):
        # replace values in xarray DataArray using locations specified by indexing_dict
        array_data = da.data
        indexing_dict_bynum = {}
        for axis,dim in enumerate(da.dims):
            if dim in indexing_dict.keys():
                indexing_dict_bynum = {**indexing_dict_bynum,**{axis:indexing_dict[dim]}}
        ndims = len(array_data.shape)
        indexing_list = [':']*ndims
        for axis in indexing_dict_bynum.keys():
            indexing_list[axis] = indexing_dict_bynum[axis]
        indexing_str = ",".join(indexing_list)

        # using exec isn't ideal, but this works for both NumPy and Dask arrays
        exec('array_data['+indexing_str+'] = replace_values')

        return da

    # u flux padding
    for tile in range(0,3):
        u_flux_padded = da_replace_at_indices(u_flux_padded,{'tile':str(tile),'i_g':'-1'},\
                                              u_flux.isel(tile=tile+3,i_g=0).data)
    for tile in range(3,6):
        u_flux_padded = da_replace_at_indices(u_flux_padded,{'tile':str(tile),'i_g':'-1'},\
                                              v_flux.isel(tile=12-tile,j_g=0,i=slice(None,None,-1)).data)
    u_flux_padded = da_replace_at_indices(u_flux_padded,{'tile':'6','i_g':'-1'},\
                                          u_flux.isel(tile=7,i_g=0).data)
    for tile in range(7,9):
        u_flux_padded = da_replace_at_indices(u_flux_padded,{'tile':str(tile),'i_g':'-1'},\
                                              u_flux.isel(tile=tile+1,i_g=0).data)
    for tile in range(10,12):
        u_flux_padded = da_replace_at_indices(u_flux_padded,{'tile':str(tile),'i_g':'-1'},\
                                              u_flux.isel(tile=tile+1,i_g=0).data)

    # v flux padding
    for tile in range(0,2):
        v_flux_padded = da_replace_at_indices(v_flux_padded,{'tile':str(tile),'j_g':'-1'},\
                                              v_flux.isel(tile=tile+1,j_g=0).data)
    v_flux_padded = da_replace_at_indices(v_flux_padded,{'tile':'2','j_g':'-1'},\
                                          u_flux.isel(tile=6,j=slice(None,None,-1),i_g=0).data)
    for tile in range(3,6):
        v_flux_padded = da_replace_at_indices(v_flux_padded,{'tile':str(tile),'j_g':'-1'},\
                                              v_flux.isel(tile=tile+1,j_g=0).data)
    v_flux_padded = da_replace_at_indices(v_flux_padded,{'tile':'6','j_g':'-1'},\
                                          u_flux.isel(tile=10,j=slice(None,None,-1),i_g=0).data)
    for tile in range(7,10):
        v_flux_padded = da_replace_at_indices(v_flux_padded,{'tile':str(tile),'j_g':'-1'},\
                                              v_flux.isel(tile=tile+3,j_g=0).data)
    for tile in range(10,13):
        v_flux_padded = da_replace_at_indices(v_flux_padded,{'tile':str(tile),'j_g':'-1'},\
                                              u_flux.isel(tile=12-tile,j=slice(None,None,-1),i_g=0).data)

    # take differences
    diff_u_flux = u_flux_padded.diff('i_g')
    diff_v_flux = v_flux_padded.diff('j_g')

    # include coordinates of input DataArrays and correct dimension/coordinate names
    diff_u_flux = diff_u_flux.assign_coords(u_flux.coords).rename({'i_g':'i'})
    diff_v_flux = diff_v_flux.assign_coords(v_flux.coords).rename({'j_g':'j'})

    diff_flux_vector_dict = {'X':diff_u_flux,'Y':diff_v_flux}

    return diff_flux_vector_dict
[48]:
# # This way of computing horizontal volume divergence has more concise code,
# # but it may not run optimally in some computing environments.
# # In that case, try doing the computation in "superchunks" of 12 months (1 year)
# # at a time, as coded in the next cell.


# # Difference of horizontal transports in x and y directions
# vol_flux_diff = ecco_xgcm_grid.diff_2d_vector({'X': vol_transport_x, \
#                                                'Y': vol_transport_y},\
#                                               boundary='fill')
#
# # volume flux divergence into each grid cell, m^3 / s
# vol_horiz_divergence = (vol_flux_diff['X'] + vol_flux_diff['Y'])
# # restore time coordinate to DataArray (lost in xgcm.diff_2d_vector operation)
# vol_horiz_divergence = vol_horiz_divergence.assign_coords({'time':vol_transport_x.time.data})
#
# # change in eta per unit time due to volumetric horizonal
# # convergence at each depth level: m/s
# # a positive DIVERGENCE leads to negative eta tendency
# G_vol_horiz_divergence = -vol_horiz_divergence / ecco_grid.rA
#
# # change in eta in each grid cell per unit time due to horiz. divergence: m/s
# G_vol_horiz_divergence_depth_integrated = G_vol_horiz_divergence.sum('k').compute()
#
# # calculate time-mean using the month length weights
# G_vol_horiz_divergence_depth_integrated_mean = \
#     (G_vol_horiz_divergence_depth_integrated * month_length_weights).sum('time')
[49]:
# complete horizontal divergence calculation, 12 time indices (1 year) at a time

import copy

area_masked = ecco_grid.rA.where(ecco_grid.hFacC.isel(k=0) > 0).compute()
area_masked_sum = area_masked.sum(dim=('i','j','tile')).values

G_vol_horiz_divergence_depth_integrated_globalmean_array = np.empty((vol_transport_x.sizes['time'],))
G_vol_horiz_divergence_depth_integrated_mean = xr.DataArray(\
    data=np.zeros(vol_transport_x.shape[-3:]),dims=['tile','j','i'])
size_superchunk = 12
for time_superchunk in range(int(np.ceil(vol_transport_x.sizes['time']/size_superchunk))):
    time_log = np.array([time.time()])

    curr_time_ind = np.arange(size_superchunk*time_superchunk,\
                              np.fmin(size_superchunk*(time_superchunk+1),\
                                      vol_transport_x.sizes['time']))

    # # Difference of horizontal transports in x and y directions
    # vol_flux_diff = ecco_xgcm_grid.diff_2d_vector({'X': vol_transport_x.isel(time=curr_time_ind), \
    #                                                'Y': vol_transport_y.isel(time=curr_time_ind)},\
    #                                               boundary='fill')
    vol_flux_diff = diff_2d_flux_llc90({'X': vol_transport_x.isel(time=curr_time_ind),\
                                         'Y': vol_transport_y.isel(time=curr_time_ind)})

    # volume flux divergence into each grid cell, m^3 / s
    vol_horiz_divergence = (vol_flux_diff['X'] + vol_flux_diff['Y'])

    # restore time coordinate to DataArray (lost in xgcm.diff_2d_vector operation)
    vol_horiz_divergence = vol_horiz_divergence.assign_coords({'time':vol_transport_x.time[curr_time_ind].data})

    # change in eta per unit time due to volumetric horizonal
    # convergence at each depth level: m/s
    # a positive DIVERGENCE leads to negative eta tendency
    curr_G_vol_horiz_divergence = -vol_horiz_divergence / ecco_grid.rA
    if curr_time_ind[0] == 0:
        G_vol_horiz_divergence = copy.deepcopy(curr_G_vol_horiz_divergence)
    else:
        G_vol_horiz_divergence = xr.concat((G_vol_horiz_divergence,\
                                            curr_G_vol_horiz_divergence),\
                                           dim='time',\
                                           coords='minimal',compat='override')

    # change in eta in each grid cell per unit time due to horiz. divergence: m/s
    curr_G_vol_horiz_divergence_depth_integrated = xr.DataArray(\
        data=curr_G_vol_horiz_divergence.sum('k'),dims=['time','tile','j','i'])
    if curr_time_ind[0] == 0:
        G_vol_horiz_divergence_depth_integrated = copy.deepcopy(curr_G_vol_horiz_divergence_depth_integrated)
    else:
        G_vol_horiz_divergence_depth_integrated = xr.concat((G_vol_horiz_divergence_depth_integrated,\
                                                             curr_G_vol_horiz_divergence_depth_integrated),\
                                                            dim='time',\
                                                            coords='minimal',compat='override')

    # compute global mean at each time using area weighting
    G_vol_horiz_divergence_depth_integrated_globalmean_array[curr_time_ind] = \
        (area_masked*curr_G_vol_horiz_divergence_depth_integrated).sum(dim=('i','j','tile')).values/area_masked_sum

    # calculate time-mean using the month length weights
    G_vol_horiz_divergence_depth_integrated_mean += \
        (curr_G_vol_horiz_divergence_depth_integrated*month_length_weights[curr_time_ind]).sum('time').values

    time_log = np.append(time_log,time.time())

    print(f'Completed calculation for time ind = {curr_time_ind[[0,-1]]}, '\
            +f'time elapsed = {np.diff(time_log)[0]} s')

# # create DataArray of global mean values
G_vol_horiz_divergence_depth_integrated_globalmean = xr.DataArray(\
    data=G_vol_horiz_divergence_depth_integrated_globalmean_array,dims=['time'],coords={'time':vol_transport_x.time})

# # calculate time-mean using the month length weights
# G_vol_horiz_divergence_depth_integrated_mean = (G_vol_horiz_divergence_depth_integrated*month_length_weights).sum('time')
Completed calculation for time ind = [ 0 11], time elapsed = 50.8755362033844 s
Completed calculation for time ind = [12 23], time elapsed = 53.78032994270325 s
Completed calculation for time ind = [24 35], time elapsed = 51.95507740974426 s
Completed calculation for time ind = [36 47], time elapsed = 49.97524380683899 s
Completed calculation for time ind = [48 59], time elapsed = 48.08865427970886 s
Completed calculation for time ind = [60 71], time elapsed = 51.046666622161865 s
Completed calculation for time ind = [72 83], time elapsed = 50.7204864025116 s
Completed calculation for time ind = [84 95], time elapsed = 52.394755125045776 s
Completed calculation for time ind = [ 96 107], time elapsed = 48.3212685585022 s
Completed calculation for time ind = [108 119], time elapsed = 48.82131290435791 s
Completed calculation for time ind = [120 131], time elapsed = 51.62622809410095 s
Completed calculation for time ind = [132 143], time elapsed = 48.991485834121704 s
Completed calculation for time ind = [144 155], time elapsed = 48.832687854766846 s
Completed calculation for time ind = [156 167], time elapsed = 48.470677852630615 s
Completed calculation for time ind = [168 179], time elapsed = 49.91401410102844 s
Completed calculation for time ind = [180 191], time elapsed = 52.36308455467224 s
Completed calculation for time ind = [192 203], time elapsed = 50.70799994468689 s
Completed calculation for time ind = [204 215], time elapsed = 50.598023414611816 s
Completed calculation for time ind = [216 227], time elapsed = 48.314342975616455 s
Completed calculation for time ind = [228 239], time elapsed = 50.95944428443909 s
Completed calculation for time ind = [240 251], time elapsed = 51.90694284439087 s
Completed calculation for time ind = [252 263], time elapsed = 50.1651508808136 s
Completed calculation for time ind = [264 275], time elapsed = 48.01549291610718 s
Completed calculation for time ind = [276 287], time elapsed = 47.79500150680542 s

Plot the time-mean G_{\text{horizontal flux divergence}}

[50]:
plt.figure(figsize=(20,8));
obj = ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              G_vol_horiz_divergence_depth_integrated_mean,\
                              show_colorbar=True,\
                              cmin=-1e-7, cmax=1e-7, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);
plt.title('Average $\partial \eta / \partial t$ [m/s] implied by horizontal flux divergence\n Positive = Water out',
          fontsize=20);
_images/ECCO_v4_Volume_budget_closure_80_0.png

Plot one example G_{\text{horizontal flux divergence}}

[51]:
G_vol_horiz_divergence_depth_integrated.shape
[51]:
(288, 13, 90, 90)
[52]:
curr_t_ind = 100
curr_time = ecco_monthly_mean.time.values[curr_t_ind]
curr_timestr = str(curr_time)[0:4]+'-'+str(curr_time)[5:7]

plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              G_vol_horiz_divergence_depth_integrated.isel(time=curr_t_ind),
                              show_colorbar=True,\
                              cmin=-2e-7, cmax=2e-7, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);
plt.title('$\partial \eta / \partial t$ [m/s] implied by horizontal flux divergence for '\
          +curr_timestr+'\nPositive = Water out',
          fontsize=20);
_images/ECCO_v4_Volume_budget_closure_83_0.png

Plot G_{\text{horizontal flux divergence}} at 100m depth

[53]:
plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              G_vol_horiz_divergence.isel(k=10,time=curr_t_ind),
                              show_colorbar=True,\
                              cmin=-2e-6, cmax=2e-6, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);
plt.title('$\partial \eta / \partial t$ [m/s] at 100 m implied by horizontal flux divergence for '\
          +curr_timestr+'\nPositive = Water out',
          fontsize=20);
_images/ECCO_v4_Volume_budget_closure_85_0.png

Comparison of LHS and RHS

Time mean difference of LHR and RHS

Plot the time-mean difference between the LHS and RHS of the volume budget equation.

[54]:
#LHS ETA TENDENCY
a = G_total_tendency_mean
# RHS ETA TENDENCY FROM VOL DIVERGENCE AND SURFACE FLUXES
b = G_vol_horiz_divergence_depth_integrated_mean
c = G_surf_fluxes_mean

delta = a - b - c
[55]:
plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, delta, \
                              show_colorbar=True,\
                              cmin=-1e-9, cmax=1e-9, \
                              cmap='RdBu_r', \
                              dx=map_dx,dy=map_dy);
plt.title('Residual $\partial \eta / \partial t$ [m/s]: LHS - RHS ',
          fontsize=20);
_images/ECCO_v4_Volume_budget_closure_90_0.png

The residual of the time-mean tendency terms is essentially zero. We can look at the distribution of residuals to get a little more confidence.

Histogram of residuals

[56]:
tmp = np.abs( a-b-c).values.ravel();
plt.figure(figsize=(10,3));

plt.hist(tmp[np.nonzero(tmp > 0)],np.linspace(0, .5e-10,1000));
plt.grid()
_images/ECCO_v4_Volume_budget_closure_93_0.png

Almost all residuals < 10^{-11} m/s. We can close the ETAN budget using UVELMASS, VVELMASS, WVELMASS and oceFWflx.

Note: As stated earlier, we don’t actually need oceFWflx because the surface velocity of WVELMASS is proportional to oceFWflx but we kept it so that the surface forcing term is explicit.

ETAN budget closure through time

Global average ETAN budget closure

Another way of demonstrating volume budget closure is to show the global spatially-averaged ETAN tendency terms through time

[57]:
# LHR and RHS through time
a = G_total_tendency
b = G_vol_horiz_divergence_depth_integrated
c = G_surf_fluxes
# residuals
d = a-b-c

# area_masked = ecco_grid.rA.where(ecco_grid.hFacC.isel(k=0) > 0).compute()
# area_masked_sum = area_masked.sum(dim=('i','j','tile')).values

# take area-weighted mean of these terms
tmp_a=(a*area_masked).sum(dim=('i','j','tile')).compute()/area_masked_sum
# tmp_b=(b*area_masked).sum(dim=('i','j','tile')).compute()/area_masked_sum     # can compute global mean earlier when vol horiz div is computed
tmp_b=G_vol_horiz_divergence_depth_integrated_globalmean
tmp_c=(c*area_masked).sum(dim=('i','j','tile')).compute()/area_masked_sum
# tmp_d=(d*area_masked).sum(dim=('i','j','tile'))/area_masked_sum
tmp_d = tmp_a - tmp_b - tmp_c

# result is four time series
tmp_a.dims
[57]:
('time',)
[58]:
fig, axs = plt.subplots(2, 2, figsize=(12,8))

plt.sca(axs[0,0])
tmp_a.plot(color='orange')
axs[0,0].set_title("d(eta)/d(t) total [m/s]")
plt.grid()
plt.ylim([-1e-8, 1e-8]);

plt.sca(axs[0,1])
tmp_b.plot(color='g')
axs[0,1].set_title("d(eta)/d(t) vol divergence [m/s]")
plt.grid()
plt.ylim([-1e-8, 1e-8]);

plt.sca(axs[1,0])
tmp_c.plot(color='r')
axs[1,0].set_title("d(eta)/d(t) surf fluxes [m/s]")
plt.grid()
plt.ylim([-1e-8, 1e-8]);

plt.sca(axs[1,1])
tmp_d.plot(color='b')
axs[1,1].set_title("d(eta)/d(t) LHS - RHS [m/s]")
plt.grid()
plt.subplots_adjust(hspace = .5, wspace=.2)
plt.ylim([-1e-8, 1e-8])
plt.suptitle('Global Volume Budget',fontsize=20);
_images/ECCO_v4_Volume_budget_closure_96_0.png

When averaged over the entire ocean surface the volumetric divergence has no net impact on \partial \eta / \partial t. This makes sence because horizontal flux divergence can only redistribute volume. Globally, \eta can only change via net surface fluxes.

Local ETAN budget closure

Locally we expect that volume divergence can impact \eta. This is demonstrated for a single point the Arctic.

[59]:
# close distributed client (removes overhead of distributed client with numerous small chunks)
client.close()
[60]:
%%time

# Recall, from above...
a = G_total_tendency
b = G_vol_horiz_divergence_depth_integrated
c = G_surf_fluxes
d = a-b-c

plot_tile = 6
plot_j = 40
plot_i = 29
# tmp_aa = a.isel(tile=plot_tile,j=plot_j,i=plot_i)
# tmp_bb = b.isel(tile=plot_tile,j=plot_j,i=plot_i)
# tmp_cc = c.isel(tile=plot_tile,j=plot_j,i=plot_i)
# tmp_dd = d.isel(tile=plot_tile,j=plot_j,i=plot_i)

# use superchunks again (1 year = 1 superchunk) to partition computation
# and improve performance in limited-memory
tmp_aa = np.empty((a.sizes['time'],))
tmp_bb = np.empty((b.sizes['time'],))
tmp_cc = np.empty((c.sizes['time'],))
size_superchunk = 12
for time_superchunk in range(int(np.ceil(b.sizes['time']/size_superchunk))):
    curr_time_ind = np.arange(size_superchunk*time_superchunk,\
                              np.fmin(size_superchunk*(time_superchunk+1),\
                                      b.sizes['time']))
    tmp_aa[curr_time_ind] = G_total_tendency.isel(time=curr_time_ind,\
                                                  tile=plot_tile,j=plot_j,i=plot_i).compute()
    tmp_bb[curr_time_ind] = G_vol_horiz_divergence_depth_integrated.isel(time=curr_time_ind,\
                                                                         tile=plot_tile,j=plot_j,i=plot_i).compute()
    tmp_cc[curr_time_ind] = G_surf_fluxes.isel(time=curr_time_ind,\
                                               tile=plot_tile,j=plot_j,i=plot_i).compute()
    print(f'Completed time_superchunk = {time_superchunk}')
tmp_dd = tmp_aa-tmp_bb-tmp_cc

tmp_aa = xr.DataArray(data=tmp_aa,dims=['time'],coords=a.isel(tile=plot_tile,j=plot_j,i=plot_i).coords)
tmp_bb = xr.DataArray(data=tmp_bb,dims=['time'],coords=b.isel(tile=plot_tile,j=plot_j,i=plot_i).coords)
tmp_cc = xr.DataArray(data=tmp_cc,dims=['time'],coords=c.isel(tile=plot_tile,j=plot_j,i=plot_i).coords)
tmp_dd = xr.DataArray(data=tmp_dd,dims=['time'],coords=d.isel(tile=plot_tile,j=plot_j,i=plot_i).coords)

fig, axs = plt.subplots(2, 2, figsize=(12,8))
plt.sca(axs[0,0])
tmp_aa.plot(color='orange')
axs[0,0].set_title("d(eta)/d(t) total [m/s]")
plt.ticklabel_format(axis='y', style='sci', useMathText=True)
plt.ylim([-5e-7, 5e-7]);
plt.grid()

plt.sca(axs[0,1])
tmp_bb.plot(color='g')
axs[0,1].set_title("d(eta)/d(t) vol divergence [m/s]")
plt.ticklabel_format(axis='y', style='sci', useMathText=True)
plt.ylim([-5e-7, 5e-7]);
plt.grid()

plt.sca(axs[1,0])
tmp_cc.plot(color='r')
axs[1,0].set_title("d(eta)/d(t) surf fluxes [m/s]")
plt.ticklabel_format(axis='y', style='sci', useMathText=True)
plt.ylim([-5e-7, 5e-7]);
plt.grid()

plt.sca(axs[1,1])
tmp_dd.plot(color='b')
axs[1,1].set_title("d(eta)/d(t) LHS - RHS [m/s]")
plt.ticklabel_format(axis='y', style='sci', useMathText=True)
plt.ylim([-5e-7, 5e-7]);
plt.grid()

plt.subplots_adjust(hspace = .5, wspace=.2)
plt.suptitle('Volume Budget for one Arctic Ocean point',
                 fontsize=20);
Completed time_superchunk = 0
Completed time_superchunk = 1
Completed time_superchunk = 2
Completed time_superchunk = 3
Completed time_superchunk = 4
Completed time_superchunk = 5
Completed time_superchunk = 6
Completed time_superchunk = 7
Completed time_superchunk = 8
Completed time_superchunk = 9
Completed time_superchunk = 10
Completed time_superchunk = 11
Completed time_superchunk = 12
Completed time_superchunk = 13
Completed time_superchunk = 14
Completed time_superchunk = 15
Completed time_superchunk = 16
Completed time_superchunk = 17
Completed time_superchunk = 18
Completed time_superchunk = 19
Completed time_superchunk = 20
Completed time_superchunk = 21
Completed time_superchunk = 22
Completed time_superchunk = 23
CPU times: user 5min 42s, sys: 33.7 s, total: 6min 15s
Wall time: 5min 36s
_images/ECCO_v4_Volume_budget_closure_99_1.png

Indeed, the volume divergence term does contribute to \eta variations at this one point.

The seasonal cycle of surface fluxes from sea-ice growth and melt can be seen in the surface fluxes term (plotted below for just the year 1995)

[61]:
tmp_cc.sel(time='1995').plot(color='r')
plt.ticklabel_format(axis='y', style='sci', useMathText=True)
plt.ylim([-5e-7, 5e-7]);
plt.grid()
plt.title('d(eta)/d(t) from surf. fluxes in 1995 [m/s]');
_images/ECCO_v4_Volume_budget_closure_101_0.png

Predicted vs. actual \eta

As we have shown, in our Boussinesq model the only term that can change global mean model sea level anomaly \eta is net surface freshwater flux. Let us compare the time-evolution of \eta implied by surface freshwater fluxes and the actual \eta from the model output.

The predicted \eta time series is calculated by time integrating :math:`G_{surface fluxes}`.
This time series is compared against the actual \eta time series anomaly relative to the \eta(t=0).
[62]:
area_masked = ecco_grid.rA.where(ecco_grid.maskC.isel(k=0) == 1)

dETA_per_month_predicted_from_surf_fluxes = \
    ((G_surf_fluxes * area_masked).sum(dim=('i','j','tile')) / \
     area_masked.sum())*secs_per_month

ETA_predicted_by_surf_fluxes = \
    np.cumsum(dETA_per_month_predicted_from_surf_fluxes.values)

ETA_from_ETAN = \
    (ecco_monthly_snaps.ETAN * area_masked).sum(dim=('i','j','tile')) / \
    area_masked.sum().compute()

# plotting
plt.figure(figsize=(14,5));

plt.plot(dETA_per_month_predicted_from_surf_fluxes.time, \
         ETA_predicted_by_surf_fluxes,'b.')
plt.plot(ETA_from_ETAN.time.values, ETA_from_ETAN-ETA_from_ETAN.values[0],'r-')
plt.grid()
plt.ylabel('global mean $\eta$');
plt.legend(('predicted', 'actual'));
plt.title('$\eta(t)$ as predicted from net surface fluxes and model ETAN [m]',
          fontsize=20);
_images/ECCO_v4_Volume_budget_closure_103_0.png

The first predicted \eta occurs at the end of the first month (one month time integral of \partial \eta / \partial t). The first actual \eta is set to be zero.

The above plot is another way of confirming that in our Boussinesq model the only term that can change global mean model sea level anomaly \eta is net surface freshwater flux. To account for changes in global mean density we must apply the Greatbatch correction, inverse-barometer correction, and a correction term to account for the fact that sea-ice does not ‘float’ on top of the ocean but in fact displaces seawater upwards. All of these corrections are made for the term SSH, dynamic sea surface height anomaly (not shown here).

One can compare the sea level rise from mass fluxes in ECCO vs those estimated from GRACE and other datasets published by the WCRP Global Sea Level Budget Group: Global sea-level budget 1993-present, Earth Syst. Sci. Data, 10, 1551-1590, https://doi.org/10.5194/essd-10-1551-2018, 2018.

Available here. See Figure 16. https://www.earth-syst-sci-data.net/10/1551/2018/essd-10-1551-2018.pdf

[63]:
# Annual mean SL calculation must account for different lengths of each month.

# step 1. weight predicted eta by seconds in each month
tmp1=ETA_predicted_by_surf_fluxes*secs_per_month
# step 2, group records by year and sum across each year
tmp2=tmp1.groupby('time.year').sum()
# step 3, group secs per month by year and sum across each year
secs_per_year = secs_per_month.groupby('time.year').sum()

# step 3, divide time-weighted ETA by seconds per year
annual_mean_GMSL_due_to_mass_fluxes =tmp2/secs_per_year
num_years = len(annual_mean_GMSL_due_to_mass_fluxes.year.values)

plt.figure(figsize=(8,5));
# the -0.13 is to make the starting value comparable with WCRP fig 16.
plt.bar(annual_mean_GMSL_due_to_mass_fluxes.year.values[12:num_years],\
        1000*(annual_mean_GMSL_due_to_mass_fluxes.values[12:num_years]-.017), \
        color='k')
plt.grid()
plt.xticks(np.arange(2005,
                     annual_mean_GMSL_due_to_mass_fluxes.year.values[-1]+1,step=1));
plt.title('Sea level (mm) caused by mass fluxes');
_images/ECCO_v4_Volume_budget_closure_106_0.png

vs the WCRP data:

[64]:
from IPython.display import Image
Image('../figures/cazenave_ocean_mass.png')
[64]:
_images/ECCO_v4_Volume_budget_closure_108_0.png

Time-mean \partial \eta / \partial t due to net freshwater fluxes

We can easily calculate the time mean rate of global sea level rise due to net freshwater flux.

[65]:
x=(G_surf_fluxes_mean * area_masked).sum() / area_masked.sum()
total_GMSLR_due_to_mass_fluxes = x*seconds_in_entire_period
print('total global mean sea level rise due to mass fluxes m: ', \
      total_GMSLR_due_to_mass_fluxes.values)
total global mean sea level rise due to mass fluxes m:  0.06982078709392495

Dividing this by the total number of years in the analysis period gives a average rate per year:

[66]:
total_number_of_years = len(secs_per_month)/12
total_number_of_years
[66]:
24.0
[67]:
mean_rate_of_GMSLR_due_to_mass_fluxes = \
    total_GMSLR_due_to_mass_fluxes/total_number_of_years

print('mean rate of GMSLR due to mass fluxes [mm/yr] ',  \
      1000*mean_rate_of_GMSLR_due_to_mass_fluxes.values)
mean rate of GMSLR due to mass fluxes [mm/yr]  2.9091994622468724

Compare with other estimates of GMSLR due to mass fluxes.