Estimating Climate Sensitivity from Observations using Bayesian Methods

Climate sensitivity refers to increases in the earth's surface temperature due to a doubling of atmospheric CO ₂ concentration. It is a crucial quantity for predicting how fast the warming will be in the future.

In this notebook, we will estimate climate sensitivity from historical observations. Specifically, we will learn:

• How to build a simple energy balance model for Earth's climate.
  
• How to estimate key parameters in that model using Bayesian approaches.
  
• Differences between equilibrium and transient climate sensitivity.
  
• How to perform a proper model-data comparison.
  

1. A simple energy balance model for Earth's Climate ¶

The simplest way to model the earth's climate is probably to consider the earth as a single dot. Accounting for energy imbalance, climate feedback, and attendant temperature change, the following equation should be a good first-order description of Earth's temperature evolution around some equilibrium state, say the pre-industrial period,

$c_p \frac{dT'}{dt} = -\lambda T' + F'$ ............. (1)

where $T'$ and $F'$ denote, respectively, temperature and forcing anomalies.

There are two parameters in this model.

• $\lambda$ is a positive feedback parameter, which damps temperature anomalies. $\lambda$ is a combination of many processes including the Planck feedback, ice-albedo feedback, water-vapor feedback, lapse-rate feedback, etc.
• $c_p$ is the heat capacity parameter, which is often linked to the depth of the ocean mixed layer, and controls how fast temperature would respond to forcing.

At equilibirum, $\frac{dT'}{dt} = 0$, and associated temperature response is $T' = F'/\lambda$, which we call climate sensitivity.

From 1750--2020, the whole world has emitted more than 1.7 trillion tons of CO ₂ cumulatively , which has increased the atmospheric CO ₂ concentration to exceed 415 ppm and global mean temperature by about 1.2°C since the 1880s.

The figure below shows our current best estimate of historical climate forcing , $F'$, and temperature evolution, $T'$. Note that for simplicity, we do not distinguish different types of forcing, such as greenhouse gases, anthropogenic aerosols, volcanos, and the sun.

So, can we use these time series and Eq. 1 to estimate $\lambda$ and $c_p$ in order to understand climate sensitivity?

2. Solving the problem using Bayesian methods ¶

After defining the problem, we need to solve it. Specifically, we need to find combinations of $\lambda$ and $c_p$ such that, when plugged in with observed $F'$, Eq. 1 produces $T'$ that reasonably agrees with observed $T'$.

Philosophically, this fitting of a differential equation has no difference from fitting a linear model, i.e., $y = kx+b$, in the sense that we are still trying to find combinations of parameters that decrease some loss function or increase the likelihood of data. However, the problem we are dealing with does not have a closed-form analytical solution like linear regression . As a result, we need new methods.

The Bayesian method is a great candidate, with several benefits,

• It allows for estimating parameters no matter how complicated the underlying function is.
• It estimates uncertainties comprehensively and captures complicated covariance between parameters.
• We can bring our prior knowledge about parameters in the process of inference.

2.1 How do the Bayesian methods work? ¶

The essence of Bayesian methods is the Bayes rule,

$P(\theta|D) = \frac{P(D|\theta) P(\theta)}{\int P(D|\theta)P(\theta) d\theta} \propto P(D|\theta) P(\theta)$ .................... (2)

On the left-hand side, we have $P(\theta|D)$, which is the probability of parameters, $\theta$, condition on data, $D$. This quantity is what we need to calculate, which we call Posterior in Bayesian terminology.

On the right-hand side, we have $P(D|\theta)$ and $P(\theta)$. Among which,

• $P(D|\theta)$ is the probability of seeing data given a combination of parameters, which we call Likelihood ,
• $P(\theta)$ is the knowledge we have on $\theta$ before seeing any data, which we call Prior .

As a result, the process of performing Bayesian analysis is to use Likelihood to update the Prior and obtain the Posterior.

If you have no statistical background, I find the following example helpful for understanding the Bayes rule.

Suppose $\theta$ can take values from two options, $\theta_1$ and $\theta_2$, and data can also take values from two options, $D_1$ and $D_2$.

We know that given $\theta_1$, there is one pathway to generate $D_1$, and given $\theta_2$ there are two pathways to generate $D_1$. If we make an observation that gives us $D_1$, we probably know that $\theta_2$ is more likely because it takes two of the three pathways to generate $D_1$. On the other hand, $\theta_1$ is less likely but not impossible.

So a reasonable estimate of how likely each value of $\theta$ could be, $P(\theta|D_1)$, can involve weighing by the number of pathways each value contributes. That is, $P(\theta_1|D_1) = 2 / (2 + 1) = 2/3$, and $P(\theta_2|D_1) = 1 / (2 + 1) = 1/3$. Essentially, the numerator, 1 or 2, is $P(D|\theta)$, and the denominator, (1+2), is $\int P(D|\theta) d\theta$ in the Bayes rule.

Coming back to our climate example, how should we calculate $P(D|\theta)$ and $P(\theta)$?

[Data Likelihood] Calculating $P(D|\theta)$ often involves components, a process-level model and an observation model .

• The process-level model is a model that generates true values of quantities to be observed. In our climate example, it is Eq. 1 that generates $T'_t$ given $\lambda$, $c_p$, and $F'_t$. Here, the subscript $_{t}$ denotes underlying true values to be observed. This model has an analytical solution, ###
  $T_t'(t) = \frac{1}{c_p} {\int_{0}^{t}} e^{-\frac{\lambda}{c_p}(t-\tau)} F'(\tau) d\tau$.
  where $t$ in the parenthesis denotes time. You may notice that we assume that forcing measurements, $F'$, are perfect, i.e., $F'_t$ = $F'$.
  This assumption is bold because the uncertainty of aerosol forcing estimate remains high and is an active field of research.

• The obervational model is a model that generates noisy observations. This is where probability is computed. In our climate example, this model could be $T' \sim N(T'_t,\sigma)$. That is, $T'$ in individual years independently follows an identical Gaussian distribution. As we will see in the later part of the notebook, uncertainties in these temperature estimates vary with time and are correlated across years .

[Prior] Calculating $P(\theta)$ depends on the prior we choose, which we will discuss later. But for now, we assume that we know nothing about $\lambda$ and $c_p$, such that their probability density function (pdf) is flat, i.e., $P(\lambda,c_p)$ = 1. We also assume the standard deviation of observational error of temperatures, $\sigma$, which is a parameter in the observational model, is 0.2°C.

Plugging these back into the Bayes rule, the posterior distribution is,

$P(\lambda, c_p, | T', F', \sigma = 0.2) \propto P(T' | \lambda, c_p, \sigma = 0.2, F') = \mathbf{normpdf}(T'_t, 0.2, T')$.

2.2 Perform the analysis ¶

Now we know how to compute the posterior, but essentially, we need to estimate, for example, the mean and uncertainty of $\lambda$. There are generally two types of methods to analyze the posterior distribution,

• Evaluating on a grid of parameters
• Sampling from the posterior distribution

2.2.1 Gird Search ¶

Let's start with the grid search method, which requires looping over every possible combination of parameters.

However, the ocean mixed layer depth is not as deep as 700m. Suppose we know that the mixed layer depth has a prior distribution, $P(c_p) \sim N(1,0.3) $hm. Let's bring in this peice of information and redo the grid search.

The prior information influences the estimate of our posterior distribution. For example, the new estimate does not have the long tail for a mixed layer deeper than 300 m, thereby excluding $\lambda$ smaller than 1.0 W/°C from our estimates.

The bayesian estimate is, therefore, an average over prior knowledge and data. In other words, Bayesian methods can be used to pull different pieces of information together. Several classic data assimilation techniques that we are going to cover in the class, such as the Kalman Filter, are posterior distributions averaged over model simulations (prior) and observations (data likelihood).

2.2.2 Sampling form the posterior distribution ¶

One disadvantage of the grid search is the number of required computations increases exponentially with the number of parameters, which can quickly become a nightmare that is often called the curse of dimensionality . Even for the 2D case above, the lower left and upper right corner of the parameter space has a low probability, and spending time to evaluate them is inefficient. Moreover, it is not an easy task to calculate, for example, mean and standard deviation, from grid search results. These quantities can be evaluated more efficiently if we can take samples from the posterior distribution.

Markov chain Monte Carlo (MCMC) methods comprise a class of algorithms for sampling from a probability distribution. In this notebook, we are going to introduce the Metropolis–Hastings algorithm , which obtains a sequence of random samples from a probability distribution from which direct sampling is difficult.

The mathematical details and proof of the algorithm can be found on this Wikipedia page.

Here, we focus on how to implement it in the climate model example.

The animation below gives you a step-by-step visualization of what happened when we run the above block.

The black dot denotes the proposal at a new time step, the gray line shows sampled chain. Also shown are probability ratio between the proposal and the existing location and whether the proposal is accepted.

Several points to notice:

• Usually, the first half of the chain is considered warm-up and is discarded.

• For the retained samples, need to check if they are stationary. That is, the starting and the end of a chain are exploring the same region in the parameter space.

• We also need to check for convergence. That is, if we start independent chains from different locations, do they explore the same region? To perform this check, we need to run several chains starting from different locations.

After warming up, the four chains converge visually. We can also use the ratio between total variance among all the chains versus mean variance within individual chains,

and $r > 1.02$ normally suggests poor convergence for a long chain.

3. Putting things together ¶

So far, we learned,

• Building a simple energy balance model to simulate the earth's surface temperature during the 20th century.

• Bayesian methods is to use data likelihood to update prior information, thereby pooling information from different sources.

• A model often consists of a deterministic process-level model that generates the true underlying quantities to be observed and an observational model that is used to evaluate data likelihood.

• Using grid search and MCMC methods to estimate parameters of model parameters in a Bayesian framework.

Can we further improve the above estimates?

3.1 Improve the error model ¶

One possibility is that the observational error, $\sigma$, is also an unknown parameter. The observational error generally decreases with time, meaning that earlier data would provide fewer constraints on the parameters. In addition, misfits between the model and data do not only arise from observational errors, $\sigma_o$. Internal climate variability, $\sigma^2_i$, such as ENSO, also leads to interannual variations that cannot be captured by the simple climate model.

As a result, a better error model would be,

$\sigma^2 = \sigma^2_i + \sigma^2_o$,

where we use estimates of $\sigma^2_o$ from HadCRUT5 and let the MCMC determine $\sigma^2_i$, we also use a parameter to remove the temperature change that is linearly associated with the Nino3.4 index .

3.2 Adding Another Mode Representing the Effect of Deeper Ocean ¶

From this analysis, it is clear that a simple model like Eq. 1 is capable of reproducing the overall temperature evolution at the earth's surface given historical estimates of external forcing.

Given the estimated $\lambda$ parameter, we can further estimate climate sensitivity as $F'/\lambda$. The effective radiative forcing of doubling CO ₂ is roughly 3.7 W/m ² . Plugging in posterior samples of $\lambda$ gives an estimate of climate sensitivity of 2.8 ([2.5, 3.4])°C. Such an estimate is tighter and, on average, lower than the 3 ([2, 5]) °C estimate by IPCC AR6.

Why that is the case?

In our simple model, we assumed that only the top ~100m of the ocean is responding to external forcing. But if a deeper ocean is also responding, it could increase the responding time scale of the system, thereby increase climate sensitivity estimates (the long tail in the fitting in section 2.2.1).

To represent the potential effect of a deeper ocean, we extend the simple model to have more than two modes, a fast model representing seasonal mixed-layer and a slower model representing responding ocean beneath the seasonal mixed-layer. Mathematically, we are adding up the solution of two modes.

$T_t'(t) = \frac{1}{c_{p1}} {\int_{0}^{t}} e^{-\frac{\lambda_1}{c_{p1}}(t-\tau)} F'(\tau) d\tau + \frac{1}{c_{p2}} {\int_{0}^{t}} e^{-\frac{\lambda_2}{c_{p2}}(t-\tau)} F'(\tau) d\tau$,

where subscripts 1 and 2 denote the fast and the slow mode, respectively. The fast mode should have a smaller $c_p$ and a larger $\lambda$, vice versa .

Now, we can repeat our Bayesian analysis,

Here, it does not necessarily mean that this updated result from the two box is 100% correct, but including the slow mode leads to higher ECS estimate and also a longer tail. We often call estimates using only fast modes transient climate sensitivity and those including slower modes equilibrium climate sensitivity .

3.3 A Proper Model-Data Comparison ¶

Look at the figure for fitted temperature in Sec. 3.2 again, although models capture the overall increasing trend, there are also misfits. As discussed before, some of them can arise from internal variability. For example, the two peaks in 1941 and 1998 can be explained by ENSO events. decadal variability, such as the Pacific decadal oscillation may have also contributed. But there are still some periods that internal variations cannot explain. For example, during World War II (WWII), observations show large positive anomalies, especially in 1945.

Of course, the model may be too simple, but when the model and data disagree, in addition to the model's being imperfect and internal variability could contribute, it is also important to keep in mind that the so-called observations we see could contain biases and errors. Below is a story about problems in sea-surface temperature measurements during WWII.

3.4 Further notes about Bayesian Analysis and MCMC ¶

Two other topics that may be of your interest:

(1) There are packages to implement MCMC-type Bayesian analysis. A famous one is Stan , which is a C package that also has command line , R , Python , Matlab , and Julia interfaces. Julia also has its own Bayesian package, called Turing . I find Turing powerful, especially when the process-level models are complicated and need to be solved using other packages, such as a differential equation solver .

(2) In practice, the Metropolis-Hastings algorithm often rejects too many proposals and performs poorly . A more involved method, called Hamiltonian Monte Carlo ( HMC ), solves this issues by considering momentum and gradient in the parameter space. I find this video a good starting point to grab some intuitions about HMC. Knowing Automatic differentiation and Dual number would also be helpful, but these would involve more advanced content in an Optimization course. Nevertheless, if your goal is to simply use HMC, both Stan and Turing has that option.