Examing the transient and equilibrium CO2 response in the CESM

This notebook is part of The Climate Laboratory by Brian E. J. Rose, University at Albany.

I have run two sets of experiments with the CESM model:

The fully coupled model:
- pre-industrial control
- 1%/year CO $_2$ ramp scenario for 80 years
The slab ocean model:
- pre-industrial control with prescribed q-flux
- 2xCO $_2$ scenario run out to equilibrium

Our main first task is to compute the two canonical measures of climate sensitivity for this model:

Equilibrium Climate Sensitivity (ECS)
Transient Climate Response (TCR)

From the IPCC AR5 WG1 report, Chapter 9, page 817:

Equilibrium climate sensitivity (ECS) is the equilibrium change in global and annual mean surface air temperature after doubling the atmos- pheric concentration of CO2 relative to pre-industrial levels.

The transient climate response (TCR) is the change in global and annual mean surface temperature from an experiment in which the CO2 concentration is increased by 1% yr $^{–1}$ , and calculated using the difference between the start of the experiment and a 20-year period centred on the time of CO2 doubling.

First, a quick demonstration that 1%/year compounded increase reaches doubling after 70 years¶

startingamount = 1.
amount = startingamount
for n in range(70):
    amount *= 1.01
amount

2.006763368395386

TCR is always smaller than ECS due to the transient effects of ocean heat uptake.

We are going to estimate the ECS of the fully coupled model by using the equilibrium response of the Slab Ocean .

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr

Load the concatenated output from the CAM output (atmosphere)¶

casenames = {'cpl_control': 'cpl_1850_f19',
             'cpl_CO2ramp': 'cpl_CO2ramp_f19',
             'som_control': 'som_1850_f19',
             'som_2xCO2':   'som_1850_2xCO2',
            }
# The path to the THREDDS server, should work from anywhere
basepath = 'http://thredds.atmos.albany.edu:8080/thredds/dodsC/CESMA/'
# For better performance if you can access the roselab_rit filesystem (e.g. from JupyterHub)
# basepath = '/roselab_rit/cesm_archive/'
casepaths = {}
for name in casenames:
    casepaths[name] = basepath + casenames[name] + '/concatenated/'

# make a dictionary of all the CAM atmosphere output
atm = {}
for name in casenames:
    path = casepaths[name] + casenames[name] + '.cam.h0.nc'
    print('Attempting to open the dataset ', path)
    atm[name] = xr.open_dataset(path, decode_times=False)

Attempting to open the dataset  http://thredds.atmos.albany.edu:8080/thredds/dodsC/CESMA/cpl_1850_f19/concatenated/cpl_1850_f19.cam.h0.nc
Attempting to open the dataset  http://thredds.atmos.albany.edu:8080/thredds/dodsC/CESMA/cpl_CO2ramp_f19/concatenated/cpl_CO2ramp_f19.cam.h0.nc
Attempting to open the dataset  http://thredds.atmos.albany.edu:8080/thredds/dodsC/CESMA/som_1850_f19/concatenated/som_1850_f19.cam.h0.nc
Attempting to open the dataset  http://thredds.atmos.albany.edu:8080/thredds/dodsC/CESMA/som_1850_2xCO2/concatenated/som_1850_2xCO2.cam.h0.nc

A plot of the prescribed CO2 concentrations in the coupled simulations¶

days_per_year = 365
fig, ax = plt.subplots()
for name in ['cpl_control', 'cpl_CO2ramp']:
    ax.plot(atm[name].time/days_per_year, atm[name].co2vmr*1E6, label=name)
ax.set_title('CO2 volume mixing ratio (CESM coupled simulations)')
ax.set_xlabel('Years')
ax.set_ylabel('pCO2 (ppm)')
ax.grid()
ax.legend();

Issues to think about:

Why do we talk about fractional changes in CO $_2$ , such as “doubling atmospheric CO $_2$ ”, and "1%/year compounded CO $_2$ increase?
Why not instead talk about changes in absolute amounts of CO $_2$ ?

The answer is closely related to the fact that the radiative forcing associated with CO $_2$ increase is approximately logarithmic in CO $_2$ amount. So a doubling of CO $_2$ represents roughly the same radiative forcing regardless of the initial CO $_2$ concentration.

Compute and plot time series of global, annual mean near-surface air temperature in all four simulations¶

# The surface air temperature, which we will use for our sensitivity metrics
atm['cpl_control'].TREFHT

#  The area weighting needed for global averaging
gw = atm['som_control'].gw
print(gw)

<xarray.DataArray 'gw' (lat: 96)> Size: 768B
[96 values with dtype=float64]
Coordinates:
  * lat      (lat) float64 768B -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
Attributes:
    long_name:  gauss weights

def global_mean(field, weight=gw):
    '''Return the area-weighted global average of the input field'''
    return (field*weight).mean(dim=('lat','lon'))/weight.mean(dim='lat')

#  Loop through the four simulations and produce the global mean timeseries
TREFHT_global = {}
for name in casenames:
    TREFHT_global[name] = global_mean(atm[name].TREFHT)

Make some pretty timeseries plots, including an approximate running annual average¶

fig, axes = plt.subplots(2,1,figsize=(10,8))
for name in casenames:
    if 'cpl' in name:
        ax = axes[0]
        ax.set_title('Fully coupled ocean')
    else:
        ax = axes[1]
        ax.set_title('Slab ocean')
    field = TREFHT_global[name]
    field_running = field.rolling(time=12, center=True).mean()
    line = ax.plot(field.time / days_per_year, 
                   field, 
                   label=name,
                   linewidth=0.75,
                   )
    ax.plot(field_running.time / days_per_year, 
            field_running, 
            color=line[0].get_color(),
            linewidth=2,
           )
for ax in axes:
    ax.legend();
    ax.set_xlabel('Years')
    ax.set_ylabel('Temperature (K)')
    ax.grid();
    ax.set_xlim(0,100)
fig.suptitle('Global mean surface air temperature in CESM simulations', fontsize=16);

Issues to think about here include:

Why is the annual average here only approximate? (think about the calendar)
Why is there an annual cycle in the global average temperature? (planet is coldest during NH winter)
Different character of the temperature variability in the coupled vs. slab model
Much more rapid warming in the Slab Ocean Model

Now we can work on computing ECS and TCR¶

# extract the last 10 years from the slab ocean control simulation
# and the last 20 years from the coupled control
nyears_slab = 10
nyears_cpl = 20
clim_slice_slab = slice(-(nyears_slab*12),None)
clim_slice_cpl = slice(-(nyears_cpl*12),None)

# extract the last 10 years from the slab ocean control simulation
T0_slab = TREFHT_global['som_control'].isel(time=clim_slice_slab).mean(dim='time')
T0_slab

# and the last 20 years from the coupled control
T0_cpl = TREFHT_global['cpl_control'].isel(time=clim_slice_cpl).mean(dim='time')
T0_cpl

# extract the last 10 years from the slab 2xCO2 simulation
T2x_slab = TREFHT_global['som_2xCO2'].isel(time=clim_slice_slab).mean(dim='time')
T2x_slab

# extract the last 20 years from the coupled CO2 ramp simulation
T2x_cpl = TREFHT_global['cpl_CO2ramp'].isel(time=clim_slice_cpl).mean(dim='time')
T2x_cpl

ECS = T2x_slab - T0_slab
TCR = T2x_cpl - T0_cpl
print('The Equilibrium Climate Sensitivity is {:.3} K.'.format(float(ECS)))
print('The Transient Climate Response is {:.3} K.'.format(float(TCR)))

The Equilibrium Climate Sensitivity is 2.89 K.
The Transient Climate Response is 1.67 K.

Some CMIP climate sensitivity results to compare against¶

Figure and table from Flato et al. (2013)

Figure 9.43 | (a) Strengths of individual feedbacks for CMIP3 and CMIP5 models (left and right columns of symbols) for Planck (P), water vapour (WV), clouds (C), albedo (A), lapse rate (LR), combination of water vapour and lapse rate (WV+LR) and sum of all feedbacks except Planck (ALL), from Soden and Held (2006) and Vial et al. (2013), following Soden et al. (2008). CMIP5 feedbacks are derived from CMIP5 simulations for abrupt fourfold increases in CO2 concentrations (4 × CO2). (b) ECS obtained using regression techniques by Andrews et al. (2012) against ECS estimated from the ratio of CO2 ERF to the sum of all feedbacks. The CO2 ERF is one-half the 4 × CO2 forcings from Andrews et al. (2012), and the total feedback (ALL + Planck) is from Vial et al. (2013).

Figure caption reproduced from the AR5 WG1 report Flato et al., 2013

Comparing against the multi-model mean of the ECS and TCR, our model is apparently slightly less sensitive than the CMIP5 mean.

Let’s make some maps to compare spatial patterns of transient vs. equilibrium warming¶

Here is a helper function that takes a 2D lat/lon field and renders it as a nice contour map with accompanying zonal average line plot.

# The map projection capabilities come from the cartopy package. There are many possible projections
import cartopy.crs as ccrs
from cartopy.util import add_cyclic_point

def make_map(field, title="", levels=8):
    '''input field should be a 2D xarray.DataArray on a lat/lon grid.
        Make a filled contour plot of the field, and a line plot of the zonal mean
    '''
    fig = plt.figure(figsize=(14,6))
    nrows = 10; ncols = 3
    mapax = plt.subplot2grid((nrows,ncols), (0,0), colspan=ncols-1, rowspan=nrows-1, projection=ccrs.Robinson())
    barax = plt.subplot2grid((nrows,ncols), (nrows-1,0), colspan=ncols-1)
    plotax = plt.subplot2grid((nrows,ncols), (0,ncols-1), rowspan=nrows-1)
    # add cyclic point so cartopy doesn't show a white strip at zero longitude
    wrap_data, wrap_lon = add_cyclic_point(field.values, coord=field.lon, axis=field.dims.index('lon'))
    cx = mapax.contourf(wrap_lon, field.lat, wrap_data, levels=levels, transform=ccrs.PlateCarree())
    mapax.set_global(); mapax.coastlines();
    plt.colorbar(cx, cax=barax, orientation='horizontal', label='Surface air temperature (K)')
    plotax.plot(field.mean(dim='lon'), field.lat)
    plotax.set_xlabel('Zonal mean (K)')
    plotax.set_ylabel('Latitude')
    plotax.grid()
    fig.suptitle(title, fontsize=16)
    return fig, (mapax, plotax, barax), cx

# Plot a single time slice of surface air temperature just as example
fig, axes, cx = make_map(atm['cpl_control'].TREFHT.isel(time=0), 
                        title="Surface air temperature in the coupled control simulation")

Make maps of the surface air temperature anomaly due to CO2 doubling in both the slab and coupled models¶

Tmap_cpl_2x = atm['cpl_CO2ramp'].TREFHT.isel(time=clim_slice_cpl).mean(dim='time')
Tmap_cpl_control = atm['cpl_control'].TREFHT.isel(time=clim_slice_cpl).mean(dim='time')
DeltaT_cpl = Tmap_cpl_2x - Tmap_cpl_control

Tmap_som_2x = atm['som_2xCO2'].TREFHT.isel(time=clim_slice_slab).mean(dim='time')
Tmap_som_control = atm['som_control'].TREFHT.isel(time=clim_slice_slab).mean(dim='time')
DeltaT_som = Tmap_som_2x - Tmap_som_control

contours = np.arange(0., 8.5, 0.5)
fig, axes, cx = make_map(DeltaT_cpl, levels=contours,
                        title="Surface air temperature anomaly (coupled transient)")
axes[1].set_xlim(0,7)  # ensure the line plots have same axes

fig, axes,cx  = make_map(DeltaT_som, levels=contours,
                         title="Surface air temperature anomaly (equilibrium SOM)")
axes[1].set_xlim(0,7)

(0.0, 7.0)

Lots of intersting phenomena to think about here, including:

Polar amplification of surface warming
Reduction in equator-to-pole temperature gradients
Much larger polar amplification in SOM than in transient -- especially over the Southern Ocean (the delayed warming of the Southern Ocean)
North Atlantic warming hole present in transient but not in equilibrium SOM.
Land-ocean warming contrast: larger in transient, but still present in equilibrium

Appendix: for later reference, here is how you can open the other output types¶

The following will open the rest of the CESM output (land, sea ice, river routing, ocean).

These are not needed for the above homework assignment, but may be useful later on.

# # make a dictionary of all the CLM land model output
# land = {}
# for name in casenames:
#     path = casepaths[name] + casenames[name] + '.clm2.h0.nc'
#     print('Attempting to open the dataset ', path)
#     land[name] = xr.open_dataset(path)

# # make a dictionary of all the sea ice model output
# ice = {}
# for name in casenames:
#     path = casepaths[name] + casenames[name] + '.cice.h.nc'
#     print('Attempting to open the dataset ', path)
#     ice[name] = xr.open_dataset(path)

# # make a dictionary of all the river transport output
# rtm = {}
# for name in casenames:
#     path = casepaths[name] + casenames[name] + '.rtm.h0.nc'
#     print('Attempting to open the dataset ', path)
#     rtm[name] = xr.open_dataset(path)

# ocn = {}
# for name in casenames:
#     if 'cpl' in name:
#         path = casepaths[name] + casenames[name] + '.pop.h.nc'
#         print('Attempting to open the dataset ', path)
#         ocn[name] = xr.open_dataset(path)

Credits¶

This notebook is part of The Climate Laboratory, an open-source textbook developed and maintained by Brian E. J. Rose, University at Albany.

It is licensed for free and open consumption under the Creative Commons Attribution 4.0 International (CC BY 4.0) license.

Development of these notes and the climlab software is partially supported by the National Science Foundation under award AGS-1455071 to Brian Rose. Any opinions, findings, conclusions or recommendations expressed here are mine and do not necessarily reflect the views of the National Science Foundation.

References¶

Flato, G., Marotzke, J., Abiodun, B., Braconnot, P., Chou, S. C., Collins, W., Cox, P., Driouech, F., Emori, S., Eyring, V., Forest, C., Gleckler, P., Guilyardi, E., Jakob, C., Kattsov, V., Reason, C., & Rummukainen, M. (2013). Evaluation of Climate Models. In T. F. Stocker, D. Qin, G.-K. Plattner, M. Tignor, S. K. Allen, J. Boschung, A. Nauels, Y. Xia, V. Bex, & P. M. Midgley (Eds.), Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press.