Calculate Surface Ocean Heat using CESM2 LENS data

Overview¶

This notebook is an adpation of a workflow in the NCAR gallery of the Pangeo collection
This notebook illustrates how to compute surface ocean heat content using potential temperature data from CESM2 Large Ensemble Dataset (Community Earth System Model 2) hosted on NCAR’s GDEX.
This data is open access and is accessed via OSDF using the pelicanFS package and demonstrates how you can stream data from NCAR’s GDEX
Please refer to the first chapter of this cookbook to learn more about OSDF, pelican or pelicanFS

Prerequisites¶

Concepts	Importance	Notes
Intro to Intake-ESM	Necessary	Used for searching CMIP6 data
Understanding of Zarr	Helpful	Familiarity with metadata structure
Matplotlib	Helpful	Package used for plotting
PelicanFS	Necessary	The python package used to stream data in this notebook
OSDF	Helpful	OSDF is used to stream data in this notebook

Time to learn: 20 mins

Imports¶

import intake
import numpy as np
import pandas as pd
import xarray as xr
import seaborn as sns
import re
import matplotlib.pyplot as plt
import dask
from dask.distributed import LocalCluster
import pelicanfs 
import cf_units as cf

# Load Catalog URL
cat_url = 'https://osdata.gdex.ucar.edu/d010092/catalogs/d010092-osdf.json'

Set up local dask cluster¶

Before we do any computation let us first set up a local cluster using dask

cluster = LocalCluster()          
client = cluster.get_client()

# Scale the cluster
n_workers = 5
cluster.scale(n_workers)
cluster

Data Loading¶

Load CESM2 LENS data from NCAR’s GDEX¶

Load CESM2 LENS zarr data from GDEx using an intake-ESM catalog that has OSDF links
For more details regarding the dataset. See, https://gdex.ucar.edu/datasets/d010092/#

col = intake.open_esm_datastore(cat_url)
col

# Uncomment this line to see all the variables
# cesm_cat.df['variable'].values

cesm_temp = col.search(variable ='TEMP', frequency ='monthly')
cesm_temp

cesm_temp.df['path'].values

array(['osdf:///ncar-gdex/d010092/ocn/monthly/cesm2LE-historical-cmip6-TEMP.zarr',
       'osdf:///ncar-gdex/d010092/ocn/monthly/cesm2LE-ssp370-cmip6-TEMP.zarr',
       'osdf:///ncar-gdex/d010092/ocn/monthly/cesm2LE-ssp370-smbb-TEMP.zarr'],
      dtype=object)

dsets_cesm = cesm_temp.to_dataset_dict()


--> The keys in the returned dictionary of datasets are constructed as follows:
	'component.experiment.frequency.forcing_variant'

dsets_cesm.keys()

dict_keys(['ocn.historical.monthly.cmip6', 'ocn.ssp370.monthly.cmip6', 'ocn.ssp370.monthly.smbb'])

historical       = dsets_cesm['ocn.historical.monthly.cmip6']
future_smbb      = dsets_cesm['ocn.ssp370.monthly.smbb']
future_cmip6     = dsets_cesm['ocn.ssp370.monthly.cmip6']

Change units¶

orig_units = cf.Unit(historical.z_t.attrs['units'])
orig_units

Unit('centimeters')

def change_units(ds, variable_str, variable_bounds_str, target_unit_str):
    orig_units = cf.Unit(ds[variable_str].attrs['units'])
    target_units = cf.Unit(target_unit_str)
    variable_in_new_units = xr.apply_ufunc(orig_units.convert, ds[variable_bounds_str], target_units, dask='parallelized', output_dtypes=[ds[variable_bounds_str].dtype])
    return variable_in_new_units

depth_levels_in_m = change_units(historical, 'z_t', 'z_t', 'm')
hist_temp_in_degK = change_units(historical, 'TEMP', 'TEMP', 'degK')
fut_cmip6_temp_in_degK = change_units(future_cmip6, 'TEMP', 'TEMP', 'degK')
fut_smbb_temp_in_degK = change_units(future_smbb, 'TEMP', 'TEMP', 'degK')
#
hist_temp_in_degK  = hist_temp_in_degK.assign_coords(z_t=("z_t", depth_levels_in_m['z_t'].data))
hist_temp_in_degK["z_t"].attrs["units"] = "m"
hist_temp_in_degK

depth_levels_in_m.isel(z_t=slice(0, -1))

Compute depth level deltas using the difference of z_t levels

depth_level_deltas = depth_levels_in_m.isel(z_t=slice(1, None)).values - depth_levels_in_m.isel(z_t=slice(0, -1)).values
# Optionally, if you want to keep it as an xarray DataArray, re-wrap the result
depth_level_deltas = xr.DataArray(depth_level_deltas, dims=["z_t"], coords={"z_t": depth_levels_in_m.z_t.isel(z_t=slice(0, -1))})
depth_level_deltas

Compute Ocean Heat content for ocean surface¶

Ocean surface is considered to be the top 100m
The formula for this is:
$H = \rho C \int_0^z T(z) dz$
(1)

Where H is ocean heat content, the value we are trying to calculate,

$\rho$ is the density of sea water, $1026 kg/m^3$ ,

$C$ is the specific heat of sea water, $3990 J/(kg K)$ ,

$z$ is the depth limit of the calculation in meters,

and $T(z)$ is the temperature at each depth in degrees Kelvin.

def calc_ocean_heat(delta_level, temperature):
    rho = 1026 #kg/m^3
    c_p = 3990 #J/(kg K)
    weighted_temperature = delta_level * temperature
    heat = weighted_temperature.sum(dim="z_t")*rho*c_p
    return heat

# Remember that the coordinate z_t still has values in cm
hist_temp_ocean_surface = hist_temp_in_degK.where(hist_temp_in_degK['z_t'] < 1e4,drop=True)
hist_temp_ocean_surface

depth_level_deltas_surface = depth_level_deltas.where(depth_level_deltas['z_t'] <1e4, drop= True)
depth_level_deltas_surface

hist_ocean_heat = calc_ocean_heat(depth_level_deltas_surface,hist_temp_ocean_surface)
hist_ocean_heat

Plot Ocean Heat¶

%%time
# Compute annual and ensemble mean
hist_oceanheat_ann_mean = hist_ocean_heat.mean('member_id').groupby('time.year').mean()
hist_oceanheat_ann_mean

CPU times: user 419 ms, sys: 47 ms, total: 466 ms
Wall time: 454 ms

hist_oceanheat_ano = \
hist_oceanheat_ann_mean.sel(year=2014) - hist_oceanheat_ann_mean.sel(year=1850)

%%time
hist_oceanheat_ano.plot()

CPU times: user 7.89 s, sys: 1.53 s, total: 9.41 s
Wall time: 1min 24s

Indeed! The surface ocean is trapping heat as the globe warms!

Summary¶

In this notebook we used sea temperature data from the Community Earth System Model (CESM2) Large Ensemble dataset to compute surface ocean heat and convince ourselves that the surface ocean is trapping extra heat as the globe warms. We used an intake-ESM catalog backed by pelican links to stream data from NCAR’s Geoscience Data Exchange via NCAR’s OSDF origin!

What’s next?¶

In the next notebook, we will see how to load data from multiple OSDF origins into a workflow. We will stream CMIP6 model data from AWS and observational data from NCAR’s GDEX.

Resources and references¶

Original notebook on the pangeo NCAR gallery
CESM2 Large Ensemble Dataset (Community Earth System Model 2) hosted on NCAR’s GDEX.