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
- Simulation versus parameterization of heat transport
- The temperature diffusion parameterization
- Solving the temperature diffusion equation with
climlab - Parameterizing the radiation terms
- The one-dimensional diffusive energy balance model
- The annual-mean EBM
- Effects of diffusivity in the EBM
- Summary: parameter values in the diffusive EBM
In the previous lectures we have seen how heat transport by winds and ocean currents acts to COOL the tropics and WARM the poles. The observed temperature gradient is a product of both the insolation and the heat transport!
We were able to ignore this issue in our models of the global mean temperature, because the transport just moves energy around between latitude bands – does not create or destroy energy.
But if want to move beyond the global mean and create models of the equator-to-pole temperature structure, we cannot ignore heat transport. Has to be included somehow!
This leads to us the old theme of simulation versus parameterization.
Complex climate models like the CESM simulate the heat transport by solving the full equations of motion for the atmosphere (and ocean too, if coupled).
Let's revisit an animation of the global 6-hourly sea-level pressure field from our slab ocean simulation with CESM. (We first saw this [back in Lecture 4](Lecture04 -- Climate system components.ipynb))
from IPython.display import YouTubeVideo
YouTubeVideo('As85L34fKYQ') <iframe
width="400"
height="300"
src="https://www.youtube.com/embed/As85L34fKYQ"
frameborder="0"
allowfullscreen
></iframe>
All these traveling weather systems tend to move warm, moist air poleward and cold, dry air equatorward. There is thus a net poleward energy transport.
A model like this needs to simulate the weather in order to model the heat transport.
Let’s emphasize: the most important role for heat transport by winds and ocean currents is to more energy from where it’s WARM to where it’s COLD, thereby reducing the temperature gradient (equator to pole) from what it would be if the planet were in radiative-convective equilibrium everywhere with no north-south motion.
This is the basis for the parameterization of heat transport often used in simple climate models.
Discuss analogy with molecular heat conduction: metal rod with one end in the fire.
Define carefully temperature gradient dTs / dy Measures how quickly the temperature decreases as we move northward (negative in NH, positive in SH)
In any conduction or diffusion process, the flux (transport) of a quantity is always DOWN-gradient (from WARM to COLD).
So our parameterization will look like
Where
Last time we wrote down an energy budget for a thin zonal band centered at latitude
where we have written every term as an explicit function of latitude to remind ourselves that this is a local budget, unlike the zero-dimensional global budget we considered at the start of the course.
Let’s now formally introduce a parameterization that approximates the heat transport as a down-gradient diffusion process:
With
The value of
Notice that we have explicitly chosen to the use surface temperature gradient to set the heat transport. This is a convenient (and traditional) choice to make, but it is not the only possibility! We could instead tune our parameterization to some measure of the free-tropospheric temperature gradient.
Plug the parameterization into our energy budget to get
If we assume that
Let's now make the same assumption we made [back at the beginning of the course](Lecture01 -- Planetary energy budget.ipynb) when we first wrote down the zero-dimensional EBM.
Most of the heat capacity is in the oceans, so that the energy content of each column $E$ is proportional to surface temperature:
where
Now our budget becomes a PDE for the surface temperature
Notice that if we were NOT on a spherical planet and didn’t have to worry about the changing size of latitude circles, this would look something like
$$ \frac{\partial T}{\partial t} = K \frac{\partial^2 T}{\partial y^2} + \text{forcing terms} $$
with
Does equation look familiar?
This is the heat equation, one of the central equations in classical mathematical physics.
This equation describes the behavior of a diffusive system, i.e. how mixing by random molecular motion smears out the temperature.
In our case, the analogy is between the random molecular motion of a metal rod, and the net mixing / stirring effect of weather systems.
Take the integral
The global average of the last term (heat transport) must go to zero (why?)
Therefore this reduces to our familiar zero-dimensional EBM.
climlab has a pre-defined process for solving the meridional diffusion equation. Let's look at a simple example in which diffusion is the ONLY process that changes the temperature.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import climlab
from climlab import constants as const
# Disable interactive plotting (use explicit display calls to show figures)
plt.ioff()# First define an initial temperature field
# that is warm at the equator and cold at the poles
# and varies smoothly with latitude in between
from climlab.utils import legendre
sfc = climlab.domain.zonal_mean_surface(num_lat=90, water_depth=10.)
lat = sfc.lat.points
initial = 12. - 40. * legendre.P2(np.sin(np.deg2rad(lat)))
fig, ax = plt.subplots()
ax.plot(lat, initial)
ax.set_xlabel('Latitude')
ax.set_ylabel('Temperature (deg C)')
fig
## Set up the climlab diffusion process
# make a copy of initial so that it remains unmodified
Ts = climlab.Field(np.array(initial), domain=sfc)
# thermal diffusivity in W/m**2/degC
D = 0.55
# meridional diffusivity in 1/s
K = D / sfc.heat_capacity
# create the climlab diffusion process
# setting the diffusivity and a timestep of ONE MONTH
d = climlab.dynamics.MeridionalDiffusion(state=Ts, K=K,
timestep=const.seconds_per_month)
print dclimlab Process of type <class 'climlab.dynamics.diffusion.MeridionalDiffusion'>.
State variables and domain shapes:
default: (90, 1)
The subprocess tree:
top: <class 'climlab.dynamics.diffusion.MeridionalDiffusion'>
# We are going to step forward one month at a time
# and store the temperature each time
niter = 5
temp = np.zeros((Ts.size, niter+1))
temp[:, 0] = np.squeeze(Ts)
for n in range(niter):
d.step_forward()
temp[:, n+1] = np.squeeze(Ts)# Now plot the temperatures
fig,ax = plt.subplots()
ax.plot(lat, temp)
ax.set_xlabel('Latitude')
ax.set_ylabel('Temperature (deg C)')
ax.legend(range(niter+1))
fig
At each timestep, the warm temperatures get cooler (at the equator) while the cold polar temperatures get warmer!
Diffusion is acting to reduce the temperature gradient.
If we let this run a long time, what should happen??
Try it yourself and find out!
Here we have used a function called the “2nd Legendre polynomial”, defined as
where we have also set
Just turns out to be a useful mathematical description of the relatively smooth changes in things like annual-mean insolation from equator to pole.
In fact these are so useful that they are coded up in a special module within climlab:
x = np.linspace(-1,1)
fig,ax = plt.subplots()
ax.plot(x, legendre.P2(x))
ax.set_title('$P_2(x)$')
fig
Let's go back to the complete budget with our heat transport parameterization
We want to express this as a closed equation for surface temperature
First, as usual, we can write the solar term as
For now, we will assume that the planetary albedo is fixed (does not depend on temperature). Therefore the entire shortwave term
Note that the solar term is (at least in annual average) larger at equator than poles… and transport term acts to flatten out the temperatures.
Now, we almost have a model we can solve for T! Just need to express the OLR in terms of temperature.
So… what’s the link between OLR and temperature????
[ discuss ]
We spent a good chunk of the course looking at this question, and developed a model of a vertical column of air.
We are trying now to build a model of the equator-to-pole (or pole-to-pole) temperature structure.
We COULD use an array of column models, representing temperature as a function of height and latitude (and time).
But instead, we will keep things simple, one spatial dimension at a time.
Introduce the following simple parameterization:
With
Let's look at the data to find reasonable values for
import xarray as xr
ncep_url = "http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/"
ncep_air = xr.open_dataset( ncep_url + "pressure/air.mon.1981-2010.ltm.nc", decode_times=False)
ncep_Ts = xr.open_dataset( ncep_url + "surface_gauss/skt.sfc.mon.1981-2010.ltm.nc", decode_times=False)
lat_ncep = ncep_Ts.lat; lon_ncep = ncep_Ts.lon
print ncep_Ts<xarray.Dataset>
Dimensions: (lat: 94, lon: 192, nbnds: 2, time: 12)
Coordinates:
* lon (lon) float32 0.0 1.875 3.75 5.625 7.5 9.375 11.25 ...
* time (time) float64 -6.571e+05 -6.57e+05 -6.57e+05 ...
* lat (lat) float32 88.542 86.6531 84.7532 82.8508 80.9473 ...
Dimensions without coordinates: nbnds
Data variables:
climatology_bounds (time, nbnds) float64 ...
skt (time, lat, lon) float64 ...
valid_yr_count (time, lat, lon) float64 ...
Attributes:
title: 4x daily NMC reanalysis
description: Data is from NMC initialized reanalysis\n...
platform: Model
Conventions: COARDS
not_missing_threshold_percent: minimum 3% values input to have non-missi...
history: Created 2011/07/12 by doMonthLTM\nConvert...
References: http://www.esrl.noaa.gov/psd/data/gridded...
dataset_title: NCEP-NCAR Reanalysis 1
Ts_ncep_annual = ncep_Ts.skt.mean(dim=('lon','time'))ncep_ulwrf = xr.open_dataset( ncep_url + "other_gauss/ulwrf.ntat.mon.1981-2010.ltm.nc", decode_times=False)
ncep_dswrf = xr.open_dataset( ncep_url + "other_gauss/dswrf.ntat.mon.1981-2010.ltm.nc", decode_times=False)
ncep_uswrf = xr.open_dataset( ncep_url + "other_gauss/uswrf.ntat.mon.1981-2010.ltm.nc", decode_times=False)
OLR_ncep_annual = ncep_ulwrf.ulwrf.mean(dim=('lon','time'))
ASR_ncep_annual = (ncep_dswrf.dswrf - ncep_uswrf.uswrf).mean(dim=('lon','time'))from scipy.stats import linregress
slope, intercept, r_value, p_value, std_err = linregress(Ts_ncep_annual, OLR_ncep_annual)
print 'Best fit is A = %0.0f W/m2 and B = %0.1f W/m2/degC' %(intercept, slope)Best fit is A = 214 W/m2 and B = 1.6 W/m2/degC
We're going to plot the data and the best fit line, but also another line using these values:
# More standard values
A = 210.
B = 2.fig, ax1 = plt.subplots(figsize=(8,6))
ax1.plot( Ts_ncep_annual, OLR_ncep_annual, 'o' , label='data')
ax1.plot( Ts_ncep_annual, intercept + slope * Ts_ncep_annual, 'k--', label='best fit')
ax1.plot( Ts_ncep_annual, A + B * Ts_ncep_annual, 'r--', label='B=2')
ax1.set_xlabel('Surface temperature (C)', fontsize=16)
ax1.set_ylabel('OLR (W m$^{-2}$)', fontsize=16)
ax1.set_title('OLR versus surface temperature from NCEP reanalysis', fontsize=18)
ax1.legend(loc='upper left')
ax1.grid()fig
Discuss these curves...
Suggestion of at least 3 different regimes with different slopes (cold, medium, warm).
Unbiased "best fit" is actually a poor fit over all the intermediate temperatures.
The astute reader will note that... by taking the zonal average of the data before the regression, we are biasing this estimate toward cold temperatures. [WHY?]
Let's take these reference values:
Note that in the global average, recall
And so this parameterization gives
And the observed global mean is
Putting the above OLR parameterization into our budget equation gives
This is the equation for a very important and useful simple model of the climate system. It is typically referred to as the (one-dimensional) Energy Balance Model.
(although as we have seen over and over, EVERY climate model is actually an “energy balance model” of some kind)
Also for historical reasons this is often called the Budyko-Sellers model, after Budyko and Sellers who both (independently of each other) published influential papers on this subject in 1969.
Recap: parameters in this model are
- C: heat capacity in J m$^{-2}$ ºC$^{-1}$
- A: longwave emission at 0ºC in W m$^{-2}$
- B: increase in emission per degree, in W m$^{-2}$ ºC$^{-1}$
- D: horizontal (north-south) diffusivity of the climate system in W m$^{-2}$ ºC$^{-1}$
We also need to specify the albedo.
Let's go back to the NCEP Reanalysis data to see how planetary albedo actually varies as a function of latitude.
days = np.linspace(1.,50.)/50 * const.days_per_year
Qann_ncep = np.mean( climlab.solar.insolation.daily_insolation(lat_ncep, days ),axis=1)
albedo_ncep = 1 - ASR_ncep_annual / Qann_ncep
albedo_ncep_global = np.average(albedo_ncep, weights=np.cos(np.deg2rad(lat_ncep)))print 'The annual, global mean planetary albedo is %0.3f' %albedo_ncep_global
fig,ax = plt.subplots()
ax.plot(lat_ncep, albedo_ncep)
ax.grid();
ax.set_xlabel('Latitude')
ax.set_ylabel('Albedo')
figThe annual, global mean planetary albedo is 0.354
The albedo increases markedly toward the poles.
There are several reasons for this:
- surface snow and ice increase toward the poles
- Cloudiness is an important (but complicated) factor.
- Albedo increases with solar zenith angle (the angle at which the direct solar beam strikes a surface)
Like temperature and insolation, this can be approximated by a smooth function that increases with latitude:
where
In effect we are using a truncated series expansion of the full meridional structure of
We will set
# Add a new curve to the previous figure
a0 = albedo_ncep_global
a2 = 0.25
ax.plot(lat_ncep, a0 + a2 * legendre.P2(np.sin(np.deg2rad(lat_ncep))))
fig
Of course we are not fitting all the details of the observed albedo curve. But we do get the correct global mean a reasonable representation of the equator-to-pole gradient in albedo.
Suppose we take the annual mean of the planetary energy budget.
If the albedo is fixed, then the average is pretty simple. Our EBM equation is purely linear, so the change over one year is just
where
Notice that once we average over the seasonal cycle, there are no time-dependent forcing terms. The temperature will just evolve toward a steady equilibrium.
The equilibrium temperature is then the solution of this Ordinary Differential Equation (setting
You will often see this equation written in terms of the independent variable
which is 0 at the equator and
This is actually a 2nd order ODE, and actually a 2-point Boundary Value Problem for the temperature
This form can be convenient for analytical solutions. As we will see, the non-dimensional number
We will leave the time derivative in our model, because this is the most convenient way to find the equilibrium solution!
There is code available in climlab to solve the diffusive EBM.
Before looking at the details of how to set up an EBM in climlab, let's look at an animation of the adjustment of the model (its temperature and energy budget) from an isothermal initial condition.
For reference, all the code necessary to generate the animation is here in the notebook.
# Some imports needed to make and display animations
from IPython.display import HTML
from matplotlib import animation
def setup_figure():
templimits = -20,32
radlimits = -340, 340
htlimits = -6,6
latlimits = -90,90
lat_ticks = np.arange(-90,90,30)
fig, axes = plt.subplots(3,1,figsize=(8,10))
axes[0].set_ylabel('Temperature (deg C)')
axes[0].set_ylim(templimits)
axes[1].set_ylabel('Energy budget (W m$^{-2}$)')
axes[1].set_ylim(radlimits)
axes[2].set_ylabel('Heat transport (PW)')
axes[2].set_ylim(htlimits)
axes[2].set_xlabel('Latitude')
for ax in axes: ax.set_xlim(latlimits); ax.set_xticks(lat_ticks); ax.grid()
fig.suptitle('Diffusive energy balance model with annual-mean insolation', fontsize=14)
return fig, axes
def initial_figure(model):
# Make figure and axes
fig, axes = setup_figure()
# plot initial data
lines = []
lines.append(axes[0].plot(model.lat, model.Ts)[0])
lines.append(axes[1].plot(model.lat, model.ASR, 'k--', label='SW')[0])
lines.append(axes[1].plot(model.lat, -model.OLR, 'r--', label='LW')[0])
lines.append(axes[1].plot(model.lat, model.net_radiation, 'c-', label='net rad')[0])
lines.append(axes[1].plot(model.lat, model.heat_transport_convergence(), 'g--', label='dyn')[0])
lines.append(axes[1].plot(model.lat,
np.squeeze(model.net_radiation)+model.heat_transport_convergence(), 'b-', label='total')[0])
axes[1].legend(loc='upper right')
lines.append(axes[2].plot(model.lat_bounds, model.diffusive_heat_transport())[0])
lines.append(axes[0].text(60, 25, 'Day 0'))
return fig, axes, lines
def animate(day, model, lines):
model.step_forward()
# The rest of this is just updating the plot
lines[0].set_ydata(model.Ts)
lines[1].set_ydata(model.ASR)
lines[2].set_ydata(-model.OLR)
lines[3].set_ydata(model.net_radiation)
lines[4].set_ydata(model.heat_transport_convergence())
lines[5].set_ydata(np.squeeze(model.net_radiation)+model.heat_transport_convergence())
lines[6].set_ydata(model.diffusive_heat_transport())
lines[-1].set_text('Day {}'.format(int(model.time['days_elapsed'])))
return lines # A model starting from isothermal initial conditions
e = climlab.EBM_annual()
e.Ts[:] = 15. # in degrees Celsius
e.compute_diagnostics()# Plot initial data
fig, axes, lines = initial_figure(e)ani = animation.FuncAnimation(fig, animate, frames=np.arange(1, 100), fargs=(e, lines))HTML(ani.to_html5_video())Here is a simple example using the parameter values we just discussed.
For simplicity, this model will use the annual mean insolation, so the forcing is steady in time.
We haven't yet selected an appropriate value for the diffusivity
D = 0.1
model = climlab.EBM_annual(A=210, B=2, D=D, a0=0.354, a2=0.25)
print modelclimlab Process of type <class 'climlab.model.ebm.EBM_annual'>.
State variables and domain shapes:
Ts: (90, 1)
The subprocess tree:
top: <class 'climlab.model.ebm.EBM_annual'>
diffusion: <class 'climlab.dynamics.diffusion.MeridionalDiffusion'>
LW: <class 'climlab.radiation.aplusbt.AplusBT'>
albedo: <class 'climlab.surface.albedo.P2Albedo'>
insolation: <class 'climlab.radiation.insolation.AnnualMeanInsolation'>
model.param{'A': 210,
'B': 2,
'D': 0.1,
'S0': 1365.2,
'a0': 0.354,
'a2': 0.25,
'timestep': 350632.51200000005,
'water_depth': 10.0}
model.integrate_years(10)Integrating for 900 steps, 3652.422 days, or 10 years.
Total elapsed time is 10.0 years.
fig, axes = plt.subplots(1,2, figsize=(12,4))
ax = axes[0]
ax.plot(model.lat, model.Ts, label=('D = %0.1f' %D))
ax.plot(lat_ncep, Ts_ncep_annual, label='obs')
ax.set_ylabel('Temperature (degC)')
ax = axes[1]
energy_in = np.squeeze(model.ASR - model.OLR)
ax.plot(model.lat, energy_in, label=('D = %0.1f' %D))
ax.plot(lat_ncep, ASR_ncep_annual - OLR_ncep_annual, label='obs')
ax.set_ylabel('Net downwelling radiation at TOA (W m$^{-2}$)')
for ax in axes:
ax.set_xlabel('Latitude'); ax.legend(); ax.grid()
fig
def inferred_heat_transport( energy_in, lat_deg ):
'''Returns the inferred heat transport (in PW) by integrating the net energy imbalance from pole to pole.'''
from scipy import integrate
from climlab import constants as const
lat_rad = np.deg2rad( lat_deg )
return ( 1E-15 * 2 * np.math.pi * const.a**2 *
integrate.cumtrapz( np.cos(lat_rad)*energy_in,
x=lat_rad, initial=0. ) )fig, ax = plt.subplots()
ax.plot(model.lat, inferred_heat_transport(energy_in, model.lat), label=('D = %0.1f' %D))
ax.set_ylabel('Heat transport (PW)')
ax.legend(); ax.grid()
ax.set_xlabel('Latitude')
fig
The upshot: compared to observations, this model has a much too large equator-to-pole temperature gradient, and not enough poleward heat transport!
Apparently we need to increase the diffusivity to get a better fit.
- Solve the annual-mean EBM (integrate out to equilibrium) over a range of different diffusivity parameters.
- Make three plots:
- Global-mean temperature as a function of
$D$ - Equator-to-pole temperature difference
$\Delta T$ as a function of$D$ - Maximum poleward heat transport
$\mathcal{H}_{max}$ as a function of$D$
- Global-mean temperature as a function of
- Choose a value of
$D$ that gives a reasonable approximation to observations:-
$\Delta T \approx 45$ ºC -
$\mathcal{H}_{max} \approx 5.5$ PW
-
Darray = np.arange(0., 2.05, 0.05)model_list = []
Tmean_list = []
deltaT_list = []
Hmax_list = []
for D in Darray:
ebm = climlab.EBM_annual(A=210, B=2, a0=0.354, a2=0.25, D=D)
ebm.integrate_years(20., verbose=False)
Tmean = ebm.global_mean_temperature()
deltaT = np.max(ebm.Ts) - np.min(ebm.Ts)
energy_in = np.squeeze(ebm.ASR - ebm.OLR)
Htrans = inferred_heat_transport(energy_in, ebm.lat)
Hmax = np.max(Htrans)
model_list.append(ebm)
Tmean_list.append(Tmean)
deltaT_list.append(deltaT)
Hmax_list.append(Hmax)color1 = 'b'
color2 = 'r'
fig = plt.figure(figsize=(8,6))
ax1 = fig.add_subplot(111)
ax1.plot(Darray, deltaT_list, color=color1)
ax1.plot(Darray, Tmean_list, 'b--')
ax1.set_xlabel('D (W m$^{-2}$ K$^{-1}$)', fontsize=14)
ax1.set_xticks(np.arange(Darray[0], Darray[-1], 0.2))
ax1.set_ylabel('$\Delta T$ (equator to pole)', fontsize=14, color=color1)
for tl in ax1.get_yticklabels():
tl.set_color(color1)
ax2 = ax1.twinx()
ax2.plot(Darray, Hmax_list, color=color2)
ax2.set_ylabel('Maximum poleward heat transport (PW)', fontsize=14, color=color2)
for tl in ax2.get_yticklabels():
tl.set_color(color2)
ax1.set_title('Effect of diffusivity on temperature gradient and heat transport in the EBM', fontsize=16)
ax1.grid()
ax1.plot([0.6, 0.6], [0, 140], 'k-')[<matplotlib.lines.Line2D at 0x1298c0f90>]
fig
When
When
The real climate seems to lie in a sweet spot in between these limits.
It looks like our fitting criteria are met reasonably well with
Also, note that the global mean temperature (plotted in dashed blue) is completely insensitive to
Our model is defined by the following equation
with the albedo given by
We have chosen the following parameter values, which seems to give a reasonable fit to the observed annual mean temperature and energy budget:
- $ A = 210 ~ \text{W m}^{-2}$
- $ B = 2 ~ \text{W m}^{-2}~^\circ\text{C}^{-1} $
- $ a_0 = 0.354$
- $ a_2 = 0.25$
- $ D = 0.6 ~ \text{W m}^{-2}~^\circ\text{C}^{-1} $
There is one parameter left to choose: the heat capacity
[Why?]
We will instead look at seasonally varying models in the next set of notes.
%load_ext version_information
%version_information numpy, scipy, matplotlib, xarray, climlab| 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 |
| numpy | 1.11.1 |
| scipy | 0.18.1 |
| matplotlib | 2.0.0 |
| xarray | 0.9.5 |
| climlab | 0.5.6 |
| Thu May 25 11:53:40 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.