{
 "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",
    "```{prf:definition} Equilibrium climate sensitivity (ECS)\n",
    "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",
    "\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class EBM:\n",
    "    \"\"\"\n",
    "    Zero Order Energy Balance Model (EBM)\n",
    "\n",
    "    The Energy Balance Model (EBM) represents the balance between incoming solar radiation and outgoing thermal radiation.\n",
    "    It also considers the greenhouse effect caused by CO2 levels. This model can simulate the temporal evolution of temperature \n",
    "    based on various parameters like albedo, solar constant, and greenhouse effect coefficients.\n",
    "\n",
    "    Attributes:\n",
    "    - T : Temperature (in Kelvin)\n",
    "    - t : Time\n",
    "    - deltat : Time step\n",
    "    - CO2 : Carbon Dioxide function that returns CO2 levels in dependency of time t\n",
    "    - C : Heat capacity\n",
    "    - a : Greenhouse effect coefficient\n",
    "    - A : Outgoing thermal radiation constant\n",
    "    - B : Temperature sensitivity of outgoing radiation\n",
    "    - CO2_PI : Pre-industrial CO2 concentration\n",
    "    - alpha : Albedo\n",
    "    - S : Solar constant\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, T, t, deltat, CO2, C, a, A, B, CO2_PI, alpha, S):\n",
    "        self.T = np.array(T)\n",
    "        self.t = t\n",
    "        self.deltat = deltat\n",
    "        self.C = C\n",
    "        self.a = a\n",
    "        self.A = A\n",
    "        self.B = B\n",
    "        self.co2_pi = CO2_PI\n",
    "        self.alpha = alpha\n",
    "        self.S = S\n",
    "        self.co2 = CO2\n",
    "\n",
    "    def absorbed_solar_radiation(self, S, alpha):\n",
    "        return (S*(1-alpha)/4)\n",
    "\n",
    "    def outgoing_thermal_radiation(self, T, A, B):\n",
    "        return A - B*T\n",
    "\n",
    "    def greenhouse_effect(self, CO2, a=5, CO2_PI = 280):\n",
    "        return a*np.log(CO2/CO2_PI)\n",
    "\n",
    "    def tendency(self):\n",
    "        current_T = self.T[-1] if self.T.size > 1 else self.T\n",
    "        current_t = self.t[-1] if self.T.size > 1 else self.t\n",
    "        \n",
    "        return 1. / self.C * (\n",
    "            + self.absorbed_solar_radiation(S=self.S, alpha=self.alpha)\n",
    "            - self.outgoing_thermal_radiation(current_T, A=self.A, B=self.B)\n",
    "            + self.greenhouse_effect(self.co2(current_t), a=self.a, CO2_PI=self.co2_pi)\n",
    "        )\n",
    "\n",
    "    @property\n",
    "    def timestep(self):\n",
    "        new_T = self.T[-1] + self.deltat * self.tendency() if self.T.size > 1 else self.T + self.deltat * self.tendency()\n",
    "        new_t = self.t[-1] + self.deltat if self.T.size > 1 else self.t + self.deltat\n",
    "        \n",
    "        self.T = np.append(self.T, new_T)\n",
    "        self.t = np.append(self.t, new_t)\n",
    "        \n",
    "    def run(self, end_year):\n",
    "        for year in range(end_year):\n",
    "            self.timestep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the ECS function"
   ]
  },
  {
   "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": [
    "model_parameters = {\n",
    "    \"T\":14,\n",
    "    \"t\":0,\n",
    "    \"deltat\":1,\n",
    "    \"CO2\":double_CO2,\n",
    "    \"C\":51,\n",
    "    \"a\":5,\n",
    "    \"A\":221.2,\n",
    "    \"B\":-1.3,\n",
    "    \"CO2_PI\": 280,\n",
    "    \"alpha\":0.3,\n",
    "    \"S\":1368\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = EBM(**model_parameters)\n",
    "model.run(300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1)\n",
    "\n",
    "ax.plot(model.t, model.T - model.T[0], label = \"$\\Delta T (t) = T(t) - T_0$\", color = \"red\")\n",
    "ax.axhline(ecs(model.B, model.a), label = \"ECS\", color = \"darkred\", ls = \"--\")\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"
   ]
  },
  {
   "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\n"
   ]
  },
  {
   "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"
   ]
  },
  {
   "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 = 50)\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",
    "model_parameters = {\n",
    "    \"T\":14,\n",
    "    \"t\":0,\n",
    "    \"deltat\":1,\n",
    "    \"CO2\":1,\n",
    "    \"C\":51,\n",
    "    \"a\":5,\n",
    "    \"A\":221.2,\n",
    "    \"B\":-1.3,\n",
    "    \"CO2_PI\": 280,\n",
    "    \"alpha\":0.3,\n",
    "    \"S\":1368\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_parameters[\"CO2\"] = my_co2function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_model = EBM(**model_parameters)"
   ]
  },
  {
   "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"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can change values before running the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = EBM(**model_parameters)\n",
    "model.B = -2\n",
    "model.run(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.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?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}