Compute meridional heat transport

This notebook shows how to compute meridional heat transport (MHT) across any particular latitude band. Additionally, we show this for both global and basin specific cases.

Oceanographic computations:

  • use xgcm to compute masks and grab values at a particular latitude band
  • use ecco_v4_py to select a specific basin
  • compute meridional heat transport at one or more latitude bands

Python basics on display:

Note that each of these tasks can be accomplished more succinctly with ecco_v4_py functions, but are shown explicitly to illustrate these tools. Throughout, we will note the ecco_v4_py (python) and gcmfaces (MATLAB) functions which can perform these computations.

[1]:
import warnings
warnings.filterwarnings('ignore')
[2]:
import os
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import cartopy as cart
import sys

Load Model Variables

Because we’re computing transport, we want the files containing ‘UVELMASS’ and ‘VVELMASS’ for volumetric transport, and ‘ADVx_TH’, ‘ADVy_TH’ and ‘DFxE_TH’, ‘DFyE_TH’ for the advective and diffusive components of heat transport, respectively.

[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.  For example, if your ecco_v4_py
##    files are in /Users/ifenty/ECCOv4-py/ecco_v4_py, then use:

sys.path.append('/home/ifenty/ECCOv4-py')
import ecco_v4_py as ecco
[4]:
## Set top-level file directory for the ECCO NetCDF files
## =================================================================
# base_dir = '/home/username/'
base_dir = '/home/ifenty/ECCOv4-release'

## define a high-level directory for ECCO fields
ECCO_dir = base_dir + '/Release3_alt'
[5]:
## Load the model grid
grid_dir= ECCO_dir + '/nctiles_grid/'

ecco_grid = ecco.load_ecco_grid_nc(grid_dir, 'ECCOv4r3_grid.nc')

## Load one year of 2D daily data, SSH, SST, and SSS
data_dir= ECCO_dir + '/nctiles_monthly'

ecco_vars = ecco.recursive_load_ecco_var_from_years_nc(data_dir, \
                                           vars_to_load=['ADVx_TH', 'ADVy_TH', \
                                                        'DFxE_TH', 'DFyE_TH'],\
                                           years_to_load = 'all')
## Merge the ecco_grid with the ecco_vars to make the ecco_ds
ecco_ds = xr.merge((ecco_grid , ecco_vars))
loading files of  ADVx_TH
loading files of  ADVy_TH
loading files of  DFxE_TH
loading files of  DFyE_TH

Grab latitude band: 26^\circN array as an example

Here we want to grab the transport values which along the band closest represented in the model to 26^\circN. In a latitude longitude grid this could simply be done by, e.g. U.sel(lat=26). However, the LLC grid is slightly more complicated. Luckily, the functionality enabled by the xgcm Grid object makes this relatively easy.

Note that this subsection can be performed with the with the ecco_v4_py modules vector_calc and scalar_calc as follows:

from ecco_v4_py import vector_calc, scalar_calc
grid = ecco_v4_py.get_llc_grid(ds)
rapid_maskW, rapid_maskS = vector_calc.get_latitude_masks(lat_val=26,yc=ds.YC,grid=grid)
rapid_maskC = scalar_calc.get_latitude_mask(lat_val=26,yc=ds.YC,grid=grid)

One can also use the gcmfaces function gcmfaces_calc/gcmfaces_lines_zonal.m.

[6]:
# Get array of 1's at and north of latitude
lat = 26
ones = xr.ones_like(ecco_ds.YC)
dome_maskC = ones.where(ecco_ds.YC>=lat,0)
[7]:
plt.figure(figsize=(12,6))
ecco.plot_proj_to_latlon_grid(ecco_ds.XC,ecco_ds.YC,dome_maskC,
                              projection_type='robin',cmap='viridis',user_lon_0=0,show_colorbar=True);
_images/ECCO_v4_Example_MHT_11_0.png

Again, if this were a lat/lon grid we could simply take a finite difference in the meridional direction. The only grid cells with 1’s remaining would be at the southern edge of grid cells at approximately 26^\circN.

However, recall that the LLC grid has a different picture.

[8]:
maskC = ecco_ds.maskC
maskS = ecco_ds.maskS
maskW = ecco_ds.maskW
[9]:
plt.figure(figsize=(12,6))
ecco.plot_tiles(dome_maskC+maskC.isel(k=0), cmap='viridis');
<Figure size 864x432 with 0 Axes>
_images/ECCO_v4_Example_MHT_14_1.png

Recall that for tiles 7-12, the y-dimension actually runs East-West. Therefore, we want to compute a finite difference in the x-dimension on these tiles to get the latitude band at 26^\circN. For tiles 1-5, we clearly want to difference in the y-dimension. Things get more complicated on tile 6.

Here we make the xgcm Grid object which allows us to compute finite differences in simple one liners. This object understands how each of the tiles on the LLC grid connect, because we provide that information. To see under the hood, checkout the utility function get_llc_grid where these connections are defined.

[10]:
grid = ecco.get_llc_grid(ecco_ds)
[11]:
lat_maskW = grid.diff(dome_maskC,'X',boundary='fill')
lat_maskS = grid.diff(dome_maskC,'Y',boundary='fill')
[12]:
plt.figure(figsize=(12,6))
ecco.plot_tiles(lat_maskW+maskW.isel(k=0), cmap='viridis');
<Figure size 864x432 with 0 Axes>
_images/ECCO_v4_Example_MHT_19_1.png
[13]:
plt.figure(figsize=(12,6))
ecco.plot_tiles(lat_maskS+maskS.isel(k=0), cmap='viridis',show_colorbar=True);
<Figure size 864x432 with 0 Axes>
_images/ECCO_v4_Example_MHT_20_1.png

Select the Atlantic ocean basin for RAPID-MOCHA MHT

Now that we have 26^\circN we just need to select the Atlantic. This can be done with the ecco_v4_py.get_basin module, specifically ecco_v4_py.get_basin.get_basin_mask. Note that this function requires a mask as an input, and then returns the values within a particular ocean basin. Therefore, provide the function with ds['maskC'] for a mask at tracer points, ds['maskW'] for west grid cell edges, etc.

Note: this mirrors the gcmfaces functionality ecco_v4/v4_basin.m.

Here we just want the Atlantic ocean, but lets look at what options we have …

[14]:
print(ecco.get_available_basin_names())
['pac', 'atl', 'ind', 'arct', 'bering', 'southChina', 'mexico', 'okhotsk', 'hudson', 'med', 'java', 'north', 'japan', 'timor', 'eastChina', 'red', 'gulf', 'baffin', 'gin', 'barents']

Notice that, for instance, ‘mexico’ exists separately from the Atlantic (‘atl’). This is to provide as fine granularity as possible (and sensible). To grab the broader Atlantic ocean basin, i.e. the one people usually refer to, use the option ‘atlExt’. Similar options exist for the Pacific and Indian ocean basins.

[15]:
atl_maskW = ecco.get_basin_mask(basin_name='atlExt',mask=maskW.isel(k=0))
atl_maskS = ecco.get_basin_mask(basin_name='atlExt',mask=maskS.isel(k=0))
loading  basins.data
data shape  (1170, 90)
dims, num_dims, llc  (1170, 90) 2 90
2 dimensions
f3 shape  (90, 90)
f5 shape  (90, 270)
2D, data_compact shape  (1170, 90)
data_tiles shape  (13, 90, 90)
data_tiles shape =  (13, 90, 90)
loading  basins.data
data shape  (1170, 90)
dims, num_dims, llc  (1170, 90) 2 90
2 dimensions
f3 shape  (90, 90)
f5 shape  (90, 270)
2D, data_compact shape  (1170, 90)
data_tiles shape  (13, 90, 90)
data_tiles shape =  (13, 90, 90)

Notice that we pass the routine a 2D mask by selecting the first depth level. This is simply to make things run faster.

[16]:
plt.figure(figsize=(12,6))
ecco.plot_proj_to_latlon_grid(ecco_ds.XC,ecco_ds.YC,atl_maskW,
                              projection_type='robin',cmap='viridis',user_lon_0=-30,show_colorbar=True);
_images/ECCO_v4_Example_MHT_28_0.png

MHT at the approximate RAPID array latitude

This can be done with the ecco_v4_py.calc_meridional_trsp module for heat, salt, and volume transport as follows:

mvt = ecco_v4_py.calc_meridional_vol_trsp(ecco_ds,lat_vals=26,basin_name='atlExt')
mht = ecco_v4_py.calc_meridional_heat_trsp(ecco_ds,lat_vals=26,basin_name='atlExt')
mst = ecco_v4_py.calc_meridional_salt_trsp(ecco_ds,lat_vals=26,basin_name='atlExt')

Additionally, one could compute the overturning streamfunction at this latitude band with ecco_v4_py.calc_meridional_stf. The inputs are the same as the functions above, see the module ecco_v4_py.calc_stf.

In MATLAB, one can use the functions:

A note about computational performance

When we originally open the dataset with all of the variables, we don’t actually load anything into memory. In fact, nothing actually happens until “the last minute”. For example, the data are only loaded once we do any computation like compute a mean value, plot something, or explicitly provide a load command for either the entire dataset or an individual DataArray within the dataset. This ‘lazy execution’ is enabled by the data structure underlying the xarray Datasets and DataArrays, the dask array.

In short, the when the data are opened, dask builds an execution task graph which it saves up to execute at the last minute. Dask also allows for parallelism, and by default runs in parallel across threads for a single core architecture. For now, this is what we will show.

Some quick tips are:

  • Don’t load the full 25 years of ECCOv4r3 output unless you’re on a machine with plenty of memory. I am doing this in the cell below because I’m on a Skylake node with 192GB. Proceed with caution before copying and pasting.
  • If you’re in this situation where you can’t load all months into memory, it’s a good idea to load a final result before plotting, in case you need to plot it many times in a row, see below…
[17]:
%%time
ecco_ds['ADVx_TH'].load();
ecco_ds['DFxE_TH'].load();
ecco_ds['ADVy_TH'].load();
ecco_ds['DFyE_TH'].load();
ecco_ds['dxG'].load();
ecco_ds['dyG'].load();
ecco_ds['drF'].load();
CPU times: user 15.8 s, sys: 52.1 s, total: 1min 7s
Wall time: 47.8 s
[17]:
<xarray.DataArray 'drF' (k: 50)>
array([ 10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.  ,  10.01,
        10.03,  10.11,  10.32,  10.8 ,  11.76,  13.42,  16.04,  19.82,
        24.85,  31.1 ,  38.42,  46.5 ,  55.  ,  63.5 ,  71.58,  78.9 ,
        85.15,  90.18,  93.96,  96.58,  98.25,  99.25, 100.01, 101.33,
       104.56, 111.33, 122.83, 139.09, 158.94, 180.83, 203.55, 226.5 ,
       249.5 , 272.5 , 295.5 , 318.5 , 341.5 , 364.5 , 387.5 , 410.5 ,
       433.5 , 456.5 ], dtype=float32)
Coordinates:
  * k        (k) int64 0 1 2 3 4 5 6 7 8 9 10 ... 40 41 42 43 44 45 46 47 48 49
    Z        (k) float32 -5.0 -15.0 -25.0 -35.0 ... -5039.25 -5461.25 -5906.25
    drF      (k) float32 10.0 10.0 10.0 10.0 10.0 ... 387.5 410.5 433.5 456.5
    PHrefC   (k) float32 49.05 147.15 245.25 ... 49435.043 53574.863 57940.312
Attributes:
    standard_name:  cell_z_size
    long_name:      cell z size
    units:          m
[18]:
%%time
trsp_x = ((ecco_ds['ADVx_TH'] + ecco_ds['DFxE_TH']) *  atl_maskW).sum('k')
trsp_y = ((ecco_ds['ADVy_TH'] + ecco_ds['DFyE_TH']) *  atl_maskS).sum('k')
CPU times: user 26.2 s, sys: 44.2 s, total: 1min 10s
Wall time: 1min 10s
[19]:
%%time
trsp_x = trsp_x * lat_maskW
trsp_y = trsp_y * lat_maskS

# Sum horizontally
trsp_x = trsp_x.sum(dim=['i_g','j','tile'])
trsp_y = trsp_y.sum(dim=['i','j_g','tile'])

# Add together to get transport at depth level
trsp_at_depth = trsp_x + trsp_y

# Sum over full depth and convert to PW
mht = (trsp_at_depth *10**-15 * 1000 * 4000)
mht.attrs['units']='PW'
CPU times: user 388 ms, sys: 788 ms, total: 1.18 s
Wall time: 1.17 s

Now that we have computed MHT, lets load the result for iterative plotting

For some reason when dask arrays are plotted, they are computed on the spot but don’t stay in memory. This takes a bit to get the hang of, but keep in mind that this allows us to scale the same code on distributed architecture, so we could use these same routines for high resolution output. This seems worthwhile!

Note that we probably don’t need this load statement if we have already loaded the underlying datasets.

[20]:
plt.figure(figsize=(12,6))
plt.plot(mht['time'],mht,'r')
plt.grid()
plt.title('Monthly Meridional Heat Transport at 26N')
plt.ylabel(f'[{mht.attrs["units"]}]')
plt.show()
_images/ECCO_v4_Example_MHT_35_0.png

Now compare global and Atlantic MHT at many latitudes

[21]:
# pick a temp directory for yourself
nc_save_dir = '/home/ifenty/tmp/'

if not os.path.isdir(nc_save_dir):
    os.makedirs(nc_save_dir)

nc_file = f'{nc_save_dir}/eccov4r3_mht.nc'
[25]:
global_lats = np.arange(-60,60,1)
mht = ecco.calc_meridional_heat_trsp(ecco_ds, lat_vals=global_lats)
mht = mht.rename({'heat_trsp':'global_heat_trsp'})
mht = mht.rename({'heat_trsp_z':'global_heat_trsp_z'})
print(' --- Done with global --- ')
 --- Done with global ---
[66]:
plt.figure(figsize=(12,6))
plt.plot(mht['lat'], mht['global_heat_trsp'].mean('time'))
plt.grid();
_images/ECCO_v4_Example_MHT_39_0.png
[27]:
basin_lats = np.arange(-35,60,1)
atl = ecco.calc_meridional_heat_trsp(ecco_ds,lat_vals=basin_lats,basin_name='atlExt')
atl = atl.rename({'heat_trsp':'atl_heat_trsp'})
atl = atl.rename({'heat_trsp_z':'atl_heat_trsp_z'})
print(' --- Done with Atlantic --- ')
loading  basins.data
data shape  (1170, 90)
dims, num_dims, llc  (1170, 90) 2 90
2 dimensions
f3 shape  (90, 90)
f5 shape  (90, 270)
2D, data_compact shape  (1170, 90)
data_tiles shape  (13, 90, 90)
data_tiles shape =  (13, 90, 90)
loading  basins.data
data shape  (1170, 90)
dims, num_dims, llc  (1170, 90) 2 90
2 dimensions
f3 shape  (90, 90)
f5 shape  (90, 270)
2D, data_compact shape  (1170, 90)
data_tiles shape  (13, 90, 90)
data_tiles shape =  (13, 90, 90)
 --- Done with Atlantic ---
[28]:
atl
[28]:
<xarray.Dataset>
Dimensions:          (k: 50, lat: 95, time: 288)
Coordinates:
  * time             (time) datetime64[ns] 1992-01-16T12:00:00 ... 2015-12-16T12:00:00
  * k                (k) int64 0 1 2 3 4 5 6 7 8 ... 41 42 43 44 45 46 47 48 49
  * lat              (lat) int64 -35 -34 -33 -32 -31 -30 ... 54 55 56 57 58 59
    drF              (k) float32 10.0 10.0 10.0 10.0 ... 387.5 410.5 433.5 456.5
    PHrefC           (k) float32 49.05 147.15 245.25 ... 53574.863 57940.312
    Z                (k) float32 -5.0 -15.0 -25.0 ... -5039.25 -5461.25 -5906.25
Data variables:
    atl_heat_trsp_z  (time, k, lat) float64 -0.2068 -0.2489 -0.2775 ... 0.0 0.0
    atl_heat_trsp    (time, lat) float64 0.1535 0.1445 0.1371 ... 0.5322 0.5043
[65]:
plt.figure(figsize=(12,6))
plt.plot(atl['lat'], atl['atl_heat_trsp'].mean('time'))
plt.grid()
_images/ECCO_v4_Example_MHT_42_0.png
[30]:
plt.figure(figsize=(12,6))
plt.plot(mht['lat'], mht['global_heat_trsp'].mean('time'))
plt.plot(atl['lat'], atl['atl_heat_trsp'].mean('time'))
plt.legend(('Global','Atlantic'))
plt.grid(linestyle='--')
plt.title(f'Meridional Heat Transport [{mht["global_heat_trsp"].attrs["units"]}]')
plt.ylabel(f'[{mht["global_heat_trsp"].attrs["units"]}]')
plt.xlabel('Latitude')
plt.show()
_images/ECCO_v4_Example_MHT_43_0.png

MHT as a function of depth

[63]:
def lat_depth_plot(mht,fld,label):
    fig = plt.figure(figsize=(12,6))

    # Set up depth coordinate
    depth = -mht['Z']
    stretch_depth = 1000

    # Set up colormap and colorbar
    cmap = 'RdBu_r'
    fld_mean = mht[fld].mean('time')
    abs_max = np.max(np.abs([fld_mean.min(),fld_mean.max()]))
    cmin = -abs_max*.1
    cmax = -cmin

    # First top 500m
    ax1 = plt.subplot(2,1,1)
    p1 = ax1.pcolormesh(mht['lat'],depth,fld_mean,cmap=cmap,vmin=cmin,vmax=cmax)
    plt.grid()

    # Handle y-axis
    ax1.invert_yaxis()
    plt.ylim([stretch_depth, 0])
    ax1.yaxis.axes.set_yticks(np.arange(stretch_depth,0,-100))
    plt.ylabel(f'Depth [{mht["Z"].attrs["units"]}]')

    # Remove top plot xtick label
    ax1.xaxis.axes.set_xticklabels([])

    # Now the rest ...
    ax2 = plt.subplot(2,1,2)
    p2 = ax2.pcolormesh(mht['lat'],depth, fld_mean, cmap=cmap,vmin=cmin,vmax=cmax)
    plt.grid()

    # Handle y-axis
    ax2.invert_yaxis()
    plt.ylim([4000, 1000])
    #yticks = np.flip(np.arange(6000,stretch_depth,-1000))
    #ax2.yaxis.axes.set_yticks(yticks)
    plt.ylabel(f'Depth [{mht["Z"].attrs["units"]}]')

    # Label  axis
    plt.xlabel('Latitude')

    # Reduce space between subplots
    fig.subplots_adjust(hspace=0.0)

    # Make a single title
    fig.suptitle(f'{label} time mean meridional heat transport',verticalalignment='top',fontsize=24)

    # Make an overyling colorbar
    fig.subplots_adjust(right=0.83)
    cbar_ax = fig.add_axes([0.87, 0.1, 0.025, 0.8])
    fig.colorbar(p2,cax=cbar_ax)
    cbar_ax.set_ylabel(f'[{mht[fld].attrs["units"]}]')

    plt.show()
[64]:
lat_depth_plot(mht,'global_heat_trsp_z','Global')
lat_depth_plot(atl,'atl_heat_trsp_z','Atlantic')
_images/ECCO_v4_Example_MHT_46_0.png
_images/ECCO_v4_Example_MHT_46_1.png

Exercise: reproduce figure from (Ganachaud and Wunsch, 2000)

[33]:
from IPython.display import Image
Image('../figures/buckley_mht.png')
[33]:
_images/ECCO_v4_Example_MHT_48_0.png

Figure from (Buckley and Marshall, 2016), which is adapted from (Ganachaud and Wunsch, 2000).

Note: to do this, you may need to pair latitude masks with those defined through the get_section_masks module. An example of this functionality is shown in the next tutorial.

References

Buckley, M. W. and Marshall, J. ( 2016), Observations, inferences, and mechanisms of the Atlantic Meridional Overturning Circulation: A review, Rev. Geophys., 54, doi:10.1002/2015RG000493.

Ganachaud, A., & Wunsch, C. (2000). Improved estimates of global ocean circulation, heat transport and mixing from hydrographic data. Nature, 408(6811), 453-7. doi:http://dx.doi.org.ezproxy.lib.utexas.edu/10.1038/35044048