2020-11-14-climate-of-the-ocean-h1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1, policy goals under uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1, policy goals under uncertainty\n",
    "A recent ground-breaking [review paper](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019RG000678) produced the most comprehensive and up-to-date estimate of the *climate feedback parameter*, which they find to be\n",
    "\n",
    "$B \\approx \\mathcal{N}(-1.3, 0.4),$\n",
    "\n",
    "i.e. our knowledge of the real value is normally distributed with a mean value $\\overline{B} = -1.3$ W/m²/K and a standard deviation $\\sigma = 0.4$ W/m²/K. These values are not very intuitive, so let us convert them into more policy-relevant numbers.\n",
    "\n",
    "**Definition:** *Equilibrium climate sensitivity (ECS)* is defined as the amount of warming $\\Delta T$ caused by a doubling of $CO_2$ (e.g. from the pre-industrial value 280 ppm to 560 ppm), at equilibrium.\n",
    "\n",
    "At equilibrium, the energy balance model equation is:\n",
    "\n",
    "$0 = \\frac{S(1 - α)}{4} - (A - BT_{eq}) + a \\ln\\left( \\frac{2\\;\\text{CO}_{2\\text{PI}}}{\\text{CO}_{2\\text{PI}}} \\right)$\n",
    "\n",
    "From this, we subtract the preindustrial energy balance, which is given by:\n",
    "\n",
    "$0 = \\frac{S(1-α)}{4} - (A - BT_{0}),$\n",
    "\n",
    "The result of this subtraction, after rearranging, is our definition of $\\text{ECS}$:\n",
    "\n",
    "$\\text{ECS} \\equiv T_{eq} - T_{0} = -\\frac{a\\ln(2)}{B}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from energy_balance_model import ebm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create ecs(B, a) function\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def double_CO2(t):\n",
    "    return 280 * 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1)\n",
    "EBM = ebm(14, 0, 1., double_CO2)\n",
    "EBM.run(300)\n",
    "ax.plot(EBM.t, EBM.T - EBM.T[0], label = \"$\\Delta T (t) = T(t) - T_0$\", color = \"red\")\n",
    "\n",
    "ax.axhline(ecs(EBM.B, EBM.a), label = \"ECS\", color = \"darkred\", ls = \"--\")\n",
    "\n",
    "ax.legend()\n",
    "ax.grid()\n",
    "ax.set_title(\"Transient response to instant doubling of CO$_2$\")\n",
    "ax.set_ylabel(\"temperature [°C]\")\n",
    "ax.set_xlabel(\"years after doubling\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot above provides an example of an \"abrupt 2xCO$_2$\" experiment, a classic experimental treatment method in climate modelling which is used in practice to estimate ECS for a particular model (Note: in complicated climate models the values of the parameters $a$ and $B$ are not specified a priori, but emerge as outputs for the simulation).\n",
    "\n",
    "The simulation begins at the preindustrial equilibrium, i.e. a temperature °C is in balance with the pre-industrial CO$_2$ concentration of 280 ppm until CO$_2$ is abruptly doubled from 280 ppm to 560 ppm. The climate responds by rapidly warming, and after a few hundred years approaches the equilibrium climate sensitivity value, by definition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a graph to visualize ECS as a function of B (B should be on the x axis)\n",
    "# calculate the range from -2 to -0.1 with 0.1 as a step size\n",
    "# Note use plt.scatter for plotting and \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question:\n",
    "\n",
    "(1) What does it mean for a climate system to have a more negative value of $B$? Explain why we call $B$ the climate feedback parameter.\n",
    "\n",
    "Answer:\n",
    "\n",
    "(2) What happens when $B$ is greater than or equal to zero?\n",
    "\n",
    "Answer:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 1.2 - _Doubling CO\n",
    "\n",
    "To compute ECS, we doubled the $CO_2$ in our atmosphere. This factor 2 is not entirely arbitrary: without substantial effort to reduce $CO_2$ emissions, we are expected to **at least** double the $CO_2$ in our atmosphere by 2100. \n",
    "\n",
    "Right now, our $CO_2$ concentration is 415 ppm -- 1.482 times the pre-industrial value of 280 ppm from 1850. \n",
    "\n",
    "The $CO_2$ concentrations in the _future_ depend on human action. There are several models for future concentrations, which are formed by assuming different _policy scenarios_. A baseline model is RCP8.5 - a \"worst-case\" high-emissions scenario. In our notebook, this model is given as a function of ``t``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CO2_RCP85(t):\n",
    "    return 280 * (1+ ((t-1850)/220)**3 * np.maximum(1., np.exp(((t-1850)-170)/100)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.arange(1850, 2100)\n",
    "plt.ylabel(\"CO$_2$ concentration [ppm]\")\n",
    "plt.plot(t, CO2_RCP85(t));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question:\n",
    "\n",
    "In what year are we expected to have doubled the $CO_2$ concentration, under policy scenario RCP8.5?\n",
    "\n",
    "Hint: the function \n",
    "```python\n",
    "np.where()\n",
    "``` \n",
    "might be useful"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Enter your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1.3 - _Uncertainty in B_\n",
    "\n",
    "The climate feedback parameter ``B`` is not something that we can control– it is an emergent property of the global climate system. Unfortunately, ``B`` is also difficult to quantify empirically (the relevant processes are difficult or impossible to observe directly), so there remains uncertainty as to its exact value.\n",
    "\n",
    "A value of ``B`` close to zero means that an increase in $CO_2$ concentrations will have a larger impact on global warming, and that more action is needed to stay below a maximum temperature. In answering such policy-related question, we need to take the uncertainty in ``B`` into account. In this exercise, we will do so using a Monte Carlo simulation: we generate a sample of values for ``B``, and use these values in our analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a probability distribution for for $B_{avg}$ above. Plot a histogram. \n",
    "\n",
    "Hint: use the functions\n",
    "\n",
    "```python\n",
    "np.random.normal() # with 50000 samples\n",
    "```\n",
    "and plot with \n",
    "```python\n",
    "plt.hist()\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = 0.4\n",
    "b_avg = -1.3\n",
    "\n",
    "samples = # Enter code here\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "    # plot here\n",
    "plt.xlabel(\"B [W/m²/K]\")\n",
    "plt.ylabel(\"samples\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a probability distribution for the ECS based on the probability distribution function for $B$ above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values =  # your code here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values = np.where((values < -20) | (values > 20) , np.nan, values) # drop outlier\n",
    "plt.hist(values, bins = 20)\n",
    "plt.xlim([0, 20])\n",
    "plt.xlabel(\"Temperature [°C]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like the ECS distribution is **not normally distributed**, even though $B$ is. \n",
    "\n",
    "Question: How does $\\overline{\\text{ECS}(B)}$ compare to $\\text{ECS}(\\overline{B})$? What is the probability that $\\text{ECS}(B)$ lies above $\\text{ECS}(\\overline{B})$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question: Does accounting for uncertainty in feedbacks make our expectation of global warming better (less implied warming) or worse (more implied warming)?\n",
    "\n",
    "Answer:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1.5 - _Running the model_\n",
    "\n",
    "In the lecture notebook we introduced a class `ebm` (_energy balance model_), which contains:\n",
    "- the parameters of our climate simulation (`C`, `a`, `A`, `B`, `CO2_PI`, `alpha`, `S`, see details below)\n",
    "- a function `CO2`, which maps a time `t` to the concentrations at that year. For example, we use the function `t -> 280` to simulate a model with concentrations fixed at 280 ppm.\n",
    "\n",
    "`ebm` also contains the simulation results, in two arrays:\n",
    "- `T` is the array of tempartures (°C, `Float64`).\n",
    "- `t` is the array of timestamps (years, `Float64`), of the same size as `T`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can set up an instance of `ebm` like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_co2function(t):\n",
    "    # here we imply NO co2 increase\n",
    "    return 280\n",
    "\n",
    "my_model = ebm(T=14, t=0, deltat=1., CO2=my_co2function)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look into our ebm object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "attrs = vars(my_model)\n",
    "print(', \\n'.join(\"%s: %s\" % item for item in attrs.items()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What function do we have?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(my_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the model\n",
    "\n",
    "EBM = ebm(14, 1850, 1, my_co2function)\n",
    "EBM.run(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, look inside `simulated_model` and notice that `T` and `t` have accumulated the simulation results.\n",
    "\n",
    "In this simulation, we used `T0 = 14` and `CO2 = 280`, which is why `T` is constant during our simulation. These parameters are the default, pre-industrial values, and our model is based on this equilibrium.\n",
    "\n",
    "`Question`: Run a simulation starting at 1850 with policy scenario RCP8.5, and plot the computed temperature graph. What is the global temperature at 2100?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CO2_RCP85(t):\n",
    "    return 280 * (1+ ((t-1850)/220)**3 * np.maximum(1., np.exp(((t-1850)-170)/100)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Run the model here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can change values before running the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EBM = ebm(15, 1850, 1, my_co2function)\n",
    "EBM.B = -2\n",
    "EBM.run(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EBM.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1.6 - _Application to policy relevant questions_ (BONUS)\n",
    "\n",
    "We talked about two _emissions scenarios_: RCP2.6 (strong mitigation - controlled CO2 concentrations) and RCP8.5 (no mitigation - high CO2 concentrations). These are given by the following functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CO2_RCP26(t):\n",
    "    return 280 * (1+ ((t-1850)/220)**3 * np.minimum(1., np.exp(-((t-1850)-170)/100)))\n",
    "def CO2_RCP85(t):\n",
    "    return 280 * (1+ ((t-1850)/220)**3 * np.maximum(1., np.exp(((t-1850)-170)/100)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are interested in how the **uncertainty in our input** $B$ (the climate feedback paramter) *propagates* through our model to determine the **uncertainty in our output** $T(t)$, for a given emissions scenario. The goal of this exercise is to answer the following by using *Monte Carlo Simulation* for *uncertainty propagation*:\n",
    "\n",
    "> What is the probability that we see more than 2°C of warming by 2100 under the low-emissions scenario RCP2.6? What about under the high-emissions scenario RCP8.5?"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}