Gravitational interactions among the many bodies in the Solar System continuously induce small perturbations in the Earth's orbit around the Sun, which in turn cause subtle changes in the distribution of solar radiation by latitude and season over geologic time. In the early 20th Century, Serbian mathematician Milutin Milankovitch proposed that these variations might be related to the timing of the Ice Ages. Because these changes are pseudo-cyclical, they are known as Milankovitch Cycles. For climate, the most import variations are changes in the eccentricity of the ellipse that describes the Earth's orbit, the tilt of the Earth's axis, and the precession of the equinoxes) which affects the seasonal distribution of radiation.
These changes can be predicted (and hindcast) by numerical integration of the multibody equations of celestial mechanics. Jacques Laskar and his colleagues at l'Observatoire de Paris has performed these calculations for 250 million years into the past and for 20 million years into the future. The data and plots here are derived from the output of their calculations, which they have generously made available on their website. The full model, calculations, results, and analysis is described in their paper (Laskar et al, 2011; available from Astronomy and Astrophysics via Open Access ).
The daily mean top-of-the atmosphere incoming solar radiation (insolation) by latitude and season is determined from the orbital parameters using well-known formulae based on spherical trigonometry (see derivation on wikipedia). For all the details, please see the “website code” tab to the right.
Some excellent web-based educational resources on Milankovitch forcing of climate are:
This website is controlled using the R package “shiny.” This one was a lot of work!
There are 10 components:
Scroll down or click links in the list above to read all about it!
laskar.R:
laskar <- function(start.time=-1000, end.time=20){
# Read timeseries of Earth's orbital eccentracity, obliquity, and precession
# from a big binary file based on the data of Laskar et al (2004)
# Subset the data between start.time and end.time (thousands of years from now)
# then calculate top-of-atmosphere daily mean insolation at 65N June solstice
# SOURCE: http://www.imcce.fr/fr/presentation/equipes/ASD/insola/earth/earth.html
load('data/orbit.LA2004.rda')
data <- subset(orbit, kiloyear>= start.time & kiloyear<= end.time)
insolation <- Qday.65N.solstice(eccentricity=data$eccentricity,
obliquity=data$obliquity,
perihelion=data$perihelion)
return(data.frame(time=data$kiloyear, eccentricity=data$eccentricity,
obliquity=data$obliquity, perihelion=data$perihelion,
insolation=insolation))
}
Qday.65N.solstice.R:
Qday.65N.solstice <- function(S0=1367, eccentricity=0.0167,
obliquity=0.4091, perihelion=1.7963){
# Given Earth's orbital parameters (eccentricity, obliquity, & perhihelion),
# Calculate the daily mean top-of-atmosphere insolation at 65N June solstice
# For a derivation see http://en.wikipedia.org/wiki/Insolation
latitude <- 65. * pi / 180.
Q.day <- S0*(1+eccentricity*sin(perihelion+pi))^2 *sin(latitude)*sin(obliquity)
return(Q.day)
}
plot.timeseries.R:
plot.timeseries <- function(orbit){
# Make a "ladder" plot of orbital parameters over time
# attach links names in the list orbit to variables in the script
attach(orbit)
# Remember original plotting parameters and then set them the way we want
orig.par <- par(no.readonly=TRUE)
par(mfrow=c(4,1), mar=c(0,4,2,1), oma=c(5,1,0,0), cex.lab=1.3, cex.axis=1.3)
# Top plot is eccentricity of Earth's orbit
plot(time,eccentricity,type='l',lwd=2,col='red',main='Eccentricity',
xaxt='n', xlab='',ylab='')
grid(col='black', ny=0)
abline(v=0, lwd=2)
# Second plot is obliquity of rotational axis
plot(time,180*obliquity/pi,type='l',lwd=2,col='forestgreen',main='Obliquity',
xaxt='n', ylab='Degrees',xlab='')
grid(col='black', ny=0)
abline(v=0, lwd=2)
# Third plot is the precession of the equinox around the orbit
plot(time,eccentricity*sin(perihelion),type='l',lwd=2,col='dodgerblue',
main=expression(bold(paste('Precession Index: ',e(1+sin(bar(omega)))))),
xaxt='n', xlab='',ylab='')
grid(col='black', ny=0)
abline(v=0, lwd=2)
# Bottom plot is the resulting changes in summer solar radiation at 65N
plot(time,insolation,type='l',lwd=2,col='goldenrod1',
main='Mean Daily Insolation at 65N on Summer Solstice',
ylab='Watts/m2', sub='kiloyears')
grid(col='black', ny=0)
abline(v=0, lwd=2)
# Add an overall x-axis to the page
mtext(side=1, line=2.5, 'kiloyears from present')
# reset plotting parameters to original values
par(orig.par)
# release variable names
detach(orbit)
}
global.insol.R:
global.insol <- function(S0=1367, ecc=0.0167, eps=0.4091, perih=1.7963,
nLats=181, nLons=181){
#######################
# Given Earth's orbital parameters (eccentricity, obliquity, and precession),
# calculate daily mean top-of-atmosphere incoming solar radiation
# for each latitude and time of year, then return the results as a 2D array
# (longitude here refers to time of year ...
# it's the "longitude" of the Earth in it's orbit around the Sun:
# Lon=0 is the March Equinox, Lon=90 is June Solstice
# Lon=180 is September Equinox, Lon=270 is December Solstice)
# Most of this function is just prep, the real calculation happens
# just before the return value is assigned: Qday as a function of orbit
# For details and a derivation see http://en.wikipedia.org/wiki/Insolation
#######################
# Set up arrays to hold lat, lon, and insolation for each point
lat <- replicate(nLats,NA)
lon <- replicate(nLons,NA)
Qday <- matrix(nrow=nLats, ncol=nLons)
# Loop over all the points
for (j in 1:nLats) for (i in 1:nLons) {
# latitude and longitude for this [i,j] pair
lat[j] <- (-90 + (j-1) * 180/(nLats-1)) * pi / 180
lon[i] <- (-180 + (i-1) * 360/(nLons-1)) * pi / 180
# Declination for this time of year (longitude)
# declination is the latitude where the Sun is directly overhead
dec <- eps * sin(lon[i])
# Hour angle at sunrise
arg <- tan(lat[j])*tan(dec)
if (arg > 1) {
h0 <- pi # perpetual summer daylight
} else if (arg < -1) {
h0 = 0} # perpetual winter darkness
else {
h0 <- acos(-arg) # normal sunrise/sunset
}
# Earth-Sun distance (ratio of mean to actual)
R0.div.RE <- 1 + ecc * cos(lon[i] - perih)
# Daily mean insolation for this lat and time
Qday[i,j] <- S0/pi * (R0.div.RE)^2 *
(h0 * sin(lat[j]) * sin(dec) + cos(lat[j]) * cos(dec) * sin(h0))
}
return(insol=list(lat=lat*180/pi, lon=lon*180/pi, Qday=Qday))
}
plot.global.insol.R:
plot.global.insol <- function(insol, plot.title='Reference time', diff=FALSE){
# Make a color-filled contour plot of global TOA insolation by lat & season
# optionally use a nice difference color bar
library(fields) # Need this to make the image.plot
# x-coordinate is longitude of Earth in orbit relative to March Equinox
# This is essentially the date in the year, so label with month names
x.labels =c('SEP','OCT','NOV','DEC','JAN','FEB','MAR',
'APR','MAY','JUN','JUL','AUG','SEP')
# y-coordinate is latitude on Earth
y.labels = c('SP','60 S','30 S','EQ','30N','60 N','NP')
# Optionally center colors on zero for difference plots
if (diff) {
z.limits <- c(-max(abs(insol$Qday)),max(abs(insol$Qday)))
} else {
z.limits <- c(min(insol$Qday),max(insol$Qday))
}
# Draw an image using colors to represent the insolation
image.plot(insol$lon,insol$lat,insol$Qday, xlim=c(-180,180), ylim=c(-90,90),
zlim=z.limits, xlab='', ylab='Latitude', axes=FALSE)
# Add contours lines for insolation
contour(insol$lon,insol$lat,insol$Qday, axes=FALSE, add=TRUE)
# Add nicely-labeled axes for months and latitude
axis(1, at=seq(-180,180,30), labels=x.labels)
axis(2, at=seq(-90,90,30), labels=y.labels)
box(lwd=2)
# Add titles
title(plot.title,line=2)
title(expression(W ~~ m^-2), line=1)
}
compare.two.times.R:
compare.two.times <- function(orbit, time.1, time.2) {
# Given a timeseries of the Earth's orbital parameters, select two times
# (specified in input as time.1 and time.2, thousands of years from now)
# For each specified time, calculate and plot the distribution of
# insolation at the top of the atmoshere by latitude and season
# Make a nice color-filled plot of geography of sunlight at time.1
# and subtract to show differences in distribution of insolation at time.2
# Extract orbital parameters for time.1
first <- subset(orbit, time == time.1)
# Extract orbital parameters for time.2
second <- subset(orbit, time == time.2)
# Calculate geography of top-of-atmosphere sunshine for time.1
first.insol <- global.insol(ecc=first$eccentricity,
eps=first$obliquity,
perih=first$perihelion)
# Calculate geography of top-of-atmosphere sunshine for time.2
second.insol <- global.insol(ecc=second$eccentricity,
eps=second$obliquity,
perih=second$perihelion)
# Subtract insolation distribution for time.2 from time.1
insol.diff <- list(lat=first.insol$lat,
lon=first.insol$lon,
Qday=first.insol$Qday - second.insol$Qday)
# Set graphics parameters to show a vertical stack of two plots on the page
orig.par <- par(no.readonly=TRUE)
par(mfrow=c(2,1))
# Plot the distribution of insolation at time.1 by latitude & season
plot.global.insol(first.insol, plot.title=paste(time.1,'kyr'))
# Plot the distribution of insolation differences at time.1
plot.global.insol(insol.diff, diff=TRUE,
plot.title=paste(time.1,'kyr minus',time.2,'kyr'))
# Reset graphics parameters to their original values
par(orig.par)
}
sensitivity.Qday.R:
sensitivity.Qday <- function(eccentricity, obliquity, perihelion) {
# Given the Earth's orbital parameters, calculate and plot the distribution of
# insolation at the top of the atmoshere by latitude and season, and compare
# to the distribution in the present day
# Calculate present-day distribution of insolation by latitude & season
present.day <- global.insol()
# Calculate perturbed distribution of insolation by latitude & season
perturbation <- global.insol(ecc=eccentricity,
eps=obliquity,
perih=perihelion)
# Subtract perturbed from present-day insolation by latitude & season
insol.diff <- list(lat=present.day$lat, lon=present.day$lon,
Qday=present.day$Qday - perturbation$Qday)
# Set up graphics parameters for a column of two plots stacked on the page
orig.par <- par(no.readonly=TRUE)
par(mfrow=c(2,1))
# Draw the plot of present-day insolation
plot.global.insol(present.day, plot.title='Present Day Conditions')
# Draw the plot of differences between perturbed and present conditions
plot.global.insol(insol.diff, diff=TRUE,
plot.title='Differences with Perturbed Orbital Parameters')
# Label the lower plot with the perturbed orbital parameters
mtext(substitute(e==esub ~","~ epsilon==osub^o ~","~ omega==psub^o,
list(esub=eccentricity,
osub=obliquity * 180/pi,
psub=perihelion* 180/pi)),
side=1, line=3, cex=2)
# Reset graphics parameters to original values
par(orig.par)
}
write.orbital.data.R:
write.orbital.data <- function(file, orbital.data) {
sink(file)
cat('Data from http://www.imcce.fr/fr/presentation/equipes/ASD/insola/earth/earth.html\n')
cat('Five columns for each time as follows:\n')
cat(' time: thousands of years from present\n')
cat(' eccentricity: (aphelion-perihelion)/(aphelion+perihelion)\n')
cat(' obliquity: radians\n')
cat(' perihelion: angle between equinox & perihelion, radians\n')
cat(' insolation: daily mean TOA solar flux at 65N summer solstice, W/m2\n')
sink()
write.table(orbital.data, file=file,row.names=F,col.names=T,
quote=FALSE, append=TRUE, sep="\t")
}
ui.R:
library(shiny)
library(markdown)
# Define user interface for Milankovitch Data Viewer
shinyUI(pageWithSidebar(
# Application title
headerPanel("Milankovitch Orbital Data Viewer"),
# Sidebar on the left with controls for the user to specify times
sidebarPanel(
# Display an image that shows the definition of each orbital parameter
img(src='orbital.jpg', align='left'),
# Allow the user to choose the start time and duration of the orbital data
conditionalPanel(condition="input.tabset == 'tab1'", wellPanel(
br(),
h4('Select start and end of timeseries'),
sliderInput("start.year", "Beginning of Timeseries (kiloyears from now)",
min = -20000, max = 10000, value = -1000, step= 10),
sliderInput("duration", "Duration: Thousands of Years",
min = 1, max = 10000, value = 1050, step= 10)
)),
# Choose two times to compare geography of insolation between them
conditionalPanel(condition="input.tabset == 'tab2'", wellPanel(
br(),
h4('Compare Distribution of Insolation at Two Different Times'),
sliderInput("first.time", "First time (kiloyears from now)",
min = -1000, max = 1000, value = 0, step= 10),
sliderInput("second.time", "Second time (kyr from now)",
min = -1000, max = 1000, value = -100, step= 10)
)),
# Choose arbitrary orbital parameters and compare insolation to present
conditionalPanel(condition="input.tabset == 'tab3'", wellPanel(
br(),
h4('Sensitivity of Insolation to Different Orbital Parameters'),
h5('Select orbital parameters, compare to present-day'),
sliderInput('eccentricity', 'Eccentricity',
min=0, max=0.060, value=0.0167, step=0.001),
sliderInput('obliquity', 'Obliquity (degrees)',
min=22, max=24.5, value=23.44, step=0.01),
sliderInput('long.perihelion', 'Longtidue of Perihelion (omega)',
min=-180, max=180, value=70.3, step=0.1)
))
),
# Show a tabset in the main panel of the browser that displays model output
mainPanel(
tabsetPanel(id='tabset',
# First tab shows plots of model output for both control and modified atmosphere
# Also displays explanatory text, links to sources, and offers data download
tabPanel("Timeseries", value='tab1',
h4('Milankovitch Orbital Parameters over Time'),
# Explanatory text about the data and sources
h5("Based on the calculations of:"),
helpText("J. Laskar et al (2004) A long term numerical solution for the
insolation quantities of the Earth. Astronomy and Astrophysics
428, 261-285, DOI: 10.1051/0004-6361:20041335"),
h5('Go to ',a('Celestial Mechanics at Observatoire de Paris',
href='http://www.imcce.fr/fr/presentation/equipes/ASD/insola/earth/earth.html',
target='_blank')),
h5('Download paper', a('Laskar et al (2011)',
href='http://www.aanda.org/articles/aa/pdf/2011/08/aa16836-11.pdf',
target='_blank')),
# Offer to download orbital data through user's web browser
downloadButton('downloadData',label='download orbital data'),
# Make the nice 4-panel tiemseries plot
plotOutput('timeseries'),
tags$style(type="text/css", ".tab-content { overflow: visible; }")),
# Second tab shows global/seasonal distribution of insolation at reference time
tabPanel("Compare", value='tab2',
h4('Global and Seasonal Distribution of Incoming Solar Radiation'),
h5('Compare distribution of daily mean radiation at two different times:'),
plotOutput('compare.times'),
tags$style(type="text/css", ".tab-content { overflow: visible; }")),
# Third tab compares insolation distribution between ref time and perturbations
tabPanel("Sensitivity", value='tab3',
h4('Global/Seasonal Distribution of Differences'),
h5('Difference: Perturbed Orbit - Reference Time:'),
plotOutput('sensitivity'),
tags$style(type="text/css", ".tab-content { overflow: visible; }")),
# Fourth tab displays a brief model description
tabPanel("Description", value='tab4',
includeMarkdown('doc/model.description.md')),
# Fifth tab displays website code
tabPanel("Website code", value='tab5',
includeMarkdown('doc/website.code.md'))
)
)
))
server.R:
library(shiny)
# Source required R scripts
source('model/laskar.R')
source('model/write.orbital.data.R')
source('model/plot.timeseries.R')
source('model/Qday.65N.solstice.R')
source('model/global.insol.R')
source('model/plot.global.insol.R')
source('model/compare.two.times.R')
source('model/sensitivity.Qday.R')
shinyServer(function(input, output) {
# Capture orbital timeseries for later display
orbital.data <- reactive({
laskar(start.time=input$start.year,
end.time=input$start.year+input$duration)
})
# 4-panel plot of orbital timeseries including insolation at 65N
output$timeseries <- renderPlot(
plot.timeseries(orbital.data()), height=800, width=600)
# Download orbital data as a text file from web browser
output$downloadData <- downloadHandler(
filename = 'milankovitch.data.txt',
content = function(file) {
write.orbital.data(file, orbital.data())
})
# Contour plot of insolation differences for given orbital parameters
output$compare.times <- renderPlot(
compare.two.times(orbital.data(),input$first.time, input$second.time),
height=950,width=600)
# Explore sensitivity to orbital parameters:
output$sensitivity <- renderPlot(
sensitivity.Qday(input$eccentricity, input$obliquity * pi/180,
input$long.perihelion * pi/180),
height=950, width=600)
})