Orbit & astronomical theory of ice ages
Gerrit Lohmann
Exercise “Earth orbital variations”
Insolation calculation
Paleoclimatic reconstructions help to discover the natural variability of the climate system over times scales ranging from years to hundreds of thousands of years. They are fundamental in climate research, especially now, because they provide a unique set of data to validate models over climatic situations largely different from those of the last 150 years. The climatic situations of the last century are indeed available in great detail, but with a very poor diversity. Among the different modes of climatic variations, the glacial‐interglacial cycles have the advantage that they provide examples of extreme climates and that their primary astronomical cause is now pretty well known. The Astronomical Theory of paleoclimates aims indeed to explain these climatic variations occurring with quasi‐periodicities situated between tens and hundreds of thousands of years. Such variations are recorded in deep‐sea sediments, in ice sheets and in continental archives. The origin of these quasi‐cycles lies in the astronomically driven changes of the latitudinal and seasonal distribution of the energy that the Earth receives from the Sun. Milutin Milankovitch extensively published about this theory between 1920 and 1950, but the relationship between the astronomical parameters, insolation and climate, had already been suggested at the beginning of the nineteenth century. There were several authors who have contributed to what is now named (often improperly) the Milankovitch astronomical theory.
Other references:
Berger A. and Loutre M.F. (1991). Insolation values for the climate of the last 10 million years. Quaternary Science Reviews, 10(4), 297-317.
Berger A. (1978). Long-term variations of daily insolation and Quaternary climatic changes. Journal of Atmospheric Science, 35(12), 2362-2367.
Calculate insolation parameters
### Insolation, converted and adapted from Huybers Code, based on Berger 1991
daily_insolation_param<-function(lat,day,ecc,obliquity,long_perh,day_type=1)
{
# Usage:
# Fsw = daily_insolation(lat,day,eccentricity,obliquity,long_perh)
#
# Optional inputs/outputs:
# [Fsw, ecc, obliquity, long_perh] = daily_insolation(kyear,lat,day,day_type)
#
# Description:
# Computes daily average insolation as a function of day and latitude at
# any point during the past 5 million years.
#
# Inputs:
# kyear: Thousands of years before present (0 to 5000).
# lat: Latitude in degrees (-90 to 90).
# day: Indicator of time of year; calendar day by default.
# day_type: Convention for specifying time of year (+/- 1,2) [optional].
# day_type=1 (default): day input is calendar day (1-365.24), where day 1
# is January first. The calendar is referenced to the vernal equinox
# which always occurs at day 80.
# day_type=2: day input is solar longitude (0-360 degrees). Solar
# longitude is the angle of the Earth's orbit measured from spring
# equinox (21 March). Note that calendar days and solar longitude are
# not linearly related because, by Kepler's Second Law, Earth's
# angular velocity varies according to its distance from the sun.
# Output:
# Fsw = Daily average solar radiation in W/m^2.
# Can also output orbital parameters.
#
# This script contains orbital parameter data for the past 50000 years
# from Berger and Loutre (1991).
#
# Detailed description of calculation:
# Values for eccentricity, obliquity, and longitude of perihelion for the
# past 5 Myr are taken from Berger and Loutre 1991 (data from
# ncdc.noaa.gov). If using calendar days, solar longitude is found using an
# approximate solution to the differential equation representing conservation
# of angular momentum (Kepler's Second Law). Given the orbital parameters
# and solar longitude, daily average insolation is calculated exactly
# following Berger 1978.
#
# References:
# Berger A. and Loutre M.F. (1991). Insolation values for the climate of
# the last 10 million years. Quaternary Science Reviews, 10(4), 297-317.
# Berger A. (1978). Long-term variations of daily insolation and
# Quaternary climatic changes. Journal of Atmospheric Science, 35(12),
# 2362-2367.
#
# Authors:
# Ian Eisenman and Peter Huybers, Harvard University, August 2006
# eisenman@fas.harvard.edu
# This file is available online at
# http://deas.harvard.edu/~eisenman/downloads
# Translated into R by Thomas Laepple
# Suggested citation:
# P. Huybers and I. Eisenman, 2006. Integrated summer insolation
# calculations. NOAA/NCDC Paleoclimatology Program Data
# Contribution #2006-079.
#
# === Get orbital parameters ===
epsilon<-obliquity*pi/180
omega<-long_perh*pi/180
# === Calculate insolation ===
lat=lat*pi/180 # latitude
# lambda (or solar longitude) is the angular distance along Earth's orbit measured from spring equinox (21 March)
if (day_type ==1) # calendar days
{
# estimate lambda from calendar day using an approximation from Berger 1978 section 3
delta_lambda_m <- (day-80)*2*pi/365.2422; ## das w‰re lambda bei gleich langen Tagen
beta<-(1-ecc^2)^(1/2)
lambda_m0<-(-2)*( (1/2*ecc+1/8*ecc^3)*(1+beta)*sin(-omega)-1/4*ecc^2*(1/2+beta)*sin(-2*omega)+1/8*ecc^3*(1/3+beta)*(sin(-3*omega)))
lambda_m<-lambda_m0+delta_lambda_m
lambda<-lambda_m+(2*ecc-1/4*ecc^3)*sin(lambda_m-omega)+(5/4)*ecc^2*sin(2*(lambda_m-omega))+(13/12)*ecc^3*sin(3*(lambda_m-omega))
#Wie kann ich das aufhaengen !? Eventuell einfach ueber sommer solstice verschieben ? d.h. zu 270°
} else if (day_type==2) #solar longitude (1-360)
{lambda=day*2*pi/360 # lambda=0 for spring equinox
} else stop('Error: invalid day_type')
So<-1365 # solar constant (W/m^2)
delta<-asin(sin(epsilon)*sin(lambda)) # declination of the sun
Ho<-acos(-tan(lat)*tan(delta)) # hour angle at sunrise/sunset
# no sunrise or no sunset: Berger 1978 eqn (8),(9)
Ho[(abs(lat) >= (pi/2 - abs(delta)) ) & ( lat*delta > 0 )] <- pi
Ho[( abs(lat) >= (pi/2 - abs(delta)) ) & ( lat*delta <= 0 )] <- 0
# Insolation: Berger 1978 eq (10)
Fsw=So/pi*(1+ecc*cos(lambda-omega))^2 /(1-ecc^2)^2 * ( Ho*sin(lat)*sin(delta) + cos(lat)*cos(delta)*sin(Ho))
return(list(Fsw=Fsw,ecc=ecc,obliquity=obliquity,long_perh=long_perh,lambda=lambda/2/pi*360 ))
}calculate annual insolation
annual_insolation<-function(kyear,lat)
{
result<-vector()
for (i in 1:length(kyear)) result[i]<-mean(daily_insolation(kyear[i],lat,1:365)$Fsw)
return(result)
}calculate daily insolation
daily_insolation<-function(kyear,lat,day,day_type=1,fast=T)
{
# Usage:
# Fsw = daily_insolation(kyear,lat,day)
#
# Optional inputs/outputs:
# [Fsw, ecc, obliquity, long_perh] = daily_insolation(kyear,lat,day,day_type)
#
# Description:
# Computes daily average insolation as a function of day and latitude at
# any point during the past 5 million years.
#
# Inputs:
# kyear: Thousands of years before present (0 to 5000).
# lat: Latitude in degrees (-90 to 90).
# day: Indicator of time of year; calendar day by default.
# day_type: Convention for specifying time of year (+/- 1,2) [optional].
# day_type=1 (default): day input is calendar day (1-365.24), where day 1
# is January first. The calendar is referenced to the vernal equinox
# which always occurs at day 80.
# day_type=2: day input is solar longitude (0-360 degrees). Solar
# longitude is the angle of the Earth's orbit measured from spring
# equinox (21 March). Note that calendar days and solar longitude are
# not linearly related because, by Kepler's Second Law, Earth's
# angular velocity varies according to its distance from the sun.
# Output:
# Fsw = Daily average solar radiation in W/m^2.
# Can also output orbital parameters.
#
# This script contains orbital parameter data for the past 50000 years
# from Berger and Loutre (1991).
#
# Detailed description of calculation:
# Values for eccentricity, obliquity, and longitude of perihelion for the
# past 5 Myr are taken from Berger and Loutre 1991 (data from
# ncdc.noaa.gov). If using calendar days, solar longitude is found using an
# approximate solution to the differential equation representing conservation
# of angular momentum (Kepler's Second Law). Given the orbital parameters
# and solar longitude, daily average insolation is calculated exactly
# following Berger 1978.
#
# References:
# Berger A. and Loutre M.F. (1991). Insolation values for the climate of
# the last 10 million years. Quaternary Science Reviews, 10(4), 297-317.
# Berger A. (1978). Long-term variations of daily insolation and
# Quaternary climatic changes. Journal of Atmospheric Science, 35(12),
# 2362-2367.
#
# Authors:
# Ian Eisenman and Peter Huybers, Harvard University, August 2006
# eisenman@fas.harvard.edu
# This file is available online at
# http://deas.harvard.edu/~eisenman/downloads
# Translated into R by Thomas Laepple
# Suggested citation:
# P. Huybers and I. Eisenman, 2006. Integrated summer insolation
# calculations. NOAA/NCDC Paleoclimatology Program Data
# Contribution #2006-079.
#
# === Get orbital parameters ===
if (fast) temp<-orbital_parameters_fast(kyear) else temp<-orbital_parameters(kyear)
ecc<-temp$ecc
epsilon<-temp$epsilon
omega<-temp$omega
# For output of orbital parameters
obliquity<-epsilon*180/pi
long_perh<-omega*180/pi
# === Calculate insolation ===
lat=lat*pi/180 # latitude
# lambda (or solar longitude) is the angular distance along Earth's orbit measured from spring equinox (21 March)
if (day_type ==1) # calendar days
{
# estimate lambda from calendar day using an approximation from Berger 1978 section 3
delta_lambda_m <- (day-80)*2*pi/365.2422;
beta<-(1-ecc^2)^(1/2)
lambda_m0<-(-2)*( (1/2*ecc+1/8*ecc^3)*(1+beta)*sin(-omega)-1/4*ecc^2*(1/2+beta)*sin(-2*omega)+1/8*ecc^3*(1/3+beta)*(sin(-3*omega)))
lambda_m<-lambda_m0+delta_lambda_m
lambda<-lambda_m+(2*ecc-1/4*ecc^3)*sin(lambda_m-omega)+(5/4)*ecc^2*sin(2*(lambda_m-omega))+(13/12)*ecc^3*sin(3*(lambda_m-omega))
} else if (day_type==2) #solar longitude (1-360)
{lambda=day*2*pi/360 # lambda=0 for spring equinox
} else stop('Error: invalid day_type')
So<-1365 # solar constant (W/m^2)
delta<-asin(sin(epsilon)*sin(lambda)) # declination of the sun
Ho<-acos(-tan(lat)*tan(delta)) # hour angle at sunrise/sunset
# no sunrise or no sunset: Berger 1978 eqn (8),(9)
Ho[(abs(lat) >= (pi/2 - abs(delta)) ) & ( lat*delta > 0 )] <- pi
Ho[( abs(lat) >= (pi/2 - abs(delta)) ) & ( lat*delta <= 0 )] <- 0
# Insolation: Berger 1978 eq (10)
Fsw=So/pi*(1+ecc*cos(lambda-omega))^2 /(1-ecc^2)^2 * ( Ho*sin(lat)*sin(delta) + cos(lat)*cos(delta)*sin(Ho))
return(list(Fsw=Fsw,ecc=ecc,obliquity=obliquity,long_perh=long_perh,lambda=lambda/2/pi*360 ))
}orbital_parameters_fast<-function(kyear)
{
if ( length(orbital_global$ecc) != 50001) halt("orbital parameters not initalized")
return(list(ecc=orbital_global$ecc[kyear*10+1],epsilon=orbital_global$epsilon[kyear*10+1],omega=orbital_global$omega[kyear*10+1]))
}##############
orbital_parameters<-function(kyear,FILEDATA="ins_data.txt")
{
# === Load orbital parameters (given each kyr for 0-5Mya) ===
# Load the matrix contains data from Berger and Loutre (1991),
# downloaded as ORBIT91 from ncdc.noaa.gov
m<-read.table(FILEDATA)
kyear0<- -1*m[,1] # kyears before present for data (kyear0>=0);
ecc0<-m[,2] # eccentricity
# add 180 degrees to omega (see lambda definition, Berger 1978 Appendix)
omega0<-m[,3]+180; # longitude of perihelion (precession angle)
omega0=unwrap(omega0*pi/180)*180/pi # remove discontinuities (360 degree jumps)
epsilon0=m[,4] # obliquity angle
# Interpolate to requested dates
ecc=spline(x=kyear0,y=ecc0,xout=kyear)$y
omega=spline(x=kyear0,y=omega0,xout=kyear)$y * pi/180
epsilon=spline(x=kyear0,y=epsilon0,xout=kyear)$y * pi / 180
return(list(ecc=ecc,epsilon=epsilon,omega=omega))
}#Q = unwrap(P) corrects the radian phase angles in array P by adding multiples of ±2pi when absolute jumps between consecutive array elements are greater than pi radians.
# based on http://ccrma.stanford.edu/~jos/sasp/Matlab_listing_unwrap_m.html
unwrap<-function(p)
{
N = length(p)
up = rep(0,N)
pm1 = p[1]
up[1] = pm1
po = 0
thr = pi
pi2 = 2*pi
for (i in 2:N)
{
cp <- p[i] + po
dp <- cp-pm1
pm1 <- cp
if (dp>=thr)
{
while (dp>=thr)
{
po <- po - pi2
dp <- dp - pi2
}
}
if (dp<=((-1)*thr))
{
while (dp<=thr)
{
po <- po + pi2
dp <- dp + pi2
}
}
cp = p[i] + po
pm1 = cp
up[i] = cp
}
return(up)
}tlag<-function(data,ilag)
{
temp<-rep(data,3)
return(temp[366:730-ilag])
}### Cover functions to return the daily insolation, either using a classical calendar (aligned with March21) or using an alignment with Dec21 summer solstice
ins.march21<-function(kyear,LAT)
{
return(daily_insolation(kyear,LAT,1:365)$Fsw)
}
ins.dec21<-function(kyear,LAT)
{
r<-daily_insolation(kyear,LAT,1:365) #Eigentlich passt 356 besser
shift<-355-which.min(abs(r$lambda-270)) #dann entprechend hinschieben
r<-tlag(r$Fsw,shift)
return(r)
}ins.dec21.param<-function(ecc,obliquity,long_perh,LAT)
{
r<-daily_insolation_param(LAT,1:365,ecc,obliquity,long_perh) #Eigentlich passt 356 besser
shift<-355-which.min(abs(r$lambda-270)) #dann entprechend hinschieben1
r<-tlag(r$Fsw,shift)
return(r)
}orbital_global<-orbital_parameters((0:50000)/10,paste("ins_data.txt",sep=""))Plot of insolation for 5.000 kys
#Validation
#plot(june.65N)
june.65N.new<-vector()
for (i in 1:5000)
{
june.65N.new[i]<-daily_insolation(kyear=i,lat=65,day=172)$Fsw
}
plot(-(1:5000),june.65N.new,col="red",type='l')
wave <- june.65N.new
wave[june.65N.new>510] <- 510
plot(-(1:5000),wave,type="l"); title("overflowed, non-linear wave"); abline(h=mean(june.65N.new),lty=3)
Plot of orbital parameters
ecc.new<-vector()
obliquity.new<-vector()
lambda.new<-vector()
time<-vector()
for (i in 1:5000)
{
ecc.new[i]<-daily_insolation(kyear=i,lat=65,day=172)$ecc
obliquity.new[i]<-daily_insolation(kyear=i,lat=65,day=172)$obliquity
lambda.new[i]<-daily_insolation(kyear=i,lat=65,day=172)$lambda
time[i] <- -i
}
plot(time,ecc.new,col="red",type='l')
plot(time,obliquity.new,col="blue",type='l')
plot(time,lambda.new,col="black",type='l')
xx <-spectrum(obliquity.new,spans=5,main="spans=5"); abline(v=0:5/100,h=0:1000/1000,lty=3);
xxx<-spec.ar(obliquity.new, order = 1)
plot(xx)
lines(xxx$freq,xxx$spec)
xx <-spectrum(wave,spans=5,main="spans=5")
xxx<-spec.ar(wave, order = 1)
plot(xx,col="green"); abline(v=0:5/100,h=0:1000/1000,lty=3);
lines(xxx$freq,xxx$spec)