Milankovitch Orbital Data Viewer

Milankovitch Orbital Parameters over Time

Based on the calculations of:
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

Global/Seasonal Distribution of Differences

Difference: Perturbed Orbit - Reference Time:

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:

How this website works (including all the code!)

This website is controlled using the R package “shiny.” This one was a lot of work!

There are 10 components:

• A script (laskar.R) that reads the orbitral data from a file
• A function (Qday.65N.solstice.R) that calculates daily mean insolation at 65N from the given orbital parameters
• A script plot.timeseries.R that draws the nice plots on the browser
• A function global.insol.R that computes daily mean top-of-atmosphere insolation by latitude & season from Earth's orbital parameters (eccentricity, obliquity, and longitude of perihelion)
• A script plot.global.insol that plots the distribution of insolation by latitude & season
• A script compare.two.times that compares the distributon of insolation at two different times chosen by the user
• A script sensitivity.Qday that lets the user explore the sensitivity of the distribution of insolation to the Earth's orbital parameters
• a function write.orbital.data that write a text file of orbital data and lets the user download it through the web browser
• The user interface (ui.R) that controls the sliders and passes inputs to the model
• The “server” (server.R) that responds to user interface events in the web browser

Read and subset orbital data between user-specified times

``````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

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))
}
``````

Calculate top-of-atmosphere insolation at 65N at June Solstice

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)
}
``````

Make a 4-panel “ladder plot” of orbital parameters over time

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

}
``````

Calculate distribution of insolation by latitude and season

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 distribution of insolation by latitude and season

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

# 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)

title(plot.title,line=2)
title(expression(W ~~ m^-2), line=1)

}
``````

Compare distribution of insolation at two specified times

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
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)
}
``````

Calculate sensitivity of insolation by latitude & season to orbital parameters

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
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 to a text files & download via web browser

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('  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")
}
``````

This is the user interface that actually controls this website

ui.R:

``````library(shiny)
library(markdown)

# Define user interface for Milankovitch Data Viewer

shinyUI(pageWithSidebar(

#  Application title

# 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
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')),
href='http://www.aanda.org/articles/aa/pdf/2011/08/aa16836-11.pdf',
target='_blank')),

# 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'))
)
)
))
``````

This is the server function that responds to user input

server.R:

``````library(shiny)

# Source required R scripts
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({
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)

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)

})
``````