Exercise 1, policy goals under uncertainty#

Exercise 1, policy goals under uncertainty#

A recent ground-breaking review paper produced the most comprehensive and up-to-date estimate of the climate feedback parameter, which they find to be

\(B \approx \mathcal{N}(-1.3, 0.4),\)

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.

Definition 1 (Equilibrium climate sensitivity (ECS))

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.

At equilibrium, the energy balance model equation is:

\[0 = \frac{S(1 - α)}{4} - (A - BT_{eq}) + a \ln\left( \frac{2\;\text{CO}_{2\text{PI}}}{\text{CO}_{2\text{PI}}} \right)\]

From this, we subtract the preindustrial energy balance, which is given by:

\[0 = \frac{S(1-α)}{4} - (A - BT_{0}),\]

The result of this subtraction, after rearranging, is our definition of \(\text{ECS}\):

\[\text{ECS} \equiv T_{eq} - T_{0} = -\frac{a\ln(2)}{B}\]
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt

class EBM:
    """
    Zero Order Energy Balance Model (EBM)

    The Energy Balance Model (EBM) represents the balance between incoming solar radiation and outgoing thermal radiation.
    It also considers the greenhouse effect caused by CO2 levels. This model can simulate the temporal evolution of temperature 
    based on various parameters like albedo, solar constant, and greenhouse effect coefficients.

    Attributes:
    - T : Temperature (in Kelvin)
    - t : Time
    - deltat : Time step
    - CO2 : Carbon Dioxide function that returns CO2 levels in dependency of time t
    - C : Heat capacity
    - a : Greenhouse effect coefficient
    - A : Outgoing thermal radiation constant
    - B : Temperature sensitivity of outgoing radiation
    - CO2_PI : Pre-industrial CO2 concentration
    - alpha : Albedo
    - S : Solar constant
    """

    def __init__(self, T, t, deltat, CO2, C, a, A, B, CO2_PI, alpha, S):
        self.T = np.array(T)
        self.t = t
        self.deltat = deltat
        self.C = C
        self.a = a
        self.A = A
        self.B = B
        self.co2_pi = CO2_PI
        self.alpha = alpha
        self.S = S
        self.co2 = CO2

    def absorbed_solar_radiation(self, S, alpha):
        return (S*(1-alpha)/4)

    def outgoing_thermal_radiation(self, T, A, B):
        return A - B*T

    def greenhouse_effect(self, CO2, a=5, CO2_PI = 280):
        return a*np.log(CO2/CO2_PI)

    def tendency(self):
        current_T = self.T[-1] if self.T.size > 1 else self.T
        current_t = self.t[-1] if self.T.size > 1 else self.t
        
        return 1. / self.C * (
            + self.absorbed_solar_radiation(S=self.S, alpha=self.alpha)
            - self.outgoing_thermal_radiation(current_T, A=self.A, B=self.B)
            + self.greenhouse_effect(self.co2(current_t), a=self.a, CO2_PI=self.co2_pi)
        )

    @property
    def timestep(self):
        new_T = self.T[-1] + self.deltat * self.tendency() if self.T.size > 1 else self.T + self.deltat * self.tendency()
        new_t = self.t[-1] + self.deltat if self.T.size > 1 else self.t + self.deltat
        
        self.T = np.append(self.T, new_T)
        self.t = np.append(self.t, new_t)
        
    def run(self, end_year):
        for year in range(end_year):
            self.timestep

# Define the ECS function

def double_CO2(t):
    return 280 * 2

model_parameters = {
    "T":14,
    "t":0,
    "deltat":1,
    "CO2":double_CO2,
    "C":51,
    "a":5,
    "A":221.2,
    "B":-1.3,
    "CO2_PI": 280,
    "alpha":0.3,
    "S":1368
}

model = EBM(**model_parameters)
model.run(300)

f, ax = plt.subplots(1)

ax.plot(model.t, model.T - model.T[0], label = "$\Delta T (t) = T(t) - T_0$", color = "red")
ax.axhline(ecs(model.B, model.a), label = "ECS", color = "darkred", ls = "--")
ax.legend()
ax.grid()
ax.set_title("Transient response to instant doubling of CO$_2$")
ax.set_ylabel("temperature [°C]")
ax.set_xlabel("years after doubling")

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).

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.

# Create a graph to visualize ECS as a function of B (B should be on the x axis)
# calculate the range from -2 to -0.1 with 0.1 as a step size
# Note use plt.scatter for plotting and

Question:

(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.

Answer:

(2) What happens when \(B\) is greater than or equal to zero?

Answer:

Exercise 1.2 - _Doubling CO#

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.

Right now, our \(CO_2\) concentration is 415 ppm – 1.482 times the pre-industrial value of 280 ppm from 1850.

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.

def CO2_RCP85(t):
    return 280 * (1+ ((t-1850)/220)**3 * np.maximum(1., np.exp(((t-1850)-170)/100)))

t = np.arange(1850, 2100)
plt.ylabel("CO$_2$ concentration [ppm]")
plt.plot(t, CO2_RCP85(t));

Question:

In what year are we expected to have doubled the \(CO_2\) concentration, under policy scenario RCP8.5?

Hint: the function

np.where()

might be useful

# Enter your code here

Answer:

Exercise 1.3 - Uncertainty in B#

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.

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.

Generate a probability distribution for for \(B_{avg}\) above. Plot a histogram.

Hint: use the functions

np.random.normal() # with 50000 samples

and plot with

plt.hist()

sigma = 0.4
b_avg = -1.3

samples = # Enter code here

# plot here
plt.xlabel("B [W/m²/K]")
plt.ylabel("samples")

Generate a probability distribution for the ECS based on the probability distribution function for \(B\) above.

values =  # your code here

values = np.where((values < -20) | (values > 20) , np.nan, values) # drop outlier
plt.hist(values, bins = 50)
plt.xlim([0, 20])
plt.xlabel("Temperature [°C]")

It looks like the ECS distribution is not normally distributed, even though \(B\) is.

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})\)?

# your code here

Question: Does accounting for uncertainty in feedbacks make our expectation of global warming better (less implied warming) or worse (more implied warming)?

Answer:

Exercise 1.5 - Running the model#

In the lecture notebook we introduced a class EBM (energy balance model), which contains:

  • the parameters of our climate simulation (C, a, A, B, CO2_PI, alpha, S, see details below)

  • 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.

EBM also contains the simulation results, in two arrays:

  • T is the array of tempartures (°C, Float64).

  • t is the array of timestamps (years, Float64), of the same size as T.

You can set up an instance of EBM like so:

def my_co2function(t):
    # here we imply NO co2 increase
    return 280

model_parameters = {
    "T":14,
    "t":0,
    "deltat":1,
    "CO2":1,
    "C":51,
    "a":5,
    "A":221.2,
    "B":-1.3,
    "CO2_PI": 280,
    "alpha":0.3,
    "S":1368
}

model_parameters["CO2"] = my_co2function

my_model = EBM(**model_parameters)

Let’s look into our ebm object

attrs = vars(my_model)
print(', \n'.join("%s: %s" % item for item in attrs.items()))

What function do we have?

help(my_model)

# Run the model

Again, look inside simulated_model and notice that T and t have accumulated the simulation results.

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.

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?

def CO2_RCP85(t):
    return 280 * (1+ ((t-1850)/220)**3 * np.maximum(1., np.exp(((t-1850)-170)/100)))

## Run the model here

We can change values before running the model.

model = EBM(**model_parameters)
model.B = -2
model.run(10)

model.T

Exercise 1.6 - Application to policy relevant questions (BONUS)#

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

def CO2_RCP26(t):
    return 280 * (1+ ((t-1850)/220)**3 * np.minimum(1., np.exp(-((t-1850)-170)/100)))
def CO2_RCP85(t):
    return 280 * (1+ ((t-1850)/220)**3 * np.maximum(1., np.exp(((t-1850)-170)/100)))

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:

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?