Feedbacks with ClimKern

Overview

While it is good to know how to manually calculate radiative feedbacks, it is also time-consuming to constantly recreate the code. Here, we will use a new Python package called ClimKern that offers functions to calculate radiative feedbacks from climate model or reanalysis output. The advantages of using ClimKern go beyond making it simpler—it standardizes the methods and assumptions that go into these sometimes complicated calculations. We will use CMIP6 output in this notebook. You will learn how to:

Download CMIP6 data from Pangeo’s space on Google Cloud.
Obtain a radiative kernel from the ClimKern repository.
Calculate the surface albedo, temperature, and water vapor feedbacks.
Compute the cloud feedbacks using the two main methods: residual and adjustment.
Quantify stratospheric feedbacks.

Prerequisites

Concepts	Importance	Notes
Intro to Xarray	Necessary
Intro to Cartopy	Necessary
Intro to Matplotlib	Helpful
Loading CMIP6 Data with Intake-ESM	Helpful

Time to learn: 20 minutes

Imports

import climkern as ck
import intake
import matplotlib.pyplot as plt
import s3fs
import fsspec
import xarray as xr
import glob
import importlib.util
import os
import cartopy.crs as ccrs

%matplotlib inline
plt.rcParams["figure.dpi"] = 100

Prepare the CMIP6 Data

To start, we will need some CMIP6 data to calculate feedbacks. We will use the preindustrial control and 4×CO\(_2\) experiments from just one model (CESM2).

# Make a list of variables and experiments we need
var_list = ["rsds","rsus","ta","ts","ps","hus"]
exp_list = ["piControl","abrupt-4xCO2"]

# Specify data location, open it
cat_url = "https://storage.googleapis.com/cmip6/pangeo-cmip6.json"
col = intake.open_esm_datastore(cat_url)

# Create a catalog of matching simulations
cat = col.search(experiment_id=exp_list,source_id="CESM2",variable_id=var_list,
                table_id="Amon")

# Convert to a dictionary of Xarray Datasets
ds_dict = cat.to_dataset_dict(zarr_kwargs={'consolidated': True})

--> The keys in the returned dictionary of datasets are constructed as follows:
	'activity_id.institution_id.source_id.experiment_id.table_id.grid_label'

100.00% [2/2 00:04<00:00]

The data, especially the preindustrial control simulation that will serve as our control climate, is huge. We are going to only use the last 50 years of the control and last 30 years of the abrupt 4×CO\(_2\) simulation. There are a few extra coordinates and/or dimensions we don’t need, hence the squeeze().

# Our control simulation
ctrl = ds_dict["CMIP.NCAR.CESM2.piControl.Amon.gn"].isel(
    time=slice(-600,None)).squeeze().chunk({"plev":1,"time":-1})

# The increase CO2 aka "perturbed" simulation
pert = ds_dict["CMIP.NCAR.CESM2.abrupt-4xCO2.Amon.gn"].isel(
    time=slice(-360,None)).squeeze().chunk({"plev":1,"time":-1})

To save a little diskspace because ClimKern still isn’t dask compatible yet, let’s make our control climatology ahead of time.

ctrl = ctrl.groupby(ctrl.time.dt.month).mean(dim="time").rename({"month":"time"})
pert = pert.groupby(pert.time.dt.month).mean(dim="time").rename({"month":"time"})

Feedback calculations

Calculate surface temperature change

ClimKern outputs feedbacks as radiative perturbations at the TOA (because there are multiple ways to quantify feedbacks). One of the most common ways is to normalize by the change in the surface temperature, which puts values in “true” feedback units of W/m\(^2\)/K. For that, we start by calculating the change in the surface temperature in our CESM2 experiment.

# Take difference between pert and ctrl ts fields
Δts = pert.ts - ctrl.ts

Surface albedo feedback

Perhaps the simplest of the radiative feedbacks to calculate is the surface albedo feedback because it requires only two variables and is a 3-dimensional field (no height coordinate). The variables we will use are as follows:

Variable name	Description	Units
rsus	The upwelling shortwave radiation at the surface	W/m\(^2\)
rsds	The downwelling shortwave radiation at the surface	W/m\(^2\)

The calc_alb_feedback() function is straightforward. Be sure to include our kernel name in the arguments.

alb = ck.calc_alb_feedback(
    ctrl.rsus,ctrl.rsds,pert.rsus,pert.rsds,kern="ERA5"
)

Let’s see what the resulting Dataarray looks like.

alb

<xarray.DataArray (time: 12, lat: 192, lon: 288)> Size: 5MB
dask.array<mul, shape=(12, 192, 288), dtype=float64, chunksize=(1, 192, 288), chunktype=numpy.ndarray>
Coordinates:
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB 0.0 1.25 2.5 3.75 ... 356.2 357.5 358.8
    member_id       <U8 32B 'r1i1p1f1'
    dcpp_init_year  float64 8B nan
  * time            (time) int64 96B 1 2 3 4 5 6 7 8 9 10 11 12

Finally, let’s plot the time average surface albedo feedback on a map, first as just the perturbation.

# Set a few variables
proj = ccrs.Robinson()

# Create figure and axis
fig, ax = plt.subplots(subplot_kw=dict(projection=proj),
                       layout="constrained")

# Plot surface albedo feedback & add coastlines
alb.mean(dim="time").plot(ax=ax,transform=ccrs.PlateCarree(),
                         cbar_kwargs={"orientation": "horizontal", "shrink": 0.7,
                                      "label":"W/m$^2$"})
ax.coastlines()

# Add a title
ax.set_title("Time-averaged Surface Albedo Feedback")

plt.show()

DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip

../../_images/55ae980841622d34d344679f565cb2953d7a299c4b30b066492db05486f148f3.png

Now normalized by the surface temperature change…

# Set a few variables
proj = ccrs.Robinson()

# Create figure and axis
fig, ax = plt.subplots(subplot_kw=dict(projection=proj),
                       layout="constrained")

# Plot surface albedo feedback & add coastlines
(alb.mean(dim="time")/Δts.mean(dim="time")).plot(ax=ax,transform=ccrs.PlateCarree(),
                         cbar_kwargs={"orientation": "horizontal", "shrink": 0.7,
                                      "label":"W/m$^2$/K"})
ax.coastlines()

# Add a title
ax.set_title("Time-averaged Surface Albedo Feedback")

plt.show()

../../_images/3d27f8ed3e6f2c9c3f6c05742792bd1819d32d1659e1bdde60ca29c72620057d.png

As one might expect, the greatest radiative change at the TOA from the surface albedo feedback is near the poles where there is melting sea ice and snow.

Let’s end by calculating the global average. ClimKern comes with a spat_avg() function.

# Output the time-averaged global surface albedo feedback
output = "The global average surface albedo feedback is {value:.2f} W/m\u00b2/K."
print(output.format(value=float(ck.spat_avg((alb.mean(dim="time")/Δts.mean(dim="time")
                                         ).expand_dims("time")))))

The global average surface albedo feedback is 0.23 W/m²/K.

Planck and lapse rate feedbacks

As a reminder, the total tropospheric temperature feedback is often decomposed into a distict Planck and lapse rate feedback. ClimKern has a single function that produces both feedbacks simultaneously called calc_T_feedbacks(). The variables we need this time are:

Variable name	Description	Units	Notes
ta	The 4-dimensional air temperature	K
ts	The 3-dimensional surface skin temperature	K
ps	The 3-dimensional surface pressure	Pa, mb, or hPa
trop_p	The 3-dimensional tropopause height	Pa, mb, or hPa	Optional

CMIP6 output rarely contains the tropopause height, so we will create a default tropopause that decreases linearly with the cosine of latitude from 100 hPa at the Equator to 300 hPa at the poles.

# Use ClimKern's "hidden" make_tropo function
pert["trop_p"] = ck.util.make_tropo(pert.ps)
pert.trop_p.attrs["units"] = "Pa"

lr,pl = ck.calc_T_feedbacks(
    ctrl.ta,ctrl.ts,ctrl.ps,pert.ta,pert.ts,pert.ps,pert_trop=pert.trop_p,kern="ERA5"
)

Simple. Let’s see what it looks like.

# Set a few variables
proj = ccrs.Robinson()

# Create figure and axis
fig, ax = plt.subplots(subplot_kw=dict(projection=proj),
                       layout="constrained")

# Plot surface albedo feedback & add coastlines
lr.mean(dim="time").plot(ax=ax,transform=ccrs.PlateCarree(),
                         cbar_kwargs={"orientation": "horizontal", "shrink": 0.7,
                                      "label":"W/m$^2$"})
ax.coastlines()

# Add a title
ax.set_title("Time-averaged Lapse Rate Feedback")

plt.show()

../../_images/284c158d0f6208d14b7e8da3aa975dc418e878e2d91512ea015495282540ab5f.png

# Set a few variables
proj = ccrs.Robinson()

# Create figure and axis
fig, ax = plt.subplots(subplot_kw=dict(projection=proj),
                       layout="constrained")

# Plot surface albedo feedback & add coastlines
(lr.mean(dim="time")/Δts.mean(dim="time")).plot(ax=ax,transform=ccrs.PlateCarree(),
                         cbar_kwargs={"orientation": "horizontal", "shrink": 0.7,
                                      "label":"W/m$^2$/K"})
ax.coastlines()

# Add a title
ax.set_title("Time-averaged Lapse Rate Feedback")

plt.show()

../../_images/a531bb84244f4153f01594b038c5bf5b720873161a61c893bb8e9236fa5ad5d5.png

# Set a few variables
proj = ccrs.Robinson()

# Create figure and axis
fig, ax = plt.subplots(subplot_kw=dict(projection=proj),
                       layout="constrained")

# Plot surface albedo feedback & add coastlines
pl.mean(dim="time").plot(ax=ax,transform=ccrs.PlateCarree(),cmap="Blues_r",
                         cbar_kwargs={"orientation": "horizontal", "shrink": 0.7,
                                      "label":"W/m$^2$"})
ax.coastlines()

# Add a title
ax.set_title("Time-averaged Planck Feedback")

plt.show()

../../_images/b256ef30be8e8c63b7883d3c7d943975089d4f19692c4ec4b2b99fd5435d5947.png

# Set a few variables
proj = ccrs.Robinson()

# Create figure and axis
fig, ax = plt.subplots(subplot_kw=dict(projection=proj),
                       layout="constrained")

# Plot surface albedo feedback & add coastlines
(pl.mean(dim="time")/Δts.mean(dim="time")).plot(
    ax=ax,transform=ccrs.PlateCarree(),cmap="Blues_r",
    cbar_kwargs={"orientation": "horizontal", "shrink": 0.7,"label":"W/m$^2$"}
)
ax.coastlines()

# Add a title
ax.set_title("Time-averaged Planck Feedback")

plt.show()

../../_images/a5b3d8a3c64c590d85f07da85f6e4e5857c56d7a4f0e6e2f827e0b55c3f1cad5.png

These results seem reasonable. The Planck feedback is often the feedback with the largest magnitude and scales nonlinearly with temperature change. The lapse rate feedback changes sign depending on the vertical structure of the atmosphere is usually positive (warming) at the poles and negative (cooling) in the tropics.

Let’s get those global averages again!

# Output the time-averaged global temperature feedbacks
output = ("The global average lapse rate and Planck feedbacks are"+
          " {lr:.2f} and {pl:.2f} W/m\u00b2/K, respectively.")
print(output.format(lr=float(ck.spat_avg((lr.mean(dim="time")/Δts.mean(dim="time")
                                         ).expand_dims("time"))),
                   pl=float(ck.spat_avg((pl.mean(dim="time")/Δts.mean(dim="time")
                                        ).expand_dims("time")))))

The global average lapse rate and Planck feedbacks are -0.78 and -3.33 W/m²/K, respectively.

Water vapor feedbacks

Here we use ClimKern to calculate the longwave and shortwave water vapor feedabcks. ClimKern has a single function that produces both called calc_q_feedbacks(). This function requires both the specific humidity and air temperature because ClimKern has to normalize the water vapor kernels to make them compatible with the specific humidity.

Variable name	Description	Units	Notes
hus	The 4-dimensional specific humidity	\(\frac{g}{g}\) or \(\frac{g}{kg}\)
ta	The 4-dimensional air temperature	K
ps	The 3-dimensional surface pressure	Pa, mb, or hPa
trop_p	The 3-dimensional tropopause height	Pa, mb, or hPa	Optional

There is one more argument we can provide called method, which is simply a number 1-4. This corresponds to whether to use a fractional approximation for logarithms, use the logarithms explicitly, or use the linear change in water vapor concentration. See here for complete descriptions. We will use method 3 which is to use the natural logarithm of specific humidity without doing a fractional approximation.

We will use the same tropopause height as before.

# Use ClimKern's water vapor feedback function
qlw,qsw = ck.calc_q_feedbacks(
    ctrl.hus,ctrl.ta,ctrl.ps,pert.hus,pert.ps,pert_trop=pert.trop_p,kern="ERA5",
    method=3
)

Let’s take a look, starting with the longwave.

# Set a few variables
proj = ccrs.Robinson()

# Create figure and axis
fig, ax = plt.subplots(subplot_kw=dict(projection=proj),
                       layout="constrained")

# Plot the longwave water vapor feedback and add coastlines
qlw.mean(dim="time").plot(ax=ax,transform=ccrs.PlateCarree(),cmap="Reds",
                         cbar_kwargs={"orientation": "horizontal", "shrink": 0.7,
                                      "label":"W/m$^2$"})
ax.coastlines()

# Add a title
ax.set_title("Time-averaged LW Water Vapor Feedback")

plt.show()

../../_images/757b0823c38b67fe5bfcdc20e9204bdddf962bc8ab0588f66796d7059da63adb.png

# Set a few variables
proj = ccrs.Robinson()

# Create figure and axis
fig, ax = plt.subplots(subplot_kw=dict(projection=proj),
                       layout="constrained")

# Plot the longwave water vapor feedback and add coastlines
(qlw.mean(dim="time")/Δts.mean(dim="time")).plot(
    ax=ax,transform=ccrs.PlateCarree(),cmap="Reds",
    cbar_kwargs={"orientation": "horizontal", "shrink": 0.7,
                 "label":"W/m$^2$/K"}
)
ax.coastlines()

# Add a title
ax.set_title("Time-averaged LW Water Vapor Feedback")

plt.show()

../../_images/0e3b88660b575f7eadfa8fa9c27578d39f83ec566b4006f193c472f24a94b84b.png

And now the shortwave…

# Set a few variables
proj = ccrs.Robinson()

# Create figure and axis
fig, ax = plt.subplots(subplot_kw=dict(projection=proj),
                       layout="constrained")

# Plot the longwave water vapor feedback and add coastlines
qsw.mean(dim="time").plot(ax=ax,transform=ccrs.PlateCarree(),cmap="Reds",
                         cbar_kwargs={"orientation": "horizontal", "shrink": 0.7,
                                      "label":"W/m$^2$"})
ax.coastlines()

# Add a title
ax.set_title("Time-averaged SW Water Vapor Feedback")

plt.show()

../../_images/ffe650dc615470dcaa076eb8c252bc4f794569a96c3c52c7f11cfb1f05b4036e.png

# Set a few variables
proj = ccrs.Robinson()

# Create figure and axis
fig, ax = plt.subplots(subplot_kw=dict(projection=proj),
                       layout="constrained")

# Plot the longwave water vapor feedback and add coastlines
(qsw.mean(dim="time")/Δts.mean(dim="time")).plot(
    ax=ax,transform=ccrs.PlateCarree(),cmap="Reds",
                         cbar_kwargs={"orientation": "horizontal", "shrink": 0.7,
                                      "label":"W/m$^2$/K"}
)
ax.coastlines()

# Add a title
ax.set_title("Time-averaged SW Water Vapor Feedback")

plt.show()

../../_images/6b36c16b0245a63d64e1cb2d8f1c60074524b65bdcb77c0d1950d708cf9e5234.png

# Output the time-averaged global surface albedo feedback
output = ("The global average longwave and shortwave water vapor feedbacks are"+
          " {lw:.2f} and {sw:.2f} W/m\u00b2/K, respectively.")
print(output.format(lw=float(ck.spat_avg((qlw.mean(dim="time")/Δts.mean(dim="time")
                                         ).expand_dims("time"))),
                   sw=float(ck.spat_avg((qsw.mean(dim="time")/Δts.mean(dim="time")
                                        ).expand_dims("time")))))

The global average longwave and shortwave water vapor feedbacks are 1.59 and 0.19 W/m²/K, respectively.

Cloud feedbacks

Download radiative forcing

It’s time for cloud feedbacks. However, like anything to do with clouds, there are some complications. You ideally need the effective radiative forcing (ERF), which is a measure of radiative imbalance imposed by CO\(_2\) in the absence of surface temperature change Smith et al., 2020. Thankfully, CESM2 radiative forcing for our scenarios is available on the Google cloud.

# Specify data location, open it
cat_url = "https://storage.googleapis.com/cmip6/pangeo-cmip6.json"
col = intake.open_esm_datastore(cat_url)

# # Create a catalog of matching simulations
cat = col.search(activity_id="RFMIP",source_id="CESM2",experiment_id=["piClim-4xCO2",
                                                                      "piClim-control"],
                variable_id=["rlut","rsut","rsdt"])

# Convert to a dictionary of Xarray Datasets
ds_dict = cat.to_dataset_dict(zarr_kwargs={'consolidated': True})

# Assign Datasets to variables
rf_ctrl = ds_dict["RFMIP.NCAR.CESM2.piClim-control.Amon.gn"]
rf_pert = ds_dict["RFMIP.NCAR.CESM2.piClim-4xCO2.Amon.gn"]

--> The keys in the returned dictionary of datasets are constructed as follows:
	'activity_id.institution_id.source_id.experiment_id.table_id.grid_label'

100.00% [2/2 00:01<00:00]

Info

Your relevant information here!

Feel free to copy this around and edit or play around with yourself. Some other admonitions you can put in:

Success

We got this done after all!

Warning

Be careful!

Danger

Scary stuff be here.

We also suggest checking out Jupyter Book’s brief demonstration on adding cell tags to your cells in Jupyter Notebook, Lab, or manually. Using these cell tags can allow you to customize how your code content is displayed and even demonstrate errors without altogether crashing our loyal army of machines!