Tutorial 6: Ocean Heat Content
Contents
Tutorial 6: Ocean Heat Content#
Week 1, Day 2: Ocean and Atmospheric Reanalysis
Content creators: Aurora Basinski
Content reviewers: Katrina Dobson, Danika Gupta, Maria Gonzalez, Will Gregory, Nahid Hasan, Sherry Mi, Beatriz Cosenza Muralles, Jenna Pearson, Chi Zhang, Ohad Zivan
Content editors: Brodie, Pearson, Jenna Pearson, Chi Zhang, Ohad Zivan
Production editors: Wesley Banfield, Jenna Pearson, Chi Zhang, Ohad Zivan
Our 2023 Sponsors: NASA TOPS and Google DeepMind
Tutorial Objectives#
In this tutorial, you will explore the ocean’s vast heat capacity, which has a significant impact on the climate system.
The ocean has a heat capacity that is approximately 1000 times greater than the entire atmosphere, due to the relatively large mass and specific heat capacity of water. This means that the ocean has a significant impact on Earth’s thermal equilibrium state. Ocean heat uptake and ocean carbon uptake mitigates the effect of anthropogenic climate change by absorbing roughly 90% of the excess heat and 25% of human-emitted CO2. As you will see in this tutorial, ocean heat uptake largely occurs in the upper ocean as it is the region in contact with the atmosphere.
The ocean’s high heat capacity also facilitates meridional ocean heat transport from the equator to the poles, which acts in addition to the meridional atmospheric heat transport.
Through this tutorial, you will explore the spatial distribution of heat in the ocean and how the ocean’s heat content is changing over time. To do this, you will utilize the Estimating the Circulation and Climate of the Ocean (ECCO) dataset.
Setup#
# installations ( uncomment and run this cell ONLY when using google colab or kaggle )
# !pip install seaborn
# !pip install cmocean
# !pip install cartopy
# !pip install geoviews
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import seaborn as sns
import cartopy as cart
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cmocean
import os
import pooch
import tempfile
from scipy import integrate
Helper functions#
# @title Helper functions
def pooch_load(filelocation=None, filename=None, processor=None):
shared_location = "/home/jovyan/shared/Data/tutorials/W1D2_StateoftheClimateOceanandAtmosphereReanalysis" # this is different for each day
user_temp_cache = tempfile.gettempdir()
if os.path.exists(os.path.join(shared_location, filename)):
file = os.path.join(shared_location, filename)
else:
file = pooch.retrieve(
filelocation,
known_hash=None,
fname=os.path.join(user_temp_cache, filename),
processor=processor,
)
return file
Section 1: Vertical Distribution of Heat Within the Ocean#
First, let’s load ECCO ocean temperature (THETA). This datasets contains the annual mean temperature over the period of 1992 to 2016.
# import preprocessed ECCO data. This data is full depth temperature data over 1992 to 2016 (annual mean)
# this file takes about 5 minutes to load
filename_theta_annual = "theta_annual_mean.nc"
url_theta_annual = "https://osf.io/c8wqt/download"
theta_annual = xr.open_dataset(pooch_load(url_theta_annual, filename_theta_annual))
theta_annual = theta_annual.THETA
theta_annual = theta_annual.where(theta_annual != 0) # make land points equal to NaN
theta_annual
Downloading data from 'https://osf.io/c8wqt/download' to file '/tmp/theta_annual_mean.nc'.
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[6], line 6
3 filename_theta_annual = "theta_annual_mean.nc"
4 url_theta_annual = "https://osf.io/c8wqt/download"
----> 6 theta_annual = xr.open_dataset(pooch_load(url_theta_annual, filename_theta_annual))
7 theta_annual = theta_annual.THETA
8 theta_annual = theta_annual.where(theta_annual != 0) # make land points equal to NaN
Cell In[3], line 10, in pooch_load(filelocation, filename, processor)
8 file = os.path.join(shared_location, filename)
9 else:
---> 10 file = pooch.retrieve(
11 filelocation,
12 known_hash=None,
13 fname=os.path.join(user_temp_cache, filename),
14 processor=processor,
15 )
17 return file
File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/pooch/core.py:239, in retrieve(url, known_hash, fname, path, processor, downloader, progressbar)
236 if downloader is None:
237 downloader = choose_downloader(url, progressbar=progressbar)
--> 239 stream_download(url, full_path, known_hash, downloader, pooch=None)
241 if known_hash is None:
242 get_logger().info(
243 "SHA256 hash of downloaded file: %s\n"
244 "Use this value as the 'known_hash' argument of 'pooch.retrieve'"
(...)
247 file_hash(str(full_path)),
248 )
File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/pooch/core.py:803, in stream_download(url, fname, known_hash, downloader, pooch, retry_if_failed)
799 try:
800 # Stream the file to a temporary so that we can safely check its
801 # hash before overwriting the original.
802 with temporary_file(path=str(fname.parent)) as tmp:
--> 803 downloader(url, tmp, pooch)
804 hash_matches(tmp, known_hash, strict=True, source=str(fname.name))
805 shutil.move(tmp, str(fname))
File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/pooch/downloaders.py:226, in HTTPDownloader.__call__(self, url, output_file, pooch, check_only)
224 progress = self.progressbar
225 progress.total = total
--> 226 for chunk in content:
227 if chunk:
228 output_file.write(chunk)
File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/requests/models.py:816, in Response.iter_content.<locals>.generate()
814 if hasattr(self.raw, "stream"):
815 try:
--> 816 yield from self.raw.stream(chunk_size, decode_content=True)
817 except ProtocolError as e:
818 raise ChunkedEncodingError(e)
File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/urllib3/response.py:628, in HTTPResponse.stream(self, amt, decode_content)
626 else:
627 while not is_fp_closed(self._fp):
--> 628 data = self.read(amt=amt, decode_content=decode_content)
630 if data:
631 yield data
File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/urllib3/response.py:567, in HTTPResponse.read(self, amt, decode_content, cache_content)
564 fp_closed = getattr(self._fp, "closed", False)
566 with self._error_catcher():
--> 567 data = self._fp_read(amt) if not fp_closed else b""
568 if amt is None:
569 flush_decoder = True
File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/urllib3/response.py:533, in HTTPResponse._fp_read(self, amt)
530 return buffer.getvalue()
531 else:
532 # StringIO doesn't like amt=None
--> 533 return self._fp.read(amt) if amt is not None else self._fp.read()
File ~/miniconda3/envs/climatematch/lib/python3.10/http/client.py:466, in HTTPResponse.read(self, amt)
463 if self.length is not None and amt > self.length:
464 # clip the read to the "end of response"
465 amt = self.length
--> 466 s = self.fp.read(amt)
467 if not s and amt:
468 # Ideally, we would raise IncompleteRead if the content-length
469 # wasn't satisfied, but it might break compatibility.
470 self._close_conn()
File ~/miniconda3/envs/climatematch/lib/python3.10/socket.py:705, in SocketIO.readinto(self, b)
703 while True:
704 try:
--> 705 return self._sock.recv_into(b)
706 except timeout:
707 self._timeout_occurred = True
File ~/miniconda3/envs/climatematch/lib/python3.10/ssl.py:1274, in SSLSocket.recv_into(self, buffer, nbytes, flags)
1270 if flags != 0:
1271 raise ValueError(
1272 "non-zero flags not allowed in calls to recv_into() on %s" %
1273 self.__class__)
-> 1274 return self.read(nbytes, buffer)
1275 else:
1276 return super().recv_into(buffer, nbytes, flags)
File ~/miniconda3/envs/climatematch/lib/python3.10/ssl.py:1130, in SSLSocket.read(self, len, buffer)
1128 try:
1129 if buffer is not None:
-> 1130 return self._sslobj.read(len, buffer)
1131 else:
1132 return self._sslobj.read(len)
KeyboardInterrupt:
# take the mean over the period 1992 to 1994
subset_theta = theta_annual.sel(year=slice("1992", "1994")).mean("year")
subset_theta
# plot a zonal mean slice of this data. we take a mean just in longitudes by dividing the dA coordinate by the
fig, ax = plt.subplots(figsize=(8, 6))
(
((subset_theta - 273.15) * subset_theta.dX).sum("longitude")
/ (subset_theta.dX.sum("longitude"))
).plot(ax=ax)
ax.set_title("Global zonal mean of temperature (C)")
Temperatures are warmest near the ocean’s surface and close to the Equator, which suggests that heat is not distributed evenly within the ocean. In this tutorial you will explore the spatial distribution of heat storage in the ocean (i.e., ocean heat content), and how this distribution is changing over time.
Heat content is typically measured in Joules, which is equivalent to the units kg\(*\)m\(^2\)/s\(^2\). To determine the heat content within a specific depth range of the global ocean, i.e., between depth \(z_1\) and the surface (height of 0), you can calculate a volume integral over the dimensions \(x,y,\) and \(z\). This integral can be written as: $\(\iiint_{-z_1}^0 c_p \cdot \rho_0 \cdot \theta(x,y,z) dz dA\)$
Here, \(dA\) represents the area integral over the \(x\) and \(y\) (lat, lon) coordinates. \(\rho_0\) is the reference density in units of kg/m\(^3\) and \(c_p\) is specific heat capacity in units of J/(kg\(*\)K)
theta_area_int = (
(subset_theta * subset_theta.dA).sum("latitude").sum("longitude")
) # we take an area integral first at each depth level
rho = 1026 # kg/m^3
c_p = 3990 # J/(kg K)
sns.set_style(style="whitegrid")
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(
-subset_theta.Zu, (rho * c_p * theta_area_int * subset_theta.dZ).cumsum() / 10**27
)
ax.set_xlabel("Depth (m)")
ax.set_ylabel("Heat content above this depth ($10^6$ ZJ)")
ax.set_title("Global Ocean Heat Content above each depth")
You can see that much of the ocean’s heat is concentrated in the upper ocean (where the line is steep), with less heat stored in the deepest ocean regions (where the line plateaus). At first glance, this seems consistent with the zonal mean plot you plotted earlier in the tutorial, where the upper ocean tends to be warmer than deeper waters. However, in the integral equation above, \(\theta\) is not the only depth-dependent term. The global ocean area (\(A\)) also varies with depth, with the area of the global ocean decreasing with depth until only a few deep trenches contain water at the greatest ocean depths.
Let’s explore whether the ocean heat content plot we just created is driven by temperature variations or global ocean area variations with depth. One way to do this is to calculate and plot an integral of the global ocean area between each depth and the surface (i.e., the volume of the ocean above a each depth): \(Volume(z) = \iiint_{-z_1}^0 dz dA\).
If the volume as a function of depth looks similar to the heat content plot above, it would suggest that the smaller heat content of the deeper ocean (i.e., the plateau at large depths) is caused by the relatively small volume of water contained at these depths, rather than the vertical variations in temperature.
area_of_ocean = (
(subset_theta * subset_theta.dA / subset_theta).sum("latitude").sum("longitude")
) # we take an area integral first at each depth level
sns.set_style(style="whitegrid")
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(-subset_theta.Zu, (area_of_ocean * subset_theta.dZ).cumsum())
ax.set_xlabel("Depth (m)")
ax.set_ylabel("Volume of the global ocean above this depth (m$^3$)")
ax.set_title("Global ocean volume above each depth")
Question 1#
Based on the last two plots, are depth-variations in ocean heat content primarily due to vertical changes in the temperature or area of the ocean?
Section 2: Changes in Ocean Heat Content#
In this section you will examine how the total (i.e., full-depth) heat content of the ocean has changed over time. You will also explore heat content changes within two specific upper-ocean layers: one above 700 m depth and another above 2000 m depth\(^*\). By analyzing these near surface layers, you will identify whether changes in the ocean’s heat content are evenly distributed through the water column.
\(^*\)Note: technically the grid of the data means you will be looking above 677 m and 1997 m respectively
# this cell may take a while to run!
rho = 1026 # kg/m^3
c_p = 3990 # J/(kg K)
global_heat, years_to_plot, global_heat_upper2000, global_heat_upper700 = [], [], [], []
for year in theta_annual.year:
years_to_plot.append(int(year))
subset_theta_year = theta_annual.sel(year=int(year))
global_heat.append(
float(
rho
* c_p
* (subset_theta_year * subset_theta_year.dZ * subset_theta_year.dA)
.sum("Z")
.sum("latitude")
.sum("longitude")
)
)
global_heat_upper2000.append(
float(
rho
* c_p
* (
(
subset_theta_year.where(subset_theta_year.Zu > -2000)
* subset_theta_year.dZ
* subset_theta_year.dA
)
.sum("Z")
.sum("latitude")
.sum("longitude")
)
)
)
global_heat_upper700.append(
float(
rho
* c_p
* (
(
subset_theta_year.where(subset_theta_year.Zu > -700)
* subset_theta_year.dZ
* subset_theta_year.dA
)
.sum("Z")
.sum("latitude")
.sum("longitude")
)
)
)
# we now have lists, and list don't support math operations (-)
# we also divide the values by 10**21 to make them easier to read.
heat_anom_fulldepth = [
(heat - global_heat[0]) / 10**21 for heat in global_heat
] # remove year 1992
heat_anom_upper2000 = [
(heat - global_heat_upper2000[0]) / 10**21 for heat in global_heat_upper2000
] # remove year 1992
heat_anom_upper700 = [
(heat - global_heat_upper700[0]) / 10**21 for heat in global_heat_upper700
] # remove year 1992
heat_anom_upper2000_700 = [
a - b for a, b in zip(heat_anom_upper2000, heat_anom_upper700)
] # difference series between 2000 to 700
heat_anom_upperfulldepth_2000 = [
a - b for a, b in zip(heat_anom_fulldepth, heat_anom_upper2000)
] # difference series between fulldepth to 2000
fig, ax = plt.subplots(figsize=(8, 6))
sns.set_style(style="whitegrid")
ax.plot(years_to_plot, heat_anom_fulldepth, "k--")
ax.plot(years_to_plot, heat_anom_upper700)
ax.plot(years_to_plot, heat_anom_upper2000_700)
ax.plot(years_to_plot, heat_anom_upperfulldepth_2000)
ax.set_xlabel("Time")
ax.set_ylabel("Heat content change (ZJ)")
ax.legend(
[
"Full depth",
"Surface to 700 meters depth",
"700 to 2000 meters depth",
"Below 2000 meters depth",
]
)
ax.set_title("Change in heat over time")
Questions 2#
The volume of the ocean model in the reanalysis product does not change over time. Thus the changes in ocean heat content that you just calculated are caused by changes in the ocean’s temperature. Most of the ocean’s warming (heat gain) has been within the upper ocean (shallower than 700 m). The deeper ocean has also warmed, but not as substantially as near-surface waters.
Based on this graph, what percentage of the ocean’s heat gain since 1992 is contained within the top 2000 meters?
Section 3: Spatial Distribution of Ocean Heat Content#
You just saw that the ocean heat increase is concentrated near the ocean surface. Now you will explore where that heat is stored as a function of latitude and longitude. You can do this by creating a global map of ocean heat content in the upper 700 m of the ocean - which is essentially the same integral as above without the horizontal area integral: \(\int_{-700m}^0 c_p\rho_0\theta(x,y,z) dz\)
# first let's plot where heat is stored in the mean
fig, ax = plt.subplots(
subplot_kw={"projection": ccrs.PlateCarree()}, figsize=(11, 12), dpi=100
) # this is from cartopy https://rabernat.github.io/research_computing_2018/maps-with-cartopy.html
p = (
(
(rho * c_p * subset_theta.where(-subset_theta.Zu < 700) * subset_theta.dZ).sum(
"Z"
)
)
).plot(
vmin=7e11,
vmax=8.15e11,
cmap=cmocean.cm.thermal,
cbar_kwargs={
"shrink": 0.75,
"orientation": "horizontal",
"extend": "both",
"pad": 0.05,
"label": "",
},
ax=ax,
)
ax.coastlines(color="grey", lw=0.5)
ax.set_xticks([-180, -120, -60, 0, 60, 120, 180], crs=ccrs.PlateCarree())
ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.add_feature(cart.feature.LAND, zorder=100, edgecolor="k")
ax.set_title(
"Ocean Heat Content of top 700 m per unit area, mean of 1992 to 1994, J m$^{-2}$"
)
fig.tight_layout()
The lower latitude ocean contains more heat than the higher latitudes. This finding is consistent with your previous plot of warmer waters near the Equator during this tutorial.
Note: the color scale of this figure was chosen to emphasize latitudinal differences in ocean heat conent. As a result, some regions with shallow water depth display as black on the plot due to their relatively low column-integrated ocean heat content (\(<7 \times 10^{11} J m^{-2}\)). These black regions do not have zero ocean heat content.
Now let’s explore the spatial pattern of (full-depth) ocean heat content rate of change between 1992 and 2016.
# we already defined an object that's the mean over years 1992 to 1994 (subset_theta)
# now define an object that's the mean over 2014 to 2016
subset_theta_future = theta_annual.sel(year=slice("2014", "2016")).mean("year")
length_of_time_period = 24 * 60 * 60 * 365 * (2015 - 1993)
full_depth_heat_content_change = (
rho * c_p * subset_theta_future * subset_theta_future.dZ
).sum("Z") - (rho * c_p * subset_theta * subset_theta.dZ).sum("Z")
upper_700m_heat_content_change = (
rho
* c_p
* subset_theta_future.where(-subset_theta.Zu < 700)
* subset_theta_future.dZ
).sum("Z") - (
rho * c_p * subset_theta.where(-subset_theta.Zu < 700) * subset_theta.dZ
).sum(
"Z"
)
fig, ax = plt.subplots(
1, 2, subplot_kw={"projection": ccrs.PlateCarree()}, figsize=(11, 12), dpi=100
) # this is from cartopy https://rabernat.github.io/research_computing_2018/maps-with-cartopy.html
(full_depth_heat_content_change / length_of_time_period).plot(
ax=ax[0],
vmin=-10,
vmax=10,
cmap=cmocean.cm.balance,
cbar_kwargs={
"shrink": 0.75,
"orientation": "horizontal",
"extend": "both",
"pad": 0.05,
"label": "",
},
)
(upper_700m_heat_content_change / length_of_time_period).plot(
ax=ax[1],
vmin=-10,
vmax=10,
cmap=cmocean.cm.balance,
cbar_kwargs={
"shrink": 0.75,
"orientation": "horizontal",
"extend": "both",
"pad": 0.05,
"label": "",
},
)
ax[0].coastlines(color="grey", lw=0.5)
ax[0].set_xticks([-180, -120, -60, 0, 60, 120, 180], crs=ccrs.PlateCarree())
ax[0].set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax[0].add_feature(cart.feature.LAND, zorder=100, edgecolor="k")
ax[0].set_title(
"Rate of change in full-depth ocean heat content \n (2014-2016 minus 1992-1994); (J m$^{-2}$)/year"
)
fig.tight_layout
ax[1].coastlines(color="grey", lw=0.5)
ax[1].set_xticks([-180, -120, -60, 0, 60, 120, 180], crs=ccrs.PlateCarree())
ax[1].set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
ax[0].set_ylabel("Latitude")
ax[1].set_ylabel("")
ax[0].set_xlabel("Longitude")
ax[1].set_xlabel("Longitude")
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax[1].add_feature(cart.feature.LAND, zorder=100, edgecolor="k")
ax[1].set_title(
"Rate of change in upper 700 m ocean heat content \n (2014-2016 minus 1992-1994); (J m$^{-2}$)/year"
)
fig.tight_layout
With these plots, you demonstrated that ocean heat gain is not evenly distributed across space. When comparing the two plots, you once again see that the upper ocean contains a large fraction of the warming (recall that Equatorial regions contribute more to the global mean than high-latitude regions becasue of their relatively large area).
Summary#
In this tutorial, you have quantified the spatial patterns and temporal changes of the ocean’s heat content. You showed that the upper layers of the ocean contain most of the ocean’s heat content, due to their relatively large area (and hence volume) compared to the deepest ocean layers. These upper layers also experience a disproportionately large fraction of the ocean warming that has been observed in recent decades. You also found that heat content distribution varies by latitude and longitude, and is typically greater in the lower latitudes, and the ocean’s heat gain over time is not uniformly distributed across different oceanic regions.