Climate Dynamics#

!pip install climlab
!wget https://clark.physics.illinois.edu/temp.data

Collecting climlab

  Downloading climlab-0.8.2.tar.gz (3.1 MB)
?25l     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/3.1 MB ? eta -:--:--
     ━━━━━╸━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.5/3.1 MB 13.9 MB/s eta 0:00:01

     ━━━━━━━━━━━━━━━━━╺━━━━━━━━━━━━━━━━━━━━━━ 1.3/3.1 MB 19.3 MB/s eta 0:00:01

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━ 2.4/3.1 MB 23.7 MB/s eta 0:00:01
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 3.1/3.1 MB 23.6 MB/s eta 0:00:00
?25h

  Installing build dependencies ... ?25l-

 \

 |

 /

 -

 done

?25h  Getting requirements to build wheel ... ?25ldone

?25h  Preparing metadata (pyproject.toml) ... ?25l-

 done
?25hRequirement already satisfied: numpy in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from climlab) (1.24.4)

Collecting xarray (from climlab)
  Downloading xarray-2023.1.0-py3-none-any.whl (973 kB)
?25l     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/973.1 kB ? eta -:--:--
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 973.1/973.1 kB 41.9 MB/s eta 0:00:00
?25h

Collecting scipy (from climlab)
  Downloading scipy-1.10.1-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (34.5 MB)
?25l     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/34.5 MB ? eta -:--:--

     ━━╺━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 2.1/34.5 MB 63.3 MB/s eta 0:00:01

     ━━━━━╸━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 4.9/34.5 MB 72.3 MB/s eta 0:00:01
     ━━━━━━━━╸━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 7.7/34.5 MB 74.8 MB/s eta 0:00:01

     ━━━━━━━━━━━━━━╺━━━━━━━━━━━━━━━━━━━━━━━━━ 12.1/34.5 MB 95.2 MB/s eta 0:00:01

     ━━━━━━━━━━━━━━━━━━╸━━━━━━━━━━━━━━━━━━━━ 16.5/34.5 MB 123.8 MB/s eta 0:00:01
     ━━━━━━━━━━━━━━━━━━━━━━━╸━━━━━━━━━━━━━━━ 21.2/34.5 MB 130.5 MB/s eta 0:00:01

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━ 25.9/34.5 MB 131.4 MB/s eta 0:00:01

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸━━━━ 30.7/34.5 MB 135.0 MB/s eta 0:00:01
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 34.5/34.5 MB 147.1 MB/s eta 0:00:01

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 34.5/34.5 MB 147.1 MB/s eta 0:00:01

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 34.5/34.5 MB 69.9 MB/s eta 0:00:00
?25h

Collecting pandas>=1.3 (from xarray->climlab)
  Downloading pandas-2.0.3-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (18 kB)
Requirement already satisfied: packaging>=21.3 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from xarray->climlab) (23.2)

Requirement already satisfied: python-dateutil>=2.8.2 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from pandas>=1.3->xarray->climlab) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from pandas>=1.3->xarray->climlab) (2023.3.post1)

Collecting tzdata>=2022.1 (from pandas>=1.3->xarray->climlab)
  Downloading tzdata-2023.3-py2.py3-none-any.whl (341 kB)
?25l     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/341.8 kB ? eta -:--:--
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 341.8/341.8 kB 66.5 MB/s eta 0:00:00
?25h

Requirement already satisfied: six>=1.5 in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from python-dateutil>=2.8.2->pandas>=1.3->xarray->climlab) (1.16.0)
Downloading pandas-2.0.3-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (12.4 MB)
?25l   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/12.4 MB ? eta -:--:--

   ━━━━━━━━━━━━━━━━━━╺━━━━━━━━━━━━━━━━━━━━━ 5.7/12.4 MB 169.9 MB/s eta 0:00:01
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸━━ 11.7/12.4 MB 170.5 MB/s eta 0:00:01

   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 12.4/12.4 MB 112.2 MB/s eta 0:00:00
?25h

Building wheels for collected packages: climlab

  Building wheel for climlab (pyproject.toml) ... ?25l-

 done
?25h  Created wheel for climlab: filename=climlab-0.8.2-py3-none-any.whl size=171440 sha256=348388e6f52f9b7e1f4e4b2d209fb686e1beb5f84aa2f06b2298dcde02dc1af7
  Stored in directory: /home/runner/.cache/pip/wheels/bd/07/56/ff3ce2e43c85c96acebbda20e5d72cdb0506dbee6ffad49364
Successfully built climlab

Installing collected packages: tzdata, scipy, pandas, xarray, climlab

Successfully installed climlab-0.8.2 pandas-2.0.3 scipy-1.10.1 tzdata-2023.3 xarray-2023.1.0

[notice] A new release of pip is available: 23.0.1 -> 23.3.1
[notice] To update, run: pip install --upgrade pip

--2023-12-01 18:07:05--  https://clark.physics.illinois.edu/temp.data
Resolving clark.physics.illinois.edu (clark.physics.illinois.edu)... 

18.220.149.166
Connecting to clark.physics.illinois.edu (clark.physics.illinois.edu)|18.220.149.166|:443... connected.

HTTP request sent, awaiting response... 200 OK
Length: 17595 (17K)
Saving to: ‘temp.data’


temp.data             0%[                    ]       0  --.-KB/s               
temp.data           100%[===================>]  17.18K  --.-KB/s    in 0.01s   

2023-12-01 18:07:05 (1.16 MB/s) - ‘temp.data’ saved [17595/17595]

import numpy as np
import matplotlib.pyplot as plt
import random
from climlab import constants as const
from climlab.solar.insolation import daily_insolation
from climlab.utils import legendre
import climlab
import scipy
def resetMe(keepList=[]):
    ll=%who_ls
    keepList=keepList+['resetMe','np','plt','random']
    for iiii in keepList:
        if iiii in ll:
            ll.remove(iiii)
    for iiii in ll:
        jjjj="^"+iiii+"$"
        %reset_selective -f {jjjj}
    ll=%who_ls
    return
tempData=np.loadtxt("temp.data")
tempData[:,1]=tempData[:,1]+14.2
import datetime;datetime.datetime.now()

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[2], line 4
      2 import matplotlib.pyplot as plt
      3 import random
----> 4 from climlab import constants as const
      5 from climlab.solar.insolation import daily_insolation
      6 from climlab.utils import legendre

File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/climlab/__init__.py:20
     18 from .utils import constants, thermo, legendre
     19 # some more useful shorcuts
---> 20 from .model.column import GreyRadiationModel, RadiativeConvectiveModel, BandRCModel
     21 from .model.ebm import EBM, EBM_annual, EBM_seasonal
     22 from .domain.field import Field, global_mean

File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/climlab/model/column.py:43
     41 from climlab.domain.initial import column_state
     42 from climlab.domain.field import Field
---> 43 from climlab.radiation.insolation import FixedInsolation
     44 from climlab.radiation.greygas import GreyGas, GreyGasSW
     45 from climlab.convection.convadj import ConvectiveAdjustment

File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/climlab/radiation/__init__.py:15
     13 from .water_vapor import ManabeWaterVapor
     14 #from radiation import Radiation, Radiation_SW, Radiation_LW
---> 15 from .cam3 import CAM3, CAM3_LW, CAM3_SW
     16 from .rrtm import RRTMG, RRTMG_LW, RRTMG_SW

File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/climlab/radiation/cam3.py:39
     37 from climlab import constants as const
     38 from climlab.utils.thermo import vmr_to_mmr
---> 39 from climlab.radiation.radiation import _Radiation_SW, _Radiation_LW
     40 import warnings
     42 #  Wrapping these imports in try/except to avoid failures during documentation building on readthedocs

File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/climlab/radiation/radiation.py:86
     84 from climlab import constants as const
     85 from climlab.domain.field import Field
---> 86 import os, warnings, pooch
     87 import xarray as xr
     88 from climlab.utils import _datapath_http

ModuleNotFoundError: No module named 'pooch'

Exercise 1: Climate and Global Warming#

a. The temperature of the dark airless Earth#

Our goal is to start understanding something about the climate. We first want two functions:

def IncomingRadiationPerSquareMeter():
  ...
  ...
  return ...

which returns the incoming radiation to the Earth and

def OutgoingRadiationPerSquareMeter():
  ...
  ...
  return ...

which returns the outgoing radiation.

We are going to iteratively improve upon these functions to add more realism eventually ending up with a one-dimensional representation of the Earth.

Let us begin, though, with the simplest thing we can do. The sun is outputting a certain amount of Energy per second. By the time it has spread out to the orbit of the Earth, it is a flux of

  • \(S=1367.6\) Watts/\(\textrm{m}^2\)

From the Sun’s point of view, the earth is a disc (circle) that is of area \(\pi R^2\) where \(R\) is the radius of the Earth.

From this, we can compute the amount of radiation that is hitting Earth. That radiation, though, is split out over the entire surface of the Earth - an area of \(4\pi R^2\).

Using this information you should be able to write a function IncomingRadiationPerSquareMeter() which specifies how much radiation hits the Earth. Notice that this doesn’t depend on the radius of the Earth \(R\).

Answer (start)
## Answer here
Answer (end)

To get an estimate of the outgoing radiation, we can use the Stefan-Boltzmann law. This law tell us that a blackbody (like the Earth) radiates at \(\sigma T^4\) where

  • \(\sigma = 5.67 \times 10^{-8}\) W/ (\(m^2\) \(K^4\))

We can now right the function OutgoingRadiationPerSquareMeter(T) which takes the temperature \(T\).

Answer (start)
### Answer here
Answer (end)

Now we need to find the temperature at which the incoming radiation and the outgoing radiation match.

Do this in two ways. First graph the two functions you’ve written as a function of temperature and determine at temperature \(T\) the two match.

Secondly, let’s just call some optimization. Define a cost function diff(T) which gives the absolute value of the difference between the incoming and outgoing radiation. Let’s then just minimize this by calling opt=scipy.optimize.minimize(diff,[10],method='COBYLA',bounds=[(0,300)])

Report the number you find and compare this with the actual temperature of the Earth which is approximately 14 degrees Celsius (pre-industrial) or 15 degrees Celsius (now).

Answer (start)
# Answer here
Answer (end)

b. The temperature of Earth with air and clouds#

It turns out that we didn’t find the correct temperature. We have ignored two important effects.

Albedo

Some fraction of the radiation that is incoming to the Earth simply gets reflected back immediately from things like clouds and ice. This fraction is approximately \(\alpha=0.3\). The incoming radiation per square meter is then instead \(\frac{S}{4}( 1-\alpha)\).

Greenhouse Gases

In practice, some fraction of the radiation that would typically be radiated back into space is capture by the greenhouse gases in the atmosphere.

We can represent this by including an emissivity term for the atmosphere writing instead that the radiation emitted is \(\epsilon \sigma T^4\). We don’t have any good way to calculate \(\epsilon\) from first principles at the moment so instead let’s back it out. Let’s instead say that the pre-industrial temperature is 14 Kelvin. Calculate using your code what the emissivity have to be?

Answer (start)
### Answer Here
Answer (end)

c. Linearizing around the equilibrium#

We will find it convenient to linearize are outgoing radiation around zero degrees celsius. This should give us

\[\epsilon \sigma (T_0 + T_c)^4 \approx \epsilon \sigma(T_0^4 + 4 T_0^3 T_c) = A + BT_c\]

where \(T_0=273.15\) and \(T_c\) is the temperature in Celsius.

This gives us \(A=195.94\) and \(B=2.87\).

Change your outgoing radiation function to use \(A+BT_c\). Your function should take \(T\) in Kelvin but needs to convert it to Celsius to work here.

Again, optimize for the equilibrium temperature. You’ll actually find that we get a temperature \(T=15.11878309\). This is off because the linearization wasn’t great. We can fix this by instead renormalizing these parameters instead using

  • \(A=221.2\)

  • \(B=1.3\)

Check that these new parameters give a much closer equilbrium temperature (for the pre-industrial average) of 13.94

Answer (start)
### Answer Here
Answer (start)

We will be building up a dictionary with the relevant parameters in it - so far

params = {
    'alpha' : 0.3,
    'A' : 221.2  ,
    'B' : 1.3
    }

d. Global Warming#

We know that putting carbon dioxide in the atmosphere increases the temperature of the Earth. We would like to quantify this. It is known that Carbox dioxide essentially contributes to decreasing the radiation which is emitted from the Earth in a way which is logarithmic to its concentration. In other words, it corrects our emitted radiation by a term $\( r \log \left(\frac{CO_{2;current}}{CO_{2;old}}\right) \)\( where \)r=-5$ is a proportionality constant.

Below is the Carbon dioxide that is in the air as a function of year. Plot it:

CO2=np.array([[1751.590620307266, 277.1900826446281],
[1777.9960515400046, 278.95316804407713],
[1790.6841547040233, 279.8347107438017],
[1797.1997437701493, 280.27548209366387],
[1806.8001722200638, 281.1570247933884],
[1815.0283007970431, 282.03856749311296],
[1822.9161897361043, 282.4793388429752],
[1832.8625284845687, 282.92011019283746],
[1845.5506316485873, 283.80165289256195],
[1859.9484389931429, 285.56473829201104],
[1872.9711111344473, 287.7685950413223],
[1891.8062104235146, 293.0578512396694],
[1906.8674850620098, 298.3471074380165],
[1917.128545475548, 303.1955922865013],
[1926.3717223056487, 306.2809917355372],
[1930.8288615623721, 306.72176308539946],
[1936.3010490722172, 309.366391184573],
[1940.7440116773605, 312.0110192837465],
[1946.573450807019, 312.45179063360877],
[1951.3708297016603, 313.3333333333333],
[1955.8194629674358, 315.09641873278235],
[1959.9221859346617, 317.30027548209364],
[1963.6789986033375, 319.9449035812672],
[1967.0927363037795, 322.58953168044076],
[1971.1784472891093, 327.4380165289256],
[1973.903199722768, 330.52341597796146],
[1976.9625211337122, 334.93112947658403],
[1980.6994864901762, 340.6611570247934],
[1982.735253657051, 344.1873278236915],
[1985.4571707603936, 347.71349862258955],
[1987.8303422348704, 352.1212121212121],
[1990.8896636458148, 356.5289256198347],
[1994.2892246946772, 361.37741046831957],
[1996.322156531236, 365.3443526170799],
[1999.3786426118645, 370.1928374655647],
[2003.7762399319522, 379.88980716253445],
[2006.4896510443468, 384.73829201101927],
[2010.5498440568329, 393.55371900826447],
[2012.2397009251579, 397.5206611570248],
[2014.9361000556567, 405.0137741046832],
[2016.9690318922155, 408.9807162534435],
[2019.6654310227143, 416.4738292011019],
[2020.6833146061517, 418.23691460055096]])

### Plot Here

Now, modify your code to include this additional forcing term into the emitted radiation.

Compute the equilibrium temperature at the value of carbon dioxide for each year and graph this equilibrium temperature against the year. Plot it on top of the actual measured temperature which is in tempData[:,1] (where the uear is in tempData[:,0]). Notice the close agreement between the predicted temperature based on the CO2 and our 0D EBM vs the actual temperature (actually the prediction is probably closer then it should be expected given the spherical cow nature of our approximation)

params = {
    'alpha' : 0.3,
    'A' : 221.2  ,
    'B' : 1.3,
    'r': 5
    'Current_CO2' : ...
    }
Answer (start)
### Answer HERE
Answer (end)

Now, in this model, let’s compute how much more carbon dioxide we can add to the atmosphere before we can reach two degrees of temperature increase from the 1850 baseline. If we assume that you continue to increase the CO2 in the atmosphere at the same rate that we have been in the last 20 years, what year do we reach a two degree warming.

Answer (start)
### Answer Here
Answer (end)

Exercise 2: Snowball Earth#

a. Climate Dynamics#

So far we’ve mainly figured out what the equilibrium temperature of the Earth has to be - i.e. the Sun shines on the Earth for a long period of time and eventually it heats the Earth up so that the radiation the Earth is radiating back to space is equivalent to the radiation that is being absorbed.

In this exercise, we are going to start by dealing with the word eventually. We’d like to understand how long it takes the temperature to approach to equilibrium temperature - this, of course, depends on what we mean by heating up the Earth to that temperature.

It shouldn’t surprise us that we to answer this question, we should start by writing down a differential equation that tells us \(\partial T / \partial t\) where \(T\) is the temperature and \(t\) is the time.

\[ C \frac{\partial T}{\partial t} = \frac{S}{4}(1-\alpha) - (A+B T) \]

(remember in \(BT\) we need that \(T\) is in Celsius).

The r.h.s. of this is just the radiation into the Earth and the radiation out of Earth that we were working on in exercise 1. Previously, we were just setting this r.h.s. to zero - i.e. we wanted the temperature at which the change in temperature over time was zero.

But when it’s not zero, it’s telling us how much the temperature changes in a unit of time. If there is more radiation coming in then going out, we expect that the temperature of the planet should heat up.

The \(C\) is the heat capacity of the planet - it’s telling us if we put a certain amount of radiation in per second how long does that take to actually heat up the Earth. If the Earth were a big chunk of iron it would heat up pretty quickly (think a pot). On the other hand, if it was a big chunk of granite, it would heat up quite slowly. It turns out that Earth (at least the relevant part that’s heating up) is essentially a big bag of water - i.e. a big ocean. It’s a little subtle to figure out which part of the ocean is heating up but a reasonable answer here turns out to be the top 10 meters. In such a case, the heat capacity is going to be \(C=4.18130000e+07\).

Go ahead and set up a first order integrator for this system. Use a time step \(dt\) of 1 month (written in seconds!). You can use const.seconds_per_month to get the correct number of seconds.

Start the temperature of the Earth at 295 degrees Kelvin and allow that temperature to evolve in time. Plot the temperature as a function of time (in years).

We want to observe three things:

  • What is the equilibrium temperature?

  • What is the functional form for approaching the equilibrium temperature?

  • How long does it take to get there?

Try a similar thing starting from a temperature which is less 270 degrees Kelvin.

Answer (start)
### Answer Here
Answer (end)

You should have observed the following things:

  • The planet eventually equilibrates to the temperature that you predicted in exercise 1. This is the fixed point of that equation (i.e. if you start there you stay there).

  • The approach to that equilibrium is exponential. You actually could have deduced this even analytically by looking at the form of the differential equation. This tells you that if you are out of equilibrium that you end up getting pushed back towards the equilibrium temperature relatively quickly (both from above and below). The equilibrium temperature is a stable fixed point - if you push away from it, you flow back to it.

  • The time scale to fix a 5 degree difference from equilibrium is on the order of a decade.

It is interesting to ask how long it would take to relax back to the equilibrium temperature if you start 100 degrees Celsius away.

b. Snowball Earth#

We have seen that the Earth is very stable and largely independent of the initial temperature. But we are missing some important qualitiative physics. In particular, once water gets cold enough it freezes into ice. This ice is actually much more reflective of the Sun and so the albedo is therefore much higher. This tells us that we should have the albedo (and hence the radiation entering the planet) depend on temperature.

  • If the temperature \(T<-10\) degrees celsius, then we should have the albedo be \(\alpha_i = 0.5\)

  • If the temperature \(T>10\) degrees celsius, then we should have the albedo be \(\alpha_0 = 0.3\)

  • In the intermediate regime \(-10 < T < 10\) we can adjust it continuously as \(\alpha_i + (\alpha_0-\alpha_i)(T+10)/20\)

Write a function GetAlpha(T) which returns the temperature dependent albedo and plot the albedo for temperatures from \(-20 < T < 20\) degrees Celsius

Answer (start)
### Answer Here
Answer (end)

Now we want to incorporate our new temperature dependent albedo in our code.

What we want to do is loop over a large number of initial temepratures from \(T=-50\) to \(T=100\) degrees Celsius and then examine what the temperature is as a function of time.

  • Plot temperature as a function of time for mamy different initial temperatures (all plotted as a different line).

  • Secondly, for each initial temperature, figure out what the final fixed point temperature is. Plot the fixed point temperature on the y-axis as a function of the initial temperature.

Answer (start)
### Answer Here
Answer (end)

What you should find is two stable fix points (and one unstable fixed point). The final temperature that the Earth evolves too depends on the initial temperature. This is coming from negative feedack due to the ice. Once you get cold enough, you end up with enough ice that causes it to even get colder. And then once you evolve to a sufficiently cold value it is quite hard to recover.

There is some geological evidence that either all (or a very large fraction) of the Earth was covered by glacial ice. This is called the snowball Earth.

This model gives a mechanism for how such a thing might happen. Consider a scenario where you have a temperature fluctuation for various reasons for a decade; maybe there is anomolous volcanic activity which puts dust in the atmosphere temporarily raising the albedo or maybe there is a change in the solar insolation because of some orbital oscillation. Within a decade you might be in a situation which is very hard to recover from.

c. Hysteresis#

This phenomena is essentially one of Hysteresis. Let’s see this explicitly. Let’s start with a warm planet (10 degrees Celsisus) and the current solar insolation. Now let’s do the following:

  • for step in range(0,10):

    • Wait 5 years

    • Plot the final temperature vs the current solar insolation

    • Decrease the value of Q=S/4 by 5 W/m^2 (i.e S-=20)

Then do the same thing but in reverse increasing the solar insolation by 5 W/m^2 at each step. This time do it for 30 steps.

Both of these line should be on the same plot. Notice that there is a hystersis. Even after the solar insolation has gotten back to where it was, you are still not back to your original temperature. In fact, there’s essentially no way to get back to your original temperature this way because you need to way overshoot before you melt the planet.

Answer (start)
### Answer Here
Answer (end)

Exercise 3: One-dimensional EBM#

We want to implement a more accurate climate model. Previously, we had worked on a zero-dimensional model. To add some realism to this model, we are going to add an additional dimension splitting up the world into different latititudes. This is the most important effect that we are missing as these different latitudes get significantly different sunlight.

The energy balance model that we will work with will be

\[ C \frac{dT_s(\phi)}{dt} = (1-\alpha(\phi;T)) Q(\phi) - (A+B T_s(\phi)) + \frac{D}{\cos⁡\phi } \frac{\partial }{\partial \phi} \left( \cos⁡\phi ~ \frac{\partial T}{\partial \phi} \right)\]

a. Latitude dependence of solar isolation and albedo#

The first thing you should notice is that there are now a lot of different terms that depend on the latitude.

To do that you need to discretize the Earth into a number of different latitude zones. You could just discretize your own latitudes but because we need to sync our functions with yours, we will give you the latitude zones that you are going to use. You can get them by calling the following function:

def GetLatitudes():
  sfc = climlab.domain.zonal_mean_surface(num_lat=180, water_depth=10.)
  lat = sfc.lat.points
  return lat

The first thing that depends on latitude which is important is the amount of radiation from the Sun is being deposited on different parts of the Earth. We know that the equator gets significantly more sunlight then the poles. In principle, we should adjust how much each pole gets depending on the season (summer vs winter). This is an interesting extension you can do but we will start by just assigning an average amount of Sun that touches the different parts of the Earth.

You can get this by calling the following function

def GetInsolation(param):
  model1 = climlab.EBM_annual(name='EBM with interactive ice line',
                            num_lat=180,timestep=const.seconds_per_month,
                            **param)
  Q= np.array([i[0] for i in model1.insolation])
  return Q

which will return the amount of watts per square meter of each latitude zone on the Earth.

We will actually need to send this a set of parameters (which tells us things like the values of differnet coefficients. Here are the parameters we will use (we will explain them in due time))

params = {
    'D':0.6,
    'A':210,
    'B':2,
    'a0':0.354,
    'a2':0.25,
    'ai':0.5,
    'Tf':-10.,
    'CO2':283.,
    'C':4.18130000e+07,
    'r':5,
    'timeStep':const.seconds_per_month
    }

  • Go ahead and plot this as as a function of latitude

  • Measure the average radiation on the Earth and verify that it is approximately 1370/4 (the number we were using before). You need to be careful when you do this. When averaging over a sphere, you need to weight by the volume at that latitude (i.e. \(\cos \phi\)). It will behoove you to write a function which takes a array over the latitudes and averages them correctly (we will need this for the average temperature of the planet as well).

Answer (start)
### Answer Here
Answer (end)

The next thing that we need to deal with is the latitude dependence of the albedo. As we saw in the 0D model, the albedo will change its reflectivity when things freeze. We are going to do something a little simpler and just have the albedo be \(\alpha_i = 0.5\) (see parameters above) when the temperature \(t_f<-10\) Celsius.

But there is also an effect of the albedo that is coming from the latitude. In particular, because the Sun is shining at more of an angle closer to the poles.

If the temperature is greater then \(-10 < t_f\) then we should have an Albedo of \(\alpha_0 + a_2 P_2(\sin \phi)\) where \(P_2\) is the second legendre polynomial and can be gotten by doing legendre.P2(x).

We’ve given you some initial conditions which you can get by doing the following

def GetInitial(lat):
  initial = np.array(12. - 40. * legendre.P2(np.sin(np.deg2rad(lat))))+273.15
  return initial
initialTemperature=GetInitial(lat)

Go ahead and write a function which generates the Albedo when given an array of temperatures (and the parameters). Then plot the initial temperature and (separately) what the Albedo looks like for the initialTemperature. Also compute the average of the initial temperature. Your albedo should have some values at 0.5.

Answer (start)
### Answer Here
Answer (end)

b. A diffusion-less planet#

We are now in a position to deal with the incoming and outgoing radiation (but not yet the diffusion). This is implementing the same energy balance that you had previously except that now you are going to simultaneously deal with all the latititudes. Write a function to time evolve this for 100 years.

  • Plot the average temperature vs time (getting a point once a year)

  • Plot the final temperature profile

  • The final temperature difference between the equator and the pole (i.e. just the max - min difference)

You should be running with time steps of 1 month.

Answer (start)
### Answer Here
Answer (end)

You should get a pretty crazy result in this case. Essentially every latitude is doing its own thing (there are no independent because there is no talking between latitudes). Thereroe, you’re going to see some pretty sharp jumps in the temperature as you change latitudes. In addition, there is a large (~117 degrees) difference between the poles and the equator.

c. Diffusion#

We are missing an important piece of the physics once we’ve moved to 1D. The temperature at different latitudes can interact with each other. Temperature can move both by diffusion and advection. Here we are going to concern ourself primarily with diffusion. We’ve seen diffusion before in this class but typically as \(\partial^2/ \partial x^2\). In our case though we are on a sphere and need to do something different essentially writing our gradient in spherical coordinates - see the last term of our energy balance model.

We’d like to go ahead and implement the diffusion just like we’ve done in previous cases using a finite difference stencil. Unfortunately, this won’t work here (we run into numerical instabilities at the poles). What we need to do is use an implicit solver of our differential equations. This is out of scope for this assignment, so instead we will give you a function you can call

def Diffuse(T,params):
  sfc = climlab.domain.zonal_mean_surface(num_lat=180, water_depth=10.)
  Ts = climlab.Field(np.array(T), domain=sfc)
  d = climlab.dynamics.MeridionalHeatDiffusion(name='Diffusion',state=Ts, D=params['D'],timestep=params['timeStep'])
  d.step_forward()
  T=np.squeeze(Ts)
  return T

by going

T=Diffuse(T,params)

This will defuse through a time-step of \(\delta T\). This said, there is an extra 30 points extra credit if you can right your own diffuse function which does the diffusion correctly without using libraries. You will almost certainly need implicit solvers to do this.

Go ahead and put everything together now. You should be calling the diffusion every time-step as well.

  • Plot the average temperature vs time (getting a point once a year)

  • Plot the final temperature profile

  • The final temperature difference between the equator and the pole (i.e. just the max - min difference)

  • The ice line: the latitude at which the temperature becomes less then -10.

The diffusion constant \(D\) has been tuned to give a reasonable difference between the poles and the Earth.

Answer (start)
### Answer Here
Answer (end)

You don’t have to do this but you can also experiment with various snowball Earth scenarios with this model.