Happy Perihelion

When better than perihelion to think about orbital forcing?

I’ve wanted extract some insolation data for the late Quaternary to compare with some proxies as part of a rebuttal of some nonsense written elsewhere. Two insolation datasets are available from the Paleoclimatology World Data Centre.

The data from Berger (1978) is in a rather challenging format (one of the banes of working with palaeoclimate data). The data from Berger and Loutre (1991) are in a much simpler format, but only has insolation values for a reduced set of latitudes and only for June and December. The two data sets are essentially identical for the last 800 kyr – differences only become apparent on longer time scales.

Here is a function to read the Berger (1978) data into R.

read.berger78<-function(file){
 x<-read.table(file, comment.char = "d")#the comment.char argument makes it skip the rows beginning with d
 colnames(x)<-c("age.ky", "latitude", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
 x[,-(1:2)]<-x[,-(1:2)]*.4843 #convert to Watts/m2
 x$age.ky<-x$age.ky*-1 #invert age scale so +ve is past
 x
}

insol<-read.berger78("d:/bein1.dat.txt")#100 kyr in each file. 
#Use rbind() to combine multiple files.

Now we can plot the insolation for a chosen month and latitude with

jun60<-insol[insol$latitude==60,c("age.ky", "Jun")]
plot(jun60, type="l", col=2, ylab="Insolation, W/m2", xlab="Date ky BP")
title(main = "June Insolation at 60 °N")

So far so good. But I also what to plot a figure showing the difference in insolation between two time slices by month and latitude. Here I’m using filled.contour in a wrapper that extracts the insolation for two time periods and finds the difference (if given just one time period, it will plot that).

isol.contour<-function(x,t1=0, t2){
  d1<-insol[x$age.ky==t1, -(1:2)]
  if(!missing(t2)){
    d2<-x[x$age.ky==t2, -(1:2)]
    dat<-d2-d1
    zlim=range(dat)
    zlim=range(c(zlim, -zlim))  #ensure scale is symmetrical
    main=paste(t2, "kyr BP -", t1, "kyr BP insolation")
  }else{
    dat<-d1                                                          
    zlim=range(dat)
    main=paste(t1, "kyr BP insolation")
  }
  dat<-dat[nrow(dat):1,]
  dat<-t(as.matrix(dat))  

  cols<- colorRampPalette(c("blue","grey80", "red"))
  filled.contour(x=1:12, y=seq(-90,90,10), z=dat, zlim=zlim, color.palette=cols,
   plot.axes={
                axis(side=1, at=1:12, label=names(x)[-(1:2)], las=2)
                axis(2, at=c(-60,0, 60))
              },
    plot.title={
                title(xlab="Months", ylab="Latitude",main=main) 
              },
    key.title={
                mtext(side=4, text=expression("Wm"^-2), line=1.5, las=0)
              }
              
   )
}

x11(5,5);par(mar=c(3,3,1.5,1), mgp=c(1.6,.5,0), tcl=-.3)
isol.contour(x=insol,t1=0, t2=10)

The plot is very pretty.

10 kyr BP – 0 kyr BP insolation

But not exactly what I had expected. This is what I had anticipated.

10 kyr BP – 0 kyr BP. Apologies for the rainbow scale.

I checked the code is correct (it’s a subtraction – cannot go far wrong), and then gone back to re-read Berger (1978).

The issue, I think, is in how months are defined. Berger (1978) gives two options. Both start the year on the spring equinox. The first option places the mid-point of each month at 30° intervals around the Earth’s orbit of the Sun (12 months * 30°= 360°). However, as the Earth’s orbit is not circular, following Kepler’s second law, the Earth does not move at a constant angular velocity. When the Earth is closest to the Sun, perihelion, it moves fastest, covering 30° in less than a month. The second option calculates the position of the Earth every month and finds the insolation at that position in the orbit.

That Berger (1978) provides the insolation calculated with the first option (the “mid-month” option) implies it is the more useful. I’m not sure why, and most plot I can find seem to use the second option. Can anybody shine some light on this for me? Or do I need to understand Laskar et al (2004)?

Advertisement