examples/calibrated_constrained_ensemble.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# A calibrated, constrained ensemble\n",
    "\n",
    "`fair`, like every other simple or complex climate model, is naive. It will produce projections for whatever emissions/concentrations/forcing scenario you ask it to produce projections for. It is up to the user to determine whether these projections are useful and sensible.\n",
    "\n",
    "We have [developed a set of parameter calibrations](https://github.com/chrisroadmap/fair-calibrate) that reproduce both observed climate change since pre-industrial and assessed climate metrics such as the equilibrium climate sensitivity from the IPCC Sixth Assessement Report. These are described in this paper in review:\n",
    "\n",
    "Smith, C., Cummins, D. P., Fredriksen, H.-B., Nicholls, Z., Meinshausen, M., Allen, M., Jenkins, S., Leach, N., Mathison, C., and Partanen, A.-I.: fair-calibrate v1.4.1: calibration, constraining and validation of the FaIR simple climate model for reliable future climate projections, EGUsphere [preprint], https://doi.org/10.5194/egusphere-2024-708, 2024.\n",
    "\n",
    "**Note**: if you are reading this tutorial online and want to reproduce the results, you will need a few additional files. Grab these from https://github.com/OMS-NetZero/FAIR/tree/master/examples/data/calibrated_constrained_ensemble and put them in a new folder relative to this notebook (`./data/calibrated_constrained_ensemble/`). This does not apply if you are running this notebook from Binder or have cloned it from GitHub - it should run out of the box.\n",
    "\n",
    "The calibrations will be continually updated, as new data for surface temperature, ocean heat content, external forcing and emissions become available. For now, we have calibration version where emissions and climate constraints are updated to 2022, and assessments of emergent climate metrics are from the IPCC AR6 WG1 Chapter 7. We use emissions data from a variety of sources (Global Carbon Project, PRIMAP-Hist, CEDS, GFED), and harmonize SSP scenarios to ensure that the projections (which originally started in 2015) have a smooth transition when recent emissions are taken into account.\n",
    "\n",
    "As described in the Smith et al. (2024) paper, a two-step constraining process is produced. The first step ensures that historical simulations match observed climate change to a root-mean-square error of less than 0.17°C. The second step simultaneously distribution-fits to the following assessed ranges:\n",
    "\n",
    "- equilibrium climate sensitivity (ECS), very likely range 2-5°C, best estimate 3°C\n",
    "- transient climate response (TCR), very likely range 1.2-2.4°C, best estimate 1.8°C\n",
    "- global mean surface temperature change 1850-1900 to 2003-2022, very likely range 0.87-1.13°C, best estimate 1.03°C\n",
    "- effective radiative forcing from aerosol-radiation interactions 1750 to 2005-2014, very likely range -0.6 to 0 W/m², best estimate -0.3 W/m²\n",
    "- effective radiative forcing from aerosol-cloud interactions 1750 to 2005-2014, very likely range -1.7 to -0.3 W/m², best estimate -1.0 W/m²\n",
    "- effective radiative forcing from aerosols 1750 to 2005-2014, very likely range -2.0 to -0.6 W/m², best estimate -1.3 W/m²\n",
    "- earth energy uptake change 1971 to 2020, very likely range 358-573 ZJ, best estimate 465 ZJ\n",
    "- CO2 concentrations in 2022, very likely range 416.2-417.8 ppm, best estimate 417.0 ppm\n",
    "\n",
    "841 posterior ensemble members are produced from an initial prior of 1.6 million.\n",
    "\n",
    "There are many, many, many different calibration and constraining possibilities, and it depends on your purposes as to what is appropriate. If you care about the carbon cycle, you might want to constrain on TCRE and ZEC in addition, or instead of, some of the other constraints above. Not all constraints are necessarily internally consistent, and there will be some tradeoff; it is impossible to hit the above ranges perfectly. As more constraints are added, this gets harder, or will require larger prior sample sizes.\n",
    "\n",
    "<a href=\"https://doi.org/10.5281/zenodo.7694879\"><img src=\"https://zenodo.org/badge/DOI/10.5281/zenodo.8399112.svg\" align=\"left\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "## 0. Get required imports\n",
    "\n",
    "[pooch](https://www.fatiando.org/pooch/latest/) is a useful package that allows downloads of external datasets to your cache, meaning that you don't have to include them in Git repositories (particularly troublesome for large files) or `.gitignore` them (difficult for exact reproduciblity, and easy to forget and accidently commit a large file)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as pl\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from fair import FAIR\n",
    "from fair.interface import fill, initialise\n",
    "from fair.io import read_properties"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "## 1. Create FaIR instance\n",
    "\n",
    "To reproduce an AR6-like run, we want to allow methane lifetime to be affected by all its relevant chemical precursors (NOx, VOCs, etc) so we set the `ch4_method` flag to `Thornhill2021` (see https://docs.fairmodel.net/en/latest/api_reference.html#fair.FAIR for all of the options for initialising `FAIR`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = FAIR(ch4_method=\"Thornhill2021\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "## 2. Define time horizon\n",
    "\n",
    "A lot of analysis uses 2100 as the time horizon, but 2300 is an interesting end point to see the effects of long-term climate change. We'll set 2300 as the last time bound."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "f.define_time(1750, 2300, 1)  # start, end, step"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "## 3. Define scenarios\n",
    "\n",
    "The new emissions file reader introduced in `fair` v2.2 makes it easy to define your own emissions files, and name the scenarios how you like. Therefore, you are not limited to using SSPs or any other existing scenarios.\n",
    "\n",
    "In this example, we will use the first draft of the indicative scenario extensions proposed for ScenarioMIP for CMIP7 (https://github.com/chrisroadmap/explore-extensions). **Note:** these are draft scenarios and will not be the final ones used for CMIP7, which will be produced by integrated assessment models, so please don't use them naively in your own work - they are provided here as an example of how to use `fair` to read in custom scenarios!\n",
    "\n",
    "We invite you to inspect the format of the emissions file at `/data/calibrated_constrained_ensemble/extensions_1750-2500.csv`. You will note that the file format is similar to the `IamDataFrame` of [`pyam`](https://pyam-iamc.readthedocs.io/en/stable/index.html), with two key exceptions:\n",
    "\n",
    "- the `model` entry is optional. You can provide it, but it will be ignored (as will any other metadata column).\n",
    "- the `scenario` entry must be unique for every scenario.\n",
    "\n",
    "Expanding the second point above, in an `IamDataFrame` we may have the same `scenario` (e.g. `SSP3-Baseline`) run in different integrated assessment models (`model` could be, for example, `MESSAGE-GLOBIOM 1.0` or `REMIND-MAgPIE 4.2`). In `fair`, if you want to distinguish similar scenarios run by different IAMs in the same emissions data file, then you would want to modify the `scenario` column:\n",
    "\n",
    "- `MESSAGE-GLOBIOM 1.0___SSP3-Baseline`\n",
    "- `REMIND-MAgPIE 4.2___SSP3-Baseline`\n",
    "- ...\n",
    "\n",
    "(I use a triple underscore - you could use any separator you like, as long as it is not a string that that is present in any model or scenario name)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "scenarios = [\n",
    "    \"high-extension\", \n",
    "    \"high-overshoot\",\n",
    "    \"medium-overshoot\", \n",
    "    \"medium-extension\", \n",
    "    \"low\", \n",
    "    \"verylow\",\n",
    "    \"verylow-overshoot\", \n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "f.define_scenarios(scenarios)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "## 4. Define configs\n",
    "\n",
    "The constrained dataset contains 841 ensemble members, and 86 parameters that define the climate response of `fair`. The parameters pertain to `climate_configs` and `species_configs` that produce a wide range of climate responses. These values are given in the CSV file at `data/calibrated_constrained_ensemble/calibrated_constrained_parameters_calibration1.4.1.csv`.\n",
    "\n",
    "This file contains a header row and 841 additional rows. Each row corresponds to a parameter set used to probablistically run `fair` with. The first column is a label that we use to refer to that particular parameter set (here, is a number between 0 and 1599999, and refers to the original prior ensemble).\n",
    "\n",
    "The column headers refer to specific parameters that we are varying in `fair`. These are automatically mapped to the appropriate parameter within `fair` (so getting the names exactly correct is important).\n",
    "\n",
    "We sample from the 11 `climate_configs` parameters that define the [stochastic three-layer energy balance model](https://journals.ametsoc.org/view/journals/clim/33/18/jcliD190589.xml), plus a random seed, and two columns that tell the model if we want to use the seed and if stochastic response should be turned on (both boolean values).\n",
    "\n",
    "The other 74 parameters are `species_configs` and override default values of `species_configs` within `fair` (an example being the parameters defining the sensitivity of the carbon cycle feedbacks). Since every species has about 30 configs attached, there's well over a thousand potential parameters that could be modified in `fair`. Outside of the 74 parameters sampled, changing from default values would make little difference, would not be relevant to a particular species, or not be sensible to change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "fair_params_1_4_1_file = 'data/calibrated_constrained_ensemble/calibrated_constrained_parameters_calibration1.4.1.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "Here, we are not actually defining any configs, but we are telling `fair` what the labels of each parameter set are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_configs = pd.read_csv(fair_params_1_4_1_file, index_col=0)\n",
    "configs = df_configs.index  # this is used as a label for the \"config\" axis\n",
    "f.define_configs(configs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "configs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_configs.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "metadata": {},
   "source": [
    "## 5. Define species and properties\n",
    "\n",
    "We will use FaIR's default list of 61 species. They are often run with default properties that are included in the model code. However, as part of the v1.4.1 calibration, some defaults are modified, such as the sensitivity of chemical precursors to methane lifetime. Rather than manually overriding this by setting `species_configs`, it is cleaner to modify the defaults in the CSV file that is read in to define the `species` and `properties`. \n",
    "\n",
    "In fact, as this only reads in and defines `species` and `properties` (not `species_configs`), the default (no `filename`) argument could be used here, but it is efficient to put species, properties and configs in the same file, and to use the same file to read in `properties` and `species_configs`.\n",
    "\n",
    "If you're following along at home, feel free to insert a new cell after this one and inspect what the `species` and `properties` actually are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "fair_species_configs_1_4_1_file = 'data/calibrated_constrained_ensemble/species_configs_properties_calibration1.4.1.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "species, properties = read_properties(filename=fair_species_configs_1_4_1_file)\n",
    "f.define_species(species, properties)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19",
   "metadata": {},
   "source": [
    "## 6. Modify run options\n",
    "\n",
    "Not necessary, as we made all of our choices on initialisation (step 1)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {},
   "source": [
    "## 7. Create input and output xarrays\n",
    "\n",
    "If this runs without error, the problem is consistently and completely set up: we then just need to add data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "f.allocate()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22",
   "metadata": {},
   "source": [
    "## 8. Fill in data\n",
    "\n",
    "### 8a. emissions, solar forcing, and volcanic forcing\n",
    "\n",
    "We can use the new (in v2.2) convenience function `fill_from_csv()` to fill in the emissions from the emissions file that we created offline. Remember that not all `species` are things that take emissions, so if you see some NaN entries below, don't panic.\n",
    "\n",
    "There are two species defined - `solar` and `volcanic` - that take offline forcing time series, so they also need to be defined in a file and read in using `fill_from_csv()`. The file structure is similar to the emissions file (and we recommend that you inspect it) - but remember that forcing is defined on timebounds rather than timepoints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23",
   "metadata": {},
   "outputs": [],
   "source": [
    "f.fill_from_csv(\n",
    "    emissions_file='data/calibrated_constrained_ensemble/extensions_1750-2500.csv',\n",
    "    forcing_file='data/calibrated_constrained_ensemble/volcanic_solar.csv',\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [],
   "source": [
    "f.emissions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "f.forcing.sel(specie=\"Volcanic\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26",
   "metadata": {},
   "source": [
    "There's one slight adjustment we need to make - in order to ensure that the solar and volcanic scale factors are picked up, we have to manually adjust the forcing time series. In future, we hope to make this a little more automatic. See https://github.com/OMS-NetZero/FAIR/issues/126."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [],
   "source": [
    "fill(\n",
    "    f.forcing,\n",
    "    f.forcing.sel(specie=\"Volcanic\") * df_configs[\"forcing_scale[Volcanic]\"].values.squeeze(),\n",
    "    specie=\"Volcanic\",\n",
    ")\n",
    "fill(\n",
    "    f.forcing,\n",
    "    f.forcing.sel(specie=\"Solar\") * df_configs[\"forcing_scale[Solar]\"].values.squeeze(),\n",
    "    specie=\"Solar\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28",
   "metadata": {},
   "outputs": [],
   "source": [
    "pl.plot(f.timebounds, f.forcing.loc[dict(specie=\"Solar\", scenario=\"medium-extension\")]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29",
   "metadata": {},
   "source": [
    "### 8b. Fill in climate_configs and species_configs\n",
    "\n",
    "The new convenience methods in v2.2 make this very easy indeed. First we fill in the default values from the `species_configs` file, and then we use our 86 parameter set for 841 ensemble members to change all of the parameters that are pertinent to the key model responses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30",
   "metadata": {},
   "outputs": [],
   "source": [
    "f.fill_species_configs(fair_species_configs_1_4_1_file)\n",
    "f.override_defaults(fair_params_1_4_1_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31",
   "metadata": {},
   "source": [
    "### 8c. Initial conditions\n",
    "\n",
    "It's important these are defined, as they are NaN by default, and it's likely you'll run into problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32",
   "metadata": {},
   "outputs": [],
   "source": [
    "initialise(f.concentration, f.species_configs[\"baseline_concentration\"])\n",
    "initialise(f.forcing, 0)\n",
    "initialise(f.temperature, 0)\n",
    "initialise(f.cumulative_emissions, 0)\n",
    "initialise(f.airborne_emissions, 0)\n",
    "initialise(f.ocean_heat_content_change, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33",
   "metadata": {},
   "source": [
    "## 9. Run\n",
    "\n",
    "We have a total of 7 scenarios and 841 ensemble members for 550 years and 61 species. This can be a little memory constrained on some systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34",
   "metadata": {},
   "outputs": [],
   "source": [
    "f.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35",
   "metadata": {},
   "source": [
    "## 10. Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36",
   "metadata": {},
   "outputs": [],
   "source": [
    "fancy_titles = {\n",
    "    'high-extension': 'High extension',\n",
    "    'high-overshoot': 'High overshoot',\n",
    "    'medium-extension': 'Medium extension',\n",
    "    'medium-overshoot': 'Medium overshoot',\n",
    "    'low': 'Low',\n",
    "    'verylow': 'Very low',\n",
    "    'verylow-overshoot': 'Very low overshoot',\n",
    "}\n",
    "\n",
    "colors = {\n",
    "    'high-extension': '#800000',\n",
    "    'high-overshoot': '#ff0000',\n",
    "    'medium-extension': '#c87820',\n",
    "    'medium-overshoot': '#d3a640',\n",
    "    'low': '#098740',\n",
    "    'verylow': '#0080d0',\n",
    "    'verylow-overshoot': '#100060',\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37",
   "metadata": {},
   "source": [
    "### Temperature anomaly\n",
    "\n",
    "We define an anomaly baseline of 1850-1900. This is 51 complete years. As FaIR temperature anomalies are on `timebounds`, we take mid-year temperatures as averages of the bounding `timebounds`; so, 1850.5 is an average of 1850.0 and 1851.0. It means we take an average period of 1850-1901 timebounds with 0.5 weights for 1850 and 1901 and 1.0 weights for other `timebounds`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38",
   "metadata": {},
   "outputs": [],
   "source": [
    "weights_51yr = np.ones(52)\n",
    "weights_51yr[0] = 0.5\n",
    "weights_51yr[-1] = 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = pl.subplots(2, 4, figsize=(12, 6))\n",
    "\n",
    "for i, scenario in enumerate(scenarios):\n",
    "    for pp in ((0, 100), (5, 95), (16, 84)):\n",
    "        ax[i // 4, i % 4].fill_between(\n",
    "            f.timebounds,\n",
    "            np.percentile(\n",
    "                f.temperature.loc[dict(scenario=scenario, layer=0)]\n",
    "                - np.average(\n",
    "                    f.temperature.loc[\n",
    "                        dict(scenario=scenario, timebounds=np.arange(1850, 1902), layer=0)\n",
    "                    ],\n",
    "                    weights=weights_51yr,\n",
    "                    axis=0\n",
    "                ),\n",
    "                pp[0],\n",
    "                axis=1,\n",
    "            ),\n",
    "            np.percentile(\n",
    "                f.temperature.loc[dict(scenario=scenario, layer=0)]\n",
    "                - np.average(\n",
    "                    f.temperature.loc[\n",
    "                        dict(scenario=scenario, timebounds=np.arange(1850, 1902), layer=0)\n",
    "                    ],\n",
    "                    weights=weights_51yr,\n",
    "                    axis=0\n",
    "                ),\n",
    "                pp[1],\n",
    "                axis=1,\n",
    "            ),\n",
    "            color=colors[scenarios[i]],\n",
    "            alpha=0.2,\n",
    "            lw=0\n",
    "        )\n",
    "\n",
    "    ax[i // 4, i % 4].plot(\n",
    "        f.timebounds,\n",
    "        np.median(\n",
    "            f.temperature.loc[dict(scenario=scenario, layer=0)]\n",
    "            - np.average(\n",
    "                f.temperature.loc[\n",
    "                    dict(scenario=scenario, timebounds=np.arange(1850, 1902), layer=0)\n",
    "                ],\n",
    "                weights=weights_51yr,\n",
    "                axis=0\n",
    "            ),\n",
    "            axis=1,\n",
    "        ),\n",
    "        color=colors[scenarios[i]],\n",
    "    )\n",
    "    ax[i // 4, i % 4].set_xlim(1850, 2300)\n",
    "    ax[i // 4, i % 4].set_ylim(-1, 10)\n",
    "    ax[i // 4, i % 4].axhline(0, color=\"k\", ls=\":\", lw=0.5)\n",
    "    ax[i // 4, i % 4].set_title(fancy_titles[scenarios[i]])\n",
    "\n",
    "pl.suptitle(\"Temperature anomalies\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40",
   "metadata": {},
   "source": [
    "### CO2 concentrations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = pl.subplots(2, 4, figsize=(12, 6))\n",
    "\n",
    "for i, scenario in enumerate(scenarios):\n",
    "    for pp in ((0, 100), (5, 95), (16, 84)):\n",
    "        ax[i // 4, i % 4].fill_between(\n",
    "            f.timebounds,\n",
    "            np.percentile(\n",
    "                f.concentration.loc[dict(scenario=scenario, specie='CO2')],\n",
    "                pp[0],\n",
    "                axis=1,\n",
    "            ),\n",
    "            np.percentile(\n",
    "                f.concentration.loc[dict(scenario=scenario, specie='CO2')],\n",
    "                pp[1],\n",
    "                axis=1,\n",
    "            ),\n",
    "            color=colors[scenarios[i]],\n",
    "            alpha=0.2,\n",
    "            lw=0\n",
    "        )\n",
    "\n",
    "    ax[i // 4, i % 4].plot(\n",
    "        f.timebounds,\n",
    "        np.median(\n",
    "            f.concentration.loc[dict(scenario=scenario, specie='CO2')],\n",
    "            axis=1,\n",
    "        ),\n",
    "        color=colors[scenarios[i]],\n",
    "    )\n",
    "    ax[i // 4, i % 4].set_xlim(1850, 2300)\n",
    "    ax[i // 4, i % 4].set_ylim(0, 2000)\n",
    "    ax[i // 4, i % 4].axhline(0, color=\"k\", ls=\":\", lw=0.5)\n",
    "    ax[i // 4, i % 4].set_title(fancy_titles[scenarios[i]])\n",
    "\n",
    "pl.suptitle(\"CO$_2$ concentration\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42",
   "metadata": {},
   "source": [
    "### Total effective radiative forcing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = pl.subplots(2, 4, figsize=(12, 6))\n",
    "\n",
    "for i, scenario in enumerate(scenarios):\n",
    "    for pp in ((0, 100), (5, 95), (16, 84)):\n",
    "        ax[i // 4, i % 4].fill_between(\n",
    "            f.timebounds,\n",
    "            np.percentile(\n",
    "                f.forcing_sum.loc[dict(scenario=scenario)],\n",
    "                pp[0],\n",
    "                axis=1,\n",
    "            ),\n",
    "            np.percentile(\n",
    "                f.forcing_sum.loc[dict(scenario=scenario)],\n",
    "                pp[1],\n",
    "                axis=1,\n",
    "            ),\n",
    "            color=colors[scenarios[i]],\n",
    "            alpha=0.2,\n",
    "            lw=0\n",
    "        )\n",
    "\n",
    "    ax[i // 4, i % 4].plot(\n",
    "        f.timebounds,\n",
    "        np.median(\n",
    "            f.forcing_sum.loc[dict(scenario=scenario)],\n",
    "            axis=1,\n",
    "        ),\n",
    "        color=colors[scenarios[i]],\n",
    "    )\n",
    "    ax[i // 4, i % 4].set_xlim(1850, 2300)\n",
    "    ax[i // 4, i % 4].set_ylim(-2, 15)\n",
    "    ax[i // 4, i % 4].axhline(0, color=\"k\", ls=\":\", lw=0.5)\n",
    "    ax[i // 4, i % 4].set_title(fancy_titles[scenarios[i]])\n",
    "\n",
    "pl.suptitle(\"Effective radiative forcing\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44",
   "metadata": {},
   "source": [
    "### CO2 airborne fraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = pl.subplots(2, 4, figsize=(12, 6))\n",
    "\n",
    "for i, scenario in enumerate(scenarios):\n",
    "    for pp in ((0, 100), (5, 95), (16, 84)):\n",
    "        ax[i // 4, i % 4].fill_between(\n",
    "            f.timebounds,\n",
    "            np.percentile(\n",
    "                f.airborne_fraction.loc[dict(scenario=scenario, specie='CO2')],\n",
    "                pp[0],\n",
    "                axis=1,\n",
    "            ),\n",
    "            np.percentile(\n",
    "                f.airborne_fraction.loc[dict(scenario=scenario, specie='CO2')],\n",
    "                pp[1],\n",
    "                axis=1,\n",
    "            ),\n",
    "            color=colors[scenarios[i]],\n",
    "            alpha=0.2,\n",
    "            lw=0\n",
    "        )\n",
    "\n",
    "    ax[i // 4, i % 4].plot(\n",
    "        f.timebounds,\n",
    "        np.median(\n",
    "            f.airborne_fraction.loc[dict(scenario=scenario, specie='CO2')],\n",
    "            axis=1,\n",
    "        ),\n",
    "        color=colors[scenarios[i]],\n",
    "    )\n",
    "    ax[i // 4, i % 4].set_xlim(1850, 2300)\n",
    "    ax[i // 4, i % 4].set_ylim(0, 1)\n",
    "    ax[i // 4, i % 4].axhline(0, color=\"k\", ls=\":\", lw=0.5)\n",
    "    ax[i // 4, i % 4].set_title(fancy_titles[scenarios[i]])\n",
    "\n",
    "pl.suptitle(\"CO$_2$ airborne fraction\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46",
   "metadata": {},
   "source": [
    "### Earth's energy uptake"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = pl.subplots(2, 4, figsize=(12, 6))\n",
    "\n",
    "for i, scenario in enumerate(scenarios):\n",
    "    for pp in ((0, 100), (5, 95), (16, 84)):\n",
    "        ax[i // 4, i % 4].fill_between(\n",
    "            f.timebounds,\n",
    "            np.percentile(\n",
    "                f.ocean_heat_content_change.loc[dict(scenario=scenario)],\n",
    "                pp[0],\n",
    "                axis=1,\n",
    "            ),\n",
    "            np.percentile(\n",
    "                f.ocean_heat_content_change.loc[dict(scenario=scenario)],\n",
    "                pp[1],\n",
    "                axis=1,\n",
    "            ),\n",
    "            color=colors[scenarios[i]],\n",
    "            alpha=0.2,\n",
    "            lw=0\n",
    "        )\n",
    "\n",
    "    ax[i // 4, i % 4].plot(\n",
    "        f.timebounds,\n",
    "        np.median(\n",
    "            f.ocean_heat_content_change.loc[dict(scenario=scenario)],\n",
    "            axis=1,\n",
    "        ),\n",
    "        color=colors[scenarios[i]],\n",
    "    )\n",
    "    ax[i // 4, i % 4].set_xlim(1850, 2300)\n",
    "    ax[i // 4, i % 4].set_ylim(0, 1e25)\n",
    "    ax[i // 4, i % 4].axhline(0, color=\"k\", ls=\":\", lw=0.5)\n",
    "    ax[i // 4, i % 4].set_title(fancy_titles[scenarios[i]])\n",
    "\n",
    "pl.suptitle(\"Earth energy uptake\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}