Climate Dynamics

Climate Dynamics#

Show code cell content Hide code cell content

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.10.17/x64/lib/python3.10/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.10.17/x64/lib/python3.10/site-packages/climlab/model/__init__.py:28
      1 '''
      2 This package contains ready-made models that can be run "off-the-shelf".
      3 
   (...)
     26 radiative-convective column model from individual components.
     27 '''
---> 28 from .column import GreyRadiationModel, RadiativeConvectiveModel, BandRCModel
     29 from .ebm import EBM, EBM_annual, EBM_seasonal

File /opt/hostedtoolcache/Python/3.10.17/x64/lib/python3.10/site-packages/climlab/model/column.py:39
     37 import numpy as np
     38 from climlab import constants as const
---> 39 from climlab.process import TimeDependentProcess
     40 from climlab.domain import column_state, Field
     41 from climlab.radiation import (FixedInsolation, GreyGas, GreyGasSW,
     42         ThreeBandSW, FourBandLW, ManabeWaterVapor)

File /opt/hostedtoolcache/Python/3.10.17/x64/lib/python3.10/site-packages/climlab/process/__init__.py:2
      1 '''The base classes for all climlab processes.'''
----> 2 from .process import Process, process_like, get_axes
      3 from .time_dependent_process import TimeDependentProcess, couple
      4 from .implicit import ImplicitProcess

File /opt/hostedtoolcache/Python/3.10.17/x64/lib/python3.10/site-packages/climlab/process/process.py:69
     67 import time, copy
     68 import numpy as np
---> 69 from climlab.domain.field import Field
     70 from climlab.domain.domain import _Domain, zonal_mean_surface
     71 from climlab.utils import walk

File /opt/hostedtoolcache/Python/3.10.17/x64/lib/python3.10/site-packages/climlab/domain/__init__.py:7
      4 __all__ = ['axis', 'domain', 'field', 'initial', 'xarray']
      6 from climlab.domain.domain import single_column, zonal_mean_surface, surface_2D, zonal_mean_column, box_model_domain
----> 7 from climlab.domain.initial import column_state, surface_state
      8 from climlab.domain.field import Field, global_mean
      9 from climlab.domain.axis import Axis

File /opt/hostedtoolcache/Python/3.10.17/x64/lib/python3.10/site-packages/climlab/domain/initial.py:4
      2 import numpy as np
      3 from climlab.domain import domain
----> 4 from climlab.domain.field import Field
      5 from climlab.utils.attrdict import AttrDict
      6 from climlab.utils import legendre

File /opt/hostedtoolcache/Python/3.10.17/x64/lib/python3.10/site-packages/climlab/domain/field.py:9
      1 #  Trying a new data model for state variables and domains:
      2 #  Create a new sub-class of numpy.ndarray
      3 #  that has as an attribute the domain itself
   (...)
      6 #
      7 # http://docs.scipy.org/doc/numpy/user/basics.subclassing.html
      8 import numpy as np
----> 9 from climlab.domain.xarray import Field_to_xarray
     12 class Field(np.ndarray):
     13     """Custom class for climlab gridded quantities, called Field.
     14 
     15     This class behaves exactly like :py:class:`numpy.ndarray`
   (...)
     83 
     84     """

File /opt/hostedtoolcache/Python/3.10.17/x64/lib/python3.10/site-packages/climlab/domain/xarray.py:3
      1 from builtins import str
      2 from builtins import object
----> 3 from xarray import Dataset, DataArray
      4 import warnings
      7 def Field_to_xarray(field):

ModuleNotFoundError: No module named 'xarray'

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 here

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 write the function OutgoingRadiationPerSquareMeter(T) which computes the outgoing radiation given the temperature $T$.

### Answer here

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

Do this in two ways.

🦉 Graph the two functions you’ve written as a function of temperature and determine at temperature $T$ the two match.

🦉 Define a cost function diff(T) which gives the absolute value of the difference between the incoming and outgoing radiation. 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 here

Q: Was your predicted temperature too high or low compared to the observed temperature? Which did you expect and why?

A:

b. The temperature of Earth with air and clouds#

We have ignored two important effects in the previous analysis.

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 blocked by the greenhouse gases in the atmosphere. We call the fraction that escapes $\epsilon$, which we do not have a good way to measure directly. We can represent this by including an emissivity term for the atmosphere so that the radiation emitted is $\epsilon \sigma T^4$.

The reason that the outgoing radiation is blocked separately from the incoming radiation is that incoming radiation is near-infrared and visible (it’s the sunlight!), while the outgoing radiation is further in the infrared. This difference in frequencies causes different absorption/reflection by the atomosphere.

🦉 Plot the earth’s temperature versus $\epsilon$, including the effect of Earth’s albedo. Note that $\epsilon$=1 means no greenhouse effect, and $\epsilon=0$ means a complete greenhouse effect (no light ). You should not include $\epsilon=0$.

🦉 Find the value of $\epsilon$ that corresponds to the pre-industrial average temperature of 14 degrees Celsius, in a similar way as you found the temperature above.

### Answer Here

c. A more accurate emission model#

To further improve our model, we need to increase the accuracy of our emission model. The Stefan-Boltzmann law from above applies to a perfect blackbody, which the Earth is not. From experiments, we can measure that near the Earth’s temperature, $$ P_{out} = A + B T_{c}, $$ where

$A=221.2$ W/m²
$B=1.3$ W/m² K
$T_c$ is the temperature in Celcius.

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

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

If you were to linearize the Stefan-Boltzmann law around 0 Celcius, you would get $$\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.

Plugging in our numbers from previous parts of the exercis, we get $A=195.94$ W/m² and $B=2.87$ W/m² K. This is close to the measured values, which is reassuring.

If you used these naiive values, you’d find that we get a temperature $T=15.11878309$.

### Answer Here

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. The effects of changing carbon dioxide#

Changing carbon dioxide changes $\epsilon$ (and therefore $A$ and $B$); however, typically this is measured in terms of ‘forcing,’ which is how much more energy from some reference that is reflected back to the Earth’s surface. This is accounted for by adding an additional term that depends on a reference amount of carbon dioxide, resulting in a modified output power of

\[ P_{out} = A + BT_c + r \log \left(\frac{CO_{2;current}}{CO_{2;old}}\right) \]

where $r=-5$ W/m² is a proportionality constant, and CO₂ is the concentration of CO₂ in the atmosphere. You can see that if the CO₂ concentration increases, the power out decreases.

We will use the oldest pre-industrial data (around 1751) as the reference CO₂ concentration.

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

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 year is in tempData[:,0]). Notice the close agreement between the predicted temperature based on the CO₂ 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 HERE

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 since 2000 (do a linear fit), in what year do we reach a two degree warming, measured from the pre-industrial average of 14°C?

### Answer Here

Exercise 2: Snowball Earth#

a. Climate Dynamics#

So far we’ve focused on the equilibrium temperature of the Earth. It turns out that dynamics can be an important part of the Earth’s temperature, since there can be feedback loops. For this part of the exercise, we are going to assume that nothing is pumping tons of CO₂ into the air.

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 to find equilibrium. $C$ is the heat capacity per meter squared of the planet – heat capacity is how much energy it takes to change the temperature by 1 Kelvin. To a good approximation, the Earth is just water, so we use that heat capacity. The question is how deep do we go, and the answer is about 10 meters. From the specific heat of water, we can do some unit analysis and get that the heat capacity is $C=4.1813 \times 10^7 \frac{\textrm{J}}{\textrm{K}\cdot \textrm{m}^2}$.

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?

Do the same thing starting from a temperature which is less 270 degrees Kelvin.

Do the same if you start 100 degrees Celsius too hot.

### Answer Here

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.

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 than water or land 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$° Celsius

### Answer Here

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

Loop over at least 10 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 Here

What you should find is two stable fixed points separated by an unstable fixed point (which you might not have seen). The final temperature that the Earth evolves to depends on the initial temperature. This history-dependent effect is a result of a feedback loop–lower temperatures means more ice, and more ice means less absorbed radiation from the Sun.

There is some geological evidence that this happened in the past; either all (or a very large fraction) of the Earth was covered by glacial ice – snowball Earth. It is possible to switch between these fixed points through large scale events, such as volcanic activity, meteor strikes, solar insolation due to changes in the Sun’s output or orbital oscillation. It’s particularly striking how quickly very similar initial conditions result in very different output.

c. Hysteresis#

Because the Earth takes some time (of order 5 years) to respond to changes, it exhibits hysteresis, which refers to a lag in physical properties compared to the stimulus. 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:

Start S = 1367.6 J/m²

for step in range(0,10):
- Propagate the Earth’s temperature forward by 5 years.
- Decrease S by 20 W/m²
for step in range(0,40):
- Propagate the Earth’s temperature forward by 5 years.
- Increase S by 5 W/m²
Now keeping S constant, propagate the temperature forward by 50 more years.

Plot the temperature versus year and irradiation versus year over the 300 year span of your simulation. Notice that there is a hystersis. Even after the solar insolation has gotten back to where it was originally, the temperature has not recovered. The only way to get back to the original high temperature is to melt the ice somehow.

### Answer Here

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)\]

where $\phi$ is the latitude (the equator is 0 degrees and the south and north poles are at -90 and +90 respectively)

Be careful here because climlab works with degrees and most numpy functions work with radians.

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, mainly the solar iradiation, temperature, and there is a diffusion part that transports heat between latitudes.

To treat this, 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 important thing that depends on latitude is the amount of radiation from the Sun is being deposited on different parts of the Earth. The equator gets significantly more sunlight then the poles.

The package climlab has a function that computes this conveniently. You can get the insolation per 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)
  model1.compute()
  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 use the following parameters to model Earth

params = {
    'D':0.6, #Diffusion constant D
    'A':210, # Emission parameters
    'B':2,
    'a0':0.354, #base albedo 
    'a2':0.25,  # How albedo changes with latitude
    'ai':0.5,  # albedo of ice
    'Tf':-10., # Temperature at which ice is completely formed.
    'CO2':283., #amount of CO2 in the atmosphere
    'C':4.18130000e+07, #Heat capacity per m^2 in the 
    'r':5, # forcing constant for CO2
    'timeStep':const.seconds_per_month
    }

Plot the insolation 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 by $\cos \phi$. It will be useful 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 Here

We will treat the latitude dependence of the albedo next. Let’s use the following simple model:

If $T<T_f$, the albedo is $\alpha_i = 0.5$.
Otherwise, set the albedo to $\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 will also start with initial conditions as follows.

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

Plot the initial temperature. Compute the average of the initial temperature and verify that it is reasonable for the Earth.

Write a function which generates the Albedo when given an array of temperatures (and the parameters). For the initial temperature, check that your albedo has some values that are exactly 0.5 which correspond to ice-covered regions. Check that the ice-covered regions occur where you expect.

### Answer Here

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 Here

You should get a pretty crazy result in this case. Essentially every latitude is doing its own thing; they are independent because the diffusion transports heat between the latitudes and we have ignored it so far. Therefore, you’re going to see some pretty sharp jumps in the temperature as you change latitudes. In addition, there is a large (~117° C) difference in temperature 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. Heat 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 params['timeStep'].

Add the diffusion to your timestepping from part b as an additional step. Then

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 than -10.

The diffusion constant $D$ has been tuned to give a reasonable temperature difference between the poles and the Earth.

### Answer Here