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

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

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

Advertisements

About richard telford

Ecologist with interests in quantitative methods and palaeoenvironments
This entry was posted in climate, R and tagged , . Bookmark the permalink.

5 Responses to Happy Perihelion

  1. Nick Stokes says:

    I looked through your code and ran it – seems perfectly OK to me. And I can’t see than the definition of month can change much. As you say, Berger’s first is strange, but it can only squeeze up the month axis a bit. I don’t think it can produce the sharp polar events of the rainbow plot.

    Laskar says the insolation calcs are classical and refers to Berger 1978, and says he’s only presenting for clarity.

    • I think what’s happening at the pole is that spring arrives a few days earlier 10 kyr BP, and the difference between no sun and 24 hour sun is huge.
      I can see why climate modellers would want to use the exact daily insolation values (rainbow plot), but I’m not particularly interested in massive anomalies caused by spring arriving a few days early, except in that they indicate a short winter.

  2. Magma says:

    I can reproduce your second figure using Ian Eisenman and Peter Huybers’ MATLAB code to calculate paleoinsolation values (http://eisenman.ucsd.edu/code.html).
    But as for the old bein1.dat data from NCDC, I can see some obvious problems right off the start. For one thing, insolation at both the north and south poles has seven (7) months with a value of 0 and only five with non-zero values. Something of a red flag, no? The reason for the mid-month calculation rather than the monthly mean (or daily values) would have been to reduce the number of calculations by a factor of 30, an obsolete consideration these days. I suggest using the recalculated data using the original algorithms as coded by Eisenman & Huybers, unless there is a specific reason to use the historical data, which by now are of largely unknown provenance and contain multiple inconsistencies with the recalculated values.
    I suspect formatting will be ugly, but here are the values in cal/cm2/day from the bein1.dat file and from the recalculated mid-month value (15th of each month except for February where I used the 14th).

    bein1.dat file from NCDC, 0 ka BP
    90 0 0 0 553 944 1082 938 546 0 0 0 0
    80 0 0 156 545 930 1065 923 538 154 0 0 0
    70 0 87 308 604 887 1016 881 597 304 86 0 0
    60 92 229 450 700 904 984 898 692 444 226 91 50
    50 228 374 579 786 939 994 933 776 571 370 227 177
    40 374 513 690 852 961 996 954 842 680 507 372 322
    30 518 639 780 895 961 980 954 884 769 631 514 468
    20 651 747 846 912 938 940 931 901 835 738 646 609
    10 768 834 887 903 890 878 884 892 875 824 763 737
    0 866 896 901 867 819 794 813 857 888 886 861 848
    -10 942 933 887 807 726 690 721 797 875 922 935 937
    -20 992 942 846 723 615 570 611 714 835 931 986 1004
    -30 1017 924 780 618 489 439 486 611 769 913 1010 1046
    -40 1017 880 690 496 354 301 351 490 680 870 1010 1064
    -50 994 812 579 362 216 166 214 358 571 802 987 1062
    -60 957 724 450 221 86 47 86 219 444 715 950 1050
    -70 939 624 308 84 0 0 0 83 304 617 932 1085
    -80 984 563 156 0 0 0 0 0 154 556 977 1137
    -90 999 572 0 0 0 0 0 0 0 565 992 1155

    recalculated values (converted to cal/cm2/day from W/m2), 0 ka BP
    90 0 0 0 464 894 1080 1005 678 161 0 0 0
    80 0 0 104 458 881 1064 990 667 242 11 0 0
    70 0 62 257 551 844 1015 945 683 383 139 34 0
    60 78 199 403 659 878 984 938 759 515 285 160 64
    50 212 345 538 754 921 996 962 829 632 427 304 195
    40 359 487 657 829 949 998 975 882 731 559 447 341
    30 504 618 756 880 956 982 968 911 808 675 581 487
    20 640 732 831 906 938 943 938 915 860 771 701 625
    10 762 826 882 906 896 881 883 894 886 845 802 749
    0 865 896 906 878 829 797 806 847 886 894 880 856
    -10 946 940 903 825 740 694 709 776 858 916 935 941
    -20 1002 958 872 747 632 574 594 684 805 910 962 1002
    -30 1032 948 814 648 509 442 465 572 727 878 964 1038
    -40 1038 912 732 530 374 304 329 446 628 820 939 1049
    -50 1021 852 628 398 236 168 192 310 509 739 890 1039
    -60 992 772 505 257 103 48 68 172 375 637 823 1018
    -70 992 683 367 115 4 0 0 46 231 521 753 1036
    -80 1040 651 220 1 0 0 0 0 83 412 754 1085
    -90 1056 661 118 0 0 0 0 0 0 407 766 1102

    • Magma says:

      Yes. Formatting was particularly ugly. I suppose I should have double spaced between paragraphs too.

      • Thanks. I guess seven mid-months with zero insolation is correct but misleading. Mean monthly would be more interpretable. I guess it makes most difference at the poles near the equinox.

        Thanks for pointing to me to the Eisenman & Huybers code. I had read the paper but not found the code. Since everybody seems to be basing their calculations on Berger and Loutre (1991), I guess the results from Eisenman & Huybers will be the same as those from Laskar et al. The latter has a friendly web interface (& fortran), rather than matlab which I don’t think I have installed any more.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s