Brian E. J. Rose, University at Albany
First, some thoughts on modeling from xkcd
Let's look again at the observations:
Last class we introduced a very simple model for the OLR or Outgoing Longwave Radiation to space:
$$ \text{OLR} = \tau \sigma T_s^4 $$where $\tau$ is the transmissivity of the atmosphere, a number less than 1 that represents the greenhouse effect of Earth's atmosphere.
We also tuned this model to the observations by choosing $ \tau \approx 0.61$.
More precisely:
OLRobserved = 238.5 # in W/m2
sigma = 5.67E-8 # S-B constant
Tsobserved = 288. # global average surface temperature
tau = OLRobserved / sigma / Tsobserved**4 # solve for tuned value of transmissivity
print(tau)
Let's now deal with the shortwave (solar) side of the energy budget.
Let's define a few terms.
From the observations, the area-averaged incoming solar radiation or insolation is 341.3 W m$^{-2}$.
Let's denote this quantity by $Q$.
Q = 341.3 # the insolation
Some of the incoming radiation is not absorbed at all but simply reflected back to space. Let's call this quantity $F_{reflected}$
From observations we have:
Freflected = 101.9 # reflected shortwave flux in W/m2
The planetary albedo is the fraction of $Q$ that is reflected.
We will denote the planetary albedo by $\alpha$.
From the observations:
alpha = Freflected / Q
print(alpha)
That is, about 30% of the incoming radiation is reflected back to space.
The Absorbed Shortwave Radiation or ASR is the part of the incoming sunlight that is not reflected back to space, i.e. that part that is absorbed somewhere within the Earth system.
Mathematically we write
$$ \text{ASR} = Q - F_{reflected} = (1-\alpha) Q $$From the observations:
ASRobserved = Q - Freflected
print(ASRobserved)
As we noted last time, this number is just slightly greater than the observed OLR of 238.5 W m$^{-2}$.
This is one of the central concepts in climate modeling.
The Earth system is in energy balance when energy in = energy out, i.e. when
$$ \text{ASR} = \text{OLR} $$We want to know:
With our simple greenhouse model, we can get an exact solution for the equilibrium temperature.
First, write down our statement of energy balance:
$$ (1-\alpha) Q = \tau \sigma T_s^4 $$Rearrange to solve for $T_s$:
$$ T_s^4 = \frac{(1-\alpha) Q}{\tau \sigma} $$and take the fourth root, denoting our equilibrium temperature as $T_{eq}$:
$$ T_{eq} = \left( \frac{(1-\alpha) Q}{\tau \sigma} \right)^\frac{1}{4} $$Plugging the observed values back in, we compute:
# define a reusable function!
def equilibrium_temperature(alpha,Q,tau):
return ((1-alpha)*Q/(tau*sigma))**(1/4)
Teq_observed = equilibrium_temperature(alpha,Q,tau)
print(Teq_observed)
And this equilibrium temperature is just slightly warmer than 288 K. Why?
Suppose that, due to global warming (changes in atmospheric composition and subsequent changes in cloudiness):
What is the new equilibrium temperature?
For this very simple model, we can work out the answer exactly:
Teq_new = equilibrium_temperature(0.32,Q,0.57)
# an example of formatted print output, limiting to two or one decimal places
print('The new equilibrium temperature is {:.2f} K.'.format(Teq_new))
print('The equilibrium temperature increased by about {:.1f} K.'.format(Teq_new-Teq_observed))
Most climate models are more complicated mathematically, and solving directly for the equilibrium temperature will not be possible!
Instead, we will be able to use the model to calculate the terms in the energy budget (ASR and OLR).
The above exercise shows us that if some properties of the climate system change in such a way that the equilibrium temperature goes up, then the Earth system receives more energy from the sun than it is losing to space. The system is no longer in energy balance.
The temperature must then increase to get back into balance. The increase will not happen all at once! It will take time for energy to accumulate in the climate system. We want to model this time-dependent adjustment of the system.
In fact almost all climate models are time-dependent, meaning the model calculates time derivatives (rates of change) of climate variables.
We will write the total energy budget of the Earth system as
$$ C \frac{dT_s}{dt} = \text{ASR} - \text{OLR} $$where
By adopting this equation, we are assuming that the energy content of the Earth system (atmosphere, ocean, ice, etc.) is proportional to surface temperature.
Important things to think about:
For our purposes here we are going to use a value of C equivalent to heating 100 meters of water:
c_w = 4E3 # Specific heat of water in J/kg/K
rho_w = 1E3 # Density of water in kg/m3
H = 100. # Depth of water in m
C = c_w * rho_w * H # Heat capacity of the model
print('The effective heat capacity is {:.1e} J/m2/K'.format(C))
Recall that the derivative is the instantaneous rate of change. It is defined as
$$ \frac{dT}{dt} = \lim_{\Delta t\rightarrow 0} \frac{\Delta T}{\Delta t}$$So we write our model as
$$ C \frac{\Delta T}{\Delta t} \approx \text{ASR} - \text{OLR}$$where $\Delta T$ is the change in temperature predicted by our model over a short time interval $\Delta t$.
We can now use this to make a prediction:
Given a current temperature $T_1$ at time $t_1$, what is the temperature $T_2$ at a future time $t_2$?
We can write
$$ \Delta T = T_2-T_1 $$$$ \Delta t = t_2-t_1 $$and so our model says
$$ C \frac{T_2-T_1}{\Delta t} = \text{ASR} - \text{OLR} $$Which we can rearrange to solve for the future temperature:
$$ T_2 = T_1 + \frac{\Delta t}{C} \left( \text{ASR} - \text{OLR}(T_1) \right) $$We now have a formula with which to make our prediction!
Notice that we have written the OLR as a function of temperature. We will use the current temperature $T_1$ to compute the OLR, and use that OLR to determine the future temperature.
The quantity $\Delta t$ is called a timestep. It is the smallest time interval represented in our model.
Here we're going to use a timestep of 1 year:
dt = 60. * 60. * 24. * 365. # one year expressed in seconds
# Try a single timestep, assuming we have working functions for ASR and OLR
T1 = 288.
T2 = T1 + dt / C * ( ASR(alpha) - OLR(T1, tau=0.57) )
print(T2)
What happened? Why?
Try another timestep
T1 = T2
T2 = T1 + dt / C * ( ASR(alpha) - OLR(T1, tau=0.57) )
print(T2)
Warmed up again, but by a smaller amount.
But this is tedious typing. Time to define a function to make things easier and more reliable:
def step_forward(T):
return T + dt / C * ( ASR(alpha) - OLR(T, tau=0.57) )
Try it out with an arbitrary temperature:
step_forward(300.)
Notice that our function calls other functions and variables we have already defined.
This is both very useful and occasionally confusing.
Now let's really harness the power of the computer by making a loop (and storing values in arrays):
import numpy as np
numsteps = 20
Tsteps = np.zeros(numsteps+1)
Years = np.zeros(numsteps+1)
Tsteps[0] = 288.
for n in range(numsteps):
Years[n+1] = n+1
Tsteps[n+1] = step_forward( Tsteps[n] )
print(Tsteps)
What did we just do?
for
statement executes a statement (or series of statements) a specified number of times (a loop!)¶Now let's draw a picture of our result!
# a special instruction for the Jupyter notebook
# Display all plots inline in the notebook
%matplotlib inline
# import the plotting package
import matplotlib.pyplot as plt
plt.plot( Years, Tsteps )
plt.xlabel('Years')
plt.ylabel('Global mean temperature (K)');
Note how the temperature adjusts smoothly toward the equilibrium temperature, that is, the temperature at which ASR = OLR.
If the planetary energy budget is out of balance, the temperature must change so that the OLR gets closer to the ASR!
The adjustment is actually an exponential decay process: The rate of adjustment slows as the temperature approaches equilibrium.
The temperature gets very very close to equilibrium but never reaches it exactly.
This is actually not native Python, but uses a special graphics library called matplotlib
.
Just about all of our notebooks will start with this:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt