Exercise Sheet

Exercise sheet for Lecture Climate Dynamics

Prof. Dr. Gerrit Lohmann

This course is part of the introductory course on mathematical modeling for climate and environment at the University of Bayreuth, held and organized by Prof. Dr. Vadym Aizinger.

First some basics on numerical schemes:

Exercise I: Euler integrator

Play with the parameter h which is the step size in time units !
Can you find a crtical h when the numerical scheme does not work anymore?
Which value of h shall be taken in order to have an accurate solution (close to the analytical one) ?

#demonstration of Euler method in 1st order ODE: dy/dt=A*y

#the function dy/dt<-f(y,A,t)
f<-function(y,A,t) 
{
    return(A*y)
}

#constants
A<-  -0.02  #growth / decay rate
T<-  500    #integration time in time units
h<-  10     #step size in time units
Y0<- 8   #inital value

n<-T/h          #number of time steps (time / timestep)
t<-(0:(n-1))*h  #create a vector of discrete timesteps 
y<-vector()     #define an empty vector for the state variable y(t)
y[1]<-Y0        #assign initial value


#integration loop
for (i in 1:(n-1))
{
    y[i+1]<-y[i]+h*f(y[i],A,t[i])   #Euler forward: y[t+h]<-y[t]+h*A*y[t]
}


plot(t,y,type="p") #plot the result against time

#additionaly plot the analytical solution in red
lines(t,Y0*exp(A*t),col="red",type="l")

Exercise II: Radiocarbon as dating tool: Exponential processes of decay

We look at a substance (anorganic, organic or living), measured by its mass or number of individuals, to which nothing is added and which decays or dies in the process of time, as in the following examples:

Trees absorb carbon dioxide from the air and store the contained carbon in their wood. A certain part of this carbon is radioactive (C-14). When a tree is felled, no more external carbon enters the wood, and the contained C-14 decays. How much of it is left after 10, 100, 1000 years? Knowing the answer, we can reciprocally conclude the wood’s age from the measured C-14 content (C-14 dating method). A simple calculation:

\[ ^{14}C(t) = ^{14}C(t_0) \cdot exp[-k (t-t_0)] \] \[ \ln [^{14}C(t)/^{14}C(t_0)] = -k (t-t_0) \\ D t = (t - t_0) = -1/k \, \, \ln [^{14}C(t)/^{14}C(t_0)] = T_h / \ln2 \cdot {\ln [^{14}C(t_0)/^{14}C(t)]}\] where $k:$ radioactive decay constant of $^{14}C$,

1/k = mean life time = $T_{h}/ \ln2$

and half life time $T_{h}$.

$^{14}C$ is continually produced in the atomosphere and it is incorporated into all living creatures (via CO$_2$ and food). The amount of $^{14}C$ in a living creature is constant. When the organism dies, the amount of $^{14}C$ starts to decrease, because no new $^{14}C$ from the atmosphere is being incorporated. By counting the remaining radioactivity, it is possible to tell how old an item is.

Question:

A piece of fossilized wood has a carbon 14 radioactivity that is 1/4 that of new wood. The halflife of carbon14 is $ T_{1/2} = 5730 $ years. How old is the fossilized wood?

1 x 5730 = 5730 years
2 x 5730 = 11,460 years
3 x 5730 = 17,190 years
0.25 x 5730 = 1432 years

Calculate the number of C-14 atoms after 10,000 years with initial #10,000 !

The homogeneous Poisson point process can be defined as a counting process, a type of stochastic process, which can be denoted as $\{N(t),t\geq 0\}.$ A counting process represents the total number of events that have happened up to and including time t.

A counting process is a homogeneous Poisson counting process with rate $\lambda >0$ if it has the following three properties:

N(0)=0; has independent increments; and the number of events in any interval of length t is a Poisson random variable with parameter (or mean) $\lambda t.$

The last property implies:

\[ \mathbb {E} [N(t)]=\lambda t. \] In other words, the probability of the random variable N(t) being equal to n is given by:

\[ \mathbb {P} \{N(t)=n\}={\frac {(\lambda t)^{n}}{n!}}e^{-\lambda t}. \] The Poisson counting process can also be defined by stating that the time differences between events of the counting process are exponential variables with mean $1/\lambda.$ The time differences between the events or arrivals are known as interarrival or interoccurence times.

Some applications:

Exercise III: Interactive climate model including sea ice

What will happen if the CO2 content in the atmosphere is doubled? Radiative forcing= 4 Wm-2
What will happen if the factor γ is 1% higher/lower in the long-wave radiation?
Describe the effect if the diffusivity is enhanced by a factor of 2!
The coalbedo of sea ice can vary between 0.3 and 0.4. Describe the effect when varying the value!
Write down the numerical scheme (time stepping etc. from the source code) !
Show the evolution at one specific latitude and discuss it!

Exercise IV: Interactive Daisy World model

Run Daisy World with the the Daisies - black and white. Notice for what range of luminosity the daisies manage to control the planet temperature.
Do you think control would be better if you set the low and high albedos (for black and white) to a wider spread?
You will notice that the living area (“total daisies”) doesn’t exceed 70%. The deathrate is set to 0.3, which may explain the living percentage being no more than 0.7. What does the deathrate do to the daisies’ ability to control their environment’s temperature? To the species mix?
What is the effect with no daisies for the global temperature ? This can be studied by increasing the death rate to 1.
What is the influence of the optimal growth temperature? Please vary the numbers and describe the consequence !
Daisyworld Model quizz

PS: here is the Python code

Exercise V: Random Systems

Simulate the velocity evolution of one particle which is determined by the following stochastic dx/dt = -bx + kW(t)
What happens if you change the timestep ?
Simulate the ensemble of multiple particles and plot the time evolution of the v-distribution and compare with the diffusion equation !
Test the ergodic theorem: time average is equal ensemble average !
Have a look at the more complex 2D diffusion equation program. When are the two colors are mixed?

Linear stochastic equation: Hasselmann climate model

Multiple particles of the linear stochastic equation

Simulation of 2D diffusion

# simulation of 2D diffusion   

#two simplifications: deltaT is hardcoded to one
#the reflection on the boundaries is simplified 
#if a particle gets out of the boundaries, the direction gets reversed



NX<-100    #horizontal and vertical resolution
NY<-100

vmax<-0.0 #maximal velocity of initial velocity distribution
beta<-0.1    #strength of friction term
K<-0.2    #strength of random term

N<-1000 #number of particles

v.particles<-matrix(1,2,N)       #matrix 2xN, v of particles[1,] are the X, particles[2,] the Y coordinates
particles<-matrix(1,3,N)         #matrix 2xN, particles[1,] are the X, particles[2,] the Y coordinates, particles[3,] the color

v.particles[1,]<-runif(NX,-vmax,vmax)   #initalize the system with uniform distributed velocities
v.particles[2,]<-runif(NX,-vmax,vmax)

particles[1,1:500]<-30  #half of the particles start at point 30/30 
particles[2,1:500]<-30
particles[3,1:500]<-1     #and has the color 1

particles[1,500:1000]<-70  #the other half at point 70/70 and has the color 100
particles[2,500:1000]<-70
particles[3,500:1000]<-100

T<-200                 #number of timesteps to run

for (i in 1:T)
{
    v.particles[1,]<-v.particles[1,]-v.particles[1,]*beta+K*rnorm(N)   #DGL for velocity
    v.particles[2,]<-v.particles[2,]-v.particles[2,]*beta+K*rnorm(N)
    
    particles[1,]<-particles[1,]+v.particles[1,]  #brownian motion for each particle
    particles[2,]<-particles[2,]+v.particles[2,]

    #rough implementation of the reflection and boundaries...
    #is not exact and therefore loses some particles
    
    index.collisionX<-((particles[1,]>NX) | (particles[1,]<1))
    v.particles[1,index.collisionX]<-(-1)*v.particles[1,index.collisionX]

    index.collisionY<-((particles[2,]>NY) | (particles[2,]<1))
    v.particles[2,index.collisionY]<-(-1)*v.particles[2,index.collisionY]

    plot(particles[1,],xlim=c(1,NX),ylim=c(1,NY),particles[2,],col=particles[3,],type="p")

}

Exercise VI: Vegetation Dynamics and the African humid period

Something for home to think about it:

Describe the vegetation dynamics for the Holocene (8000 years ago)! Which evidences?
Which feedbacks are acting in the system ?
Study Feynman’ study
Is it possible to generate a green Sahara and Sinai under present conditions?
Is this an option solving political conflicts ?

Exercise VII: Leapfrog scheme

studying Rossby and Kelvin waves in the shallow water equations

Study the Leapfrog scheme in the following R code !
Write down the scheme for an arbritray 1d system!
Can you find a crtical h when the numerical scheme does not work anymore?

#This is just a definition of a function to plot vectorplots
 {
    u <- par("usr") 
    p <- par("pin") 
    c(p[1]/(u[2] - u[1]), p[2]/(u[4] - u[3]))
 }

quiver<-function(lon,lat,u,v,scale=1,length=0.05,maxspeed=200, ...) 
 # from http://tolstoy.newcastle.edu.au/R/help/01c/2711.html

  {
    ypos <- lat[col(u)] 
    xpos <- lon[row(u)]

    speed <- sqrt(u*u+v*v) 

    u <- u*scale/maxspeed 
    v <- v*scale/maxspeed 
    

    matplot(xpos,ypos,type="p",cex=0,xlab="lon",ylab="lat", ...) 
    arrows(xpos,ypos,xpos+u,ypos+v,length=length*min(par.uin())) 
  }



#Program starts here

#Shallow water 2D,cyclic boundary conditions + Coriolis term
nn<- 50
ni<- 2*nn+1   #number of gridcells in one direction
nt<-10000 #number of timesteps



#The physical constants
g<-0.1 #low gravity, 0.1 m/s^2
dx<-1e5 #gridcell 10km 
dy<-1e5
dt<-500 #timstep 1000 second
H<-1e3 #1km depth
Omega<-1e-4


#define three index vectors.. the middle one , one shifted one cell to the left, and one to the right 
#(including the periodic boundary conditions)
ia.0<-1:ni
ia.m1<-c(ni,1:(ni-1))
ia.p1<-c(2:ni,1)
u<-matrix(0,ni,ni) #speed at each point
v<-matrix(0,ni,ni) #speed at each point
h<-matrix(0,ni,ni) #pertubation at each point
f<-matrix(0,ni,ni) #pertubation at each point

lat<-c(-nn:nn)*90/nn
weight<-sin(lat*pi/180)
lon<-c(-nn:nn)*180/nn

f<-rep(weight*2*Omega,each=ni)
dim(f)<-c(ni,ni)
#filled.contour(f)

u.new<-u
h.new<-h
v.new<-v

#Inital condition: One smooth blobs at each side of the "equator"(sin)
idit=nn/5*2
inix=ni-idit-1
iniy=ni-2*idit-1
endx=ni-1
endy=ni-1
endy2=2*idit+1
h[inix:endx,iniy:endy]<-sin(0:20/2*pi/10)%*%t(sin(0:40/2*pi/20))
h[inix:endx,1:endy2]<-sin(0:20/2*pi/10)%*%t(sin(0:40/2*pi/20))

#equator to study the Kelvin wave:
ii=idit+1
iy=nn-10
iy2=nn+10
h[1:ii,iy:iy2]<- -sin(0:20/2*pi/10)%*%t(sin(0:20/2*pi/10))


#1st step euler forward
u.new[ia.0,ia.0]<-u[ia.0,ia.0]-g*dt/2/dx*(h[ia.p1,ia.0]-h[ia.m1,ia.0]) 
v.new[ia.0,ia.0]<-v[ia.0,ia.0]-g*dt/2/dy*(h[ia.0,ia.p1]-h[ia.0,ia.m1]) 
h.new[ia.0,ia.0]<-h[ia.0,ia.0]-H*dt/2*((u[ia.p1,ia.0]-u[ia.m1,ia.0])/dx + (v[ia.0,ia.p1]-v[ia.0,ia.m1])/dy)  #(du/dx + dv/dy)

#Divide the screen in two parts
#par(mfcol=c(2,1))

#Leapfrog from the third step on
for (n in 3:(nt-1))
{
    
    u.old<-u
    v.old<-v
    h.old<-h
    h<-h.new
    u<-u.new
    v<-v.new

    u.new[ia.0,ia.0]<-u.old[ia.0,ia.0]-g*dt/dx*(h[ia.p1,ia.0]-h[ia.m1,ia.0])+dt*f*v
    v.new[ia.0,ia.0]<-v.old[ia.0,ia.0]-g*dt/dy*(h[ia.0,ia.p1]-h[ia.0,ia.m1])-dt*f*u
    h.new[ia.0,ia.0]<-h.old[ia.0,ia.0]-H*dt*((u[ia.p1,ia.0]-u[ia.m1,ia.0])/dx + (v[ia.0,ia.p1]-v[ia.0,ia.m1])/dy)  #(du/dx + dv/dy)
    
    
    #plot every 50th image
    if ((n %% 50) == 0) {
        #quiver(lon,lat,u,v,scale=200,maxspeed=1.5,length=3)
        image(lon,lat,h,zlim=c(-1,1),col=rainbow(200)) #this time as color coated
        #persp(h/3, theta = 0, phi = 40, scale = FALSE, ltheta = -120, shade = 0.6, border = NA, box = FALSE,zlim=c(-0.3,0.3))
        
    }

}

well done