Let the surface albedo be larger wherever the temperature is below some threshold $T_f$:
$$ \alpha\left(\phi, T(\phi) \right) = \left\{\begin{array}{ccc} \alpha_0 + \alpha_2 P_2(\sin\phi) & ~ & T(\phi) > T_f \\ a_i & ~ & T(\phi) \le T_f \\ \end{array} \right. $$where $P_2(\sin\phi) = \frac{1}{2}\left( 3\left(\sin\phi\right)^2 - 1 \right) $ is called the second Legendre Polynomial (just a mathematically convenient description of a smooth variation between the equator and pole).
In [ ]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import climlab
In [ ]:
# for convenience, set up a dictionary with our reference parameters
param = {'D':0.55, 'A':210, 'B':2, 'a0':0.3, 'a2':0.078, 'ai':0.62, 'Tf':-10.}
model1 = climlab.EBM_annual( num_lat=180, D=0.55, A=210., B=2., Tf=-10., a0=0.3, a2=0.078, ai=0.62)
print model1
Because we provided a parameter ai for the icy albedo, our model now contains several sub-processes contained within the process called albedo. Together these implement the step-function formula above.
The process called iceline simply looks for grid cells with temperature below $T_f$.
In [ ]:
print model1.param
In [ ]:
# A python shortcut... we can use the dictionary to pass lots of input arguments simultaneously:
# same thing as before, but written differently:
model1 = climlab.EBM_annual( num_lat=180, **param)
print model1
In [ ]:
def ebm_plot(e, return_fig=False):
templimits = -60,32
radlimits = -340, 340
htlimits = -6,6
latlimits = -90,90
lat_ticks = np.arange(-90,90,30)
fig = plt.figure(figsize=(8,12))
ax1 = fig.add_subplot(3,1,1)
ax1.plot(e.lat, e.Ts)
ax1.set_ylim(templimits)
ax1.set_ylabel('Temperature (deg C)')
ax2 = fig.add_subplot(3,1,2)
ax2.plot(e.lat, e.ASR, 'k--', label='SW' )
ax2.plot(e.lat, -e.OLR, 'r--', label='LW' )
ax2.plot(e.lat, e.net_radiation, 'c-', label='net rad' )
ax2.plot(e.lat, e.heat_transport_convergence(), 'g--', label='dyn' )
ax2.plot(e.lat, e.net_radiation.squeeze() + e.heat_transport_convergence(), 'b-', label='total' )
ax2.set_ylim(radlimits)
ax2.set_ylabel('Energy budget (W m$^{-2}$)')
ax2.legend()
ax3 = fig.add_subplot(3,1,3)
ax3.plot(e.lat_bounds, e.heat_transport() )
ax3.set_ylim(htlimits)
ax3.set_ylabel('Heat transport (PW)')
for ax in [ax1, ax2, ax3]:
ax.set_xlabel('Latitude')
ax.set_xlim(latlimits)
ax.set_xticks(lat_ticks)
ax.grid()
if return_fig:
return fig
In [ ]:
# Integrate out to equilibrium.
model1.integrate_years(5)
# Check for energy balance
print climlab.global_mean(model1.net_radiation)
f = ebm_plot(model1)
In [ ]:
# There is a diagnostic that tells us the current location of the ice edge:
model1.icelat
In [ ]:
deltaA = 4.
# This is a very handy way to "clone" an existing model:
model2 = climlab.process_like(model1)
# Now change the longwave parameter:
model2.subprocess['LW'].A = param['A'] - deltaA
# and integrate out to equilibrium again
model2.integrate_years(5, verbose=False)
plt.plot(model1.lat, model1.Ts)
plt.plot(model2.lat, model2.Ts)
The warming is polar-amplified: more warming at the poles than elsewhere.
Why?
Also, the current ice line is now:
In [ ]:
model2.icelat
There is no ice left!
Let's do some more greenhouse warming:
In [ ]:
model3 = climlab.process_like(model1)
model3.subprocess['LW'].A = param['A'] - 2*deltaA
model3.integrate_years(5, verbose=False)
plt.plot(model1.lat, model1.Ts)
plt.plot(model2.lat, model2.Ts)
plt.plot(model3.lat, model3.Ts)
plt.xlim(-90, 90)
plt.grid()
In the ice-free regime, there is no polar-amplified warming. A uniform radiative forcing produces a uniform warming.
In [ ]:
m = climlab.EBM_annual( num_lat=180, **param )
# The current (default) solar constant, corresponding to present-day conditions:
m.subprocess.insolation.S0
In [ ]:
# First, get to equilibrium
m.integrate_years(5.)
# Check for energy balance
print climlab.global_mean(m.net_radiation)
In [ ]:
m.icelat
In [ ]:
# Now make the solar constant smaller:
m.subprocess.insolation.S0 = 1300.
In [ ]:
# Integrate to new equilibrium
m.integrate_years(10.)
# Check for energy balance
print climlab.global_mean(m.net_radiation)
In [ ]:
m.icelat
In [ ]:
ebm_plot(m)
A much colder climate! The ice line is sitting at 54º. The heat transport shows that the atmosphere is moving lots of energy across the ice line, trying hard to compensate for the strong radiative cooling everywhere poleward of the ice line.
In [ ]:
# Now make the solar constant smaller:
m.subprocess.insolation.S0 = 1200.
# First, get to equilibrium
m.integrate_years(5.)
# Check for energy balance
print climlab.global_mean(m.net_radiation)
ebm_plot(m)
Something very different happened! Where is the ice line now?
In [ ]:
# Now make the solar constant smaller:
m.subprocess.insolation.S0 = 1365.2
# First, get to equilibrium
m.integrate_years(5.)
# Check for energy balance
print climlab.global_mean(m.net_radiation)
In [ ]:
ebm_plot(m)
Is this the same climate we started with?
To melt all the ice and get out of this so-called "Snowball Earth", we need to raise the solar constant much higher.
You will look at this in your homework assignment.
In [ ]: