An Introduction to Climlab

Author of this notebook: Maybritt Schillinger

Welcome! This is a little notebook to help you try some things with the Python package Climlab. The documentation can be found here: https://climlab.readthedocs.io/en/latest/intro.html

Some setups

First, here are some setups and helper functions. No need for you to bother with this. These are just for importing the necessary libraries and defining our own class for an EBM model that we will use later.

In [1]:
from __future__ import division
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import climlab
from builtins import next
from climlab.process.diagnostic import DiagnosticProcess
from climlab.domain.field import Field, global_mean

class tanalbedo(DiagnosticProcess):

    def __init__(self, **kwargs):
        super(tanalbedo, self).__init__(**kwargs)

        self.add_diagnostic('albedo')
        Ts = self.state['Ts']
        self._compute_fixed()

    def _compute_fixed(self):
        Ts = self.state['Ts']
        try:
            lon, lat = np.meshgrid(self.lon, self.lat)
        except:
            lat = self.lat
        phi = lat

        try:
            albedo=np.zeros(len(phi));
            albedo=0.42-0.20*np.tanh(0.052*(Ts-3))

        except:
            albedo = np.zeros_like(phi)

        dom = next(iter(self.domains.values()))
        self.albedo = Field(albedo, domain=dom)

    def _compute(self):
        self._compute_fixed()
        return {}


class simple_EBM_seasonal(climlab.model.ebm.EBM_seasonal):
    
    def __init__(self, orbit={'ecc': 0.0167643, 'long_peri': 280.32687, 'obliquity': 23.459277}, albedo='constant', Tf=-10.0, CO2=300, **kwargs):
        super().__init__(**kwargs)
        self.remove_subprocess('insolation')
        self.remove_subprocess('SW')
        self.remove_subprocess('albedo')
        self.remove_subprocess('LW')

        # catch model domain for subprocess creation
        surface = self.domains['Ts']
        # define new insolation and SW process
        insolation = climlab.radiation.DailyInsolation(domains=surface, orb = orbit, **self.param)
        
        # now define albedo process
        if albedo=='constant':
            alb = climlab.surface.albedo.ConstantAlbedo(domains=surface, **self.param)
        elif albedo=='step':
            alb = climlab.surface.albedo.StepFunctionAlbedo(state=self.state, Tf=Tf, **self.param)
        elif albedo=='tanh':
            alb = tanalbedo(state=self.state, **self.param)
        else:
            print('The given albedo option is not known. Will put in constant albedo.')
            alb = climlab.surface.albedo.ConstantAlbedo(domains=surface, **self.param)
        
        SW = climlab.radiation.absorbed_shorwave.SimpleAbsorbedShortwave(insolation=insolation.insolation, state = self.state, albedo = alb.albedo, **self.param)
        LW = climlab.radiation.aplusbt.AplusBT_CO2(CO2=CO2, state=self.state, **self.param)
        self.add_subprocess('LW', LW)
        self.add_subprocess('insolation', insolation)
        self.add_subprocess('SW', SW)
        self.add_subprocess('albedo', alb)
        self.compute_diagnostics()
    
    def global_annual_mean(self):
        num_steps_per_year = int(self.time['num_steps_per_year'])
        mean_year = np.empty(num_steps_per_year)
        for m in range(num_steps_per_year):
            self.step_forward()
            mean_year[m] = self.global_mean_temperature()
        gam = np.mean(mean_year)
        return gam
        
    def plot_temp_rad(self):
        # creating plot figure
        fig = plt.figure(figsize=(15,5))

        # Temperature plot
        ax1 = fig.add_subplot(121)
        ax1.plot(self.lat,self.Ts)

        ax1.set_xticks([-90,-60,-30,0,30,60,90])
        ax1.set_xlim([-90,90])
        ax1.set_title('Surface Temperature', fontsize=14)
        ax1.set_ylabel('(degC)', fontsize=12)
        ax1.grid()

        # Net Radiation plot
        ax3 = fig.add_subplot(122, sharex = ax1)
        ax3.plot(self.lat, self.OLR, label='OLR',
                                               color='cyan')
        ax3.plot(self.lat, self.ASR, label='ASR',
                                               color='magenta')
        ax3.plot(self.lat, self.ASR-self.OLR, 
                                               label='net radiation',
                                               color='red')

        ax3.set_title('Net Radiation', fontsize=14)
        ax3.set_ylabel('(W/m$^2$)', fontsize=12)
        ax3.legend(loc='best')
        ax3.grid()
        plt.show()
        
    def plot_timeave(self):
        # integrate for one year
        self.integrate_years(1.)
        
        # creating plot figure
        fig = plt.figure(figsize=(15,5))

        # Temperature plot
        ax1 = fig.add_subplot(121)
        ax1.plot(self.lat,self.timeave['Ts'])

        ax1.set_xticks([-90,-60,-30,0,30,60,90])
        ax1.set_xlim([-90,90])
        ax1.set_title('Annual mean Surface Temperature', fontsize=14)
        ax1.set_ylabel('(degC)', fontsize=12)
        ax1.grid()

        # Net Radiation plot
        ax3 = fig.add_subplot(122, sharex = ax1)
        ax3.plot(self.lat, self.timeave['OLR'], label='Annual mean OLR',
                                               color='cyan')
        ax3.plot(self.lat, self.timeave['ASR'], label='Annual mean ASR',
                                               color='magenta')
        ax3.plot(self.lat, self.timeave['ASR']-self.timeave['OLR'], 
                                               label='Annual mean net radiation',
                                               color='red')

        ax3.set_title('Annual mean Net Radiation', fontsize=14)
        ax3.set_ylabel('(W/m$^2$)', fontsize=12)
        ax3.legend(loc='best')
        ax3.grid()
        plt.show()
    
    def plot_albedo(self):
        fig = plt.figure(figsize=(7.5,5))

        # Temperature plot
        ax2 = fig.add_subplot(111)
        ax2.plot(self.lat,self.albedo)

        ax2.set_title('Albedo', fontsize=14)
        ax2.set_xlabel('latitude', fontsize=10)
        ax2.set_ylabel('', fontsize=12)

        ax2.set_xticks([-90,-60,-30,0,30,60,90])
        ax2.set_xlim([-90,90])
        ax2.set_ylim([0,1])
        ax2.grid()
    
        plt.show()
        
    def plot_temp_year(self):

        num_steps_per_year = int(self.time['num_steps_per_year'])
        Tyear = np.empty((self.lat.size, num_steps_per_year))
        for m in range(num_steps_per_year):
            self.step_forward()
            Tyear[:,m] = np.squeeze(self.Ts)
        Tmin=np.min(Tyear)
        Tmax=np.max(Tyear)
        
        fig = plt.figure(figsize=(5,5))
        ax = fig.add_subplot(111)
        
        factor = 365. / num_steps_per_year
        cax = ax.contourf(factor * np.arange(num_steps_per_year),
                      self.lat, Tyear[:,:], 
                      cmap=plt.cm.seismic, vmin=Tmin, vmax=Tmax)
        cbar1 = plt.colorbar(cax)
        ax.set_title('Temperatures throughout the year', fontsize=14)
        ax.set_xlabel('Days of year', fontsize=14)
        ax.set_ylabel('Latitude', fontsize=14)
        
    def plot_ice_year(self):
        
        if 'Tf' in self.subprocess['albedo'].param.keys():
            Tf = self.subprocess['albedo'].param['Tf']
        else:
            print('No ice considered in this model. Can not plot.')
            return

        num_steps_per_year = int(self.time['num_steps_per_year'])
        ice_year = np.empty((self.lat.size, num_steps_per_year))
        for m in range(num_steps_per_year):
            self.step_forward()
            ice_year[:,m] = np.where(np.squeeze(self.Ts) <= Tf, 0, 1)
        
        fig = plt.figure(figsize=(5,5))
        ax = fig.add_subplot(111)
        
        factor = 365. / num_steps_per_year
        cax = ax.contourf(factor * np.arange(num_steps_per_year),
                      self.lat, ice_year[:,:], 
                      cmap=plt.cm.seismic, vmin=0, vmax=1, levels=2)
        cbar1 = plt.colorbar(cax)
        ax.set_title('Ice throughout the year', fontsize=14)
        ax.set_xlabel('Days of year', fontsize=14)
        ax.set_ylabel('Latitude', fontsize=14)

Quick Overview

Climlab is a flexible package, where you can define your own climate models by coupling together different processes. But it has certain types of models implemented already. Setting up a climate model in climlab can be as easy as this:

In [2]:
ebm = climlab.EBM()

Let's print what we just created:

In [5]:
print(ebm)

climlab Process of type <class 'climlab.model.ebm.EBM'>. 
State variables and domain shapes: 
  Ts: (90, 1) 
The subprocess tree: 
Untitled: <class 'climlab.model.ebm.EBM'>
   LW: <class 'climlab.radiation.aplusbt.AplusBT'>
   insolation: <class 'climlab.radiation.insolation.P2Insolation'>
   albedo: <class 'climlab.surface.albedo.StepFunctionAlbedo'>
      iceline: <class 'climlab.surface.albedo.Iceline'>
      warm_albedo: <class 'climlab.surface.albedo.P2Albedo'>
      cold_albedo: <class 'climlab.surface.albedo.ConstantAlbedo'>
   SW: <class 'climlab.radiation.absorbed_shorwave.SimpleAbsorbedShortwave'>
   diffusion: <class 'climlab.dynamics.meridional_heat_diffusion.MeridionalHeatDiffusion'>

We see that the mdel has one state variable: Ts, that is the surface temperature. The domain where it is defined on has shape (90, 1), that means it has 90 points in latitudinal direction and 1 point in the depth-direction, i.e. downwards. So this is a 1D model now.

Important note: The model is assuming we just have a planet of water of uniform depth. We can get the water depth by ebm.param, which also gives us some other parameters, for example S0, the solar constant.

In [6]:
ebm.param
Out[6]:
{'timestep': 350632.51200000005,
 'S0': 1365.2,
 's2': -0.48,
 'A': 210.0,
 'B': 2.0,
 'D': 0.555,
 'Tf': -10.0,
 'water_depth': 10.0,
 'a0': 0.3,
 'a2': 0.078,
 'ai': 0.62}

The model consists of subprocesses. For example, it has a process 'LW' which is computing the outgoing longwave radiation or 'albedo' which is calculating the planetary albedo at each latitude, or 'diffusion' for heat transport. Furthermore, there is a process called 'insolation'. Currently, it is just approximating the annual mean insolation at each latitude by a polynomial which only depends on the latitude but not on the time. So there are no seasonal variations yet. This is the process we will adjust first:

Seasonal EBM

Let us first look at a seasonal EBM. Climlab has a class EBM_seasonal implemented, which we adjust for our purposes. We first change the insolation process as well as the process for absorbed shortwave radiation since that builds up on the insolation. Also, to simplify, we assume a constant albedo value (0.33). At the beginning we defined a class that does this for you, as well as defining some methods plotting. You don't need to worry about the details but if you are curious you can check out the source code at the top of this notebook.

You can now create a simple seasonal EBM by calling:

In [4]:
ebm_seasonal = simple_EBM_seasonal()

Again, print the model:

In [5]:
print(ebm_seasonal)

climlab Process of type <class '__main__.simple_EBM_seasonal'>. 
State variables and domain shapes: 
  Ts: (90, 1) 
The subprocess tree: 
Untitled: <class '__main__.simple_EBM_seasonal'>
   diffusion: <class 'climlab.dynamics.meridional_heat_diffusion.MeridionalHeatDiffusion'>
   LW: <class 'climlab.radiation.aplusbt.AplusBT_CO2'>
   insolation: <class 'climlab.radiation.insolation.DailyInsolation'>
   SW: <class 'climlab.radiation.absorbed_shorwave.SimpleAbsorbedShortwave'>
   albedo: <class 'climlab.surface.albedo.ConstantAlbedo'>

By default, this model was given a temperature profile, which we can plot using our own new method. We can also look at the radiations: OLR = outgoing longwave radiation. ASR = absorbed shortwave radiation. Call a plot function as follows:

In [13]:
ebm_seasonal.plot_temp_rad()

Observations: The OLR (in the right graph) is symmetric around the equator just like the temperature. This makes sense: in this model, OLR is modeled as a linear function of temperature. Note that the model starts at January 1, so it is winter on the northern hemisphere. At the north pole, there is no incoming sunlight.

Question: Where is the earth gaining energy and where is it loosing energy?

So far, this is just the initial state at the time of model creation. Now we let the model involve in time. We can step forward, for example for 30 days.

Question: What do you expect? Where will it get colder and where warmer?

To try this, run the following cell to see the output. In order to run that cell, the packages need to be loaded and the variables defined. In Jupyter Lab, there is an option Run -> Run All Above Selected Cell. In Jupyter Notebook, this is Cell -> Run All Above. Through this, the above cells are run and all the setups are done. After completing this, you can make your way downwards through the notebook by running the selected cell and advancing.

In [14]:
ebm_seasonal.integrate_days(30.)
ebm_seasonal.plot_temp_rad()

Integrating for 7 steps, 30.000000000000004 days, or 0.08213727767492365 years.
Total elapsed time is 16.366666666666255 years.

Task: Step forward further in time through additional calls of ebm_seasonal.integrate_days(some number of days) and check those profiles at the equinoxes and the summer and winter solstices.

Hints:

  • Note that the model time is now already at January 30. Whenever you call one of the methods that step further in time, the current time will be saved. Should you loose track of how many days are already elapsed, you can call ebm_seasonal.time to print out some features. (Currently we have 90 timesteps per year, so this might lead to some non-integer numbers for the elapsed days and days of the year.)
  • Remember that the model is initialised with some temperature distribution which is different from the equilibrium distribution that will be reached after some years.
  • You can also call the method integrate_years() to get closer to that equilibrium, for example ebm_seasonal.integrate_years(5.0)
  • Our own class has a method to plot temperature distributions throughout the year: plot_temp_year(), for example ebm_seasonal.plot_temp_year(). It will step forward for 365 days and and make a contour plot.
In [16]:
# your code here
# suggestion

# start again at Jan 1
ebm_seasonal = simple_EBM_seasonal()
# first integrate forward some years to get closer to equilibrium
ebm_seasonal.integrate_years(5.0)
# look at March equinox
ebm_seasonal.integrate_days(31+28+20)
print('March 20')
ebm_seasonal.plot_temp_rad()
# look at summer solstice (June 20)
ebm_seasonal.integrate_days(31+30+31)
print('June 20')
ebm_seasonal.plot_temp_rad()
# September equinox (Sept 20)
ebm_seasonal.integrate_days(30+31+30)
print('Sept 20')
ebm_seasonal.plot_temp_rad()
# winter solstice (Dec 20)
ebm_seasonal.integrate_days(30+31+30)
print('Dec 20')
ebm_seasonal.plot_temp_rad()

Integrating for 450 steps, 1826.2110000000002 days, or 5.0 years.
Total elapsed time is 5.000000000000044 years.
Integrating for 19 steps, 79.0 days, or 0.21629483121063228 years.
Total elapsed time is 5.211111111111158 years.
March 20

Integrating for 22 steps, 91.99999999999999 days, or 0.2518876515364325 years.
Total elapsed time is 5.4555555555556055 years.
June 20

Integrating for 22 steps, 91.0 days, or 0.24914974228060174 years.
Total elapsed time is 5.700000000000048 years.
Sept 20

Integrating for 22 steps, 91.0 days, or 0.24914974228060174 years.
Total elapsed time is 5.944444444444482 years.
Dec 20
In [19]:
# check where we are at during the year
ebm_seasonal.time
Out[19]:
{'timestep': 350632.51200000005,
 'num_steps_per_year': 90.0,
 'day_of_year_index': 85,
 'steps': 625,
 'days_elapsed': 2536.404166666665,
 'years_elapsed': 6,
 'days_of_year': array([  0.        ,   4.05824667,   8.11649333,  12.17474   ,
         16.23298667,  20.29123333,  24.34948   ,  28.40772667,
         32.46597333,  36.52422   ,  40.58246667,  44.64071333,
         48.69896   ,  52.75720667,  56.81545333,  60.8737    ,
         64.93194667,  68.99019333,  73.04844   ,  77.10668667,
         81.16493333,  85.22318   ,  89.28142667,  93.33967333,
         97.39792   , 101.45616667, 105.51441333, 109.57266   ,
        113.63090667, 117.68915333, 121.7474    , 125.80564667,
        129.86389333, 133.92214   , 137.98038667, 142.03863333,
        146.09688   , 150.15512667, 154.21337333, 158.27162   ,
        162.32986667, 166.38811333, 170.44636   , 174.50460667,
        178.56285333, 182.6211    , 186.67934667, 190.73759333,
        194.79584   , 198.85408667, 202.91233333, 206.97058   ,
        211.02882667, 215.08707333, 219.14532   , 223.20356667,
        227.26181333, 231.32006   , 235.37830667, 239.43655333,
        243.4948    , 247.55304667, 251.61129333, 255.66954   ,
        259.72778667, 263.78603333, 267.84428   , 271.90252667,
        275.96077333, 280.01902   , 284.07726667, 288.13551333,
        292.19376   , 296.25200667, 300.31025333, 304.3685    ,
        308.42674667, 312.48499333, 316.54324   , 320.60148667,
        324.65973333, 328.71798   , 332.77622667, 336.83447333,
        340.89272   , 344.95096667, 349.00921333, 353.06746   ,
        357.12570667, 361.18395333]),
 'active_now': True}
In [20]:
# we need five further steps to get to the start of the year again
for i in range(5):
    ebm_seasonal.step_forward()
ebm_seasonal.time
Out[20]:
{'timestep': 350632.51200000005,
 'num_steps_per_year': 90.0,
 'day_of_year_index': 0,
 'steps': 630,
 'days_elapsed': 2556.6953999999973,
 'years_elapsed': 7,
 'days_of_year': array([  0.        ,   4.05824667,   8.11649333,  12.17474   ,
         16.23298667,  20.29123333,  24.34948   ,  28.40772667,
         32.46597333,  36.52422   ,  40.58246667,  44.64071333,
         48.69896   ,  52.75720667,  56.81545333,  60.8737    ,
         64.93194667,  68.99019333,  73.04844   ,  77.10668667,
         81.16493333,  85.22318   ,  89.28142667,  93.33967333,
         97.39792   , 101.45616667, 105.51441333, 109.57266   ,
        113.63090667, 117.68915333, 121.7474    , 125.80564667,
        129.86389333, 133.92214   , 137.98038667, 142.03863333,
        146.09688   , 150.15512667, 154.21337333, 158.27162   ,
        162.32986667, 166.38811333, 170.44636   , 174.50460667,
        178.56285333, 182.6211    , 186.67934667, 190.73759333,
        194.79584   , 198.85408667, 202.91233333, 206.97058   ,
        211.02882667, 215.08707333, 219.14532   , 223.20356667,
        227.26181333, 231.32006   , 235.37830667, 239.43655333,
        243.4948    , 247.55304667, 251.61129333, 255.66954   ,
        259.72778667, 263.78603333, 267.84428   , 271.90252667,
        275.96077333, 280.01902   , 284.07726667, 288.13551333,
        292.19376   , 296.25200667, 300.31025333, 304.3685    ,
        308.42674667, 312.48499333, 316.54324   , 320.60148667,
        324.65973333, 328.71798   , 332.77622667, 336.83447333,
        340.89272   , 344.95096667, 349.00921333, 353.06746   ,
        357.12570667, 361.18395333]),
 'active_now': True}
In [21]:
# maybe now try to look at distribution throughout the year
ebm_seasonal.plot_temp_year()

Varying Orbit Parameters

Next, we want to learn more about the influence of different orbit paramters. The class simple_EBM_seasonal() can take a parameter orbit for initialisation. This should be a dictionary of the form orbit = {'ecc': value, 'long_peri': value, 'obliquity': value} and by default it is set to the current values: orbit = {'ecc': 0.0167643, 'long_peri': 280.32687, 'obliquity': 23.459277}.

In [44]:
# feel free to adjust these parameters later
orbit = {'ecc': 0.0167643, 'long_peri': 280.32687, 'obliquity': 23.459277}
ebm2 = simple_EBM_seasonal(orbit)

# step forward in time for some years so model is in equilibrium
ebm2.integrate_years(10.)

Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.

Task: You can now adjust orbit parameters and see how this influences the temperature profile and the radiations throughout the year.

Hints:

  • Don't forget to rerun the above cell so the variable ebm2 is defined (and rerun it every time you change the orbit parameters).
  • One handy feature of climlab process 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 (and many other diagnostics) saved in the dictionary timeave.

  • The method plot_timeave() of our own class simple_EBM_seasonal() uses this: it calls integrate_years(1.) once and then plots the quantities from timeave. Example use:

      ebm2.plot_timeave().

  • Climlab classes also have a method for computing the global mean temperature: global_mean_temperature(). But this only refers to the one point in time where the model is currently at.

  • Hence, our own class includes another method global_annual_mean() which steps forward in time and computes the global mean at each step, later computing the mean of these values to get the annual mean. Example usage:

      print(ebm2.gobal_annual_mean())
In [24]:
# more code here
# step forward in time, use plotting functions or functions to compute the global, annual mean

# suggestions
# make some very tilted example
print('An example with obl = 90')
orbit = {'ecc': 0.0167643, 'long_peri': 280.32687, 'obliquity': 90}
ebm3 = simple_EBM_seasonal(orbit)
ebm3.integrate_years(10.)
ebm3.plot_timeave()
ebm3.plot_temp_year()
print('The global annual mean is {} \n'.format(ebm3.global_annual_mean()))

# what about no tilt at all?
print('An example with obl = 0')
orbit = {'ecc': 0.0167643, 'long_peri': 280.32687, 'obliquity': 0}
ebm4 = simple_EBM_seasonal(orbit)
ebm4.integrate_years(10.)
ebm4.plot_timeave()
ebm4.plot_temp_year()
print('The global annual mean is {}'.format(ebm4.global_annual_mean()))

An example with obl = 90
Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.
Integrating for 90 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 10.999999999999819 years.

The global annual mean is 11.082443957471053 

An example with obl = 0
Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.
Integrating for 90 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 10.999999999999819 years.

The global annual mean is 11.075053175391721
In [27]:
# Now let's vary the eccentricity - what happens with a very eccentric orbit?
print('An example with ecc = 0.2')
orbit = {'ecc': 0.2, 'long_peri': 280.32687, 'obliquity': 23}
ebm5 = simple_EBM_seasonal(orbit)
ebm5.integrate_years(10.)
ebm5.plot_timeave()
ebm5.plot_temp_year()
print('The global annual mean is {}\n'.format(ebm5.global_annual_mean()))

# And on a circular orbit?
print('An example with ecc = 0')
orbit = {'ecc': 0, 'long_peri': 280.32687, 'obliquity': 23}
ebm6 = simple_EBM_seasonal(orbit)
ebm6.integrate_years(10.)
ebm6.plot_timeave()
ebm6.plot_temp_year()
print('The global annual mean is {}'.format(ebm6.global_annual_mean()))

An example with ecc = 0.2
Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.
Integrating for 90 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 10.999999999999819 years.

The global annual mean is 13.471688370209217

An example with ecc = 0
Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.
Integrating for 90 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 10.999999999999819 years.

The global annual mean is 11.061493390088481
In [32]:
# lastly, change the longitude of the perihelion
print('An example with long_peri = 0 and still an eccentric orbit')
orbit = {'ecc': 0.2, 'long_peri': 0, 'obliquity': 23}
ebm7 = simple_EBM_seasonal(orbit)
ebm7.integrate_years(10.)
ebm7.plot_timeave()
ebm7.plot_temp_year()
plt.show() # this line ensures that the plots and print statements are in the correct order
ebm7.integrate_days(31+28+20)
print('Looking at March equinox')
ebm7.plot_temp_rad()
print('The global annual mean is {}\n'.format(ebm7.global_annual_mean()))

# And on a circular orbit?
print('An example with long_peri = 180 and still an eccentric orbit')
orbit = {'ecc': 0.2, 'long_peri': 180, 'obliquity': 23}
ebm8 = simple_EBM_seasonal(orbit)
ebm8.integrate_years(10.)
ebm8.plot_timeave()
ebm8.plot_temp_year()
plt.show()  # this line ensures that the plots and print statements are in the correct order
ebm8.integrate_days(31+28+20)
print('Looking at March equinox')
ebm8.plot_temp_rad()
print('The global annual mean is {}'.format(ebm8.global_annual_mean()))

An example with long_peri = 0 and still an eccentric orbit
Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.
Integrating for 90 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 10.999999999999819 years.

Integrating for 19 steps, 79.0 days, or 0.21629483121063228 years.
Total elapsed time is 12.21111111111088 years.
Looking at March equinox

The global annual mean is 13.471706372410656

An example with long_peri = 180 and still an eccentric orbit
Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.
Integrating for 90 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 10.999999999999819 years.

Integrating for 19 steps, 79.0 days, or 0.21629483121063228 years.
Total elapsed time is 12.21111111111088 years.
Looking at March equinox

The global annual mean is 13.4717062378678

Changing the Albedo

So far, we've made a simplifying assumption: The albedo value was set to a constant. This can be seen by using the method plot_albedo() of our class to plot the albedo.

In [34]:
ebm = simple_EBM_seasonal()
ebm.plot_albedo()

Now we will change this. Climlab implements a step function for the albedo value. At each timestep, the albedo is set depending on the temperature:

  • If the temperature is below some threshold value (by default -10.0°C), the model assumes that there is ice and sets the albedo constant to 0.62.
  • If the temperature is above this threshold, there is assumed to be no ice and the albedo is approximated by a degree-two-polynomial.

Note that changes between ice and no-ice are instantaneous, i.e. there is no melting process or built-up time for ice sheets. The step function albedo can be included by setting albedo='step' in the initialisation of our class. This gives the following picture for our initial temperature curve:

In [33]:
ebm_stepalb = simple_EBM_seasonal(albedo='step')
ebm_stepalb.plot_albedo()

Task: Try to repeat some of the above analysis for this new albedo function. How do the temperature distributions change throughout the year?

Hints:

  • plot_ice_year() is a simple method to plot the ice distribution throughout the year. For example:

      ebm_stepalb.plot_ice_year()

  • The threshold value specifying when the albedo is set to the ice-value can be changed by giving the argument Tf to the constructor. For example

      ebm_stepalb = simple_EBM_seasonal(albedo='step', Tf=-5.0).

  • Instead of using the step function albedo, you can also set albedo=tanh for an implementation of the formula

      albedo = 0.42 - 0.20 * np.tanh(0.052 * (Ts-3)).
In [37]:
# your code here
# step forward in time, use plotting functions, look at mean temperatures etc.
# maybe change the Tf parameter and see how this influences the ice and the temperatures

# suggestion
# old seasonal with constant albedo for reference
ebm_constalb = simple_EBM_seasonal(albedo='constant')
print('Constant albedo')
ebm_constalb.integrate_years(10.)
ebm_constalb.plot_timeave()
ebm_constalb.plot_temp_year()
plt.show()  # this line ensures that the plots and print statements are in the correct order

print()
print('Step Function albedo')
ebm_stepalb.integrate_years(10.)
ebm_stepalb.plot_timeave()
ebm_stepalb.plot_temp_year()
ebm_stepalb.plot_ice_year()

Constant albedo
Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.
Integrating for 90 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 10.999999999999819 years.

Step Function albedo
Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 34.000000000001414 years.
Integrating for 90 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 35.0000000000016 years.
In [38]:
# some example with more ice

ebm_ice = simple_EBM_seasonal(albedo='step', Tf = -3.0)
ebm_ice.integrate_years(10.)
ebm_ice.plot_timeave()
ebm_ice.plot_temp_year()
ebm_ice.plot_ice_year()

Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.
Integrating for 90 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 10.999999999999819 years.

Investigating the Effect of CO2 concentration

So far, we have played around with orbit parameters and albedo. Next, we will adjust the levels of CO2 in our atmosphere. The CO2 concentration (in ppm) can be changed by adjusting the parameter CO2, for example:

ebm_co2 = simple_EBM_seasonal(albedo='step', CO2=411)

to set it to the 2019 value.

Aside: The seasonal EBM in Climlab does NOT implement a complex radiation model. The relationship between temperature and outgoing longwave radiation is simply parametrised by

OLR = A + B * T,

where A and B depend on the CO2 concentration. If you want to use a more realistic radiation model, you will have to work with Climlab's column models. However, these do not have any horizontal heat transport.

Task: Adjust the CO2 value to todays value - or even higher to represent the future. Possibly integrate forward in time for some years to bring the model to equilibrium. See how the temperatures / ice / radiations change. Feel free to try different albedo parameters also.

Further suggestion: It could also be interesting to calculate a climate sensitivity, i.e. the temperature change from doubling the CO2 value in the atmosphere. You could define two models:

ebm_co2_low = simple_EBM_seasonal(albedo='step', CO2=280)
ebm_co2_high = simple_EBM_seasonal(albedo='step', CO2=560)

Then for each of them, call the integrate_years() method for a few years, maybe around 10, and compare the global annual mean with

ebm_co2_low.global_annual_mean()
ebm_co2_high.global_annual_mean().

How much warmer did it get?

In [41]:
# your code here
# feel free to change the CO2 value or define several models with different concentrations and compare

# co2 concentration of today
print('Co2 concentration of today')
ebm_co2 = simple_EBM_seasonal(albedo='step', CO2=411)
ebm_co2.integrate_years(10.)
ebm_co2.plot_timeave()
ebm_co2.plot_temp_year()
ebm_co2.plot_ice_year()
plt.show()
print('The global annual mean is {}'.format(ebm_co2.global_annual_mean()))

print()
print('Co2 in 2100 in RCP 8.5 (BAU)')
# co2 concentration in 2100 with business as usual scenario
ebm_co2 = simple_EBM_seasonal(albedo='step', CO2=935)
ebm_co2.integrate_years(10.)
ebm_co2.plot_timeave()
ebm_co2.plot_temp_year()
ebm_co2.plot_ice_year()
plt.show()
print('The global annual mean is {}'.format(ebm_co2.global_annual_mean()))

Co2 concentration of today
Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.
Integrating for 90 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 10.999999999999819 years.

The global annual mean is 14.551232715545078

Co2 in 2100 in RCP 8.5 (BAU)
Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.
Integrating for 90 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 10.999999999999819 years.

The global annual mean is 16.73907993835788

Recreating past times

Here is an overview of parameters for some past times:

Simulation Name Time CO2 value Eccentricity Obliquity Precession
Pre-Industrial 1800 284.6 0.0167643 23.459277 280.32687
Future 560 0.0167643 23.459277 280.32687
Last Interglacial 127 k 278 0.038231 24.2441 49.097
Mid-Holocene Warm Period 6 k 280 0.018682 24.1048 180.918

Task: You could try to feed those parameters into the model and see how temperatures were like in the past.

In [44]:
# your code here

# look at Last Interglacial
orbit = {'ecc': 0.038231, 'long_peri': 49.097, 'obliquity': 24.2441}
ebm = simple_EBM_seasonal(orbit, CO2=278, albedo='step')
ebm.integrate_years(10.)
ebm.plot_timeave()
ebm.plot_temp_year()
ebm.plot_ice_year()

Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.
Integrating for 90 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 10.999999999999819 years.

Influence of water depth

As mentioned above, in Climlab we just have a planet of water of uniform water depth. The water depth is used to calculate the heat capacity. It is an interesting experiment to change the value and see how this effects seasonal variation. You can try this by specifying the parameter water_depth when setting up the model, for example

ebm_deep = simple_EBM_seasonal(albedo='step', water_depth=100).

Task: Define different models with different water depths and plot the seasonal changes with plot_timeave() or plot_ice_year() (the latter only if you used albedo='step'). Compare.

In [42]:
# your code here

ebm_deep = simple_EBM_seasonal(albedo='step', water_depth=100)
ebm_deep.integrate_years(10.)
ebm_deep.plot_temp_year()

Integrating for 900 steps, 3652.4220000000005 days, or 10.0 years.
Total elapsed time is 9.999999999999863 years.