Compute MOC along the approximate OSNAP array from ECCO¶

Here we compute volumetric transport along the OSNAP lines in depth space, which can be compared to recent observations.

Here we show:

how to get masks denoting the great circle arc between two points in space
how to compute the transport or streamfunction across this section
a comparison to observations

Datasets¶

These are the ShortNames of the ECCO datasets needed to run this tutorial:

ECCO_L4_GEOMETRY_LLC0090GRID_V4R4
ECCO_L4_OCEAN_3D_VOLUME_FLUX_LLC0090GRID_MONTHLY_V4R4 (1992-2017)
ECCO_L4_OCEAN_3D_TEMPERATURE_FLUX_LLC0090GRID_MONTHLY_V4R4 (1992-2017)

Make sure you have the ecco_access library in your Python path before you run this tutorial; it will handle the access to the datasets (either download or in-cloud direct access from S3, depending on the access mode you specify).

[1]:

import warnings
warnings.filterwarnings('ignore')

[2]:

import os
import sys
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import glob
import copy
import cartopy as cart
import cartopy.crs as ccrs
notebook_path = os.getcwd()

import ecco_access as ea

# indicate mode of access
# options are:
# 'download': direct download from internet to your local machine
# 'download_ifspace': like download, but only proceeds
#                     if your machine have sufficient storage
# 's3_open': access datasets in-cloud from an AWS instance
# 's3_open_fsspec': use jsons generated with fsspec and
#                   kerchunk libraries to speed up in-cloud access
# 's3_get': direct download from S3 in-cloud to an AWS instance
# 's3_get_ifspace': like s3_get, but only proceeds if your instance
#                   has sufficient storage
access_mode = 'download_ifspace'

[3]:

# setting up a dask LocalCluster
from dask.distributed import Client
from dask.distributed import LocalCluster
cluster = LocalCluster()
client = Client(cluster)
client

[3]:

Client

Client-f2fb873a-8cd3-11ef-b071-ea833bf36b53

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

Cluster Info

LocalCluster

b4a35307

Dashboard: http://127.0.0.1:36525/status	Workers: 4
Total threads: 4	Total memory: 15.27 GiB
Status: running	Using processes: True

Scheduler Info

Scheduler

Scheduler-13421e3f-b34f-41ee-98d2-9f9abefe87cd

Comm: tcp://127.0.0.1:45601	Workers: 4
Dashboard: http://127.0.0.1:36525/status	Total threads: 4
Started: Just now	Total memory: 15.27 GiB

Workers

Worker: 0

Comm: tcp://127.0.0.1:34881	Total threads: 1
Dashboard: http://127.0.0.1:35355/status	Memory: 3.82 GiB
Nanny: tcp://127.0.0.1:38397
Local directory: /tmp/dask-scratch-space/worker-qefry5p7

Worker: 1

Comm: tcp://127.0.0.1:35259	Total threads: 1
Dashboard: http://127.0.0.1:46137/status	Memory: 3.82 GiB
Nanny: tcp://127.0.0.1:35563
Local directory: /tmp/dask-scratch-space/worker-f075goyn

Worker: 2

Comm: tcp://127.0.0.1:40867	Total threads: 1
Dashboard: http://127.0.0.1:39767/status	Memory: 3.82 GiB
Nanny: tcp://127.0.0.1:44229
Local directory: /tmp/dask-scratch-space/worker-1gdae72j

Worker: 3

Comm: tcp://127.0.0.1:41935	Total threads: 1
Dashboard: http://127.0.0.1:33797/status	Memory: 3.82 GiB
Nanny: tcp://127.0.0.1:36595
Local directory: /tmp/dask-scratch-space/worker-9kr4ntby

[4]:

## 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

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

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

import ecco_v4_py as ecco

[5]:

## Set top-level file directory for the ECCO NetCDF files
## =================================================================

## currently set to ~/Downloads/ECCO_V4r4_PODAAC
ECCO_dir = join(user_home_dir,'ECCO_V4r4_PODAAC')

# # for access_mode = 's3_open_fsspec', need to specify the root directory
# # containing the jsons
# jsons_root_dir = join('/efs_ecco','mzz-jsons')

[6]:

## Access datasets needed for this tutorial

ShortNames_list = ["ECCO_L4_GEOMETRY_LLC0090GRID_V4R4",\
                   "ECCO_L4_OCEAN_3D_VOLUME_FLUX_LLC0090GRID_MONTHLY_V4R4",\
                   "ECCO_L4_OCEAN_3D_TEMPERATURE_FLUX_LLC0090GRID_MONTHLY_V4R4"]
StartDate = '1992-01'
EndDate = '2017-12'
ds_dict = ea.ecco_podaac_to_xrdataset(ShortNames_list,\
                                            StartDate=StartDate,EndDate=EndDate,\
                                            mode=access_mode,\
                                            download_root_dir=ECCO_dir,\
                                            max_avail_frac=0.5)

[7]:

## Load the model grid
ecco_grid = ds_dict[ShortNames_list[0]]

## Load ECCO variables
ecco_vars_uvw = ds_dict[ShortNames_list[1]]
ecco_vars_adv_th = ds_dict[ShortNames_list[2]]

ecco_vars = xr.merge((ecco_vars_adv_th[['ADVx_TH','ADVy_TH']],ecco_vars_uvw[['UVELMASS','VVELMASS']]))

## Merge the ecco_grid with the ecco_vars to make the ecco_ds
ds = xr.merge((ecco_grid , ecco_vars))

Load O-SNAP data¶

These data are available at https://www.o-snap.org/observations/data/ and are published in

Lozier, M. S., Li, F., Bacon, S., Bahr, F., Bower, A. S., Cunningham, S. A., … Zhao, J. (2019). A sea change in our view of overturning in the subpolar North Atlantic. Science, 363(6426), 516–521. https://doi.org/10.1126/science.aau6592

Download these files from the O-SNAP data page:

OSNAP_MOC_MHT_MFT_TimeSeries_201408_202006_2023.nc : https://repository.gatech.edu/bitstreams/e039e311-dd2e-4511-a525-c2fcfb3be85a/download OSNAP_Streamfunction_201408_202006_2023.nc : https://repository.gatech.edu/bitstreams/5edf4cba-a28f-40a6-a4da-24d7436a42ab/download

and put them in a path:

[8]:

# change this path as needed based on where you have stored the O-SNAP data files above
osnap_data_dir = join(user_home_dir,'Downloads','OSNAP')

[9]:

obs1 = xr.open_dataset(join(osnap_data_dir,'OSNAP_MOC_MHT_MFT_TimeSeries_201408_202006_2023.nc'))
obs2 = xr.open_dataset(join(osnap_data_dir,'OSNAP_Streamfunction_201408_202006_2023.nc'))
obs = xr.merge((obs1,obs2))
print(obs)

<xarray.Dataset> Size: 834kB
Dimensions:       (TIME: 71, LEVEL: 481)
Coordinates:
  * TIME          (TIME) datetime64[ns] 568B 2014-08-01T12:00:00 ... 2020-06-...
  * LEVEL         (LEVEL) float64 4kB 23.3 23.31 23.32 ... 28.08 28.09 28.1
Data variables: (12/21)
    MOC_ALL       (TIME) float64 568B ...
    MOC_ALL_ERR   (TIME) float64 568B ...
    MOC_EAST      (TIME) float64 568B ...
    MOC_EAST_ERR  (TIME) float64 568B ...
    MOC_WEST      (TIME) float64 568B ...
    MOC_WEST_ERR  (TIME) float64 568B ...
    ...            ...
    MFT_EAST_ERR  (TIME) float64 568B ...
    MFT_WEST      (TIME) float64 568B ...
    MFT_WEST_ERR  (TIME) float64 568B ...
    T_ALL         (LEVEL, TIME) float64 273kB ...
    T_EAST        (LEVEL, TIME) float64 273kB ...
    T_WEST        (LEVEL, TIME) float64 273kB ...
Attributes: (12/13)
    title:                    OSNAP MOC MHT MFT time series (2014-2020)
    project:                  OSNAP
    contributor_name:         Yao Fu, M. Susan Lozier, Tiago Carrilho Biló, A...
    contributor_role:         data design, collection and/or processing
    contributor_institution:  Georgia Institute of Technology, USA; National ...
    publisher_name:           M. Susan Lozier; Yao Fu
    ...                       ...
    data_assembly_center:     Georgia Institute of Technology
    references:               http://www.o-snap.org; http://www.o-snap.org/data/
    citation:                 OSNAP data were collected and made freely avail...
    date_created:             2023-04-24 11:33:07
    time_coverage_start:      2014-08-1 00:00:00
    time_coverage_end:        2020-06-30 00:00:00

Define the OSNAP lines¶

We define the OSNAP lines roughly by point pairs in [longitude, latitude] space. The lines are then computed via the function ecco_v4_py.get_section_line_masks as the great circle arc between these two points.

See below for similar MATLAB functions

[10]:

pt1_east = [-44, 60]
pt2_east = [-5, 56]

pt1_west = [-56, 51]
pt2_west = [-45, 60]

[11]:

maskC_east, maskW_east, maskS_east = ecco.get_section_line_masks(pt1_east,pt2_east,ds)
maskC_west, maskW_west, maskS_west = ecco.get_section_line_masks(pt1_west,pt2_west,ds)
maskC_tot = (maskC_east+maskC_west).where(maskC_east+maskC_west==1,0)
maskW_tot = (maskW_east+maskW_west).where(np.abs(maskW_east)+np.abs(maskW_west)==1,0)
maskS_tot = (maskS_east+maskS_west).where(np.abs(maskS_east)+np.abs(maskS_west)==1,0)

[12]:

plt.figure(figsize=(12,6))
ecco.plot_proj_to_latlon_grid(ds.XC,ds.YC,maskC_tot,cmap='viridis',projection_type='robin',user_lon_0=0);

[13]:

plt.figure(figsize=(8,8))
# use dx=.1, dy=.1 so that plot shows the osnap array as a thin
# line.  remember, plot_proj_to_latlon_grid first interpolates
# the model grid to lat-lon with grid spacing as dx, dy
P = ecco.plot_proj_to_latlon_grid(ds.XC, ds.YC, maskC_tot, \
                                cmap='viridis',\
                                projection_type='PlateCarree',\
                                lat_lim=45,dx=.1,dy=.1)
#ax.add_feature(cart.feature.LAND,facecolor='0.7',zorder=2)
P[1].set_extent([-80, 20, 40, 80], crs=ccrs.PlateCarree())
plt.show()

We have defined many commonly used sections in oceanography so that users can access these easily. The available sections are shown below. This allows one to, e.g. compute volumetric transport across the Drake Passage as follows:

drake_vol = ecco_v4_py.calc_section_vol_trsp(ds,section_name='Drake Passage')

Similarly, we can do the same with calc_section_heat_trsp and calc_section_salt_trsp.

One can also get these pre-defined section masks as follows:

pt1,pt2 = ecco_v4_py.get_section_endpoints('Drake Passage')
maskC, maskW, maskS = ecco_v4_py.get_section_line_masks(ds,pt1,pt2)

Finally, one can see similar functions in MATLAB:

define general section masks: gcmfaces_calc/gcmfaces_lines_transp.m
see pre-defined section endpoints: gcmfaces_calc/gcmfaces_lines_pairs.m

[14]:

ecco.get_available_sections()

[14]:

['Bering Strait',
 'Gibraltar',
 'Florida Strait',
 'Florida Strait W1',
 'Florida Strait S1',
 'Florida Strait E1',
 'Florida Strait E2',
 'Florida Strait E3',
 'Florida Strait E4',
 'Davis Strait',
 'Denmark Strait',
 'Iceland Faroe',
 'Faroe Scotland',
 'Scotland Norway',
 'Fram Strait',
 'Barents Sea',
 'Labrador Greenland',
 'Hudson Strait',
 'English Channel',
 'Newfoundland Iberia',
 'Drake Passage',
 'Indonesia W1',
 'Indonesia W2',
 'Indonesia W3',
 'Indonesia W4',
 'Australia Antarctica',
 'Madagascar Channel',
 'Madagascar Antarctica',
 'South Africa Antarctica']

Compute the overturning streamfunction in depth space¶

The function calc_section_stf computes the overturning streamfunction across the plane normal to the section denoted by the west and south masks. It is also possible to compute the overturning streamfunction at a particular latitude band, as is done to compare to the RAPID array, for instance. Please see the function calc_meridional_stf to do this, which is also in ecco_v4_py.calc_stf.

Note that we can also compute the volumetric, heat, or salt transport across these sections using the first three functions defined in ecco_v4_py.calc_section_trsp.

In MATLAB, we can compute meridional overturning streamfunctions with gcmfaces_calc/calc_overturn.m. Section transports can be computed with gcmfaces_calc/calc_transports.m once the masks are defined.

[15]:

%%time

# compute streamfunctions, looping through individual time coordinates (this works more consistently in limited-memory environments)
time_chunksize = 12
for time_chunk in range(int(np.ceil(ds.sizes['time']/time_chunksize))):
    curr_t_isel = np.arange(time_chunk*time_chunksize,np.fmin((time_chunk+1)*time_chunksize,ds.sizes['time']))
    curr_osnap_z_stf_east = ecco.calc_section_stf(ds.isel(time=curr_t_isel),\
                                                  pt1=pt1_east, \
                                                  pt2=pt2_east,\
                                                  section_name='OSNAP East Overturning Streamfunction').compute()

    curr_osnap_z_stf_west = ecco.calc_section_stf(ds.isel(time=curr_t_isel), \
                                                  pt1=pt1_west, \
                                                  pt2=pt2_west,\
                                                  section_name='OSNAP West Overturning Streamfunction').compute()

    curr_osnap_z_stf_tot = ecco.calc_section_stf(ds.isel(time=curr_t_isel),\
                                                 maskW=maskW_tot, \
                                                 maskS=maskS_tot,\
                                                 section_name='OSNAP Total Overturning Streamfunction').compute()

    if time_chunk == 0:
        osnap_z_stf_east = copy.deepcopy(curr_osnap_z_stf_east)
        osnap_z_stf_west = copy.deepcopy(curr_osnap_z_stf_west)
        osnap_z_stf_tot = copy.deepcopy(curr_osnap_z_stf_tot)
    else:
        osnap_z_stf_east = xr.concat((osnap_z_stf_east,curr_osnap_z_stf_east),dim='time',\
                                     compat='override',data_vars='minimal',coords='minimal')
        osnap_z_stf_west = xr.concat((osnap_z_stf_west,curr_osnap_z_stf_west),dim='time',\
                                     compat='override',data_vars='minimal',coords='minimal')
        osnap_z_stf_tot = xr.concat((osnap_z_stf_tot,curr_osnap_z_stf_tot),dim='time',\
                                     compat='override',data_vars='minimal',coords='minimal')

    print('Processed time indices',curr_t_isel[0],'-',curr_t_isel[-1],'( out of',ds.sizes['time'],')')

Processed time indices 0 - 11 ( out of 312 )
Processed time indices 12 - 23 ( out of 312 )
Processed time indices 24 - 35 ( out of 312 )
Processed time indices 36 - 47 ( out of 312 )
Processed time indices 48 - 59 ( out of 312 )
Processed time indices 60 - 71 ( out of 312 )
Processed time indices 72 - 83 ( out of 312 )
Processed time indices 84 - 95 ( out of 312 )
Processed time indices 96 - 107 ( out of 312 )
Processed time indices 108 - 119 ( out of 312 )
Processed time indices 120 - 131 ( out of 312 )
Processed time indices 132 - 143 ( out of 312 )
Processed time indices 144 - 155 ( out of 312 )
Processed time indices 156 - 167 ( out of 312 )
Processed time indices 168 - 179 ( out of 312 )
Processed time indices 180 - 191 ( out of 312 )
Processed time indices 192 - 203 ( out of 312 )
Processed time indices 204 - 215 ( out of 312 )
Processed time indices 216 - 227 ( out of 312 )
Processed time indices 228 - 239 ( out of 312 )
Processed time indices 240 - 251 ( out of 312 )
Processed time indices 252 - 263 ( out of 312 )
Processed time indices 264 - 275 ( out of 312 )
Processed time indices 276 - 287 ( out of 312 )
Processed time indices 288 - 299 ( out of 312 )
Processed time indices 300 - 311 ( out of 312 )
CPU times: user 53.7 s, sys: 4.29 s, total: 58 s
Wall time: 3min 50s

[16]:

def osnap_depth_stf_vs_time(stf_ds,label):
    fig = plt.figure(figsize=(18,6))

    # Time evolving
    plt.subplot(1,4,(1,3))
    time_edge_extrap = np.hstack((stf_ds['time'].values[0] - (0.5*np.diff(stf_ds['time'].values[0:2])),\
                                  stf_ds['time'].values[:-1] + (0.5*np.diff(stf_ds['time'].values)),\
                                  stf_ds['time'].values[-1] + (0.5*np.diff(stf_ds['time'].values[-2:]))))
    Z_edge_extrap = np.hstack((np.array([0]),\
                               stf_ds['Z'].values[:-1] + (0.5*np.diff(stf_ds['Z'].values)),\
                               np.array([-6134.5])))
    plt.pcolormesh(time_edge_extrap,Z_edge_extrap,stf_ds['psi_moc'].T)
    plt.title('ECCOv4r4\nOverturning streamfunction across OSNAP %s [Sv]' % label)
    plt.ylabel('Depth [m]')
    plt.xlabel('Month')
    plt.xticks(rotation=45)
    cb = plt.colorbar()
    cb.set_label('[Sv]')

    plt.subplot(1,4,4)
    plt.plot(stf_ds['psi_moc'].mean('time'),stf_ds['Z'])
    plt.title('ECCOv4r4\nTime mean streamfunction, OSNAP %s' % label)
    plt.ylabel('Depth [m]')
    plt.xlabel('[Sv]')
    plt.grid()
    plt.show()

[17]:

osnap_depth_stf_vs_time(osnap_z_stf_east,'East')
osnap_depth_stf_vs_time(osnap_z_stf_west,'West')
osnap_depth_stf_vs_time(osnap_z_stf_tot,'Total')

[18]:

osnap_z_ov_east = osnap_z_stf_east['psi_moc'].max(dim='k')
osnap_z_ov_west = osnap_z_stf_west['psi_moc'].max(dim='k')
osnap_z_ov_tot = osnap_z_stf_tot['psi_moc'].max(dim='k')

[19]:

fig = plt.figure(figsize=(12,6))
plt.plot(osnap_z_ov_east['time'],osnap_z_ov_east)
plt.plot(osnap_z_ov_west['time'],osnap_z_ov_west)
plt.plot(osnap_z_ov_tot['time'],osnap_z_ov_tot)
plt.title('ECCOv4r4\nOverturning strength across OSNAP lines [Sv]')
plt.ylabel('[Sv]')
plt.legend((('East, avg. = %.2f $\pm$ %.2f [Sv]' %
             (osnap_z_ov_east.mean('time').values,osnap_z_ov_east.std('time').values)),
            ('West, avg. = %.2f $\pm$ %.2f [Sv]' %
             (osnap_z_ov_west.mean('time').values,osnap_z_ov_west.std('time').values)),
            ('Total, avg. = %.2f $\pm$ %.2f [Sv]' %
             (osnap_z_ov_tot.mean('time').values,osnap_z_ov_tot.std('time').values))))
plt.grid()
plt.show()

Compare to Observations¶

[20]:

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

var_list = ['MOC_EAST','MOC_WEST','MOC_ALL']
err_list = ['MOC_EAST_ERR','MOC_WEST_ERR','MOC_ALL_ERR']
obs_handles = []
for var, err in zip(var_list,err_list):
    if 'MOC_ALL' in var:
        clr = 'k'
    elif 'MOC_EAST' in var:
        clr = 'r'
    elif 'MOC_WEST' in var:
        clr = 'b'


    curr_hand = plt.plot(obs['TIME'],obs[var], color=clr)
    obs_handles.append(curr_hand[0])

    plt.fill_between(obs['TIME'].values,
                 obs[var]-obs[err],
                 obs[var]+obs[err],
                 color=[0.8,0.8,0.8])

# plot ECCO v4r4 equivalent
ecco_hand0 = plt.plot(osnap_z_ov_east['time'],osnap_z_ov_east, 'r--')
ecco_hand1 = plt.plot(osnap_z_ov_west['time'],osnap_z_ov_west, 'b--',)
ecco_hand2 = plt.plot(osnap_z_ov_tot['time'],osnap_z_ov_tot, 'k--')
ecco_handles = [ecco_hand0[0],ecco_hand1[0],ecco_hand2[0]]
plt.xlim(np.array(['2014-08','2018-01']).astype('datetime64[ns]'))

plt.ylabel('[%s]' % obs['MOC_ALL'].units)
plt.grid()
plt.title('OSNAP MOC Derived from observations [%s]' % obs['MOC_ALL'].units)
plt.legend(handles=(obs_handles + ecco_handles),labels=(\
    ('East MOC, (Lozier et al, 2019)\nTime mean: %.2f $\pm$ %.2f %s' %
            (obs['MOC_EAST'].mean('TIME'),obs['MOC_EAST'].std('TIME'),obs['MOC_EAST'].units)),
    ('West MOC, (Lozier et al, 2019)\nTime mean: %.2f $\pm$ %.2f %s' %
            (obs['MOC_WEST'].mean('TIME'),obs['MOC_WEST'].std('TIME'),obs['MOC_WEST'].units)),
    ('Total MOC, (Lozier et al, 2019)\nTime mean: %.2f $\pm$ %.2f %s' %
            (obs['MOC_ALL'].mean('TIME'),obs['MOC_ALL'].std('TIME'),obs['MOC_ALL'].units)),
    ('East MOC, ECCOv4r4\nTime mean: %.2f $\pm$ %.2f %s' %
            (osnap_z_ov_east.mean('time'),osnap_z_ov_east.std('time'), \
             osnap_z_stf_east['psi_moc'].attrs['units'])),
    ('West MOC, ECCOv4r4\nTime mean: %.2f $\pm$ %.2f %s' %
            (osnap_z_ov_west.mean('time'),osnap_z_ov_west.std('time'), \
             osnap_z_stf_west['psi_moc'].attrs['units'])),
    ('Total MOC, ECCOv4r4\nTime mean: %.2f $\pm$ %.2f %s' %
            (osnap_z_ov_tot.mean('time'),osnap_z_ov_tot.std('time'), \
             osnap_z_stf_tot['psi_moc'].attrs['units']))),
    loc='center left', bbox_to_anchor=(1, 0.5))

plt.show()