The seasonal cycle of surface temperature

Look at the observed seasonal cycle in the NCEP reanalysis data.

Read in the necessary data from the online server.

The catalog is here: http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/catalog.html

In [1]:
import ClimateUtils as clim
import orbital
import ebm
import netCDF4 as nc

datapath = "http://ramadda.atmos.albany.edu:8080/repository/opendap/latest/Top/Users/Brian+Rose/CESM+runs/"
endstr = "/entry.das"
In [2]:
ncep_url = "http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/"
ncep_air = nc.Dataset( ncep_url + "pressure/air.mon.1981-2010.ltm.nc" )
ncep_Ts = nc.Dataset( ncep_url + "surface_gauss/skt.sfc.mon.1981-2010.ltm.nc" )
lat_ncep = ncep_Ts.variables['lat'][:]; lon_ncep = ncep_Ts.variables['lon'][:]
Ts_ncep = ncep_Ts.variables['skt'][:]
print Ts_ncep.shape
(12, 94, 192)

Load the topography file from CESM, just so we can plot the continents.

In [3]:
topo = nc.Dataset( datapath + 'som_input/USGS-gtopo30_1.9x2.5_remap_c050602.nc' + endstr )
lat_cesm = topo.variables['lat'][:]
lon_cesm = topo.variables['lon'][:]

Make two maps: one of annual mean surface temperature, another of the seasonal range (max minus min).

In [4]:
maxTs = np.max(Ts_ncep,axis=0)
minTs = np.min(Ts_ncep,axis=0)
In [5]:
fig = figure( figsize=(16,6) )

ax1 = fig.add_subplot(1,2,1)
cax1 = ax1.pcolormesh( lon_ncep, lat_ncep, mean(Ts_ncep, axis=0), cmap=plt.cm.seismic )
cbar1 = colorbar(cax1)
ax1.set_title('Annual mean surface temperature (degC)', fontsize=14 )

ax2 = fig.add_subplot(1,2,2)
cax2 = ax2.pcolormesh( lon_ncep, lat_ncep, maxTs - minTs )
cbar2 = colorbar(cax2)
ax2.set_title('Seasonal temperature range (degC)', fontsize=14)

for ax in [ax1,ax2]:
    ax.contour( lon_cesm, lat_cesm, topo.variables['LANDFRAC'], [0.5,0.5], colors='k');
    ax.axis([0, 360, -90, 90])
    ax.set_xlabel('Longitude', fontsize=14 ); ax.set_ylabel('Latitude', fontsize=14 )

Make a contour plot of the zonal mean temperature as a function of time of year

In [6]:
[-65,-55,-45,-35,-25,-15,-5,5,15,25,35,45,55,65]
Out[6]:
[-65, -55, -45, -35, -25, -15, -5, 5, 15, 25, 35, 45, 55, 65]
In [7]:
arange(-65,75,10)
Out[7]:
array([-65, -55, -45, -35, -25, -15,  -5,   5,  15,  25,  35,  45,  55,  65])
In [8]:
Tmax = 65; Tmin = -Tmax; delT = 10
clevels = arange(Tmin,Tmax+delT,delT)

fig_zonobs = figure( figsize=(10,6) )
ax = fig_zonobs.add_subplot(111)
cax = contourf( arange(12)+0.5, lat_ncep, transpose(mean( Ts_ncep, axis=2 )), levels=clevels, 
               cmap=plt.cm.seismic, vmin=Tmin, vmax=Tmax )
ax.set_xlabel('Month', fontsize=16)
ax.set_ylabel('Latitude', fontsize=16 )
cbar = colorbar(cax)
ax.set_title('Zonal mean surface temperature (degC)', fontsize=20)
show()

Exploring the amplitude of the seasonal cycle with an EBM

We are looking at the 1D (zonally averaged) energy balance model with diffusive heat transport. The equation is

\(C \frac{\partial T(\phi,t)}{\partial t} = \big(1-\alpha\big) Q(\phi,t) - \Big(A+B T(\phi,t) \Big) + \frac{K}{\cos\phi} \frac{\partial}{\partial \phi} \bigg( \cos\phi \frac{\partial T}{\partial \phi} \bigg)\)

and the code in ebm.py solves this equation numerically.

One handy feature of the EBM code: the function integrate_years() automatically calculates the time averaged temperature. So if we run it for exactly one year, we get the annual mean temperature saved in the field T_timeave.

We will look at the seasonal cycle of temperature in three different models with different heat capacities (which we express though an equivalent depth of water in meters):

In [9]:
model1 = ebm.EBM_seasonal()
model1.integrate_years(1, verbose=True)

water_depths = array([2., 10., 50.])

num_depths = water_depths.size
Tann = empty( [model1.num_points, num_depths] )
models = repeat(ebm.EBM_seasonal(), num_depths )

for n in range(num_depths):
    models[n] = ebm.EBM_seasonal()
    models[n].set_water_depth(water_depths[n])
    models[n].integrate_years(20., verbose=False )
    models[n].integrate_years(1., verbose=False)
    Tann[:,n] = models[n].T_timeave
Integrating for 90 steps or 1 years.
Total elapsed time is 1.0 years.

All models should have the same annual mean temperature:

In [10]:
plot(models[0].lat, Tann)
xlim(-90,90)
xlabel('Latitude')
ylabel('Temperature (degC)')
title('Annual mean temperature in the EBM')
legend( water_depths.astype(str) )
show()

There is no automatic function in the ebm.py code to keep track of minimum and maximum temperatures (though we might add that in the future!)

Instead we'll step through one year "by hand" and save all the temperatures.

In [11]:
Tyear = empty([model1.num_points, model1.num_steps_per_year, num_depths])
for n in range(num_depths):
    for m in range(model1.num_steps_per_year):
        models[n].step_forward()
        Tyear[:,m,n] = models[n].T

Make a figure to compare the observed zonal mean seasonal temperature cycle to what we get from the EBM with different heat capacities:

In [12]:
fig = figure( figsize=(20,10) )

ax = fig.add_subplot(2,num_depths,2)
cax = contourf( arange(12)+0.5, lat_ncep, transpose(mean( Ts_ncep, axis=2 )), levels=clevels, 
               cmap=plt.cm.seismic, vmin=Tmin, vmax=Tmax )
ax.set_xlabel('Month')
ax.set_ylabel('Latitude')
cbar = colorbar(cax)
ax.set_title('Zonal mean surface temperature - observed (degC)', fontsize=20)

for n in range(num_depths):
    ax = fig.add_subplot(2,num_depths,num_depths+n+1)
    cax = contourf( 4*arange(models[n].num_steps_per_year), models[n].lat, Tyear[:,:,n], levels=clevels, 
               cmap=plt.cm.seismic, vmin=Tmin, vmax=Tmax )
    cbar1 = colorbar(cax)
    ax.set_title('water_depth = ' + str(models[n].water_depth) + ' meters', fontsize=20 )
    ax.set_xlabel('Days of year', fontsize=14 )
    ax.set_ylabel('Latitude', fontsize=14 )

#fig.set_title('Temperature in seasonal EBM with various water depths', fontsize=14)
show()

Which one looks more realistic? Depends a bit on where you look. But overall, the observed seasonal cycle matches the 10 meter case best. The effective heat capacity governing the seasonal cycle of the zonal mean temperature is closer to 10 meters of water than to either 2 or 50 meters.

Making an animation of the EBM solutions

In [15]:
fpath = '/Users/Brian/Dropbox/PythonStuff/ebm_seasonal_frames/'

fig = figure( figsize=(20,5) )
#for m in range(model2.num_steps_per_year):
for m in range(1):
    thisday = m * models[0].timestep / clim.seconds_per_day
    Q = orbital.daily_insolation( models[n].lat, thisday )
    for n in range(num_depths):
        c1 = 'b'
        ax = fig.add_subplot(1,num_depths,n+1)
        ax.plot( models[n].lat, Tyear[:,m,n], c1 )
        ax.set_title('water_depth = ' + str(models[n].water_depth) + ' meters', fontsize=20 )
        ax.set_xlabel('Latitude', fontsize=14 )
        if n is 0:
            ax.set_ylabel('Temperature', fontsize=14, color=c1 )
        ax.set_xlim([-90,90])
        ax.set_ylim([-60,60])
        for tl in ax.get_yticklabels():
            tl.set_color(c1)
        grid()
        
        c2 = 'r'
        ax2 = ax.twinx()
        ax2.plot( models[n].lat, Q, c2)
        if n is 2:
            ax2.set_ylabel('Insolation (W m$^{-2}$)', color=c2, fontsize=14)
        for tl in ax2.get_yticklabels():
            tl.set_color(c2)
        ax2.set_xlim([-90,90])
        ax2.set_ylim([0,600])
    
    filename = fpath + 'ebm_seasonal' + str(m).zfill(4) + '.png'
    #fig.savefig(filename)
    #fig.clear()
    show()

The seasonal cycle for a planet with 90º obliquity

The EBM code uses our familiar orbital.py code to calculate insolation, and therefore it's easy to set up a model with different orbital parameters. Here is an example with very different orbital parameters: 90ยบ obliquity. We looked at the distribution of insolation by latitude and season for this type of planet in the last homework.

In [16]:
model_highobl = ebm.EBM_seasonal()
orb_highobl = orbital.lookup_parameters()
orb_highobl['obliquity'] = 90.
orb_highobl['ecc'] = 0.
print orb_highobl
model_highobl.make_insolation_array( orb = orb_highobl )
{'long_peri': array(281.37), 'ecc': 0.0, 'obliquity': 90.0}

Repeat the same procedure to calculate and store temperature throughout one year, after letting the models run out to equilibrium.

In [17]:
Tann_highobl = empty( [model_highobl.num_points, num_depths] )
models_highobl = repeat(model_highobl, num_depths )

for n in range(num_depths):
    models_highobl[n] = ebm.EBM_seasonal()
    models_highobl[n].set_water_depth(water_depths[n])
    models_highobl[n].make_insolation_array( orb = orb_highobl )
    models_highobl[n].integrate_years(40., verbose=False )
    models_highobl[n].integrate_years(1.)
    Tann_highobl[:,n] = models_highobl[n].T_timeave

Tyear_highobl = empty([model_highobl.num_points, model_highobl.num_steps_per_year, num_depths])
for n in range(num_depths):
    for m in range(model_highobl.num_steps_per_year):
        models_highobl[n].step_forward()
        Tyear_highobl[:,m,n] = models_highobl[n].T
Integrating for 90 steps or 1.0 years.
Total elapsed time is 41.0 years.
Integrating for 90 steps or 1.0 years.
Total elapsed time is 41.0 years.
Integrating for 90 steps or 1.0 years.
Total elapsed time is 41.0 years.

And plot the seasonal temperature cycle same as we did above:

In [18]:
fig = figure( figsize=(20,5) )
Tmax_highobl = 125; Tmin_highobl = -Tmax_highobl; delT_highobl = 10
clevels_highobl = arange(Tmin_highobl, Tmax_highobl+delT_highobl, delT_highobl)

for n in range(num_depths):
    ax = fig.add_subplot(1,num_depths,n+1)
    cax = contourf( 4*arange(models_highobl[n].num_steps_per_year), models_highobl[n].lat, Tyear_highobl[:,:,n], 
        levels=clevels_highobl, cmap=plt.cm.seismic, vmin=Tmin_highobl, vmax=Tmax_highobl )
    cbar1 = colorbar(cax)
    ax.set_title('water_depth = ' + str(models_highobl[n].water_depth) + ' meters', fontsize=20 )
    ax.set_xlabel('Days of year', fontsize=14 )
    ax.set_ylabel('Latitude', fontsize=14 )

show()

Note that the temperature range is much larger than for the Earth-like case above (but same contour interval, 10 degC).

Why is the temperature so uniform in the north-south direction with 50 meters of water?

To see the reason, let's plot the annual mean insolation at 90ยบ obliquity, alongside the present-day annual mean insolation:

In [19]:
lat = linspace(-90, 90, 181)
days = linspace(1.,50.)/50 * clim.days_per_year
Q_present = orbital.daily_insolation( lat, days )
Q_highobl = orbital.daily_insolation( lat, days, **orb_highobl )
Q_present_ann = mean( Q_present, axis=1 )
Q_highobl_ann = mean( Q_highobl, axis=1 )
In [20]:
plot( lat, Q_present_ann, label='Earth' )
plot( lat, Q_highobl_ann, label='90deg obliquity' )
grid()
legend(loc='lower center')
xlabel('Latitude', fontsize=14 )
ylabel('W m$^{-2}$', fontsize=14 )
title('Annual mean insolation for two different obliquities', fontsize=16)
show()

Though this is a bit misleading, because our model prescribes an increase in albedo from the equator to the pole. So the absorbed shortwave gradients look even more different.

In []: