Skip to article frontmatterSkip to article content

Feedbacks with ClimKern

ClimKern Python Package Logo

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 Python package called ClimKern Janoski et al., 2025 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:

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

Prerequisites

ConceptsImportanceNotes
Intro to XarrayNecessary
Intro to CartopyNecessary
Intro to MatplotlibHelpful
Loading CMIP6 Data with Intake-ESMHelpful
  • 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

Download the kernel

Normally, when ClimKern is installed, the user needs to download data on the command line from the Zenodo repository. However, we now have the data on Jetstream2, so we can download it from there and save it in our package repository.

First, set the URL and path to point to ClimKern. Also, specify which kernel you want.

URL = "https://js2.jetstream-cloud.org:8001/" # URL for jetstream access
path = f"pythia/ClimKern" # Location of ClimKern
kernel = "ERA5"

Next, read in the data from Jetstream2

# Read in data
# Set up access to jetstream2
fs = fsspec.filesystem("s3", anon=True, client_kwargs=dict(endpoint_url=URL))
pattern = f"s3://{path}/kernels/"+kernel+f"/TOA*.nc"

# Grab the data
files = sorted(fs.glob(pattern))
fs.invalidate_cache() # This is necessary to deal with peculiarities of objects served from jetstream2

# Open file and make it Xarray Dataset
kern = xr.open_dataset(fs.open(files[0]))

# Save path for later
path_out = files[0].split(kernel+"/",1)[1]

To save this data in ClimKern’s directory, we have to figure out where it is on the machine. After that, we will go ahead and save out the kernel as a netCDF in ClimKern’s local data directory.

# Get the package location
spec = importlib.util.find_spec("climkern")
package_dir = os.path.dirname(spec.origin)

# print(f"The package directory is: {package_dir}")

# Define the path where you want to save the netCDF file within the package directory
netcdf_path = os.path.join(package_dir,"data/kernels",kernel,path_out)

# Ensure the directory exists
os.makedirs(os.path.dirname(netcdf_path), exist_ok=True)

# Save the dataset as a netCDF file
kern.to_netcdf(netcdf_path)

Now, let’s just make sure we can retrieve the kernel.

ck.util.get_kern(kernel)
Loading...

Prepare the CMIP6 Data

To start, we will need some CMIP6 data to calculate feedbacks. We will use the preindustrial control and 4×CO2_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})
Loading...

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×CO2_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().compute()

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

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/m2^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 nameDescriptionUnits
rsusThe upwelling shortwave radiation at the surfaceW/m2^2
rsdsThe downwelling shortwave radiation at the surfaceW/m2^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
Loading...

This function allows for the quick creation of a pre-configured map, eliminating the need to repeat setup code by setting up the figure, axes, and Cartopy Robinson projection.

def base_map():
    proj= ccrs.Robinson()  #select Robinson projection

    fig, ax = plt.subplots(subplot_kw=dict(projection=proj)) #Create the figure and axis 
 
    ax.coastlines() # Add coastlines

    return fig, ax

Similar to the function above, these dictionaries encapsulate the repetitive code setup and provide a set of reusable arguments to simplify the plotting of data onto this map.

plotargs = {'transform': ccrs.PlateCarree(),
            'cbar_kwargs': {"orientation": "horizontal", "shrink": 0.7, 
                            "label":"W/m$^2$"}
            }
plotargs_2 = {'transform': ccrs.PlateCarree(),
            'cbar_kwargs': {"orientation": "horizontal", "shrink": 0.7,
                           "label":"W/m$^2$/K"}
            }

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

# Create the base map
fig, ax = base_map()

# Plot the surface albedo feedback 
alb.mean(dim="time").plot(ax=ax, **plotargs)

# Add a title
ax.set_title("Time-averaged Surface Albedo Feedback");
DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
<Figure size 640x480 with 2 Axes>

Now normalized by the surface temperature change...

# Create the base map
fig, ax = base_map()

# Plot surface albedo feedback & add coastlines
(alb.mean(dim="time")/Δts.mean(dim="time")).plot(ax=ax, **plotargs_2)

# Add a title
ax.set_title("Time-averaged Surface Albedo Feedback");
<Figure size 640x480 with 2 Axes>

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 nameDescriptionUnitsNotes
taThe 4-dimensional air temperatureK
tsThe 3-dimensional surface skin temperatureK
psThe 3-dimensional surface pressurePa, mb, or hPa
trop_pThe 3-dimensional tropopause heightPa, mb, or hPaOptional

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.

# Create the base map
fig, ax = base_map()

# Plot the surface albedo feedback 
lr.mean(dim="time").plot(ax=ax, **plotargs)

# Add a title
ax.set_title("Time-averaged Lapse Rate Feedback");
<Figure size 640x480 with 2 Axes>
# Create the base map
fig, ax = base_map()

# Plot surface albedo feedback
(lr.mean(dim="time")/Δts.mean(dim="time")).plot(ax=ax, **plotargs_2)

# Add a title
ax.set_title("Time-averaged Lapse Rate Feedback");
<Figure size 640x480 with 2 Axes>
# Create the base map
fig, ax = base_map()

# Plot surface albedo feedback
pl.mean(dim="time").plot(ax=ax, cmap="Blues_r", **plotargs)

# Add a title
ax.set_title("Time-averaged Planck Feedback");
<Figure size 640x480 with 2 Axes>
# Create the base map
fig, ax = base_map()

# Plot surface albedo feedback
(pl.mean(dim="time")/Δts.mean(dim="time")).plot(ax=ax,cmap="Blues_r", **plotargs_2)

# Add a title
ax.set_title("Time-averaged Planck Feedback");
<Figure size 640x480 with 2 Axes>

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 nameDescriptionUnitsNotes
husThe 4-dimensional specific humiditygg\frac{g}{g} or gkg\frac{g}{kg}
taThe 4-dimensional air temperatureK
psThe 3-dimensional surface pressurePa, mb, or hPa
trop_pThe 3-dimensional tropopause heightPa, mb, or hPaOptional

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 the discussion in Section 3.2 of Janoski et al. (2025), and see here in the climkern docs for how these are implemented.

We will use method=1 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=1
)

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

# Create the base map
fig, ax = base_map()

# Plot the longwave water vapor feedback
qlw.mean(dim="time").plot(ax=ax,cmap="Reds", **plotargs)

# Add a title
ax.set_title("Time-averaged LW Water Vapor Feedback");
<Figure size 640x480 with 2 Axes>
# Create the base map
fig, ax = base_map()

# Plot the longwave water vapor feedback 
(qlw.mean(dim="time")/Δts.mean(dim="time")).plot(ax=ax,cmap="Reds", **plotargs_2)
ax.coastlines()

# Add a title
ax.set_title("Time-averaged LW Water Vapor Feedback");
<Figure size 640x480 with 2 Axes>

And now the shortwave...

# Create the base map
fig, ax = base_map()

# Plot the longwave water vapor feedback and add coastlines
qsw.mean(dim="time").plot(ax=ax,cmap="Reds", **plotargs)


# Add a title
ax.set_title("Time-averaged SW Water Vapor Feedback");
<Figure size 640x480 with 2 Axes>
# Create the base map
fig, ax = base_map()

# Plot the longwave water vapor feedback and add coastlines
(qsw.mean(dim="time")/Δts.mean(dim="time")).plot(ax=ax,cmap="Reds", **plotargs_2)

# Add a title
ax.set_title("Time-averaged SW Water Vapor Feedback");
<Figure size 640x480 with 2 Axes>
# 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 CO2_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"]
Loading...

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!

References
  1. Janoski, T. P., Mitevski, I., Kramer, R. J., Previdi, M., & Polvani, L. M. (2025). ClimKern v1.2: a new Python package and kernel repository for calculating radiative feedbacks. Geoscientific Model Development, 18(10), 3065–3079. 10.5194/gmd-18-3065-2025
  2. Smith, C. J., Kramer, R. J., Myhre, G., Alterskjær, K., Collins, W., Sima, A., Boucher, O., Dufresne, J.-L., Nabat, P., Michou, M., Yukimoto, S., Cole, J., Paynter, D., Shiogama, H., O’Connor, F. M., Robertson, E., Wiltshire, A., Andrews, T., Hannay, C., … Forster, P. M. (2020). Effective radiative forcing and adjustments in CMIP6 models. Atmospheric Chemistry and Physics, 20(16), 9591–9618. 10.5194/acp-20-9591-2020
Foundations
Feedback analysis using radiative kernels