Brian E. J. Rose, University at Albany
This document uses the interactive Jupyter notebook format. The notes can be accessed in several different ways:
- The interactive notebooks are hosted on
githubat https://github.com/brian-rose/ClimateModeling_courseware - The latest versions can be viewed as static web pages rendered on nbviewer
- A complete snapshot of the notes as of May 2017 (end of spring semester) are available on Brian's website.
Also here is a legacy version from 2015.
Many of these notes make use of the climlab package, available at https://github.com/brian-rose/climlab
- A single layer atmosphere
- Introducing the two-layer grey gas model
- Tuning the grey gas model to observations
- Level of emission
- Radiative forcing in the 2-layer grey gas model
- Radiative equilibrium in the 2-layer grey gas model
- Summary
We will make our first attempt at quantifying the greenhouse effect in the simplest possible greenhouse model: a single layer of atmosphere that is able to absorb and emit longwave radiation.
from IPython.display import Image
Image('../images/1layerAtm_sketch.png')
- Atmosphere is a single layer of air at temperature
$T_a$ - Atmosphere is completely transparent to shortwave solar radiation.
- The surface absorbs shortwave radiation
$(1-\alpha) Q$ - Atmosphere is completely opaque to infrared radiation
- Both surface and atmosphere emit radiation as blackbodies (
$\sigma T_s^4, \sigma T_a^4$ ) - Atmosphere radiates equally up and down (
$\sigma T_a^4$ ) - There are no other heat transfer mechanisms
We can now use the concept of energy balance to ask what the temperature need to be in order to balance the energy budgets at the surface and the atmosphere, i.e. the radiative equilibrium temperatures.
\begin{align} \text{energy in} &= \text{energy out} \ (1-\alpha) Q + \sigma T_a^4 &= \sigma T_s^4 \ \end{align}
The presence of the atmosphere above means there is an additional source term: downwelling infrared radiation from the atmosphere.
We call this the back radiation.
\begin{align} \text{energy in} &= \text{energy out} \ \sigma T_s^4 &= A\uparrow + A\downarrow = 2 \sigma T_a^4 \ \end{align}
which means that $$ T_s = 2^\frac{1}{4} T_a \approx 1.2 T_a $$
So we have just determined that, in order to have a purely radiative equilibrium, we must have
The surface must be warmer than the atmosphere.
Now plug this into the surface equation to find
and use the definition of the emission temperature
In fact, in this model, $T_e$ is identical to the atmospheric temperature $T_a$, since all the OLR originates from this layer.
Solve for the surface temperature: $$ T_s = 2^\frac{1}{4} T_e $$
Putting in observed numbers,
This model is one small step closer to reality: surface is warmer than atmosphere, emissions to space generated in the atmosphere, atmosphere heated from below and helping to keep surface warm.
BUT our model now overpredicts the surface temperature by about 15ºC (or K).
Ideas about why?
Basically we just need to read our list of assumptions above and realize that none of them are very good approximations:
- Atmosphere absorbs some solar radiation.
- Atmosphere is NOT a perfect absorber of longwave radiation
- Absorption and emission varies strongly with wavelength (atmosphere does not behave like a blackbody).
- Emissions are not determined by a single temperature
$T_a$ but by the detailed vertical profile of air temperture. - Energy is redistributed in the vertical by a variety of dynamical transport mechanisms (e.g. convection and boundary layer turbulence).
Let's generalize the above model just a little bit to build a slighly more realistic model of longwave radiative transfer.
We will address two shortcomings of our single-layer model:
- No vertical structure
- 100% longwave opacity
Relaxing these two assumptions gives us what turns out to be a very useful prototype model for understanding how the greenhouse effect works.
- The atmosphere is transparent to shortwave radiation (still)
- Divide the atmosphere up into two layers of equal mass (the dividing line is thus at 500 hPa pressure level)
- Each layer **absorbs only a fraction
$\epsilon$ ** of whatever longwave radiation is incident upon it. - We will call the fraction
$\epsilon$ the absorptivity of the layer. - Assume
$\epsilon$ is the same in each layer
This is called the grey gas model, where grey here means the emission and absorption have no spectral dependence.
We can think of this model informally as a "leaky greenhouse".
Note that the assumption that
Out of our two most important absorbers:
- CO$_2$ is well mixed
- H$_2$O is not (mostly confined to lower troposphere due to strong temperature dependence of the saturation vapor pressure).
But we will ignore this aspect of reality for now.
In order to build our model, we need to introduce one additional piece of physics known as Kirchoff's Law:
So if a layer of atmosphere at temperature
both up and down.
Image('../images/2layerAtm_sketch.png')
- Surface temperature is
$T_s$ - Atm. temperatures are
$T_0, T_1$ where$T_0$ is closest to the surface. - absorptivity of atm layers is
$\epsilon$ - Surface emission is
$\sigma T_s^4$ - Atm emission is
$\epsilon \sigma T_0^4, \epsilon \sigma T_1^4$ (up and down) - Absorptivity = emissivity for atmospheric layers
- a fraction
$(1-\epsilon)$ of the longwave beam is transmitted through each layer
This two-layer grey gas model is simple enough that we can work out all the details algebraically. There are three temperatures to keep track of
We all know how to work these things out with pencil and paper. But it can be tedious and error-prone.
Symbolic math software lets us use the computer to automate a lot of tedious algebra.
The sympy package is a powerful open-source symbolic math library that is well-integrated into the scientific Python ecosystem.
import sympy
# Allow sympy to produce nice looking equations as output
sympy.init_printing()
# Define some symbols for mathematical quantities
# Assume all quantities are positive (which will help simplify some expressions)
epsilon, T_e, T_s, T_0, T_1, sigma = \
sympy.symbols('epsilon, T_e, T_s, T_0, T_1, sigma', positive=True)
# So far we have just defined some symbols, e.g.
T_s# We have hard-coded the assumption that the temperature is positive
sympy.ask(T_s>0)True
Let's denote the emissions from each layer as \begin{align} E_s &= \sigma T_s^4 \ E_0 &= \epsilon \sigma T_0^4 \ E_1 &= \epsilon \sigma T_1^4 \end{align}
recognizing that
# Define these operations as sympy symbols
# And display as a column vector:
E_s = sigma*T_s**4
E_0 = epsilon*sigma*T_0**4
E_1 = epsilon*sigma*T_1**4
E = sympy.Matrix([E_s, E_0, E_1])
ESince we have assumed the atmosphere is transparent to shortwave, the incident beam
# Define some new symbols for shortwave radiation
Q, alpha = sympy.symbols('Q, alpha', positive=True)
# Create a dictionary to hold our numerical values
tuned = {}
tuned[Q] = 341.3 # global mean insolation in W/m2
tuned[alpha] = 101.9/Q.subs(tuned) # observed planetary albedo
tuned[sigma] = 5.67E-8 # Stefan-Boltzmann constant in W/m2/K4
tuned
# Numerical value for emission temperature
#T_e.subs(tuned)Let
The upward flux from the surface to layer 0 is $$ U_0 = E_s $$ (just the emission from the suface).
U_0 = E_s
U_0Following this beam upward, we can write the upward flux from layer 0 to layer 1 as the sum of the transmitted component that originated below layer 0 and the new emissions from layer 0:
U_1 = (1-epsilon)*U_0 + E_0
U_1Continuing to follow the same beam, the upwelling flux above layer 1 is $$ U_2 = (1-\epsilon) U_1 + E_1 $$
U_2 = (1-epsilon) * U_1 + E_1Since there is no more atmosphere above layer 1, this upwelling flux is our Outgoing Longwave Radiation for this model:
U_2The three terms in the above expression represent the contributions to the total OLR that originate from each of the three levels.
Let's code this up explicitly for future reference:
# Define the contributions to OLR originating from each level
OLR_s = (1-epsilon)**2 *sigma*T_s**4
OLR_0 = epsilon*(1-epsilon)*sigma*T_0**4
OLR_1 = epsilon*sigma*T_1**4
OLR = OLR_s + OLR_0 + OLR_1
print 'The expression for OLR is'
OLRThe expression for OLR is
Let
fromspace = 0
D_2 = fromspaceBetween layer 1 and layer 0 the beam contains emissions from layer 1:
D_1 = (1-epsilon)*D_2 + E_1
D_1Finally between layer 0 and the surface the beam contains a transmitted component and the emissions from layer 0:
D_0 = (1-epsilon)*D_1 + E_0
D_0This
In building our new model we have introduced exactly one parameter, the absorptivity
We will tune our model so that it reproduces the observed global mean OLR given observed global mean temperatures.
To get appropriate temperatures for
First, we set
From the lapse rate plot, an average temperature for the layer between 1000 and 500 hPa is
Defining an average temperature for the layer between 500 and 0 hPa is more ambiguous because of the lapse rate reversal at the tropopause. We will choose
From the graph, this is approximately the observed global mean temperature at 275 hPa or about 10 km.
# add to our dictionary of values:
tuned[T_s] = 288.
tuned[T_0] = 275.
tuned[T_1] = 230.
tunedFrom the [observed global energy budget](Lecture01 -- Planetary energy budget.ipynb) we set
We wrote down the expression for OLR as a function of temperatures and absorptivity in our model above.
We just need to equate this to the observed value and solve a quadratic equation for
This is where the real power of the symbolic math toolkit comes in.
Subsitute in the numerical values we are interested in:
# the .subs() method for a sympy symbol means
# substitute values in the expression using the supplied dictionary
# Here we use observed values of Ts, T0, T1
OLR2 = OLR.subs(tuned)
OLR2We have a quadratic equation for
Now use the sympy.solve function to solve the quadratic:
# The sympy.solve method takes an expression equal to zero
# So in this case we subtract the tuned value of OLR from our expression
eps_solution = sympy.solve(OLR2 - 238.5, epsilon)
eps_solutionThere are two roots, but the second one is unphysical since we must have
Just for fun, here is a simple of example of filtering a list using powerful Python list comprehension syntax:
# Give me only the roots that are between zero and 1!
list_result = [eps for eps in eps_solution if 0<eps<1]
print list_result
# The result is a list with a single element.
# We need to slice the list to get just the number:
eps_tuned = list_result[0]
print eps_tuned[0.586041150248834]
0.586041150248834
We conclude that our tuned value is
This is the absorptivity that guarantees that our model reproduces the observed OLR given the observed tempertures.
tuned[epsilon] = eps_tuned
tunedEven in this very simple greenhouse model, there is no single level at which the OLR is generated.
The three terms in our formula for OLR tell us the contributions from each level.
OLRterms = sympy.Matrix([OLR_s, OLR_0, OLR_1])
OLRtermsNow evaluate these expressions for our tuned temperature and absorptivity:
OLRtuned = OLRterms.subs(tuned)
OLRtunedSo we are getting about 67 W m$^{-2}$ from the surface, 79 W m$^{-2}$ from layer 0, and 93 W m$^{-2}$ from the top layer.
In terms of fractional contributions to the total OLR, we have (limiting the output to two decimal places):
sympy.N(OLRtuned / 239., 2)Notice that the largest single contribution is coming from the top layer. This is in spite of the fact that the emissions from this layer are weak, because it is so cold.
Comparing to observations, the actual contribution to OLR from the surface is about 22 W m$^{-2}$ (or about 9% of the total), not 67 W m$^{-2}$. So we certainly don't have all the details worked out yet!
As we will see later, to really understand what sets that observed 22 W m$^{-2}$, we will need to start thinking about the spectral dependence of the longwave absorptivity.
Adding some extra greenhouse absorbers will mean that a greater fraction of incident longwave radiation is absorbed in each layer.
Thus
Suppose we have
Suppose further that this increase happens abruptly so that there is no time for the temperatures to respond to this change. We hold the temperatures fixed in the column and ask how the radiative fluxes change.
Do you expect the OLR to increase or decrease?
Let's use our two-layer leaky greenhouse model to investigate the answer.
The components of the OLR before the perturbation are
OLRtermsAfter the perturbation we have
delta_epsilon = sympy.symbols('delta_epsilon')
OLRterms_pert = OLRterms.subs(epsilon, epsilon+delta_epsilon)
OLRterms_pertLet's take the difference
deltaOLR = OLRterms_pert - OLRterms
deltaOLRTo make things simpler, we will neglect the terms in
Telling sympy to set the quadratic terms to zero gives us
deltaOLR_linear = sympy.expand(deltaOLR).subs(delta_epsilon**2, 0)
deltaOLR_linearRecall that the three terms are the contributions to the OLR from the three different levels. In this case, the changes in those contributions after adding more absorbers.
Now let's divide through by
deltaOLR_per_deltaepsilon = \
sympy.simplify(deltaOLR_linear / delta_epsilon)
deltaOLR_per_deltaepsilonNow look at the sign of each term. Recall that
THIS IS VERY IMPORTANT, SO STOP AND THINK ABOUT IT.
The contribution from the surface must decrease, while the contribution from the top layer must increase.
When we add absorbers, the average level of emission goes up!
In this model, only the longwave flux can change, so we define the radiative forcing as
(with the minus sign so that
We just worked out that whenever we add some extra absorbers, the emissions to space (on average) will originate from higher levels in the atmosphere.
What does this mean for OLR? Will it increase or decrease?
To get the answer, we just have to sum up the three contributions we wrote above:
R = -sum(deltaOLR_per_deltaepsilon)
RIs this a positive or negative number? The key point is this:
It depends on the temperatures, i.e. on the lapse rate.
Stop and think about this question:
If the surface and atmosphere are all at the same temperature, does the OLR go up or down when
Understanding this question is key to understanding how the greenhouse effect works.
We will just set
R.subs([(T_0, T_s), (T_1, T_s)])which then simplifies to
sympy.simplify(R.subs([(T_0, T_s), (T_1, T_s)]))For an isothermal atmosphere, there is no change in OLR when we add extra greenhouse absorbers. Hence, no radiative forcing and no greenhouse effect.
Why?
The level of emission still must go up. But since the temperature at the upper level is the same as everywhere else, the emissions are exactly the same.
For a more realistic example of radiative forcing due to an increase in greenhouse absorbers, we can substitute in our tuned values for temperature and
We'll express the answer in W m$^{-2}$ for a 1% increase in
The three components of the OLR change are
deltaOLR_per_deltaepsilon.subs(tuned) * 0.01And the net radiative forcing is
R.subs(tuned) * 0.01So in our example, the OLR decreases by 2.2 W m$^{-2}$, or equivalently, the radiative forcing is +2.2 W m$^{-2}$.
What we have just calculated is this:
Given the observed lapse rates, a small increase in absorbers will cause a small decrease in OLR.
The greenhouse effect thus gets stronger, and energy will begin to accumulate in the system -- which will eventually cause temperatures to increase as the system adjusts to a new equilibrium.
In the previous section we:
- made no assumptions about the processes that actually set the temperatures.
- used the model to calculate radiative fluxes, given observed temperatures.
- stressed the importance of knowing the lapse rates in order to know how an increase in emission level would affect the OLR, and thus determine the radiative forcing.
A key question in climate dynamics is therefore this:
What sets the lapse rate?
It turns out that lots of different physical processes contribute to setting the lapse rate.
Understanding how these processes acts together and how they change as the climate changes is one of the key reasons for which we need more complex climate models.
For now, we will use our prototype greenhouse model to do the most basic lapse rate calculation: the radiative equilibrium temperature.
We assume that
- the only exchange of energy between layers is longwave radiation
- equilibrium is achieved when the net radiative flux convergence in each layer is zero.
First, the net upwelling flux is just the difference between flux up and flux down:
# Upwelling and downwelling beams as matrices
U = sympy.Matrix([U_0, U_1, U_2])
D = sympy.Matrix([D_0, D_1, D_2])
# Net flux, positive up
F = U-D
F(difference between what's coming in the bottom and what's going out the top of each layer)
# define a vector of absorbed radiation -- same size as emissions
A = E.copy()
# absorbed radiation at surface
A[0] = F[0]
# Compute the convergence
for n in range(2):
A[n+1] = -(F[n+1]-F[n])
AThe only other heat source is the shortwave heating at the surface.
In matrix form, here is the system of equations to be solved:
radeq = sympy.Equality(A, sympy.Matrix([(1-alpha)*Q, 0, 0]))
radeqJust as we did for the 1-layer model, it is helpful to rewrite this system using the definition of the emission temperture
radeq2 = radeq.subs([((1-alpha)*Q, sigma*T_e**4)])
radeq2In this form we can see that we actually have a linear system of equations for a set of variables
We can solve this matrix problem to get these as functions of
# Solve for radiative equilibrium
fourthpower = sympy.solve(radeq2, [T_s**4, T_1**4, T_0**4])
fourthpowerThis produces a dictionary of solutions for the fourth power of the temperatures!
A little manipulation gets us the solutions for temperatures that we want:
# need the symbolic fourth root operation
from sympy.simplify.simplify import nthroot
fourthpower_list = [fourthpower[key] for key in [T_s**4, T_0**4, T_1**4]]
solution = sympy.Matrix([nthroot(item,4) for item in fourthpower_list])
# Display result as matrix equation!
T = sympy.Matrix([T_s, T_0, T_1])
sympy.Equality(T, solution)In more familiar notation, the radiative equilibrium solution is thus
\begin{align} T_s &= T_e \left( \frac{2+\epsilon}{2-\epsilon} \right)^{1/4} \ T_0 &= T_e \left( \frac{1+\epsilon}{2-\epsilon} \right)^{1/4} \ T_1 &= T_e \left( \frac{ 1}{2 - \epsilon} \right)^{1/4} \end{align}
Plugging in the tuned value
Tsolution = solution.subs(tuned)
# Display result as matrix equation!
sympy.Equality(T, Tsolution)Now we just need to know the Earth's emission temperature
(Which we already know is about 255 K)
# Here's how to calculate T_e from the observed values
sympy.solve(((1-alpha)*Q - sigma*T_e**4).subs(tuned), T_e)# Need to unpack the list
Te_value = sympy.solve(((1-alpha)*Q - sigma*T_e**4).subs(tuned), T_e)[0]
Te_value# Output 4 significant digits
Trad = sympy.N(Tsolution.subs([(T_e, Te_value)]), 4)
sympy.Equality(T, Trad)Compare these to the values we derived from the observed lapse rates:
sympy.Equality(T, T.subs(tuned))The radiative equilibrium solution is substantially warmer at the surface and colder in the lower troposphere than reality.
This is a very general feature of radiative equilibrium, and we will see it again very soon in this course.
-
Putting a layer of longwave absorbers above the surface keeps the surface substantially warmer, because of the backradiation from the atmosphere (greenhouse effect).
-
The grey gas model assumes that each layer absorbs and emits a fraction
$\epsilon$ of its blackbody value, independent of wavelength. -
With incomplete absorption (
$\epsilon < 1$ ), there are contributions to the OLR from every level and the surface (there is no single level of emission) -
Adding more absorbers means that contributions to the OLR from upper levels go up, while contributions from the surface go down.
-
This upward shift in the weighting of different levels is what we mean when we say the level of emission goes up.
-
The radiative forcing caused by an increase in absorbers depends on the lapse rate.
-
For an isothermal atmosphere the radiative forcing is zero and there is no greenhouse effect
-
The radiative forcing is positive for our atmosphere because tropospheric temperatures tends to decrease with height.
-
Pure radiative equilibrium produces a warm surface and cold lower troposphere.
-
This is unrealistic, and suggests that crucial heat transfer mechanisms are missing from our model.
Did we need sympy to work all this out? No, of course not. We could have solved the 3x3 matrix problems by hand. But computer algebra can be very useful and save you a lot of time and error, so it's good to invest some effort into learning how to use it.
Hopefully these notes provide a useful starting point.
You are now ready to tackle [Assignment 5](../Assignments/Assignment05 -- Radiative forcing in a grey radiation atmosphere.ipynb), where you are asked to extend this grey-gas analysis to many layers.
For more than a few layers, the analytical approach we used here is no longer very useful. You will code up a numerical solution to calculate OLR given temperatures and absorptivity, and look at how the lapse rate determines radiative forcing for a given increase in absorptivity.
%load_ext version_information
%version_information sympy| Software | Version |
|---|---|
| Python | 2.7.12 64bit [GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2336.11.00)] |
| IPython | 5.3.0 |
| OS | Darwin 16.5.0 x86_64 i386 64bit |
| sympy | 1.0 |
| Thu May 25 13:44:14 2017 EDT | |
The author of this notebook is Brian E. J. Rose, University at Albany.
It was developed in support of ATM 623: Climate Modeling, a graduate-level course in the Department of Atmospheric and Envionmental Sciences
Development of these notes and the climlab software is partially supported by the National Science Foundation under award AGS-1455071 to Brian Rose. Any opinions, findings, conclusions or recommendations expressed here are mine and do not necessarily reflect the views of the National Science Foundation.