Guest Post by Willis Eschenbach

After I’d published my previous post on the Hurst Exponent entitled A Way To Calculate Effective N, I got an email from Dan Hughes which contained a most interesting idea. He thought that it would be productive to compare the Hurst analysis of the records of weather phenomena such as temperatures and the like, with the Hurst analysis of the corresponding climate model outputs. He proposed that we take a look at the question and discuss and share our findings. Hey, what’s not to like?

Along the way, as such things happen, the topic of our discussion and investigation turned to a larger and more interesting question—`which of the various natural datasets (temperatures, rainfall, eruptions, pressures, etc) and/or the corresponding global climate model outputs are chaotic?

I must, of course, start with the obvious question … what is meant by “chaotic”? A chaotic system is a system wherein nearby initial states either converge or diverge exponentially. An example is the ocean surface. If you drop two sealed empty bottles overboard in mid-ocean, one on each side of a boat, they will drift apart over time. This separation will be slow at first, and then faster and faster as the two bottles encounter different winds and currents in different areas.

Whether a dataset is chaotic or not is generally assessed by looking at the Lyapunov exponent. This is a measure of the “stretching factor”. In our ocean example above, the stretching factor measures how fast the two points are moving apart over time. In a chaotic dataset, the stretching factor generally increases or decreases with time. In non-chaotic datasets, on the other hand, the stretching factor doesn’t vary with time. Figure 1 shows the “Lyapunov curves” of the evolution of the stretching factor with time for a variety of natural and calculated datasets.

Figure 1. Lyapunov curves for a variety of datasets. All datasets have been detrended and standardized before analysis.

Now, I’ve graphed four types of datasets above, indicated by the four colors. The first type, shown in red and mostly obscured by the blue lines, shows four different varieties of random numbers—normal, uniform, poisson, and high Hurst exponent fractional Gaussian random numbers. Basically the Lyapunov curves of the random number datasets are all plotting right on top of each other. Starting from time = 0, they climb rapidly to their maximum value and then just stay there. As we would expect from random data, there’s no trend in the stretching factor over time.

The next group, in blue, shows the Lyapunov curves for a half-dozen climate-related datasets, viz:

• HadCRUT4 Monthly Mean Surface Air Temperatures 1850-2015

• Annual Nilometer Minimum River Height 622-1284

• Stockholm Monthly Tides 1801-2001

• Central England Temperature Record Daily Maximum Temperatures 1878-2015

• Armagh Ireland Daily Mean Temperatures 1865-2001

• Annual Average Nile Flow, Cubic Metres/Sec 1870-1944

As you can see, in terms of the Lyapunov analysis, all six of these climate-related datasets (blue lines) are indistinguishable from the four random datasets (red lines), which in turn are indistinguishable from each other. None of them show any trace of chaotic behavior.

Another group of datasets, those at the bottom in gold colors, are quite different from the random and the observational datasets. They are a variety of chaotic datasets. Note that they all share a common factor—as mentioned above, over time the rate of separation (as measured by the “stretching factor”) increases. The rate of separation doesn’t just go to a certain high point and stay there like the random or climate-related datasets did. The rate of separation in chaotic datasets continues to rise over time.

Finally, there are a couple of other datasets in purple. These show observations of phenomena that are usually thought of as “forcings”. One is changes in the solar activity, with daily sunspots as a proxy for the activity, from 1880 to 2015. The other is annual aerosol optical depth from 800 to 2000, which is generally a function of volcanic action and is calculated from ice core data. Curiously, these two datasets plot somewhere in between the random observations at the top, and the chaotic datasets at the bottom. In addition, both of them show significant variation in stretching factor over time. The sunspots have a slight but significant increase. The aerosol optical depth goes down and then looks like it’s starting back up. So it seems that these two datasets are weakly chaotic.

Now, these results were a great surprise to me. I’ve long believed, without ever checking it, that the climate was chaotic … however, this analysis shows that at least those six observational datasets that I analyzed above are not chaotic in the slightest. Hey, what do I know … I was born yesterday.

Are there any climate datasets which are chaotic, even weakly chaotic? I think so. It appears that tropical ocean temperatures are weakly chaotic … but that’s a question for the next post, which will look into Dan Hughes’s idea regarding the use of the Hurst analysis to distinguish between chaotic and non-chaotic datasets.

Regards to everyone,

w.

MY USUAL REQUEST: If you disagree with someone, please quote the exact words you disagree with. That way, we can all understand the exact nature of your objection.

CODE: To calculate the Lyapunov exponent I’ve used the lyap_k function from the R package tseriesChaos. Here are the functions I used to make Figure 1:

lineout=function(col="black",testdata=testdata,line0=F,lwd=3,s=110){
  output <-lyap_k(scale(dtrendlin(testdata)), m=3, d=2, s=s, t=40, ref=100, k=4, eps=4)
  lines(output,main="",cex.main=.95,col=col,lwd=lwd)
  invisible(output)
}

plotout=function(col="black",testdata=testdata,line0=T,lwd=3,s=110){
  output <-lyap_k(scale(dtrendlin(testdata)), m=3, d=2, s=s, t=40, ref=100, k=4, eps=4)
  plot(output,main="",cex.main=.95,new=F,ylim=c(-6,1),col=col,lwd=lwd,
       ylab="Log(Stretching Factor)")
  if (line0) abline(h=c(0,1),lty="dotted")
  invisible(output)
}

dtrendlin=function(timeseries,doplot=FALSE){
  outseries=timeseries
  thelm=lm(timeseries~time(timeseries))
  thegood=which(is.finite(timeseries))
  outseries[thegood]=outseries[thegood]-thelm$fit
  if (doplot){
    plot(outseries)
  }
  outseries
}