Monday, April 13, 2009

Published Comments on Rahmstorf 2007

After Rahmstorf 2007 (R07) was published in science there were two comments published which were each critical of the results. One comment was by Torben Schmith, Søren Johansen, and Peter Thejll (Schmith et al). The other was by Simon Holgate, Svetlana Jevrejeva, Philip Woodworth, and Simon Brewer. (Holgate et al). Dr. Rahmstorf wrote a reply (RR07) to both of these comments. A year later he followed up with a technical correction to the reply. I will touch briefly on Schmith et al. But focus on Holgate et al. In this post I will show that the response to Holgate et al. was not nearly as robust as presented. My conclusion is the same as Holgate et al. "...we do not agree that simplistic projections of the nature presented in [R07] substantially contribute to our understanding of the uncertainties in the nonlinear relationships of the climate system."

Schmith essentially said that due to the fact that there is a trend in both the temperature and sea level series that it violated the "basic assumptions of the statistical methods used." Remarkably they didn't comment on the smoothing. Even more importantly they ignored the fact that the residuals show a very high level of autocorrelation further showing that the model is mis-specified. Even so they concluded that the likelihood of the model was overstated. In RR07 Dr. Rahmstorf attempted to show that even with the trends removed a good fit remained, but this hardly addresses the main points which stand. The rest of the this post will show that the point is largely moot anyway.

The point that Holgate et al. made was simple. R07 didn't test to see if it could predict withheld data by modeling on the rest of the data. There are two problems with their presentation that confused the point when Dr. Rahmstorf responded. First they never duplicated his original result, instead using an SSA algorithm of their own devising. Second they used a different method to build the models. Using their own methods they were unable to predict the second half of the data set using the first half, and of course likewise in reverse. It is true that Dr. Rahmstorf didn't publish his code until he wrote RR07, but it seems like the should have sent him an email so they could start from the same place.

(I also note that at the beginning of his response to comments Dr. Rahmstorf said he was making "the computer code used in the analysis available for use by other researchers." What he neglects to mention is that the computer code he supplied was essentially useless without ssatrend.m. And nothing in his source code indicates where you could find that. The rest of his code is just simple regression and plotting, not much use to other researchers.)

Of course in RR07 Dr. Rahmstorf went back to his original method and reported that he could predict the second half using the first half of the data, and likewise the first half using the second half. In making his response he failed to report certain important points, and more importantly he made a significant error. A year later in October 2008 Science published his "technical correction" for that error, but as I will show that technical correction also fails.

The question is does the first half of the data predict the second half. Or, from my point of view, put more simply would a model built from the first half of the data be similar to a model built from the second half.

To evaluate this I regressed a model using the first twelve of the twenty-four five year bins. In this case the intercept is the same as the full model, but the coefficient is .44 versus .34 in the model built from the full data set. (RR07 reports .42, I'm not sure why the difference but it doesn't really matter. Also it incorrectly reports this as .42mm/year/degree but it is actually .42 cm/year/degree.) Using the second twelve of the five year bins the coefficient is only .24. (RR07 doesn't report this.) So simply put a model based on the first part of the period would have predicted nearly twice the sensitivity to temperature change that was experienced in the second set. RR07 shows graphs that it claims shows that it is making useful predictions, but based on this analysis it seems like a pretty poor match to me.

I've thrown some R code into the bottom of the source file so that the reader can project 2100 sea level with the different coefficients. But suffice it to say it makes a huge difference whether the coefficient is .42 or .24.

In any event I wonder whether the reader of these two posts has noticed the error in this methodology? The problem is that prior to being put into the bins the data had already been smoothed by the ssatrend algorithm with a 15 year window. This means that some of the original data from the second period is actually influencing the smoothed data in the first period. I don't know how Dr. Rahmstorf discovered this but as I said in October 2008 he published a technical correction. In the technical correction he noted the problem;

"This is correct, but it was illustrated by an incorrect figure (Fig. 1),in which the first half of the smoothed sea-level curve (1882 to 1941) was used to predict the sea level for 1942 to 2001. Because the smoothing procedure used a 15-year time window, the smoothed sea-level curve up to 1941 effectively contains sea-level information up to 1948. When this error is corrected and only annual sea-level measurements from 1882 to 1941 are used, the obtained fitgives a sea-level slope of 0.35 mm/year per °C"

(.35 mm/year/°C should, of course, be .35cm/year/°C)

This, he helpfully pointed out, was almost equal to the full trend of the original model even more strongly demonstrating how well a model from the first period predicted the rest of the data.

He failed to note three things. First in the response he said that he trained the data using the period 1880-1940. In the technical correction this was deftly changed to 1882-1941. It turns out the result is highly sensitive to this choice. Second based on my work I have determined that he changed the window size to 10 years from 15. Third using this exact technique the coefficient in the second period is only .22 this is unreported in the technical correction and is still much different than the trend from the first period.

In the data sets supplied by RR07 the sea level data begins in 1870, but the temperature data begins in 1880. As I've said before it isn't clear why he chose to bin the data at all, since it doesn't make any difference to the result, but the way that he did the binning in R07 starts the analysis in 1882. You can only find this out from the code, as R07 says it is using the period 1880-2001. I note that it doesn't change the results of the initial analysis.

In RR07 he says; "...but using only the first half of the data set (1880 to 1940) for deriving the statistical fit." RR07 supplied the code for R07 but not the code for the analysis in the response. This doesn't really matter because as I have duplicated he got the results in RR07 by using the approach of just regressing the first twelve bins. (This effectively started in 1882)

I have duplicated the .35 figure reported in the technical correction. To do this you have to use exactly 1882-1941. You also have to reduce the window to 10 years from 15, this is unreported in the technical correction. The result is non-robust to changes in either the exact year range or the window. If you move the window one year back this lowers the first period coefficient to .26. If you move it one year forward it raises the coefficient to .41. Combining the fact that the first and second periods don't match, to the fact that trivial changes in date ranges and window selections make large changes in the results shows that this model is not well specified. It is certainly not useful for updating sea level predictions beyond the results of AR4.

As a closing note on this analysis I want to say that it may not have been intentional, but R07, RR07, and the technical note were not nearly transparent enough. Instead it seems that they were written to make a point. It also points out, once again, the need for complete code disclosure.

Code for this analysis can be found here.

Saturday, April 11, 2009

Duplicating Rahmstorf 2007

I had quite an interesting time replicating the calculations from "A Semi-Empirical Approach to Projecting Future Sea-Level Rise" (Science 1/19/2007 R07) by Stefan Rahmstorf. The general concept in the paper is very simple. He estimates a linear equation for the rate of sea level rise based on global temperature from measured data. The difficulty came in because the paper relies on an algorithm that is not generally available. To demonstrate this it seems that although there were two critical comments on his paper published in Science neither replicated the original calculations of R07. It only became possible to replicate the results when Dr. Rahmstorf published his code as part of his response to those comments. Even then you had to do a little hunting. In this post I will show that I have duplicated those results, so that my subsequent comments make sense.

R07 is quite simple in concept. His theory is that the rate of sea level rise has a linear correlation with surface temperature. The idea is that you start with a situation where sea level is stable. Then you raise the temperature causing sea level to rise until is reaches a new equilibrium. I won't explain the whole paper, but this makes perfect sense. There are, of course, multiple inputs to sea level rise, but pretty much all of them should respond to a change in surface temperature, the question is how fast and how much.

To estimate this linear equation he uses sea level data from Church and White (IPCC), as well as the GISS global temperature record. But this is where things get a little tricky. He doesn't do a simple regression on the temperature and sea level because this data is "noisy." Instead runs both the temperature and the sea level data through a process to separate the trend from the noise. This isn't explained in the main text of the paper, but it is referred to in the text below figures 2 and 3. "Both temperature andsea-level curves were smoothed by computing nonlinear trend lines, with an embedding period of 15 years (14)." He uses the word "smoothed" but this is not a typical smoothing algorithm. Use of a typical polynomial fit will not get as "good" a result, as I discovered in early attempts.

I might have discovered this more quickly if I had read all the way to the bottom of Dr. Rahmstorf's reply to comments where a link to the code and data was published. The code and data weren't linked from the original paper which is where I started. Both published comments also didn't have the benefit of this code and data which is clear from reading them, and I guess after that everyone called it a day. But I will get to this stuff in a later post.

Reference 14 from the paper is "J. C. Moore, A. Grinsted, S. Jevrejeva, Eos 86, 226 (2005)." This paper is titled "New Tools for Analyzing Time Series Relationships and Trends." It is actually a short review article of several new mathematical techniques. At first glance it certainly isn't immediately apparent that this tells us how the "smoothing" was done. But there is a section titled "Nonlinear Trends in Sea Level and Temperature." This section refers to the use of Single Spectrum Analysis (SSA) to extract a nonlinear trend. It refers to Ghil et. al. 2002. So this is at best an indirect reference.

At this point I still hadn't discovered the code at the bottom of the reply, but I did have a lead. So I looked into SSA and discovered software at UCLA that was available. To make a long story short SSA is necessary but not sufficient to duplicate the R07 results. In fact further searches, and eventually the code led me to the fact that R07 relied on a matlab function called ssatrend.m. Dr. Rahmstorf on Real Climate indicated that this was available from Aslak Grinsted.

I wrote Dr. Grinsted who wrote me back very promptly, and sent me the source code to ssatrend.m. He also commented that he had no idea how Dr. Rahmstorf had gotten a copy of it, and that he had never meant for it to be distributed. I think that he was just concerned about it being unsupported. My own view is that it is pretty strange to use some random piece of code in a published paper without making the code your own. Especially where, as in this case, the strength of ther result depends on using this specific algorithm. It isn't available from any of the usual Matlab repositories which is the first place I looked.

In any event with a little effort I found an SSA algorithm written in R on source forge. I checked its results agains the UCLA program and they were identical so I knew the foundation was good. Then with some effort, I translated ssatrend.m into R so that I could run it on an open platform.

Finally I translated the code supplied in the reply written by Dr. Rahmstorf into R so that I could at see if I was starting from the same place. This was successful.

Using the supplied input files for sea level and temperature I get a correlation coefficient (R) of .88 exactly as reported in R07. In addition the following graphs are identical to figures two and three (top) of R07. The code is here.

At the end of the day the use of this unpublished algorithm made this much harder than it should have been as the underlying concept is very simple. I have no idea how a reviewer could have evaluated whether the use of this algorithm made sense. Having said that I don't really see any problem with it, and I look forward to trying ssatrend with other data.

But now that I can get the exact results from the paper, I have a couple of comments that go beyond what has already been written in Science.

New consensus on sea level rise?

Recently because of some posts on stoat I have become interested in how the consensus is moving on projected sea level rise. There are clearly some scientists who now believe that 1M above 1990 by 2100 is likely. These include Stefan Rahmstorf, and Aslak Grinsted who have published papers based on empirical analysis to come to this conclusion. There are a number of issues with the paper by Dr. Rahmstorf some of which were published as comments in Science. I will post some additional comments which follow the publishing of the code and data. The paper by Dr. Grinsted is more interesting to me, but because of the time periods he uses to train and test his model much of the input data and assumptions are necessarily pretty fuzzy.

In AR4 the consensus estimates for sea level rise are on page 821 figure 10.33. Table 10.7 on the previous page breaks down the components. The range for the various scenarios including the error bars is roughly .2 meters to .6 meters. Many people seem to feel that they just threw up their hands at faster ice sheet discharge but table 10.33 includes figures for that at the bottom. They declined to add these into the projections because they couldn't assess the likelihood of this happening. It is important to note that these would have only added .09 to .17 meters to the high end of the range bringing the top to about .8 meters. (It would have slightly lowered the bottom of the range as well.)

So I wonder how far the consensus has actually changed in the last couple of years. Or are there simply some scientists who believe the consensus is low? Should the best estimate now be considered 1M? Or should we stay with the IPCC conclusions? Of course I don't attend the conferences with the types of experts who are called on to determine the consensus view. But I note that it doesn't seem to me that any of the lead authors of Chapter 10 are authors of these recent studies. (I could easily be wrong as cross checking that type of thing is difficult.)

Anyway I am going to write a couple of posts looking at these empirical papers. I am also looking at building an empirical model of my own to see if it will improve on the results of R07.

Saturday, March 14, 2009

Troposhpere Temps from Satellites and Surface Temps

In an earlier post I took at look at MM07 versus S09. One of the interesting results was that the rate of post secondary education (PSE) in a country was inversely proportional to the rate of temperature increase in that country. Taking a deeper look I saw that if you take the top quartile of grid points based on PSE from the analysis that the surface warmed more slowly than the troposphere, while in the rest of the grid points the surface warmed much faster than the troposphere. (This is using HadCRUT for the surface and RSS for the troposphere.) I suppose you could shrug this off to coincidence except that according to the Model E data supplied by Dr. Schmidt the troposhpere is supposed to be generally warming faster than the surface everywhere.

Over the 440 grid cells of the analysis the Model E troposphere warmed faster than the surface (.16 degrees per decade versus .14) This contrasts with the measured data from HadCRUT and RSS where the surface warmed faster than the troposphere (.27 versus .23).

Looking at the segment with high PSE we can start with the US. Now I know there has been a lot of sniping about the US temperature network, but I'm guessing that it is really pretty good. Out of the 440 grid points 52 are in the US. For the these grid points the troposhpere is warming faster than the surface (.26 versus .24). There are 85 grid points in the top quartile outside the US. For those points the troposphere is also warming faster than the surface (.21 versus .19).

Now the truth is that this data is very convenient for me to look at because it was already layed out by others, and I haven't looked at any other time periods to confirm that this isn't some kind of fluke.

But I think it is pretty interesting that in the countries that probably have the best surface temperature networks the actual measurements are in line with the theory as proposed by the results of Model E. The conclusion would be that perhaps climate scientists ought to be focused on troposhperic temperatures as measured by satellite, and reduce their dependence on ground based measurements. Switching to satellite measurement seems to be happening elsewhere with a good example being sea level rise.

I should add that they ought to be noticing this type of agreement with models and be pleased with the vindication, but I don't sense that they are. I have a theory as to why.

When the satellite temperatures were first introduced by UAH they were used by climate skeptics to show that there was no warming. In addition Drs. Christy and Spencer aligned themselves to some degree with the skeptics camp. Even though subsequent events have corrected the satellite trends and there now is an independent satelitte measurement from RSS this seems to have put satellite measurements out of favor in this area. This is particularly true for the UAH data.

In fact you can get a hint of this from S09. At one point Dr. Schmidt comments that the differing regression results he got by using RSS versus UAH might be caused by the "higher trend" in RSS. In fact this is uncited and he provides no results to back this up. I think he just assumed it was true, because a quick test of the trends show that for this set of grid cells over this period RSS and UAH have identical average trends. The point being that Dr. Schmidt believed so strongly that of course UAH would show less warming than RSS that he didn't even test the conclusion.

I think it would be quite interesting if the climate modeling community would look at their results relative to the troposhperic measurements from RSS/UAH and deemphasize the surface network. There is plenty of warming in the satellite measurements, and they may be a whole lot more accurate

Follow Up on ERA-Interim from ECMWF

After looking at this post, Ryan Maue suggested that I do a further analysis of humidity trends using the ERA Interim data set of ECMWF. This is the most up to date information, although it covers a much shorter period than the ERA 40 data I used in the earlier analysis. The results are that over the 19 year period from 1989 to 2007 there are no significant trends in specific humidity (q) in this data set. Current theory would say that q should be increasing. Thus between the ERA 40 data, and the ERA Interim data there is no confirmation of the theory, and in fact many of the trends are negative, particularly over the NH mid latitudes, but not significant.

The data for this study is from the ERA interim data set downloaded from the ECMWF servers.

The results can be found here.

Wednesday, March 11, 2009

Humidity Trends from ECMWF

There was an interesting post on Climate Audit about a paper covering trends in atmospheric humidity. The substance of the post was discussing why a paper by Gareth Paltridge and others was rejected by the journal of climate. But I was quite interested to learn that there was long term trend data on humidity.

As most people who follow the issue understand a doubling of CO2 by itself should increase atmospheric temperature about 1 degree C at equilibrium. The big question is what feedback is caused by this temperature increase. The largest feedback is caused by an increase in water vapor as the temperature increases. The theory is that relative humidity (r) should remain constant, which would mean that as temperature increases specific humidity (q) would increase. This means that the total water vapor of the atmosphere would increase causing an increased forcing.

Over time I have seen a few papers which have confirmed that relative humidity is indeed staying constant. But these papers have had a lot of caveats and have covered limited areas and time frames. So I was interested in the paper by Dr. Paltridge which used a re-analysis data set called NCEP covering a 35 year period from 1973 to 2007.

I don't have access to the paper but the summary is that in this data set r decreased over the period in some relevant regions. This would be counter to theory, and interesting.

Ryan Maue, who is a student at Florida State made several comments and here on CA pointing out issues with the Paltridge paper. Essentially the objection is that prior to 1979 the data is based on a small radiosonde network, and that the subsequent data is based on a combination of satellites and radiosondes. He felt that at the minimum the results should be compared to other re-analysis data sets. He chose ECMWF as the best. He even commented that he would take a look himself.

So I thought I would go ahead and see if I could figure out how to download the ECMWF data and do a quick analysis. It turns out to be quite possible.

The ECMWF data covers parts of 46 years, but only 44 complete years from 1958 to 2001. So I downloaded the r and q data for all grid locations for those years. The results have something for everyone I suppose.

For completeness I started out looking at the trends for the entire period and the entire globe. In this case q is only negative at the very highest altitudes above 100hPa. Below 500hPa q is positive. R on the other hand is negative above 400hPa, and positive below 925hPa with the altitudes in between not being significant. Again according to theory the theory the trend in r is supposed to be zero and the trend in q is supposed to be positive.

But I get the impression that the global figure over this time period is not the most interesting. As Mr. Maue points out most of the radiosondes are in a small band in the Northern Hemisphere. So the trends that cover that area are called out both by Dr. Paltridge and By Mr. Maue in a subsequent post.

Looking at the NH results for the entire period q is significantly negative all the way down to 700mb. It only becomes significantly positive at 925hPa. In a result I don't completely understand r is negative above 400hPa and insignificant below that. I would have thought that in a warming atmosphere that if q was negative r would have to be negative as well. In any event the negative q trend over the NH which is where the majority of the real measurements would have been made in this time period seems to be different than theory and in line with the results from NCEP. Note that I tried two definitions of the Northern mid latitudes with no change in results.

In the SH the trends are largely positive for both q and insignificant for r which would be in line with theory. And this is true for the entire mid latitude and tropic region, which has negative r only for the altitudes above 400hPa.

In summary then over the 44 year period the area that shows negative q seems to be the mid latitudes of the Northern Hemisphere. Since the areas outside of the measurement regions are computed using climate models I'm not sure of the relationship of the "real" data to areas where there were no radiosondes.

I took a separate look at the "post satellite" period. Unfortunately this is a very short time in this data set since it ends in 2001 unlike the NCEP data which goes through 2007. The trends were not particularly significant over this period.

The R code for this analysis can be found here.

The data for this post is from ERA-40 and was graciously supplied by the ECMWF data server.

Tuesday, March 3, 2009

IPCC Sea Ice Forecasts

There has been a lot of discussion in recent years of the large drop in Arctic sea ice. People who are particularly concerned about CO2 induced global warming have pointed to this drop as clear current evidence that the effects of warming are accelerating. People who are generally skeptical about CO2 induced global warming have brought up the fact that Antarctic sea ice levels appear to be increasing as a counter argument. It has always seemed logical to me that global sea ice level was probably the most complete measure in this area, but I had read several informal documents that suggested that the Arctic was the true signature test for CO2 induced global warming.

In particular earlier this year a blogger wrote a post noting that global sea ice was at the same level as 1979. This generated a response from "Cryosphere Today." The response was that while the information was correct it was referring to the global trend and not on the Arctic trend. They point out that the reduced sea ice in the Arctic was being compensated for by increased ice in the Antarctic. The response includes the following. "In the context of climate change, GLOBAL sea ice area may not be the most relevant indicator. Almost all global climate models project a decrease in the Northern Hemisphere sea ice area over the next several decades under increasing greenhouse gas scenarios. But, the same model responses of the Southern Hemisphere sea ice are less certain."

I had heard statements like this before, but I had never looked to see what the basis of the statements might be. So I went to the AR4 WG1 report. On page 770 you can find chapter 10.3.3.1 "Changes in Sea Ice Cover." In that section you can find the following. "In 20th and 21st-century simulations, antarctic sea ice cover is projected to decrease more slowly than in the Arctic (Figures 10.13c,d and 10.14)." I have copied figure 10.13 below.

The top two projections are for the Arctic, and the bottom two are for the Antarctic. The black line represents the ensemble mean for the various scenarios and models. What is obvious from looking at this picture is that while the Arctic sea ice is projected to decrease more rapidly, Antarctic sea ice is projected to decrease as well. Not only that but the decrease should have been occurring fairly steadily in the last few decades, and increasing at this point.

In that section, which it is true is quite short, I found no reference to higher uncertainty regarding Antarctic sea ice. The final paragraph of the section just refers to the general uncertainty of projections including the amount of climate change in the polar regions in general.

So what evidence did CT have for their assertion? They link to a single 2005 study discussing a simulation where increased snow fall in the Antarctic tends to preserve sea ice.

So the consensus view as established by the IPCC is that while the Arctic is expected to decline more rapidly we should see declining sea ice at both poles. This means that increased Antarctic sea ice over the last 30 years is not in line with these projections. Does this prove that increasing CO2 does not cause global warming? Not even a little. But it does mean that looking at global sea ice level is relevant, and that focusing exclusively on Arctic sea ice might reasonably be considered cherry picking.