Estimating Equilibrium Climate Sensitivity

---

Overview

Equilibrium Climate Sensitivity (ECS) is defined as change in global-mean near-surface air temperature (GMST) change due to an instantaneous doubling of CO\(_2\) concentrations and once the coupled ocean-atmosphere-sea ice system has acheived a statistical equilibrium (i.e. at the top-of-atmosphere, incoming solar shortwave radiation is balanced by reflected solar shortwave and outgoing thermal longwave radiation).

This notebook uses the “Gregory method” to approximate the ECS of CMIP6 models based on the first 150 years after an abrupt doubling of CO\(_2\) concentrations. The Gregory method extrapolates the quasi-linear relationship between GMST and radiative imbalance at the top-of-atmosphere to estimate how much warming would occur if the system were in radiative balance at the top-of-atmosphere, which is by definition the equilibrium response. In particular, we extrapolate the linear relationship that occurs between 100 and 150 years after the abrupt quadrupling.

Since the radiative forcing due to CO\(_2\) is a logarithmic function of the CO\(_2\) concentration, the GMST change from a first doubling is roughly the same as for a second doubling (to first order, we can assume feedbacks as constant), which means that the GMST change due to a quadrupling of CO\(_2\) is roughly \(\Delta T_{4\times\mathrm{CO}_2}=2\times\mathrm{ECS}\). See also Mauritsen et al. 2019 for a detailed application of the Gregory method (with modifications) for the case of one specific CMIP6 model, the MPI-M Earth System Model.

For another take on applying the Gregory method to estimate ECS, see Angeline Pendergrass’ code.

Prerequisites

Concepts	Importance	Notes
Intro to Xarray	Necessary
Understanding of NetCDF	Helpful	Familiarity with metadata structure
Dask	Helpful
Climate sensitivity	Helpful

• Time to learn: 30 minutes

---

Imports

from matplotlib import pyplot as plt
import sys
import numpy as np
import pandas as pd
import xarray as xr
import cartopy
import dask
from tqdm.autonotebook import tqdm  # Fancy progress bars for our loops!
import intake
import fsspec
from dask_gateway import Gateway
from dask.distributed import Client

%matplotlib inline
plt.rcParams['figure.figsize'] = 12, 6

/tmp/ipykernel_2669/1804016931.py:8: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  from tqdm.autonotebook import tqdm  # Fancy progress bars for our loops!

Compute Cluster

Here we use a dask cluster to parallelize our analysis.

platform = sys.platform

if (platform == 'win32'):
    import multiprocessing.popen_spawn_win32
else:
    import multiprocessing.popen_spawn_posix

Initiate the Dask client:

client = Client()
client

Client

Client-944f1b53-8cca-11ef-8a6d-6045bdb5c525

Connection method: Cluster object	Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:8787/status

Cluster Info

LocalCluster

d4235bcf

Dashboard: http://127.0.0.1:8787/status	Workers: 4
Total threads: 4	Total memory: 15.61 GiB
Status: running	Using processes: True

Scheduler Info

Scheduler

Scheduler-f9b351a0-5702-4e1d-946a-933a783f036c

Comm: tcp://127.0.0.1:40071	Workers: 4
Dashboard: http://127.0.0.1:8787/status	Total threads: 4
Started: Just now	Total memory: 15.61 GiB

Workers

Worker: 0

Comm: tcp://127.0.0.1:46373	Total threads: 1
Dashboard: http://127.0.0.1:44017/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:36377
Local directory: /tmp/dask-scratch-space/worker-lh0zh583

Worker: 1

Comm: tcp://127.0.0.1:42069	Total threads: 1
Dashboard: http://127.0.0.1:33543/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:40419
Local directory: /tmp/dask-scratch-space/worker-3o9fh9nh

Worker: 2

Comm: tcp://127.0.0.1:38899	Total threads: 1
Dashboard: http://127.0.0.1:36021/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:41763
Local directory: /tmp/dask-scratch-space/worker-ztaoxnwa

Worker: 3

Comm: tcp://127.0.0.1:42395	Total threads: 1
Dashboard: http://127.0.0.1:46269/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:42053
Local directory: /tmp/dask-scratch-space/worker-v1wijhbb

Data catalogs

This notebook uses intake-esm to ingest and organize climate model output from the fresh-off-the-supercomputers Phase 6 of the Coupled Model Intercomparison Project (CMIP6).

The file https://storage.googleapis.com/cmip6/cmip6-zarr-consolidated-stores.csv in Google Cloud Storage contains thousands of lines of metadata, each describing an individual climate model experiment’s simulated data.

For example, the first line in the .csv file contains the precipitation rate (variable_id = 'pr'), as a function of latitude, longitude, and time, in an individual climate model experiment with the BCC-ESM1 model (source_id = 'BCC-ESM1') developed by the Beijing Climate Center (institution_id = 'BCC'). The model is forced by the forcing experiment SSP370 (experiment_id = 'ssp370'), which stands for the Shared Socio-Economic Pathway 3 that results in a change in radiative forcing of \(\Delta F=7.0\) W m\(^{-2}\) from pre-industrial to 2100. This simulation was run as part of the AerChemMIP activity, which is a spin-off of the CMIP activity that focuses specifically on how aerosol chemistry affects climate.

df = pd.read_csv('https://storage.googleapis.com/cmip6/cmip6-zarr-consolidated-stores.csv')
df.head()

	activity_id	institution_id	source_id	experiment_id	member_id	table_id	variable_id	grid_label	zstore	dcpp_init_year	version
0	HighResMIP	CMCC	CMCC-CM2-HR4	highresSST-present	r1i1p1f1	Amon	ps	gn	gs://cmip6/CMIP6/HighResMIP/CMCC/CMCC-CM2-HR4/...	NaN	20170706
1	HighResMIP	CMCC	CMCC-CM2-HR4	highresSST-present	r1i1p1f1	Amon	rsds	gn	gs://cmip6/CMIP6/HighResMIP/CMCC/CMCC-CM2-HR4/...	NaN	20170706
2	HighResMIP	CMCC	CMCC-CM2-HR4	highresSST-present	r1i1p1f1	Amon	rlus	gn	gs://cmip6/CMIP6/HighResMIP/CMCC/CMCC-CM2-HR4/...	NaN	20170706
3	HighResMIP	CMCC	CMCC-CM2-HR4	highresSST-present	r1i1p1f1	Amon	rlds	gn	gs://cmip6/CMIP6/HighResMIP/CMCC/CMCC-CM2-HR4/...	NaN	20170706
4	HighResMIP	CMCC	CMCC-CM2-HR4	highresSST-present	r1i1p1f1	Amon	psl	gn	gs://cmip6/CMIP6/HighResMIP/CMCC/CMCC-CM2-HR4/...	NaN	20170706

The file pangeo-cmip6.json describes the structure of the CMIP6 metadata and is formatted so as to be read in by the intake.open_esm_datastore method, which categorizes all of the data pointers into a tiered collection. For example, this collection contains the simulated data from 28691 individual experiments, representing 48 different models from 23 different scientific institutions. There are 190 different climate variables (e.g. sea surface temperature, sea ice concentration, atmospheric winds, dissolved organic carbon in the ocean, etc.) available for 29 different forcing experiments.

Use Intake-ESM

Intake-ESM is a new package designed to make working with these data archives a bit simpler.

col = intake.open_esm_datastore("https://storage.googleapis.com/cmip6/pangeo-cmip6.json")
col

pangeo-cmip6 catalog with 7674 dataset(s) from 514818 asset(s):

	unique
activity_id	18
institution_id	36
source_id	88
experiment_id	170
member_id	657
table_id	37
variable_id	700
grid_label	10
zstore	514818
dcpp_init_year	60
version	736
derived_variable_id	0

Here, we show the various forcing experiments that climate modellers ran in these simulations.

df['experiment_id'].unique()

array(['highresSST-present', 'piControl', 'control-1950', 'hist-1950',
       'historical', 'amip', 'abrupt-4xCO2', 'abrupt-2xCO2',
       'abrupt-0p5xCO2', '1pctCO2', 'ssp585', 'esm-piControl', 'esm-hist',
       'hist-piAer', 'histSST-1950HC', 'ssp245', 'hist-1950HC', 'histSST',
       'piClim-2xVOC', 'piClim-2xNOx', 'piClim-2xdust', 'piClim-2xss',
       'piClim-histall', 'hist-piNTCF', 'histSST-piNTCF',
       'aqua-control-lwoff', 'piClim-lu', 'histSST-piO3', 'piClim-CH4',
       'piClim-NTCF', 'piClim-NOx', 'piClim-O3', 'piClim-HC',
       'faf-heat-NA0pct', 'ssp370SST-lowCH4', 'piClim-VOC',
       'ssp370-lowNTCF', 'piClim-control', 'piClim-aer', 'hist-aer',
       'faf-heat', 'faf-heat-NA50pct', 'ssp370SST-lowNTCF',
       'ssp370SST-ssp126Lu', 'ssp370SST', 'ssp370pdSST', 'histSST-piAer',
       'piClim-ghg', 'piClim-anthro', 'faf-all', 'hist-nat', 'hist-GHG',
       'ssp119', 'piClim-histnat', 'piClim-4xCO2', 'ssp370',
       'piClim-histghg', 'highresSST-future', 'esm-ssp585-ssp126Lu',
       'ssp126-ssp370Lu', 'ssp370-ssp126Lu', 'land-noLu', 'histSST-piCH4',
       'ssp126', 'esm-pi-CO2pulse', 'amip-hist', 'piClim-histaer',
       'amip-4xCO2', 'faf-water', 'faf-passiveheat', '1pctCO2-rad',
       'faf-stress', '1pctCO2-bgc', 'aqua-control', 'amip-future4K',
       'amip-p4K', 'aqua-p4K', 'amip-lwoff', 'amip-m4K', 'aqua-4xCO2',
       'amip-p4K-lwoff', 'hist-noLu', '1pctCO2-cdr',
       'land-hist-altStartYear', 'land-hist', 'omip1', 'esm-pi-cdr-pulse',
       'esm-ssp585', 'abrupt-solp4p', 'piControl-spinup', 'hist-stratO3',
       'abrupt-solm4p', 'midHolocene', 'lig127k', 'aqua-p4K-lwoff',
       'esm-piControl-spinup', 'ssp245-GHG', 'ssp245-nat',
       'dcppC-amv-neg', 'dcppC-amv-ExTrop-neg', 'dcppC-atl-control',
       'dcppC-amv-pos', 'dcppC-ipv-NexTrop-neg', 'dcppC-ipv-NexTrop-pos',
       'dcppC-atl-pacemaker', 'dcppC-amv-ExTrop-pos',
       'dcppC-amv-Trop-neg', 'dcppC-pac-control', 'dcppC-ipv-pos',
       'dcppC-pac-pacemaker', 'dcppC-ipv-neg', 'dcppC-amv-Trop-pos',
       'piClim-BC', 'piClim-2xfire', 'piClim-SO2', 'piClim-OC',
       'piClim-N2O', 'piClim-2xDMS', 'ssp460', 'ssp434', 'ssp534-over',
       'deforest-globe', 'historical-cmip5', 'hist-bgc',
       'piControl-cmip5', 'rcp26-cmip5', 'rcp45-cmip5', 'rcp85-cmip5',
       'pdSST-piArcSIC', 'pdSST-piAntSIC', 'piSST-piSIC', 'piSST-pdSIC',
       'ssp245-stratO3', 'hist-sol', 'hist-CO2', 'hist-volc',
       'hist-totalO3', 'hist-nat-cmip5', 'hist-aer-cmip5',
       'hist-GHG-cmip5', 'pdSST-futAntSIC', 'futSST-pdSIC', 'pdSST-pdSIC',
       'ssp245-aer', 'pdSST-futArcSIC', 'dcppA-hindcast', 'dcppA-assim',
       'dcppC-hindcast-noPinatubo', 'dcppC-hindcast-noElChichon',
       'dcppC-hindcast-noAgung', 'hist-resIPO', 'ssp245-cov-modgreen',
       'ssp245-cov-fossil', 'ssp245-cov-strgreen', 'ssp245-covid', 'lgm',
       'ssp585-bgc', '1pctCO2to4x-withism', '1pctCO2-4xext', 'past1000',
       'pa-futArcSIC', 'pa-pdSIC', 'historical-ext', 'pdSST-futArcSICSIT',
       'pdSST-futOkhotskSIC', 'pdSST-futBKSeasSIC', 'pa-piArcSIC',
       'pa-piAntSIC', 'pa-futAntSIC', 'pdSST-pdSICSIT'], dtype=object)

Loading Data

Intake-ESM enables loading data directly into an xarray.DataArray, a metadata-aware extension of numpy arrays. Xarray objects leverage Dask to only read data into memory as needed for any specific operation (i.e. lazy evaluation). Think of Xarray Datasets as ways of conveniently organizing large arrays of floating point numbers (e.g. climate model data) on an n-dimensional discrete grid, with important metadata such as units, variable, names, etc.

Note that data on the cloud are in Zarr format, an extension of the metadata-aware format NetCDF commonly used in the geosciences.

Intake-ESM has rules for aggegating datasets; these rules are defined in the collection-specification file.

Here, we choose the piControl experiment (in which CO\(_2\) concentrations are held fixed at a pre-industrial level of ~300 ppm) and abrupt-2xCO2 experiment (in which CO\(_2\) concentrations are instantaneously doubled from a pre-industrial control state). Since the radiative forcing of CO\(_2\) is roughly a logarithmic function of CO\(_2\) concentrations, the ECS is roughly independent of the initial CO\(_2\) concentration.

Warning

The version of this notebook in the Pangeo Gallery uses the abrupt-4xCO2 forcing experiment, but fewer abrupt-2xCO2 datasets are currently avaiable in Google Cloud Storage, which significantly reduces run time. If you want to run this notebook on your own computer with the abrupt-4xCO2 experiment instead, change co2_option in the cell below. You will also need to take half of ecs, as described in the Overview.

co2_option = 'abrupt-2xCO2'

Prepare Data

query = dict(
    experiment_id=[co2_option,'piControl'], # pick the `abrupt-2xCO2` and `piControl` forcing experiments
    table_id='Amon',                            # choose to look at atmospheric variables (A) saved at monthly resolution (mon)
    variable_id=['tas', 'rsut','rsdt','rlut'],  # choose to look at near-surface air temperature (tas) as our variable
    member_id = 'r1i1p1f1',                     # arbitrarily pick one realization for each model (i.e. just one set of initial conditions)
)

col_subset = col.search(require_all_on=["source_id"], **query)
col_subset.df.groupby("source_id")[
    ["experiment_id", "variable_id", "table_id"]
].nunique()

	experiment_id	variable_id	table_id
source_id
CESM2	2	4	1
GISS-E2-1-G	2	4	1
GISS-E2-1-H	2	4	1
GISS-E2-2-G	2	4	1
IPSL-CM6A-LR	2	4	1
MIROC6	2	4	1
MRI-ESM2-0	2	4	1

The following functions help us load and homogenize the data. We use some dask.delayed programming to open the datasets in parallel.

def drop_all_bounds(ds):
    """Drop coordinates like 'time_bounds' from datasets,
    which can lead to issues when merging."""
    drop_vars = [vname for vname in ds.coords
                 if (('_bounds') in vname ) or ('_bnds') in vname]
    return ds.drop_vars(drop_vars)

def open_dsets(df):
    """Open datasets from cloud storage and return xarray dataset."""
    dsets = [xr.open_zarr(fsspec.get_mapper(ds_url), consolidated=True)
             .pipe(drop_all_bounds)
             for ds_url in df.zstore]
    try:
        ds = xr.merge(dsets, join='exact')
        return ds
    except ValueError:
        return None

def open_delayed(df):
    """A dask.delayed wrapper around `open_dsets`.
    Allows us to open many datasets in parallel."""
    return dask.delayed(open_dsets)(df)

Create a nested dictionary of models and experiments. It will be structured like this:

{'CESM2':
  {
    'piControl': <xarray.Dataset>,
    'abrupt-2xCO2': <xarray.Dataset>
  },
  ...
}

from collections import defaultdict

dsets = defaultdict(dict)
for group, df in col_subset.df.groupby(by=['source_id', 'experiment_id']):
    dsets[group[0]][group[1]] = open_delayed(df)

Open one of the datasets directly, just to show what it looks like:

%time open_dsets(df)

CPU times: user 1.51 s, sys: 206 ms, total: 1.71 s
Wall time: 17.9 s

<xarray.Dataset> Size: 7GB
Dimensions:  (lat: 160, lon: 320, time: 8412)
Coordinates:
    height   float64 8B ...
  * lat      (lat) float64 1kB -89.14 -88.03 -86.91 -85.79 ... 86.91 88.03 89.14
  * lon      (lon) float64 3kB 0.0 1.125 2.25 3.375 ... 355.5 356.6 357.8 358.9
  * time     (time) object 67kB 1850-01-16 00:00:00 ... 2550-12-16 00:00:00
Data variables:
    tas      (time, lat, lon) float32 2GB dask.array<chunksize=(600, 160, 320), meta=np.ndarray>
    rsut     (time, lat, lon) float32 2GB dask.array<chunksize=(600, 160, 320), meta=np.ndarray>
    rsdt     (time, lat, lon) float32 2GB dask.array<chunksize=(600, 160, 320), meta=np.ndarray>
    rlut     (time, lat, lon) float32 2GB dask.array<chunksize=(600, 160, 320), meta=np.ndarray>
Attributes: (12/47)
    Conventions:            CF-1.7 CMIP-6.2
    activity_id:            CMIP
    branch_method:          standard
    branch_time_in_child:   0.0
    branch_time_in_parent:  365243.0
    cmor_version:           3.4.0
    ...                     ...
    tracking_id:            hdl:21.14100/a0bd7d8f-4785-4687-b510-fd87808051b6
    variable_id:            tas
    variant_label:          r1i1p1f1
    status:                 2019-10-25;created;by nhn2@columbia.edu
    netcdf_tracking_ids:    hdl:21.14100/a0bd7d8f-4785-4687-b510-fd87808051b6
    version_id:             v20190222

xarray.Dataset

• Dimensions:
  • lat: 160
  • lon: 320
  • time: 8412
• Coordinates: (4)
  • height
    ()
    float64
    ...

    axis :

    Z

    long_name :

    height

    positive :

    up

    standard_name :

    height

    units :

    m

    [1 values with dtype=float64]

  • lat
    (lat)
    float64
    -89.14 -88.03 ... 88.03 89.14

    axis :

    Y

    bounds :

    lat_bnds

    long_name :

    Latitude

    standard_name :

    latitude

    units :

    degrees_north

    array([-89.14152, -88.02943, -86.91077, -85.79063, -84.66992, -83.54895,
           -82.42782, -81.30659, -80.18531, -79.06398, -77.94262, -76.82124,
           -75.69984, -74.57843, -73.45701, -72.33558, -71.21414, -70.09269,
           -68.97124, -67.84978, -66.72833, -65.60686, -64.4854 , -63.36393,
           -62.24246, -61.12099, -59.99952, -58.87804, -57.75657, -56.63509,
           -55.51361, -54.39214, -53.27066, -52.14917, -51.02769, -49.90621,
           -48.78473, -47.66325, -46.54176, -45.42028, -44.29879, -43.17731,
           -42.05582, -40.93434, -39.81285, -38.69137, -37.56988, -36.44839,
           -35.32691, -34.20542, -33.08393, -31.96244, -30.84096, -29.71947,
           -28.59798, -27.47649, -26.355  , -25.23351, -24.11203, -22.99054,
           -21.86905, -20.74756, -19.62607, -18.50458, -17.38309, -16.2616 ,
           -15.14011, -14.01862, -12.89713, -11.77564, -10.65415,  -9.53266,
            -8.41117,  -7.28968,  -6.16819,  -5.0467 ,  -3.92521,  -2.80372,
            -1.68223,  -0.56074,   0.56074,   1.68223,   2.80372,   3.92521,
             5.0467 ,   6.16819,   7.28968,   8.41117,   9.53266,  10.65415,
            11.77564,  12.89713,  14.01862,  15.14011,  16.2616 ,  17.38309,
            18.50458,  19.62607,  20.74756,  21.86905,  22.99054,  24.11203,
            25.23351,  26.355  ,  27.47649,  28.59798,  29.71947,  30.84096,
            31.96244,  33.08393,  34.20542,  35.32691,  36.44839,  37.56988,
            38.69137,  39.81285,  40.93434,  42.05582,  43.17731,  44.29879,
            45.42028,  46.54176,  47.66325,  48.78473,  49.90621,  51.02769,
            52.14917,  53.27066,  54.39214,  55.51361,  56.63509,  57.75657,
            58.87804,  59.99952,  61.12099,  62.24246,  63.36393,  64.4854 ,
            65.60686,  66.72833,  67.84978,  68.97124,  70.09269,  71.21414,
            72.33558,  73.45701,  74.57843,  75.69984,  76.82124,  77.94262,
            79.06398,  80.18531,  81.30659,  82.42782,  83.54895,  84.66992,
            85.79063,  86.91077,  88.02943,  89.14152])

  • lon
    (lon)
    float64
    0.0 1.125 2.25 ... 357.8 358.9

    axis :

    X

    bounds :

    lon_bnds

    long_name :

    Longitude

    standard_name :

    longitude

    units :

    degrees_east

    array([  0.   ,   1.125,   2.25 , ..., 356.625, 357.75 , 358.875])

  • time
    (time)
    object
    1850-01-16 00:00:00 ... 2550-12-...

    array([cftime.Datetime360Day(1850, 1, 16, 0, 0, 0, 0, has_year_zero=True),
           cftime.Datetime360Day(1850, 2, 16, 0, 0, 0, 0, has_year_zero=True),
           cftime.Datetime360Day(1850, 3, 16, 0, 0, 0, 0, has_year_zero=True), ...,
           cftime.Datetime360Day(2550, 10, 16, 0, 0, 0, 0, has_year_zero=True),
           cftime.Datetime360Day(2550, 11, 16, 0, 0, 0, 0, has_year_zero=True),
           cftime.Datetime360Day(2550, 12, 16, 0, 0, 0, 0, has_year_zero=True)],
          dtype=object)

• Data variables: (4)
  • tas
    (time, lat, lon)
    float32
    dask.array<chunksize=(600, 160, 320), meta=np.ndarray>

    cell_measures :

    area: areacella

    cell_methods :

    area: time: mean

    comment :

    near-surface (usually, 2 meter) air temperature

    history :

    2019-02-20T08:13:34Z altered by CMOR: Treated scalar dimension: 'height'. 2019-02-20T08:13:34Z altered by CMOR: replaced missing value flag (-9.99e+33) with standard missing value (1e+20).

    long_name :

    Near-Surface Air Temperature

    original_name :

    TA

    standard_name :

    air_temperature

    units :

    K

    Array Chunk Bytes 1.60 GiB 117.19 MiB Shape (8412, 160, 320) (600, 160, 320) Dask graph 15 chunks in 2 graph layers Data type float32 numpy.ndarray

  • rsut
    (time, lat, lon)
    float32
    dask.array<chunksize=(600, 160, 320), meta=np.ndarray>

    cell_measures :

    area: areacella

    cell_methods :

    area: time: mean

    comment :

    at the top of the atmosphere

    history :

    2019-02-20T09:29:39Z altered by CMOR: replaced missing value flag (-9.99e+33) with standard missing value (1e+20).

    long_name :

    TOA Outgoing Shortwave Radiation

    original_name :

    USWT

    standard_name :

    toa_outgoing_shortwave_flux

    units :

    W m-2

    Array Chunk Bytes 1.60 GiB 117.19 MiB Shape (8412, 160, 320) (600, 160, 320) Dask graph 15 chunks in 2 graph layers Data type float32 numpy.ndarray

  • rsdt
    (time, lat, lon)
    float32
    dask.array<chunksize=(600, 160, 320), meta=np.ndarray>

    cell_measures :

    area: areacella

    cell_methods :

    area: time: mean

    comment :

    Shortwave radiation incident at the top of the atmosphere

    history :

    2019-02-20T09:27:14Z altered by CMOR: replaced missing value flag (-9.99e+33) with standard missing value (1e+20).

    long_name :

    TOA Incident Shortwave Radiation

    original_name :

    DSWT

    standard_name :

    toa_incoming_shortwave_flux

    units :

    W m-2

    Array Chunk Bytes 1.60 GiB 117.19 MiB Shape (8412, 160, 320) (600, 160, 320) Dask graph 15 chunks in 2 graph layers Data type float32 numpy.ndarray

  • rlut
    (time, lat, lon)
    float32
    dask.array<chunksize=(600, 160, 320), meta=np.ndarray>

    cell_measures :

    area: areacella

    cell_methods :

    area: time: mean

    comment :

    at the top of the atmosphere (to be compared with satellite measurements)

    history :

    2019-02-20T09:32:16Z altered by CMOR: replaced missing value flag (-9.99e+33) with standard missing value (1e+20).

    long_name :

    TOA Outgoing Longwave Radiation

    original_name :

    ULWT

    standard_name :

    toa_outgoing_longwave_flux

    units :

    W m-2

    Array Chunk Bytes 1.60 GiB 117.19 MiB Shape (8412, 160, 320) (600, 160, 320) Dask graph 15 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (3)
  • lat
    PandasIndex

    PandasIndex(Index([-89.14152, -88.02943, -86.91077, -85.79063, -84.66992, -83.54895,
           -82.42782, -81.30659, -80.18531, -79.06398,
           ...
            79.06398,  80.18531,  81.30659,  82.42782,  83.54895,  84.66992,
            85.79063,  86.91077,  88.02943,  89.14152],
          dtype='float64', name='lat', length=160))

  • lon
    PandasIndex

    PandasIndex(Index([    0.0,   1.125,    2.25,   3.375,     4.5,   5.625,    6.75,   7.875,
               9.0,  10.125,
           ...
            348.75, 349.875,   351.0, 352.125,  353.25, 354.375,   355.5, 356.625,
            357.75, 358.875],
          dtype='float64', name='lon', length=320))

  • time
    PandasIndex

    PandasIndex(CFTimeIndex([1850-01-16 00:00:00, 1850-02-16 00:00:00, 1850-03-16 00:00:00,
                 1850-04-16 00:00:00, 1850-05-16 00:00:00, 1850-06-16 00:00:00,
                 1850-07-16 00:00:00, 1850-08-16 00:00:00, 1850-09-16 00:00:00,
                 1850-10-16 00:00:00,
                 ...
                 2550-03-16 00:00:00, 2550-04-16 00:00:00, 2550-05-16 00:00:00,
                 2550-06-16 00:00:00, 2550-07-16 00:00:00, 2550-08-16 00:00:00,
                 2550-09-16 00:00:00, 2550-10-16 00:00:00, 2550-11-16 00:00:00,
                 2550-12-16 00:00:00],
                dtype='object', length=8412, calendar='360_day', freq='30D'))

• Attributes: (47)

  Conventions :

  CF-1.7 CMIP-6.2

  activity_id :

  CMIP

  branch_method :

  standard

  branch_time_in_child :

  0.0

  branch_time_in_parent :

  365243.0

  cmor_version :

  3.4.0

  creation_date :

  2019-02-20T08:13:34Z

  data_specs_version :

  01.00.29

  experiment :

  pre-industrial control

  experiment_id :

  piControl

  external_variables :

  areacella

  forcing_index :

  1

  frequency :

  mon

  further_info_url :

  https://furtherinfo.es-doc.org/CMIP6.MRI.MRI-ESM2-0.piControl.none.r1i1p1f1

  grid :

  native atmosphere TL159 gaussian grid (160x320 latxlon)

  grid_label :

  gn

  history :

  2019-02-20T08:13:34Z ; CMOR rewrote data to be consistent with CMIP6, CF-1.7 CMIP-6.2 and CF standards.; Output from run-Dr060_piControl_100 (sfc_avr_mon.ctl)

  initialization_index :

  1

  institution :

  Meteorological Research Institute, Tsukuba, Ibaraki 305-0052, Japan

  institution_id :

  MRI

  license :

  CMIP6 model data produced by MRI is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License (https://creativecommons.org/licenses/). Consult https://pcmdi.llnl.gov/CMIP6/TermsOfUse for terms of use governing CMIP6 output, including citation requirements and proper acknowledgment. Further information about this data, including some limitations, can be found via the further_info_url (recorded as a global attribute in this file). The data producers and data providers make no warranty, either express or implied, including, but not limited to, warranties of merchantability and fitness for a particular purpose. All liabilities arising from the supply of the information (including any liability arising in negligence) are excluded to the fullest extent permitted by law.

  mip_era :

  CMIP6

  nominal_resolution :

  100 km

  parent_activity_id :

  CMIP

  parent_experiment_id :

  piControl-spinup

  parent_mip_era :

  CMIP6

  parent_source_id :

  MRI-ESM2-0

  parent_time_units :

  days since 1850-01-01

  parent_variant_label :

  r1i1p1f1

  physics_index :

  1

  product :

  model-output

  realization_index :

  1

  realm :

  atmos

  source :

  MRI-ESM2.0 (2017): aerosol: MASINGAR mk2r4 (TL95; 192 x 96 longitude/latitude; 80 levels; top level 0.01 hPa) atmos: MRI-AGCM3.5 (TL159; 320 x 160 longitude/latitude; 80 levels; top level 0.01 hPa) atmosChem: MRI-CCM2.1 (T42; 128 x 64 longitude/latitude; 80 levels; top level 0.01 hPa) land: HAL 1.0 landIce: none ocean: MRI.COM4.4 (tripolar primarily 0.5 deg latitude/1 deg longitude with meridional refinement down to 0.3 deg within 10 degrees north and south of the equator; 360 x 364 longitude/latitude; 61 levels; top grid cell 0-2 m) ocnBgchem: MRI.COM4.4 seaIce: MRI.COM4.4

  source_id :

  MRI-ESM2-0

  source_type :

  AOGCM AER CHEM

  sub_experiment :

  none

  sub_experiment_id :

  none

  table_id :

  Amon

  table_info :

  Creation Date:(14 December 2018) MD5:b2d32d1a0d9b196411429c8895329d8f

  title :

  MRI-ESM2-0 output prepared for CMIP6

  tracking_id :

  hdl:21.14100/a0bd7d8f-4785-4687-b510-fd87808051b6

  variable_id :

  tas

  variant_label :

  r1i1p1f1

  status :

  2019-10-25;created;by nhn2@columbia.edu

  netcdf_tracking_ids :

  hdl:21.14100/a0bd7d8f-4785-4687-b510-fd87808051b6

  version_id :

  v20190222

Now use dask to do this in parallel on all of the datasets:

dsets_ = dask.compute(dict(dsets))[0]

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/coding/times.py:992: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/core/indexing.py:526: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(self.get_duck_array(), dtype=dtype)
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/coding/times.py:992: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/core/indexing.py:526: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(self.get_duck_array(), dtype=dtype)
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/coding/times.py:992: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/core/indexing.py:526: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(self.get_duck_array(), dtype=dtype)
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/coding/times.py:992: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/core/indexing.py:526: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(self.get_duck_array(), dtype=dtype)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/coding/times.py:992: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/core/indexing.py:526: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(self.get_duck_array(), dtype=dtype)
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/coding/times.py:992: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/core/indexing.py:526: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(self.get_duck_array(), dtype=dtype)
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/coding/times.py:992: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/core/indexing.py:526: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(self.get_duck_array(), dtype=dtype)
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/coding/times.py:992: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/coding/times.py:992: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/core/indexing.py:526: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(self.get_duck_array(), dtype=dtype)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/coding/times.py:992: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/core/indexing.py:526: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(self.get_duck_array(), dtype=dtype)
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/core/indexing.py:526: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(self.get_duck_array(), dtype=dtype)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/coding/times.py:992: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/core/indexing.py:526: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(self.get_duck_array(), dtype=dtype)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/coding/times.py:992: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)
/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/xarray/core/indexing.py:526: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  return np.asarray(self.get_duck_array(), dtype=dtype)

Reduce Data via Global Mean

We don’t want to load all of the raw model data into memory right away. Instead, we want to reduce the data by taking the global mean. We need to remember to weight this global mean by a factor proportional to cos(lat).

def get_lat_name(ds):
    """Figure out what is the latitude coordinate for each dataset."""
    for lat_name in ['lat', 'latitude']:
        if lat_name in ds.coords:
            return lat_name
    raise RuntimeError("Couldn't find a latitude coordinate")

def global_mean(ds):
    """Return global mean of a whole dataset."""
    lat = ds[get_lat_name(ds)]
    weight = np.cos(np.deg2rad(lat))
    weight /= weight.mean()
    other_dims = set(ds.dims) - {'time'}
    return (ds * weight).mean(other_dims)

We now apply this function, plus resampling to annual mean data, to all of the datasets. We also concatenate the experiments together into a single Dataset for each model. This is the most complex cell in the notebook. A lot is happening here.

expts = ['piControl', co2_option]
expt_da = xr.DataArray(expts, dims='experiment_id',
                       coords={'experiment_id': expts})

dsets_aligned = {}

for k, v in tqdm(dsets_.items()):
    expt_dsets = v.values()
    if any([d is None for d in expt_dsets]):
        print(f"Missing experiment for {k}")
        continue

    for ds in expt_dsets:
        ds.coords['year'] = ds.time.dt.year - ds.time.dt.year[0]

    # workaround for
    # https://github.com/pydata/xarray/issues/2237#issuecomment-620961663
    dsets_ann_mean = [v[expt].pipe(global_mean)
                             .swap_dims({'time': 'year'})
                             .drop_vars('time')
                             .coarsen(year=12).mean()
                      for expt in expts]

    # align everything with the 2xCO2 experiment
    dsets_aligned[k] = xr.concat(dsets_ann_mean, join='right',
                                 dim=expt_da)

Do the Computation

Up to this point, no computations have actually happened. Everything has been “lazy”. Now we trigger the computation to actual occur and load the global/annual mean data into memory.

dsets_aligned_ = dask.compute(dsets_aligned)[0]

Now we concatenate across models to produce one big dataset with all the required variables.

source_ids = list(dsets_aligned_.keys())
source_da = xr.DataArray(source_ids, dims='source_id',
                         coords={'source_id': source_ids})

big_ds = xr.concat([ds.reset_coords(drop=True)
                    for ds in dsets_aligned_.values()],
                   dim=source_da)
big_ds

<xarray.Dataset> Size: 114kB
Dimensions:        (source_id: 7, experiment_id: 2, year: 250)
Coordinates:
  * year           (year) float64 2kB 0.0 1.0 2.0 3.0 ... 247.0 248.0 249.0
  * experiment_id  (experiment_id) <U12 96B 'piControl' 'abrupt-2xCO2'
  * source_id      (source_id) <U12 336B 'CESM2' 'GISS-E2-1-G' ... 'MRI-ESM2-0'
Data variables:
    rsut           (source_id, experiment_id, year) float64 28kB 99.36 ... nan
    tas            (source_id, experiment_id, year) float64 28kB 287.0 ... nan
    rsdt           (source_id, experiment_id, year) float64 28kB 340.3 ... nan
    rlut           (source_id, experiment_id, year) float64 28kB 240.2 ... nan

xarray.Dataset

• Dimensions:
  • source_id: 7
  • experiment_id: 2
  • year: 250
• Coordinates: (3)
  • year
    (year)
    float64
    0.0 1.0 2.0 ... 247.0 248.0 249.0

    array([  0.,   1.,   2., ..., 247., 248., 249.])

  • experiment_id
    (experiment_id)
    <U12
    'piControl' 'abrupt-2xCO2'

    array(['piControl', 'abrupt-2xCO2'], dtype='<U12')

  • source_id
    (source_id)
    <U12
    'CESM2' ... 'MRI-ESM2-0'

    array(['CESM2', 'GISS-E2-1-G', 'GISS-E2-1-H', 'GISS-E2-2-G', 'IPSL-CM6A-LR',
           'MIROC6', 'MRI-ESM2-0'], dtype='<U12')

• Data variables: (4)
  • rsut
    (source_id, experiment_id, year)
    float64
    99.36 99.8 98.94 ... nan nan nan

    array([[[ 99.36393586,  99.80283632,  98.94191234, ...,          nan,
                      nan,          nan],
            [ 97.82037343,  97.79000894,  97.01762283, ...,          nan,
                      nan,          nan]],
    
           [[101.4390611 , 101.56296953, 101.4972957 , ...,          nan,
                      nan,          nan],
            [101.0366182 , 101.00232906, 100.84826946, ...,          nan,
                      nan,          nan]],
    
           [[100.66435273, 100.51116956, 100.89353578, ...,          nan,
                      nan,          nan],
            [100.60859904, 100.53582333, 100.60002885, ...,          nan,
                      nan,          nan]],
    
           ...,
    
           [[102.75021362, 103.10281372, 102.87644196, ...,          nan,
                      nan,          nan],
            [100.99295807, 101.50688934, 101.97058105, ...,          nan,
                      nan,          nan]],
    
           [[108.8646907 , 108.49889761, 108.90849884, ..., 108.70599882,
             108.96811199, 108.72818105],
            [108.12545905, 108.27118498, 107.91076284, ..., 107.34994386,
             107.12170546, 107.09244583]],
    
           [[ 98.00645259,  97.99298826,  98.13119018, ...,          nan,
                      nan,          nan],
            [ 96.58950119,  96.63950599,  96.21623569, ...,          nan,
                      nan,          nan]]])

  • tas
    (source_id, experiment_id, year)
    float64
    287.0 287.1 287.0 ... nan nan nan

    array([[[286.9850739 , 287.06653368, 286.98444303, ...,          nan,
                      nan,          nan],
            [287.48014449, 287.7394101 , 288.01731092, ...,          nan,
                      nan,          nan]],
    
           [[286.9343543 , 286.9199378 , 286.92800874, ...,          nan,
                      nan,          nan],
            [287.47643137, 287.64304425, 287.92338648, ...,          nan,
                      nan,          nan]],
    
           [[287.23218778, 287.35340132, 287.36071829, ...,          nan,
                      nan,          nan],
            [287.82533396, 288.24506914, 288.49300187, ...,          nan,
                      nan,          nan]],
    
           ...,
    
           [[286.09158325, 286.10717773, 286.03372192, ...,          nan,
                      nan,          nan],
            [286.59127808, 287.17440796, 287.36087036, ...,          nan,
                      nan,          nan]],
    
           [[288.31629875, 288.45178258, 288.4777738 , ..., 288.45273041,
             288.37086858, 288.48120669],
            [288.8184018 , 289.08524988, 289.31636523, ..., 290.57012689,
             290.34771036, 290.18482772]],
    
           [[286.90136665, 286.87541207, 286.83983702, ...,          nan,
                      nan,          nan],
            [287.44460171, 287.60650445, 287.87638255, ...,          nan,
                      nan,          nan]]])

  • rsdt
    (source_id, experiment_id, year)
    float64
    340.3 340.3 340.3 ... nan nan nan

    array([[[340.25427344, 340.25427344, 340.25427344, ...,          nan,
                      nan,          nan],
            [340.25427344, 340.25427344, 340.25427344, ...,          nan,
                      nan,          nan]],
    
           [[340.31299703, 340.31299703, 340.31299703, ...,          nan,
                      nan,          nan],
            [340.31299573, 340.31299573, 340.31299573, ...,          nan,
                      nan,          nan]],
    
           [[340.31299573, 340.31299573, 340.31299573, ...,          nan,
                      nan,          nan],
            [340.31299573, 340.31299573, 340.31299573, ...,          nan,
                      nan,          nan]],
    
           ...,
    
           [[341.65182495, 341.65182495, 341.66079712, ...,          nan,
                      nan,          nan],
            [341.65182495, 341.65182495, 341.66079712, ...,          nan,
                      nan,          nan]],
    
           [[340.22406601, 340.24657384, 340.24656962, ..., 340.24657516,
             340.2240696 , 340.24657389],
            [340.22406187, 340.24658263, 340.24657275, ..., 340.2465728 ,
             340.22406728, 340.24657341]],
    
           [[340.254947  , 340.254947  , 340.2640629 , ...,          nan,
                      nan,          nan],
            [340.254947  , 340.254947  , 340.2640629 , ...,          nan,
                      nan,          nan]]])

  • rlut
    (source_id, experiment_id, year)
    float64
    240.2 240.5 240.1 ... nan nan nan

    array([[[240.191604  , 240.48018287, 240.14496437, ...,          nan,
                      nan,          nan],
            [237.65958533, 238.36057267, 238.90966217, ...,          nan,
                      nan,          nan]],
    
           [[238.63641591, 238.79870786, 238.71907378, ...,          nan,
                      nan,          nan],
            [235.99094115, 236.40002371, 236.7402137 , ...,          nan,
                      nan,          nan]],
    
           [[239.20591515, 239.51517967, 239.280247  , ...,          nan,
                      nan,          nan],
            [236.23198909, 237.14592423, 237.39700028, ...,          nan,
                      nan,          nan]],
    
           ...,
    
           [[238.27838135, 238.13494873, 238.2039032 , ...,          nan,
                      nan,          nan],
            [236.03692627, 236.82038879, 237.15605164, ...,          nan,
                      nan,          nan]],
    
           [[230.1649978 , 230.32354978, 230.40610909, ..., 230.28280645,
             230.28065609, 230.35605828],
            [228.11024163, 228.6672724 , 229.18270831, ..., 231.98662039,
             231.32286213, 231.08566087]],
    
           [[241.70191145, 241.46380205, 241.5131782 , ...,          nan,
                      nan,          nan],
            [239.70625678, 240.17113091, 240.5943591 , ...,          nan,
                      nan,          nan]]])

• Indexes: (3)
  • year
    PandasIndex

    PandasIndex(Index([  0.0,   1.0,   2.0,   3.0,   4.0,   5.0,   6.0,   7.0,   8.0,   9.0,
           ...
           240.0, 241.0, 242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 249.0],
          dtype='float64', name='year', length=250))

  • experiment_id
    PandasIndex

    PandasIndex(Index(['piControl', 'abrupt-2xCO2'], dtype='object', name='experiment_id'))

  • source_id
    PandasIndex

    PandasIndex(Index(['CESM2', 'GISS-E2-1-G', 'GISS-E2-1-H', 'GISS-E2-2-G', 'IPSL-CM6A-LR',
           'MIROC6', 'MRI-ESM2-0'],
          dtype='object', name='source_id'))

• Attributes: (0)

Calculated Derived Variables

We need to calculate the net radiative imbalance, plus the anomaly of the abrupt 2xCO2 run compared to the piControl run.

big_ds['imbalance'] = big_ds['rsdt'] - big_ds['rsut'] - big_ds['rlut']

ds_mean = big_ds[['tas', 'imbalance']].sel(experiment_id='piControl').mean(dim='year')
ds_anom = big_ds[['tas', 'imbalance']] - ds_mean

# add some metadata
ds_anom.tas.attrs['long_name'] = 'Global Mean Surface Temp Anom'
ds_anom.tas.attrs['units'] = 'K'
ds_anom.imbalance.attrs['long_name'] = 'Global Mean Radiative Imbalance'
ds_anom.imbalance.attrs['units'] = 'W m$^{-2}$'

ds_anom

<xarray.Dataset> Size: 58kB
Dimensions:        (year: 250, experiment_id: 2, source_id: 7)
Coordinates:
  * year           (year) float64 2kB 0.0 1.0 2.0 3.0 ... 247.0 248.0 249.0
  * experiment_id  (experiment_id) <U12 96B 'piControl' 'abrupt-2xCO2'
  * source_id      (source_id) <U12 336B 'CESM2' 'GISS-E2-1-G' ... 'MRI-ESM2-0'
Data variables:
    tas            (source_id, experiment_id, year) float64 28kB 0.06312 ... nan
    imbalance      (source_id, experiment_id, year) float64 28kB -0.04953 ......

xarray.Dataset

• Dimensions:
  • year: 250
  • experiment_id: 2
  • source_id: 7
• Coordinates: (3)
  • year
    (year)
    float64
    0.0 1.0 2.0 ... 247.0 248.0 249.0

    array([  0.,   1.,   2., ..., 247., 248., 249.])

  • experiment_id
    (experiment_id)
    <U12
    'piControl' 'abrupt-2xCO2'

    array(['piControl', 'abrupt-2xCO2'], dtype='<U12')

  • source_id
    (source_id)
    <U12
    'CESM2' ... 'MRI-ESM2-0'

    array(['CESM2', 'GISS-E2-1-G', 'GISS-E2-1-H', 'GISS-E2-2-G', 'IPSL-CM6A-LR',
           'MIROC6', 'MRI-ESM2-0'], dtype='<U12')

• Data variables: (2)
  • tas
    (source_id, experiment_id, year)
    float64
    0.06312 0.1446 0.06249 ... nan nan

    long_name :

    Global Mean Surface Temp Anom

    units :

    K

    array([[[ 0.06312165,  0.14458142,  0.06249077, ...,         nan,
                     nan,         nan],
            [ 0.55819223,  0.81745784,  1.09535866, ...,         nan,
                     nan,         nan]],
    
           [[ 0.02001347,  0.00559697,  0.01366791, ...,         nan,
                     nan,         nan],
            [ 0.56209054,  0.72870342,  1.00904565, ...,         nan,
                     nan,         nan]],
    
           [[-0.04873387,  0.07247968,  0.07979664, ...,         nan,
                     nan,         nan],
            [ 0.54441231,  0.96414749,  1.21208022, ...,         nan,
                     nan,         nan]],
    
           ...,
    
           [[ 0.19686259,  0.21245707,  0.13900126, ...,         nan,
                     nan,         nan],
            [ 0.69655741,  1.2796873 ,  1.4661497 , ...,         nan,
                     nan,         nan]],
    
           [[-0.03490859,  0.10057523,  0.12656646, ...,  0.10152306,
              0.01966124,  0.12999935],
            [ 0.46719446,  0.73404254,  0.96515788, ...,  2.21891955,
              1.99650302,  1.83362038]],
    
           [[-0.0086573 , -0.03461188, -0.07018693, ...,         nan,
                     nan,         nan],
            [ 0.53457776,  0.6964805 ,  0.9663586 , ...,         nan,
                     nan,         nan]]])

  • imbalance
    (source_id, experiment_id, year)
    float64
    -0.04953 -0.777 0.4191 ... nan nan

    long_name :

    Global Mean Radiative Imbalance

    units :

    W m$^{-2}$

    array([[[-0.04953143, -0.77701077,  0.41913172, ...,         nan,
                     nan,         nan],
            [ 4.02604967,  3.35542682,  3.57872343, ...,         nan,
                     nan,         nan]],
    
           [[ 0.12200747, -0.16419291, -0.018885  , ...,         nan,
                     nan,         nan],
            [ 3.16992385,  2.79513042,  2.60900003, ...,         nan,
                     nan,         nan]],
    
           [[ 0.36922609,  0.21314474,  0.06571119, ...,         nan,
                     nan,         nan],
            [ 3.39890584,  2.55774641,  2.24246484, ...,         nan,
                     nan,         nan]],
    
           ...,
    
           [[-0.091017  , -0.30018448, -0.13379501, ...,         nan,
                     nan,         nan],
            [ 3.90769363,  2.61029984,  1.81991745, ...,         nan,
                     nan,         nan]],
    
           [[ 0.05491239,  0.28466133, -0.20750343, ...,  0.11830477,
             -0.1641636 ,  0.02286943],
            [ 2.84889607,  2.16866012,  2.01363648, ..., -0.22945657,
              0.64003458,  0.92900159]],
    
           [[-0.41629225, -0.16471851, -0.34318069, ...,         nan,
                     nan,         nan],
            [ 2.99631382,  2.48143489,  2.49059289, ...,         nan,
                     nan,         nan]]])

• Indexes: (3)
  • year
    PandasIndex

    PandasIndex(Index([  0.0,   1.0,   2.0,   3.0,   4.0,   5.0,   6.0,   7.0,   8.0,   9.0,
           ...
           240.0, 241.0, 242.0, 243.0, 244.0, 245.0, 246.0, 247.0, 248.0, 249.0],
          dtype='float64', name='year', length=250))

  • experiment_id
    PandasIndex

    PandasIndex(Index(['piControl', 'abrupt-2xCO2'], dtype='object', name='experiment_id'))

  • source_id
    PandasIndex

    PandasIndex(Index(['CESM2', 'GISS-E2-1-G', 'GISS-E2-1-H', 'GISS-E2-2-G', 'IPSL-CM6A-LR',
           'MIROC6', 'MRI-ESM2-0'],
          dtype='object', name='source_id'))

• Attributes: (0)

Plot Timeseries

Here we plot the global mean surface temperature for each model:

ds_anom.tas.plot.line(col='source_id', x='year', col_wrap=5)

<xarray.plot.facetgrid.FacetGrid at 0x7f2fbc87a890>

We can see that the models cover different time intervals. Let’s limit the rest of our analysis to the first 150 years.

first_150_years = slice(0, 149)
ds_anom.tas.sel(year=first_150_years).plot.line(col='source_id', x='year', col_wrap=5)

<xarray.plot.facetgrid.FacetGrid at 0x7f2fbc30fb20>

Same thing for radiative imbalance:

ds_anom.imbalance.sel(year=first_150_years).plot.line(col='source_id', x='year', col_wrap=5)

<xarray.plot.facetgrid.FacetGrid at 0x7f2fbc30fe20>

Calculate ECS

ds_abrupt = ds_anom.sel(year=first_150_years, experiment_id=co2_option).reset_coords(drop=True)

def calc_ecs(tas, imb):
    a, b = np.polyfit(tas, imb, 1)
    ecs = -1.0 * (b/a) # Change -1.0 to -0.5 if using 4xCO2
    return xr.DataArray(ecs)

ds_abrupt['ecs'] = xr.apply_ufunc(calc_ecs, ds_abrupt.tas, ds_abrupt.imbalance, vectorize=True, input_core_dims=[['year'], ['year']])
ds_abrupt.compute()

<xarray.Dataset> Size: 18kB
Dimensions:    (year: 150, source_id: 7)
Coordinates:
  * year       (year) float64 1kB 0.0 1.0 2.0 3.0 ... 146.0 147.0 148.0 149.0
  * source_id  (source_id) <U12 336B 'CESM2' 'GISS-E2-1-G' ... 'MRI-ESM2-0'
Data variables:
    tas        (source_id, year) float64 8kB 0.5582 0.8175 1.095 ... 1.646 1.616
    imbalance  (source_id, year) float64 8kB 4.026 3.355 3.579 ... 1.022 0.7031
    ecs        (source_id) float64 56B 3.526 2.659 2.932 2.501 3.827 2.235 2.552

xarray.Dataset

• Dimensions:
  • year: 150
  • source_id: 7
• Coordinates: (2)
  • year
    (year)
    float64
    0.0 1.0 2.0 ... 147.0 148.0 149.0

    array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,  11.,
            12.,  13.,  14.,  15.,  16.,  17.,  18.,  19.,  20.,  21.,  22.,  23.,
            24.,  25.,  26.,  27.,  28.,  29.,  30.,  31.,  32.,  33.,  34.,  35.,
            36.,  37.,  38.,  39.,  40.,  41.,  42.,  43.,  44.,  45.,  46.,  47.,
            48.,  49.,  50.,  51.,  52.,  53.,  54.,  55.,  56.,  57.,  58.,  59.,
            60.,  61.,  62.,  63.,  64.,  65.,  66.,  67.,  68.,  69.,  70.,  71.,
            72.,  73.,  74.,  75.,  76.,  77.,  78.,  79.,  80.,  81.,  82.,  83.,
            84.,  85.,  86.,  87.,  88.,  89.,  90.,  91.,  92.,  93.,  94.,  95.,
            96.,  97.,  98.,  99., 100., 101., 102., 103., 104., 105., 106., 107.,
           108., 109., 110., 111., 112., 113., 114., 115., 116., 117., 118., 119.,
           120., 121., 122., 123., 124., 125., 126., 127., 128., 129., 130., 131.,
           132., 133., 134., 135., 136., 137., 138., 139., 140., 141., 142., 143.,
           144., 145., 146., 147., 148., 149.])

  • source_id
    (source_id)
    <U12
    'CESM2' ... 'MRI-ESM2-0'

    array(['CESM2', 'GISS-E2-1-G', 'GISS-E2-1-H', 'GISS-E2-2-G', 'IPSL-CM6A-LR',
           'MIROC6', 'MRI-ESM2-0'], dtype='<U12')

• Data variables: (3)
  • tas
    (source_id, year)
    float64
    0.5582 0.8175 1.095 ... 1.646 1.616

    long_name :

    Global Mean Surface Temp Anom

    units :

    K

    array([[0.55819223, 0.81745784, 1.09535866, ..., 2.78579546, 2.81815608,
            2.84298753],
           [0.56209054, 0.72870342, 1.00904565, ..., 2.30772005, 2.2496876 ,
            2.14807954],
           [0.54441231, 0.96414749, 1.21208022, ..., 2.63851405, 2.42792187,
            2.34847531],
           ...,
           [0.69655741, 1.2796873 , 1.4661497 , ..., 3.09743632, 3.02495707,
            3.20421733],
           [0.46719446, 0.73404254, 0.96515788, ..., 2.01300088, 1.91896167,
            1.96757883],
           [0.53457776, 0.6964805 , 0.9663586 , ..., 1.89364471, 1.64573962,
            1.61609827]])

  • imbalance
    (source_id, year)
    float64
    4.026 3.355 3.579 ... 1.022 0.7031

    long_name :

    Global Mean Radiative Imbalance

    units :

    W m$^{-2}$

    array([[4.02604967, 3.35542682, 3.57872343, ..., 0.65531467, 1.46648058,
            1.09784156],
           [3.16992385, 2.79513042, 2.60900003, ..., 0.46905707, 0.54642323,
            0.5188314 ],
           [3.39890584, 2.55774641, 2.24246484, ..., 0.16186796, 0.48404896,
            1.12914552],
           ...,
           [3.90769363, 2.61029984, 1.81991745, ..., 0.58088089, 0.90067459,
            1.13007523],
           [2.84889607, 2.16866012, 2.01363648, ..., 1.03411792, 0.76245993,
            0.66630028],
           [2.99631382, 2.48143489, 2.49059289, ..., 0.81439704, 1.02247347,
            0.70310016]])

  • ecs
    (source_id)
    float64
    3.526 2.659 2.932 ... 2.235 2.552

    array([3.52587433, 2.65858194, 2.93219464, 2.50138386, 3.82672137,
           2.23504735, 2.55244043])

• Indexes: (2)
  • year
    PandasIndex

    PandasIndex(Index([  0.0,   1.0,   2.0,   3.0,   4.0,   5.0,   6.0,   7.0,   8.0,   9.0,
           ...
           140.0, 141.0, 142.0, 143.0, 144.0, 145.0, 146.0, 147.0, 148.0, 149.0],
          dtype='float64', name='year', length=150))

  • source_id
    PandasIndex

    PandasIndex(Index(['CESM2', 'GISS-E2-1-G', 'GISS-E2-1-H', 'GISS-E2-2-G', 'IPSL-CM6A-LR',
           'MIROC6', 'MRI-ESM2-0'],
          dtype='object', name='source_id'))

• Attributes: (0)

Also, make sure that we set a couple of the variables to be coordinates.

ds_abrupt = ds_abrupt.set_coords(['tas', 'imbalance'])

fg = ds_abrupt.plot.scatter(x='tas', y='imbalance', col='source_id', col_wrap=4, add_colorbar=False)

def calc_and_plot_ecs(x, y, **kwargs):
    x = x[~np.isnan(x)]
    y = y[~np.isnan(y)]
    a, b = np.polyfit(x, y, 1)
    ecs = -1.0 * b/a
    plt.autoscale(False)
    plt.plot([0, 10], np.polyval([a, b], [0, 10]), 'k')
    plt.text(2, 3, f'ECS = {ecs:3.2f}', fontdict={'weight': 'bold', 'size': 12})
    plt.grid()

fg.map(calc_and_plot_ecs, 'tas', 'imbalance')

<xarray.plot.facetgrid.FacetGrid at 0x7f2fb496b1c0>

ds_abrupt.ecs.plot.hist();

ds_abrupt.ecs.to_dataframe().sort_values('ecs').plot(kind='bar')

<Axes: xlabel='source_id'>

We’re at the end of the notebook, so let’s shutdown our Dask cluster.

client.shutdown()

---

Summary

In this notebook, we estimated ECS for a subset of CMIP6 models using the Gregory method.

What’s next?

We will plot global average surface air temperature for a historical run and two branching emissions scenarios.

Resources and references

• Original notebook in the Pangeo Gallery by Henri Drake and Ryan Abernathey