In [1]:
%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
import climlab
import netCDF4 as nc
So far in this class we have talked a lot about Equilibrium Climate Sensitivity: the surface warming that is necessary to bring the planetary energy budget back into balance (energy in = energy out) after a doubling of atmospheric CO2.
Although this concept is very important, it is not the only important measure of climate change, and not the only question for which we need to apply climate models.
Consider two basic facts about climate change in the real world:
We will now extend our climate model to deal with both of these issues simultaneously.
In [2]:
# Need the ozone data again for our Radiative-Convective model
ozone_filename = 'ozone_1.9x2.5_L26_2000clim_c091112.nc'
ozone = nc.Dataset(ozone_filename)
lat = ozone.variables['lat'][:]
lon = ozone.variables['lon'][:]
lev = ozone.variables['lev'][:]
O3_zon = np.mean( ozone.variables['O3'],axis=(0,3) )
O3_global = np.sum( O3_zon * np.cos(np.deg2rad(lat)), axis=1 ) / np.sum( np.cos(np.deg2rad(lat) ) )
In [3]:
#steps_per_year = 20
steps_per_year = 180
# parameter set 1 -- store in a dictionary for easy re-use
p1 = {'albedo_sfc': 0.22,
'adj_lapse_rate': 6.,
'timestep': climlab.constants.seconds_per_year/steps_per_year,
'qStrat': 5E-6,
'relative_humidity': 0.77}
# parameter set 2
p2 = {'albedo_sfc': 0.2025,
'adj_lapse_rate': 7.,
'timestep': climlab.constants.seconds_per_year/steps_per_year,
'qStrat': 2E-7,
'relative_humidity': 0.6}
# Make a list of the two parameter sets
param = [p1, p2]
# And a list of two corresponding Radiative-Convective models!
slab = []
for p in param:
model = climlab.BandRCModel(lev=lev, **p)
model.absorber_vmr['O3'] = O3_global
slab.append(model)
Run both models out to equilibrium and check their surface temperatures:
In [4]:
for model in slab:
model.integrate_converge()
print 'The equilibrium surface temperatures are:'
print 'Model 0: %0.2f K' %slab[0].Ts
print 'Model 1: %0.2f K' %slab[1].Ts
Ok so our two models (by construction) start out with nearly identical temperatures.
Now we double CO2 and calculate the Equilibrium Climate Sensitivity:
In [5]:
# We will make copies of each model and double CO2 in the copy
slab_2x = []
for model in slab:
model_2x = climlab.process_like(model)
model_2x.absorber_vmr['CO2'] *= 2.
model_2x.integrate_converge()
slab_2x.append(model_2x)
# Climate sensitivity
DeltaT = []
for n in range(len(slab)):
DeltaT.append(slab_2x[n].Ts - slab[n].Ts)
print 'The equilibrium climate sensitivities are:'
print 'Model 0: %0.2f K' %DeltaT[0]
print 'Model 1: %0.2f K' %DeltaT[1]
So Model 0 is more sensitive than Model 1. It has a larger system gain, or a more positive overall climate feedback.
This is actually due to differences in how we have parameterized the water vapor feedback in the two models. We could look at this more carefully if we wished.
In [6]:
slab[0].depth_bounds
Out[6]:
The "ocean" in these models is just a "slab" of water 1 meter deep.
That's all we need to calculate the equilibrium temperatures, but it tells us nothing about the timescales for climate change in the real world.
For this, we need a deep ocean that can exchange heat with the surface.
We are now going to build two new models. The atmosphere (radiative-convective model) will be identical to the two "slab" models we just used. But these will be coupled to a column of ocean water 2000 m deep!
We will parameterize the ocean heat uptake as a diffusive mixing process. Much like when we discussed the diffusive parameterization for atmospheric heat transport -- we are assuming that ocean dynamics result in a vertical mixing of heat from warm to cold temperatures.
The following code will set this up for us.
We will make one more assumption, just for the sake of illustration:
The more sensitive model (Model 0) is also more efficent at taking up heat into the deep ocean
In [7]:
# Create the domains
ocean_bounds = np.arange(0., 2010., 100.)
depthax = climlab.Axis(axis_type='depth', bounds=ocean_bounds)
levax = climlab.Axis(axis_type='lev', points=lev)
atm = climlab.domain.Atmosphere(axes=levax)
ocean = climlab.domain.Ocean(axes=depthax)
# Model 0 has a higher ocean heat diffusion coefficient --
# a more efficent deep ocean heat sink
ocean_diff = [5.E-4, 3.5E-4]
# List of deep ocean models
deep = []
for n in range(len(param)):
# Create the state variables
Tinitial_ocean = slab[n].Ts * np.ones(ocean.shape)
Tocean = climlab.Field(Tinitial_ocean.copy(), domain=ocean)
Tatm = climlab.Field(slab[n].Tatm.copy(), domain=atm)
# By declaring Ts to be a numpy view of the first element of the array Tocean
# Ts becomes effectively a dynamic reference to that element and is properly updated
Ts = Tocean[0:1]
atm_state = {'Tatm': Tatm, 'Ts': Ts}
model = climlab.BandRCModel(state=atm_state, **param[n])
model.set_state('Tocean', Tocean)
diff = climlab.dynamics.diffusion.Diffusion(state={'Tocean': model.Tocean},
K=ocean_diff[n], diffusion_axis='depth', **param[n])
model.add_subprocess('OHU', diff)
model.absorber_vmr['O3'] = O3_global
deep.append(model)
print deep[0]
Now consider the CO2 increase. In the real world, CO2 has been increasing every year since the beginning of industrialization. Future CO2 concentrations depend on collective choices made by human societies about how much fossil fuel to extract and burn.
We will set up a simple scenario. Suppose that CO2 increases by 1% of its existing concentration every year until it reaches 2x its initial concentration. This takes about 70 years.
After 70 years, we assume that all anthropogenic emissions, and CO2 concentration is stabilized at the 2x level.
What happens to the surface temperature?
How do the histories of surface and deep ocean temperature compare in our two models?
We are going to simulation 400 years of transient global warming in the two models.
In [8]:
# This code will generate a 'live plot' showing the transient warming in the two models.
# The figure will update itself after every 5 years of simulation
num_years = 400
years = np.arange(num_years+1)
Tsarray = []
Tocean = []
netrad = []
for n in range(len(param)):
thisTs = np.nan * np.zeros(num_years+1)
thisnetrad = np.nan * np.zeros(num_years+1)
thisTocean = np.nan * np.zeros((deep[n].Tocean.size, num_years+1))
thisTs[0] = deep[n].Ts
thisnetrad[0] = deep[n].ASR - deep[n].OLR
thisTocean[:, 0] = deep[n].Tocean
Tsarray.append(thisTs)
Tocean.append(thisTocean)
netrad.append(thisnetrad)
CO2initial = deep[0].absorber_vmr['CO2'][0]
CO2array = np.nan * np.zeros(num_years+1)
CO2array[0] = CO2initial * 1E6
colorlist = ['b', 'r']
co2color = 'k'
num_axes = len(param) + 1
fig, ax = plt.subplots(num_axes, figsize=(12,14))
# Twin the x-axis twice to make independent y-axes.
topaxes = [ax[0], ax[0].twinx(), ax[0].twinx()]
# Make some space on the right side for the extra y-axis.
fig.subplots_adjust(right=0.85)
# Move the last y-axis spine over to the right by 10% of the width of the axes
topaxes[-1].spines['right'].set_position(('axes', 1.1))
# To make the border of the right-most axis visible, we need to turn the frame
# on. This hides the other plots, however, so we need to turn its fill off.
topaxes[-1].set_frame_on(True)
topaxes[-1].patch.set_visible(False)
for n, model in enumerate(slab_2x):
topaxes[0].plot(model.Ts*np.ones_like(Tsarray[n]), '--', color=colorlist[n])
topaxes[0].set_ylabel('Surface temperature (K)')
topaxes[0].set_xlabel('Years')
topaxes[0].set_title('Transient warming scenario: 1%/year CO2 increase to doubling, followed by CO2 stabilization', fontsize=14)
topaxes[0].legend(['Model 0', 'Model 1'], loc='lower right')
topaxes[1].plot(CO2array, color=co2color)
topaxes[1].set_ylabel('CO2 (ppm)', color=co2color)
for tl in topaxes[1].get_yticklabels():
tl.set_color(co2color)
topaxes[1].set_ylim(300., 1000.)
topaxes[2].set_ylabel('Radiative forcing (W/m2)', color='b')
for tl in topaxes[2].get_yticklabels():
tl.set_color('b')
topaxes[2].set_ylim(0, 4)
contour_levels = np.arange(-1, 4.5, 0.25)
for n in range(len(param)):
cax = ax[n+1].contourf(years, deep[n].depth, Tocean[n] - Tsarray[n][0], levels=contour_levels)
ax[n+1].invert_yaxis()
ax[n+1].set_ylabel('Depth (m)')
ax[n+1].set_xlabel('Years')
fig.subplots_adjust(bottom=0.12)
cbar_ax = fig.add_axes([0.25, 0.02, 0.5, 0.03])
fig.colorbar(cax, cax=cbar_ax, orientation='horizontal');
# Increase CO2 by 1% / year for 70 years (until doubled), and then hold constant
for y in range(num_years):
if deep[0].absorber_vmr['CO2'][0] < 2 * CO2initial:
for model in deep:
model.absorber_vmr['CO2'] *= 1.01
CO2array[y+1] = deep[0].absorber_vmr['CO2'][0] * 1E6
for n, model in enumerate(deep):
model.integrate_years(1, verbose=False)
Tsarray[n][y+1] = deep[n].Ts
Tocean[n][:, y+1] = deep[n].Tocean
netrad[n][y+1] = deep[n].ASR - deep[n].OLR
# is it time to update the plots?
plot_freq = 5
if y%plot_freq == 0:
for n, model in enumerate(deep):
topaxes[0].plot(Tsarray[n], color=colorlist[n])
topaxes[2].plot(netrad[n], ':', color=colorlist[n])
for n in range(len(param)):
cax = ax[n+1].contourf(years, deep[n].depth, Tocean[n] - Tsarray[n][0], levels=contour_levels)
topaxes[1].plot(CO2array, color=co2color)
fig.canvas.draw()