Simple model that gets you in spitting distance of a reasonable estimate of Equilibrium Climate Sensitivity

Background

ECS - Equilibrium Climate Sensitivity - is the temperature change on Earth as a result of a doubling of atmospheric CO2.  That means going from 280 ppm (the level prior to the industrial revolution) past 420 (where we are now) and on to 560 ppm.  Most people probably understand why this would increase temperatures by now, but I'm going to repeat the argument briefly anyway, for completeness

Earth is heated by the sun's rays and its surface warms up.  As it does it emits energy in the infrared spectrum as shown in the bottom right of the above figure.  The total energy emitted per square metre is, according to the Stefan-Boltzmann law, $\sigma T^4$, where $\sigma = 5.67\times 10^{-8} Wm^{-2}K^{-4}$.

In equilbrium the total energy coming in from the sun must equal the total energy going out.  The energy arriving at the top of the atmosphere per square metre facing the sun, $S$ is 1360 $Wm^{-2}$.  But we can ignore the fraction of this that is reflected, known as the albedo $\alpha$, equal to about 0.3.  So we're left with total energy coming in

$$
S(1-\alpha) \pi R^2
$$

over an area of $4\pi R^2$, giving an average input from the sun, per square metre, of

$$
\frac{S(1-\alpha)}{4} \approx 238 Wm^{-2}
$$

If you set this equal to $\sigma T^4$ and solve, you get 255K which is about 33C lower than the true average for Earth, namely 288K.  The reason for the discrepancy is of course the atmosphere.  The figure above shows how, at each wavelength radiated by Earth, the different components of the atmosphere absorb between 0 and 100% of the outgoing radiation.  One of those components is CO2 and, like the others it re-radiates half up and half back down towards Earth, warming it up like a blanket.  We're going to come up with a basic model for how the temperature of Earth varies in the atmosphere and use this to estimate the ECS.

The word "equilibrium" does a lot of heavy lifting in the term ECS.  One might assume that it means you double the CO2, allow your model Earth to settle for $\infty$ years, and then measure the temperature change.  It doesn't mean that for a couple of reasons.  The first reason is that it would involve modelling long term changes that we're not good at modelling, such as feedbacks from melting permafrost (which releases greenhouse gases), and receding ice sheets (which lowers the albedo).  The second reason is that we're a selfish and short sighted species and really only care about our lives, our children's, and - at a push - our grandchildren's.  So, the "Equilibrium" in "Equilibrium Climate Sensitivity" means we allow the atmosphere to settle into a new regime, and also the sea surface (including the sea ice), but we do not wait long enough for Greenland to melt or the forests to turn to savanna.  But our model is only intended to get us to spitting distance, so none those details really matter to us.

The model

We're going to make some assumptions in our model

1. The incoming sunlight is not absorbed by the atmosphere at all (this is a bad approximation, only partially justified by the graphic above)
2. The atmosphere is stratified and static - i.e. there is no convection or latent heat transfer caused by water vapour condensing in clouds.  The only vertical energy transfer is caused by infrared emissions by the ground or the atmosphere.  (We will come back to this later on.)
3. Sunlight falls equally everywhere at all times.
   
4. There are no feedbacks, i.e. no increased water vapour caused by increases in sea surface temperature, no change in cloud cover or sea ice extent altering the albedo, other than that...
   
5. Every atmospheric layer is at equilibrium, i.e. not warming up or cooling down.
   

Yup, it's going to be pretty basic.

First off we need to introduce the concept of optical depth, $\tau$.  This is a property of a block of some material, and has the definition

$$
\tau = ln\left(\frac{F_{in}}{F_{out}}\right)
$$

where $F_{in}$ is the flux of radiative energy going in, and $F_{out}$ is the flux coming out, not including any energy emitted by the material itself due to its temperature.  This definition is not complete without specifying either a specific wavelength $\lambda$ or an intensity distribution.  From here onwards we're going to assume the mostly infrared distribution of light emitted by the surface of Earth and its atmosphere.

$\tau$ seems an odd definition but it has the advantage that if you have two materials back to back with depths $\tau_1$ and $\tau_2$, e.g. horizontal layers of the atmosphere, then the overall depth is just $\tau_1 + \tau_2$.  Neat.  Let's consider a block of material at temperature $T$, with one surface and infinite optical depth.  Each thin layer emits an amount of radiation towards the surface that is proportional to its thickness $d\tau$, let's call that amount $d\tau E$.   If this layer is at optical depth $\tau$ then the energy reaching the surface is $d\tau E e^{-\tau}$.  The total amount reaching the surface is the flux predicted by Stefan-Boltzmann, i.e.

$$
\sigma T^4 = \int_0^{\infty} E e^{-\tau} d\tau
$$

This is pretty easy to solve and gives us simply that $E=\sigma T^4$.  Now let's consider a block of material of optical depth $\tau$ instead of $\infty$.  Using our new value of $E$ we have that the flux reaching the surface must be

$$
\int_0^{\tau} \sigma T^4 e^{-\tau'} d\tau' = \sigma T^4 (1 - e^{-\tau})
$$

Now a few definitions:

• $\tau$ = the optical depth between the top of atmosphere and a given atmospheric layer, relative to the infrared spectrum emitted by Earth and its atmosphere
  
• $\tau_0$ = the total optical depth of the atmosphere (i.e. $\tau$ at $z=0$)
• $T_0$ = the temperature at $z=0$
  
• $F^+(\tau)$  = the infrared flux going up at depth $\tau$
• $F^-(\tau)$  = the infrared flux going down at depth $\tau$
• $F_{net} = F^+ - F^-$ the net infrared flux going out to space

We know that $F_{net}$ must equal $S(1-\alpha)/4$ everywhere, as the downward flux of visible light must be matched for the layers below to be in equilibrium.  We can also write equations for $F^+$, and $F^-$ using the formulas we derived earlier

$$
\begin{align}
F^+(\tau) &= \int_{\tau}^{\tau_0}\sigma T^4 e^{-(\tau'-\tau)}d\tau' + \sigma T_0^4 e^{-(\tau_0-\tau)}\\
F^-(\tau) &= \int_{0}^{\tau}\sigma T^4 e^{-(\tau-\tau')}d\tau'
\end{align}
$$

If we agree to define $T=T_0$ for $\tau > \tau_0$ we can simplify $F^+$ to  

$$
F^+(\tau) = \int_{\tau}^{\infty}\sigma T^4 e^{-(\tau'-\tau)}d\tau'
$$

Differentiating $F_{net}$ twice using simple calculus gives us

$$
\frac{d^2F_{net}}{d\tau^2} = \frac{d(-2\sigma T^4)}{d\tau} + F_{net}
$$

We already noted that $F_{net}$ is constant, so we can integrate to get

$$
F_{net}\tau + C = 2\sigma T^4
$$

The following equation then holds because it's just the same thing done to both sides

$$
\int_0^{\infty} F_{net}\tau e^{-\tau}d\tau + \int_0^{\infty} Ce^{-\tau}d\tau \\
= \int_0^{\infty} 2\sigma T^4 e^{-\tau}d\tau
$$

The RHS is just $2F^+(0)$, or alternatively $2F_{net}$.  The LHS reduces to $F_{net} + C$ and therefore $C = F_{net}$.  Now that we know $C$ we can rearrange to get an equation for $T$ valid between $0 < \tau < \tau_0$

$$
T = \sqrt[4]{F_{net}\frac{\tau+1}{2\sigma}}
$$

Knowing that $F_{net} = 238 Wm^{-2}$, and that $T_0=288K$ allows us to calculate that $\tau_0 = 2.3$.  This corresponds to ~90% of the infrared absorbed which, if you look at the line marked "Total" in the plot at the top, looks pretty reasonable.  Now let's do something incredibly handwavy and try to guess roughly how much of the infrared absorption is due to just the CO2 (I did say spitting distance).  Looking at the plot again, I'd say it's about 20%.  That means that the optical depth contributed by the current level of CO2 is

$$
ln\left(\frac{1}{0.8}\right) \approx 0.2
$$

Let's divide this by 1.5 to get the pre-industrial value and then add that to $\tau_0$ to represent a doubling of pre-industrial CO2.  This then predicts $T_0=291.7$ i.e. an ECS of 3.7C.

Visualizations

What does our curve look like then?

Our model

This is great but we'd prefer a plot of $T$ vs $z$.  We can use the ideal gas law to make an approximation that pressure is proportional to density i.e. $p = \frac{p_0}{\rho_0}\rho$. (See note 1.)  But pressure is also just the weight of the column above, i.e. $p = g \int_z^\infty \rho dz$.  These can only both be true if $\rho = \rho_0e^{-\frac{\rho_0 g z}{p_0}}$.  If we assume a homogenous atmosphere then $d\tau \propto \rho dz$ which we can solve to get $\tau = \tau_0 e^{-\frac{\rho_0 g z}{p_0}}$.  This can be substituted into our equation for $T$ to give us a picture like this

Same model, but plotted vs z

Model improvements

Of course, there's something wrong with our reasoning.  Whilst we got an ECS in the right ballpark, it's 2 or 3C  too high for an estimate that doesn't include any short term feedbacks like an increase in water vapour or loss of sea ice.  The problem is that, although CO2 is well mixed in the atmosphere, other greenhouse gases like water vapour are not.  This means that a disproportionate amount of the optical thickness is distributed in the lower part of the atmosphere.  This in turn means that, in our model, the temperature gradient in the lower part of the atmosphere would be greater than the "adiabatic lapse rate", the maximum gradient before convection occurs (which turns out to be $g/c_p$ or about $10^\circ C/km$).  So a better model would determine the height of the convection layer $z_1$, set the temperature difference $T_1-T_0$ to be proportional to that, and then work out $T_1$ using the method above.  If the CO2 is doubled this would increase $z_1$ as well as $\tau_1$.

There are other improvements too.  For example, we could use the fact that some sunlight is in fact absorbed by lighter molecules in the upper atmosphere.   You can go on tweaking for eternity, but it's nice to know we can get within an order of magnitude using a simple model.  The CMIP6 models used by the IPCC for their latest assessment report give a range of 1.83C to 5.64C for the ECS, so let's not beat ourselves up too much.

So what's the actual ECS?

In reality we have next to no idea.  Contemporary estimates range enormously from values below 1C to above 10C, but a lot of studies end up giving central 66% ranges of around 1.5-4.5C, just like in the 1979 Charney report.  Given in 2021 we're at 1.1C for only a 50% increase, and the mostly positive feedbacks haven't had time to do their work, it seems very unlikely that the lowest estimates from the past are correct.  But as for an upper bound: who knows!  And there are some, like Steffen et al. who assert that ECS isn't even particularly meaningful.  In their perspective the climate system has "limit cycles" like the current ~100kyr cycles between ice ages.  Earth entered this cycle when the Himalayas went up and the increased weathering drew a lot of CO2 out of the air, bringing temperatures low enough to make that trajectory feasible.  In this model, if temperatures rise above the highest interglacial values (i.e pre-industrial+2C) the climate could spin off into a completely different limit cycle, perhaps never to return.  Maybe this isn't backed up by hard maths, but then the hard maths used in most quantitative studies employs some seriously dodgy assumptions (like linearity).  And given the huge error bars and incompatible estimates, most studies are either unhelpful or wrong.  What we do know is that the lower bounds look bad enough to justify reorganising society around a low carbon future.

NOTES

1. This works really well, but you can still solve for $T$ in terms of $z$ if you don't want to make this approximation, it's just a bit fiddly.  See here for how I did it.
   

REFERENCES

1. Andrew Dessler, 2021, Introduction to modern climate change, 3rd ed.
2. Tyler D. Robinson and David C. Catling, 2012, AN ANALYTIC RADIATIVE–CONVECTIVE MODEL FOR PLANETARY ATMOSPHERES
3. Knutti et al., 2017, Beyond equilibrium climate sensitivity
4. Steffen et al., 2018, Trajectories of the Earth System in the Anthropocene