Archive for the ‘smooth’ Category


Rahmstorf (2009): Off the mark again (Part 12). A mathematical comedy

February 13, 2011

Here is one more post about the laughably bad PNAS “Global Sea level linked to global temperature” by Vermeer and Rahmstorf.  Will this  fount of absurdity never run dry?

Much has been said about Rahmstorf’s data smoothing techniques.  But the little gem you are about read may make your head spin.

Remember the Chao reservoir correction?  This was the correction that VR2009 applied to the Church and White sea level data to compensate for water that has been impounded in man-made reservoirs.  Never mind the fact that VR2009 paid lip service to, but did not include, a counter-correction for water that has been pumped from the aquifers and has artificially added to the sea level.  Let’s look at some details of how VR2009 handled this correction.

Here is something amazing…

VR2009 had the 2006 Church and White sea level data, which is rather noisy.  They also had the Chao reservoir correction data, which is also noisy.  They correctly saw the need to smooth the noisy data.  It seems that they could have done it one of  two ways: smooth each set separately, then  add the smoothed Chao data to the smoothed Church and White data, or add the unsmoothed Chao data to the unsmoothed Church and White data and then smooth the result.

When I reproduced VR2009’s basic algorithm, I choose the first method.  But VR2009 doubled up on smoothing the Choa reservoir correction.  They smoothed the Chao data, added it to the unsmoothed Church and White data, then smoothed the sum again.  So, the Chao data was effectively smoothed twice.

But here is the really amazing thing:  Look at the overlay of Chao’s data, VR2009’s smooth for the Chao data, and my smooth for the Chao data…

Figure 1. Chao correction to sea level rise rate with VR2009 smooth and Moriarty smooth

Wow! All I can say is “Wow!”  Can you believe how terrible the VR2009 fit for the additional sea level rise rate is?  It’s just amazingly bad! 

How did VR2009 come up with this bizarre data smooth?

In the Matlab program file that VR2009 uses to find the relationship between sea level and temperature (sealevel2.m, get copy here) they first import the unsmoothed Church and White data (church_13221.txt, get copy here) with the following code…

% load the church & white sea level data
load church_13221.txt;
seayear = church_13221(:,1);
sealevel = church_13221(:,2)/10;

Two arrays are created, one with the year, one with the sea level.  The “/10” in the last line of code converts the sea level data from mm to cm.

Then they apply their Chao reservoir correction.  Instead of importing a time series with the Chao data, they apply a function…

% Apply Chao et al (2008) reservoir correction:
if chao == ‘y’
     sealevel = sealevel + 1.65 + (3.7/3.1415)*atan2(seayear-1978,13);

So, VR2009 claims the term “1.65 + (3.7/3.1415)*atan2(seayear-1978,13)” is a representation of the Chao reservoir correction.  Figure 1, above shows the derivative of the Chao reservoir correction (which you can see as figure 3 in Chao’s Science paper).  So the derivative of VR2009’s Chao correction term should at least be close to the derivative provided in Chao’s paper.  Alas, instead it looks like the blue peak in figure 1, above. 

How did VR2009 come up with this strange correction that “fits” the Chao reservoir correction to an inverse tangent (atan2) function?  VR2009 claims to use sophisticated single spectrum analysis (SSA) to smooth its sea level and temperature data.  But their SSA code yields a numerical result, not an analytic one (that is, a time series of numbers, not a formula).  So SSA was NOT used to generate VR2009’s Chao correction term.

If you use my smooth of the Chao data as a baseline, then the VR2009 fit is about 0.2 mm too low around 1960 and about 0.3 mm too high by 1980.  By using their fit to the Chao reservoir sea level rise rate correction, they have effectively increased the sea level rise rate from 1960 to 1980 by an additional 0.5 mm per year.  They have pushed the Chao sea leve rise rate correction to later in the century which, of course, fits their general theme.

The following plot shows the 2006 Church and White sea level data with the questionable VR2009 version of the  Chao reservoir correction data and my version of the Chao reservoir correction.  At first they do not look much different.  But consider this: The VR2009 version causes the average sea level rise rate from 1950 to 1970 to be 1.66 mm/year, and for 1970 to 1990 to be  1.99 mm/year.  That’s a 16% increase.  If my version is used there is an average DECREASE in sea level rise rate, from 1.87 mm/year to 1.78 mm/year.  That is a 5% drop.  Look at figure 1, above, and ask yourself “Whose smooth of the Chao data is better?”

I will not attempt to assign motivation for this laughably bad smooth of the Chao reservoir correction data.  Suffice it to say that it is just one more in long series of blunders and bizarre consequences for VR2009.

Read more about the comedy known as the PNAS “Global Sea level linked to global temperature” by Vermeer and Rahmstorf.


Rahmstorf (2009): Off the mark again (part 8): Reproducing VR2009 results

October 14, 2010

Here is my effort to reproduce the Vermeer and Rahmstorf’s model linking sea level rise rate to temperature.  I am reproducing their results here is so I will have more credibility in subsequent posts when I criticize some of their claims about the significance of those results.  The 2009 PNAS paper by Vermeer and Rahmstorf will be referred to as VR2009 for the rest of this post.

Note that in 2007 Rahmstorf proposed the following model linking global temperature and sea level rise rate…

Rahmstorf’s model had some very serious shortcomings, so VR2009 proposed an improvement by adding another term…


H is the sea level

T is the temperature

T0 is a constant “equilibrium temperature”

t is the time

a and b are constants

The VR2009 idea is to take the measured sea level (H) and extract its time derivative (dH/dt, the sea level rise rate) and the measured temperature and extract its time derivative (dT/dt, the rate of temperature change) and insert them into equation 2.  Then the values of a, b and T0 can  be found by adjusting their values to minimize the difference between both sides of the equation. 

VR2009 claims that the model in equation 2 is superior to the equation 1 model.  They explicitly state that their ”analysis shows that the dual model [equation 2] is the preferred model, suggesting also here that the second term is physically meaningful.”  In subsequent posts I will show that this claim cannot be backed up.

The data

There are a variety of sets of sea level and temperature data.  VR2009 used the commonly quoted GISS global temperature and Church and White’s sea level data  .  Church and White’s data were further modified with Chao’s correction for artificial reservoir storage.  The temperature data and sea level data are quite noisy and their time derivatives are even noisier.  The only possible hope of using them to solve for a, b and T0 is to apply a smoothing algorithm.  The VR2009 Supporting Information states that they

“smoothed with a 15-year smoothing period (embedding dimension) using Singular Spectrum Analysis, SSA,  in order to reduce the extraneous impact of mostly short-term natural variability, in which we are not interested…The precise choice of smoothing method is not critical to our results.” 

 I have chosen the much simpler smoothing technique of applying a gaussian filter, with some minor variations.


Sea level

Figure 1 shows the Church and White sea level data overlaid with its derivative, the sea level rise rate, dH/dt.    The point of the following plot is the extremely noisy nature of the sea level rise rate data. 

Figure 1. Sea level, H, and sea level rise rate, dH/dt

Vermeer and Rahmstorf added Chao’s reservoir correction from the sea level.  This correction is supposed to represent the fraction of water that has been impounded on land in artificial reservoirs which would otherwise be adding to the sea level.  I don’t think this correction necessarily makes the sea level data any more accurate because it neglects the counteracting effect of aquifer depletion.  However, I include it  in order to reproduce Vermeer’s and Rahmstorf’s results.   The Chao data, figure 2, was acquired by digitizing the plot from Chao’s paper.

figure 2. Chao reservoir correction

Figure 3 shows the annual Church and White sea level data with the Chao reservoir correction added (blue) as well as a 15 year gaussian smoothed version (red) that is used in my implementation of the Vermeer and Rahmstorf model.

Figure 3. Church and White monthly sea level with Chao reservoir correction (blue), and 15 year FWHM Gaussian smoothed version (red).

Figure 4 is an animation of the sea level smoothing process and the resulting sea level rise rate.  I have included it to make clear the profound effect that smoothing has on the sea level rise rate.  The top plot shows the Church  and White sea level data and its smoothed version as smoothing increases from a 1 year FWHM gaussian filter to a 15 year FWHM gaussian filter.  The middle plot shows the time derivative (sea level rise rate, dH/dt) along with its smoothed version.  the bottom plot shows the same smoothed sea level rise rate data as the middle plot, but expanded for easier examination.

Notice that the sea level rise rate in the middle plot (thin white plot) looks nearly random.  Distinct features appear as smoothing (red plot) starts out but get ironed away as smoothing increases. (As an aside, here is Simon Holgate’s version of sea level rise rate from his 2007 Geophysical Research Letters paper.  Holgate sea level data was included in 2004 IPCC assessment.  Notice that in Holgate’s judgement short period features, on the order of ten years, are significant.)   At what smoothing level are real sea level rise rate features being removed? I don’t know.  I question whether it makes sense to apply this extreme level of smoothing, but did so to mimic Vermeer’s and Rahmstorf’s 15 year smoothing period.

Figure 4. Smoothing animation.


Like Vermeer and Rahmstorf, I used the GISS Global Annual Mean Surface Air Temperature Change.  I smoothed the data with a 20 year gaussian smooth for 1880 to 1894, and a 15 guassian smooth from 1895 to 2000.

Figure 5. GISS global temperature and 15 year gaussian smooth.

Turning the crank

Inserting the time derivative for H (dH/dt) into the left side of equation 1, and T and its time derivative (dT/dt) into the right side of equation 1, then finding the values of a, b,and To that yield the best fit gives the following results:

Vermeer and Rahmstorf found

a = 5.6 ± 0.5 mm/year/K

b= -49 ± 10 mm/K

To = -0.41 ± 0.03 K

I found

a = 5.6  mm/year/K

b= -52 mm/K

To = -0.42 K

So, my results were a very close match and well within Vermeer’s and Rahmstorf’s quoted uncertainties.

Figure 6 (click to enlarge) is a side by side comparison of Vermeer’s and Rahmstorf’s sea level data (left) and my sea level data (right).  Vermeer’s and Rahmstorf’s data is a copy of their plot from the bottom of their figure 3.  They say the red line is “Observations based rate of sea-level rise (with tectonic and reservoir effects removed).  The gray line is their result for equation 1, above, for their best values of a and To applied to their smoothed temperature data.   The blue line is their result for equation 2, above, for their best values of a, b,and To applied to their smoothed temperature data. 

My results use the same color scheme.  The red line is the Church and White sea level data (which already has tectonic effects removed) with the reservoir correction from Chao added.  The gray line is my result for equation 1 and the blue line is my result for equation 2.

Figure 6. Sea level data with reservoir correction (red) and the model results. Vermeer and Rahmstorf (left). ClimateSanity (right).

Figure 7 (click to enlarge) is a side by side comparison of Vermeer’s and Rahmstorf’s sea level rise rate data (left) and my sea level rise rate data (right).  Here the red curve is the input sea level rise rate with smoothing applied .  The gray and blue curves are sea level rise rates that result from equations 1 and 2, respectively.  VR2009 data is on the left, ClimateSanity data is on the right.

Figure 7. Sea level rise rate data with reservoir correction (red) and the model results. Vermeer and Rahmstorf (left). ClimateSanity (right).

Some observations.

The red curve on the right side of figure 7 is the same same as the red curve that forms in the middle and bottom plots of the animation in figure 4.  It is my version of the smoothed Church and White sea level rise rate with the Chao reservoir correction, dH/dt.  Why is it so different from VR2009’s version of the same data (the red curve on the left side of figure 7)?  Simply put, it is because of the degree of smoothing  and the smoothing technique. 

As I pointed out above, I have my doubts about the appropriateness of applying this heavy smoothing.  Perhaps it might have been just as good to fit the sea level rise data to a polynomial.  Figure 8, below, shows the sea level rise rate raw data, a 6th order polynomial fit, VR2009’s smoothed version, and my smoothed version.

Figure 8. Church and White sea level rise rate with Chao correction and various smoothing functions.

Here is a minor curiosity that leads me to doubt sea level smoothing of VR2009.  Look closely at their smoothed sea level rise rate from 1965 to 1985 and the unsmoothed Church and White data with the Chao reservoir correction, highlighted in figure 9.  This is a 20 year period where the VR2009 smoothing algorithm falls flat on its face.  I am not claiming that my smoothing technique is any better.  I am just pointing out the difficulty and uncertainty in smoothing data that is this noisy.

Figure 9.


My goal with this post was to demonstrate the I can closely reproduce VR2009’s results in terms of a, b and T0 in equation 2.  I have succeeded in that.  I have shown VR2009 is correct when they say “The precise choice of smoothing method is not critical to our results.”   My ability to reproduce their results will lend some credibility to my criticisms that I will express in a subsequent post after laying more foundational work.  I will show that their analysis which which they claim shows “that the dual model [equation 2] is the preferred model” is just plain wrong.  Rahmstorf’s first model (equation 1) was bogus and the “preferred” Vermeer and Rahmstorf model (equation 2) is no better.


My experience with Rahmstorf’s non-linear trend line

July 20, 2009

One of the original impetuses for me to start blogging was my experience with Stefan Rahmstorf concerning his 2007 paper “A Semi-Empirical Approach to Projecting Future Sea-Level Rise” (Science, 315, 2007).  I posted a several part critique on my old blogspot site, which I later ported over to this wordpress site. 

But this was only part of the story.  I  have decided to tell the rest of the story after reading “The Secret of the Rahmstorf ‘Non-Linear Trend Line’” at Steve McIntyre’s Climate Audit

Rahmstorf’s sea-level rise paper was based on plotting 120 years worth of sea -level rise rates vs each year’s corresponding global temperature.  Since both of these sets of data are quite noisy, Rahmstorf said ” Both temperature and sea-level curves were smoothed by computing nonlinear trend lines with an embedding period of 15 years.”

Rahmstorf referenced “New Tools for Analyzing Time Series Relationships and Trends” by Moore, et. al. (Eos, 86, 2005) for his nonlinear trend line smoothing technique.  This short paper refers to a variety of techniques for handling time series, including varieties of wavelet analysis and spectrum analysis.  The Moore paper invested several paragraphs on the use of Monte Carlo Single Spectrum Analysis for finding nonlinear trends in sea level and sea temperature, with the reader referred to a variety of  other papers to get the details. 

I waded hip deep into these papers  to get a handle on this new “nonlinear trend line” technique that led Rahmstorf to his startling projection of a huge sea level rise over this century.  I shouldn’t have wasted my time.  I found that I could essentially reproduce his results in an Excel spreadsheet by simply smoothing the original sea-level and temperature data with a 15 year FWHM Gaussian filter. 

Here is Rahmstorf’s sea-level rise rate vs. temperature after his nonlinear trend line smoothing and 5 year binning, followed by my sea-level rise rate vs. temperature after my 15 year FWHM Gaussian smoothing and 5 year binning, and finally, my version of the data without binning.

Rahmstorf's sea level rise vs T

Moriarty's sea level rise vs T binned

Moriarty's sea level rise vs T not binned

Rahmstorf binned his 120 data points into 24 bins containing 5 points each.   The binning was not needed to  remove noise that obfuscated his salient point – the data had already been smoothed through his non-linear trend line technique.     The binning could not be justified by claiming that it somehow made the plot of sea level rise rate vs temperature easier to read.  It actually reduced the amount of information to the reader by removing obvious real structure in the data.

I believe that Rahmstorf deliberately presented his data in a way calculated to deceive.  These are harsh words, and I say them with regret.

The only plausible reason that I can come up with for binning the 120 data points into 24 bins is because the resulting 24 points looked like they could conceivably be fit to a line without failing the laugh test.  Seeing the original 120 smoothed data points made it perfectly clear that there was not a linear relationship between the sea-level rise rate  and the temperature.  The full set of 120 data points also make it clear that when the temperature remains constant the sea level rise rate drops, in direct contradiction of one of Rahmstorf’s own working assumptions.

It turned out that Rahmstorf’s startling conclusion about extreme sea-level rise had nothing to do with any new sophisticated data analysis techniques for deriving nonlinear trend lines.  I got the same results as him using a simple spreadsheet.  Rather, his startling results came from his bogus interpretation.  Specifically, here are the three problems I identified:

1) The assumption that the time required to arrive at the new equilibrium is “on the order or millennia” is not borne out by the data. More…

2) Sea level rise rate vs. temperature is displayed in a way that erroneously implies that it is well fit to a line.  More…

3) Rahmstorf extrapolates out more than five times the measured temperature domain. More…

Rahmstorf’s code and peer review

In the midst of my wandering through a mathematical labyrinth to reproduce Rahmstorf’s results, before my simple excel spreadsheet approach, I asked Rahmstorf several questions via email.  Amazingly, he offered to send me his code, to which I happily accepted.  Here is what he said when he sent it (emphasis added by me):

From: Stefan Rahmstorf [mailto:rahmstorf@xxxxxxxxxxxx.xx]
Sent: Monday, August 20, 2007 13:20
To: Moriarty, Tom
Subject: Science paper

Dear Tom, see attached. Please report any issues you encounter, you are the first outside person to test this code.

Cheers, Stefan

Stefan Rahmstorf

So, the punchline is that although his data and results had been published a half a year before in the journal Science,  the highly regarded, unassailable, peer reviewed pinnacle of scientific research , I was “the first outside person to test his code.”

Again, I offer this harsh criticsm with regret, because Rahmstorf was, after all, kind enough to send me his code.