Nature Unleashes a Flood … of Bad Science.

Guest Post by Willis Eschenbach

Recently, Nature Magazine published a paywalled paper called “Human contribution to more-intense precipitation extremes” by Seung-Ki Min, Xuebin Zhang, Francis W. Zwiers & Gabriele C. Hegerl (hereinafter MZZH11) was published in Nature Magazine. The Supplementary Information is available here. The study makes a very strong claim to have shown that CO2 and other greenhouse gases are responsible for increasing extreme rainfall events, viz:

Here we show that human-induced increases in greenhouse gases have contributed to the observed intensification of heavy precipitation events found over approximately two-thirds of data-covered parts of Northern Hemisphere land areas.

Figure 1. Extreme 1-day rainfall. New Orleans, Katrina. Photo Source

There are two rainfall indices which are used in their analysis, called the RX1day and RX5day indices. The RX1day and RX5day indices give the maximum one-day precipitation and five-day precipitation for a given station for a given month. These individual station datasets (available here, free registration required) have been combined into a gridded dataset, called HADEX (Hadley Climate Extremes Dataset) . It is this gridded dataset that was used in the MZZH11 study.

So what’s wrong with the study? Just about everything. Let me peel the layers off it for you, one by one.

Other people have commented on a variety of problems with the study, including Roger Pielke Jr., Andy Revkin,  Judith Curry . But to begin with, I didn’t read them, I did what I always do. I went for the facts. I thrive on facts. I went to get the original data. For me, this is not the HADEX data, as that data has already been gridded. I went to the actual underlying data used to create the HADEX dataset, as cited above. Since they don’t provide a single datablock file with all of the areas (grrrr … pet peeve), I started by looking at the USA data.

And as is my habit, the first thing I do is just to look at the individual records. There are 2,661 stations in the USA database, of which some 731 contain some RX1day maximum one day rainfall data. However, as is usual with weather records of all kinds, many of these have missing data.  In addition, only 9% of the stations contain a significant trend at the 95% confidence level. Since with a 95% confidence interval (CI) we would expect 5% of the stations to exceed that in any random dataset, we’re only slightly above what would be expected in a random dataset. In addition, the number of stations available varies over time..

Now, let me repeat part of that, because it is important.

91% of the rainfall stations in the US do not show a significant trend in precipitation extremes, either up or down.

So overwhelmingly in the US there has been

No significant change in the extreme rainfall.

And as if that wasn’t enough …

Of the remaining 9% that have significant trends, 5% of the trends are probably from pure random variation.

So this means that

Only about 5% of the stations in the US show any significant change in rainfall extremes.

So when you see claims about changes in US precipitation extremes, bear in mind that they are talking about a situation where only ~ 5% of the US rainfall stations show a significant trend in extreme rainfall. The rest of the nation is not doing anything.

Now, having seen that, let’s compare that to the results shown in the study:

Figure 2. The main figure of the MZZH11 study, along with the original caption. This claims to show that the odds of extreme events have increased in the US.

Hmmmm …. so how did they get that result, when the trends of the individual station extreme precipitation show that some 95% of the stations aren’t doing anything out of the ordinary? Let me go over the stages step by step as they are laid out in the study. Then I’ll return to discuss the implications of each step.

1. The HADEX folks start with the individual records. Then, using a complex formula based on the distance and the angle from the center of the enclosing gridcell, they take a weighted station average of each month’s extreme 1-day rain values from all stations inside the gridcell. This converts the raw station data into the HADEX gridded station data.

2. Then in this study they convert each HADEX gridcell time series to a “Probability-Based Index” (PI) as follows:

Observed and simulated  annual extremes are converted to PI by fitting a separate generalized extreme value (GEV) distribution to each 49-year time series of annual extremes and replacing values with their corresponding percentiles on the fitted distribution. Model PI values are interpolated onto the HadEX grid to facilitate comparison with observations (see Methods Summary and Supplementary Information for details).

In other words, they separately fit a generalized three-parameter probability function each to gridcell time series, to get a probability distribution. The fitting is done iteratively, by repeatedly adjusting each parameter to find the best fit. Then they replace that extreme rainfall value (in millimetres per day) with the corresponding probability distribution value, which is between zero and 1.

They explain this curious transformation as follows:

Owing to the high spatial variability of precipitation and the sparseness of the observing network in many regions, estimates of area means of extreme precipitation may be uncertain; for example, for regions where the distribution of individual stations does not adequately sample the spatial variability of extreme values across the region. In order to reduce the effects of this source of uncertainty on area means, and to improve representativeness and inter-comparability, we standardized values at each grid-point before estimating large area averages by mapping extreme precipitation amounts onto a zero-to-one scale. The resulting ‘probability-based index’ (PI) equalizes the weighting given to grid-points in different locations and climatic regions in large area averages and facilitates comparison between observations and model simulations.

Hmmm … moving right along …

3. Next, they average the individual gridcells into “Northern Hemisphere”, “Northern Tropics”, etc.

4. Then the results from the models are obtained. Of course, models don’t have point observations, they already have gridcell averages. However, the model gridcells are not the same as the HADEX gridcells. So the model values have to be area-averaged onto the HADEX gridcells, and then the models averaged together.

5. Finally, they use a technique optimistically called “optimal fingerprinting”. As near as I can tell this method is unique to climate science. Here’s their description:

In this method, observed patterns are regressed onto multi-model simulated responses to external forcing (fingerprint patterns). The resulting best estimates and uncertainty ranges of the regression coefficients (or scaling factors) are analysed to determine whether the fingerprints are present in the observations. For detection, the estimated scaling factors should be positive and uncertainty ranges should exclude zero. If the uncertainty ranges also include unity, the model patterns are considered to be consistent with observations.

In other words, the “optimal fingerprint” method looks at the two distributions H0 and H1 (observational data and model results) and sees how far the distributions overlap. Here’s a graphical view of the process, from Bell, one of the developers of the technique.

Figure 2a. A graphical view of the “optimal fingerprint” technique.

As you can see, if the distributions are anything other than Gaussian (bell shaped), the method gives incorrect results. Or as Bell says (op. cit.) the optimal footprint model involves several crucial assumptions, viz:

•  It assumes the probability distribution of the model dataset and the actual dataset are Gaussian

•  It assumes the probability distribution of the model dataset and the actual have approximately the same width

While it is possible that the extreme rainfall datasets fit these criteria, until we are shown that they do fit them we don’t know if the analysis is valid. However, it seems extremely doubtful that the hemispheric averages of the probability based indexes will be normal. The MZZH11 folks haven’t thought through all of the consequences of their actions. They have fitted an extreme value distribution to standardize the gridcell time series.

This wouldn’t matter a bit, if they hadn’t then tried to use optimal fingerprinting. The problem is that the average of a PI of a number of extreme value distributions will be an extreme value distribution, not a Gaussian distribution. As you can see in Figure 2a above, for the “optimal fingerprint” method to work, the distributions have to be Gaussian. It’s not as though the method will work with other distributions but just give poorer results. Unless the data is Gaussian, the “optimal fingerprint” method is worse than useless … it is actively misleading.

It also seems doubtful that the two datasets have the same width. While I do not have access to their model dataset, you can see from Figure 1 that the distribution of the observations is wider, both regarding increases and decreases, than the distribution of the model results.

This seems extremely likely to disqualify the use of optimal fingerprinting in this particular case even by their own criteria. In either case, they need to show that the “optimal fingerprint” model is actually appropriate for this study. Or in the words of Bell, the normal distribution “should be verified for the particular choice of variables”. If they have done so there is no indication of that in the study.

I think that whole concept of using a selected group of GCMs for “optimal fingerprinting” is very shaky. While I have seen theoretical justifications for the procedure, I have not seen any indication that it has been tested against real data (not used on real data, but tested against a selected set of real data where the answer is known). The models are tuned to match the past. Because of that, if you remove any of the forcings, it’s almost a certainty that the model will not perform as well … duh, it’s a tuned model. And without knowing how or why the models are chosen, how can they say their results are solid?

OK, I said above that I would first describe the steps of their analysis. Those are the steps. Now let’s look at the implications of each step individually.

STEP ONE: We start with what underlies the very first step, which is the data. I didn’t have to look far to find that the data used to make the HADEX gridded dataset contains some really ugly errors. One station shows 48 years of August rains with a one-day maximum of 25 to 50 mm (one to two inches), and then has one August (1983) with one day when it is claimed to have rained 1016 mm (40 inches) … color me crazy, but I think that once again, as we have seen time after time, the very basic steps have been skipped. Quality doesn’t seem to be getting controlled. So … we have an unknown amount of uncertainty in the data simply due to bad individual data points. I haven’t done an analysis of how much, but a quick look revealed a dozen stations with that egregious an error in the 731 US datasets … no telling about the rest of the world.

The next data issue is “inhomogeneities” (sudden changes in volume or variability) in the data. In a Finnish study, 70% of the rainfall stations had inhomogeneities. While there are various mathematical methods used by the HADEX folks to “correct” for this, it introduces additional uncertainty into the data. I think it would be preferable to split the data at the point of the inhomogeneous change, and analyze each part as a separate station. Either way, we have an uncertainty of at least the difference in results of the two methods. In addition, the Norwegian study found that on average, the inhomogeneities tended to increase the apparent rainfall over time, introducing a spurious trend into the data.

In addition, extreme rainfall data is much harder to quality control than mean temperature data. For example, it doesn’t ever happen that the January temperature at a given station averages 40 degrees every January but one, when it averages 140 degrees. But extreme daily rainfall could easily change from 40 mm one January to an unusual rain of 140 mm. This makes for very difficult judgements as to whether a large daily reading is erroneous.

In addition, an extreme value is one single value, so if that value is incorrectly large it is not averaged out by valid data. It carries through, and is wrong for the day, the month, the year, and the decade.

Rainfall extreme data also suffers in the recording itself. If I have a weather station and I go away for the weekend, my maximum thermometer will record the maximum temperature of the two days I missed. But the rainfall gauge can only give me the average of the two days I missed … or I could record the two days as one with no rain on the other day. Either way … uncertainties.

Finally, up to somewhere around the seventies, the old rain gauges were not self emptying. This means that if the gauge were not manually emptied, it could not record an extreme rain. All of these problems with the collection of the extreme rainfall data means it is inherently less accurate than either mean or extreme temperature data.

So that’s the uncertainties in the data itself. Next we come to the first actual mathematical step, the averaging of the station data to make the HADEX gridcells. HADEX, curiously, uses the averaging method rejected by the MZZH11 folks. HADEX averages the actual rainfall extreme values, and did not create a probability-based index (PI) as in the MZZH11 study. I can make a cogent argument for either one, PI or raw data, for the average. But using a PI based average of a raw data average seems like an odd choice, which would result in unknown uncertainties. But I’m getting ahead of myself. Let me return to the gridding of the HADEX data.

Another problem increasing the uncertainty of the gridding is the extreme spatial and temporal variability of rainfall data. They are not well correlated, and as the underlying study for HADEX says (emphasis mine):

[56] The angular distance weighting (ADW) method of calculating grid point values from station data requires knowledge of the spatial correlation structure of the station data, i.e., a function that relates the magnitude of correlation to the distance between the stations. To obtain this we correlate time series for each station pairing within defined latitude bands and then average the correlations falling within each 100 km bin. To optimize computation only pairs of stations within 2000 km of each other are considered. We assume that at zero distance the correlation function is equal to one. This may not necessarily be the best assumption for the precipitation indices because of their noisy nature but it does provide a good compromise to give better gridded coverage.

Like most AGW claims, this seems reasonable on the surface. It means that stations closer to the gridbox center get weighted more than distant stations. It is based on the early observation by Hansen and Lebedeff in 1987 that year-to-year temperature changes were well correlated between nearby stations, and that correlation fell off with distance. In other words, if this year is hotter than last year in my town, it’s likely hotter than last year in a town 100 km. away. Here is their figure showing that relationship:

Figure 3. Correlation versus Inter-station Distance. Original caption says “Correlation coefficients between annual mean temperature changes for pairs of randomly selected stations having at least 50 common years in their records.”

Note that at close distances there is good correlation between annual temperature changes, and that at the latitude of the US (mostly the bottom graph in Figure 3) the correlation is greater than 50% out to around 1200 kilometres.

Being a generally suspicious type fellow, I wondered about their claim that changes in rainfall extremes could be calculated by assuming they follow the same distribution used for temperature changes. So I calculated the actual relationship between correlation and inter-station distance for the annual change in maximum one-day rainfall. Figure 4 shows that result. It is very different from temperature data, which has good correlation between nearby stations and drops off slowly with increasing distance. Extreme rainfall does not follow that pattern in the slightest.

Figure 4. Correlation of annual change in 1-day maximum rainfall versus the distance between the stations. Scatterplot shows all station pairs between all 340 mainland US stations which have at least 40 years of data per station. Red line is a 501 point Gaussian average of the data.

As you can see, there is only a slight relationship at small distances between extreme rainfall event correlation and distance between stations. There is an increase in correlation with decreasing distance as we saw with temperature, but it drops to zero very quickly. In addition, there are a significant number of negative correlations at all distances. In the temperature data shown in Figure 3, the decorrelation distance (the distance where the average correlation drops to 0.50) is on the order of 1200 km. The corresponding decorrelation distance for one-day extreme precipitation is only 40 km …

Thinking that the actual extreme values might correlate better than the annual change in the extreme values, I plotted that as well … it is almost indistinguishable from Figure 4. Either way, there is only a very short-range (less than 40 km) relation between distance and correlation for the RX1day data.

In summary, the method of weighting averages by angular distances used for gridding temperature records is supported by the Hansen/Lebedeff temperature data in Figure 3. On the other hand, the observations of extreme rainfall events in Figure 4 means that we cannot use same method for gridding of extreme rainfall data. It makes no sense, and reduces accuracy, to average data weighted by distance when the correlation doesn’t vary with anything but the shortest distances, and the standard deviation for the correlation is so large at all distances.

STEP 2: Next, they fit a generalized extreme value (GEV) probability distribution to each individual gridcell. I object very strongly to this procedure. The GEV distribution has three different parameters. Depending on how you set the three GEV dials, it will give you distributions ranging from a normal to an exponential to a Weibull distribution. Setting the dials differently for each gridcell introduces an astronomical amount of uncertainty into the results. If one gridcell is treated as a normal distribution, and the next gridcell is treated as an exponential distribution, how on earth are we supposed to compare them? I would throw out the paper based on this one problem alone.

If I decided to use their method, I would use a Zipf distribution rather than a GEV. The Zipf distribution is found in a wide range of this type of natural phenomena. One advantage of the Zipf distribution is that it only has one parameter, sigma. Well, two, but one is the size of the dataset N. Keeps you from overfitting. In addition, the idea of fitting a probability distribution to the angular-distance weighted average of raw extreme event data is … well … nuts. If you’re going to use a PI, you need to use it on the individual station records, not on some arbitrary average somewhere down the line.

STEP 3: Hemispheric and zonal averages. In addition to the easily calculable statistical error propagation in such averaging, we have the fact that in addition to statistical error each individual gridpoint has its own individual error. I don’t see any indication that they have dealt with this source of uncertainty.

STEP 4: Each model needs to have its results converted from the model grid to the HADEX grid. This, of course, gives a different amount of uncertainty to each of the HADEX gridboxes for each of the models. In addition, this uncertainty is different from the uncertainty of the corresponding observational gridbox …

There are some other model issues. The most important one is that they have not given any ex-ante criteria for selecting the models used. There are 24 models in the CMIP database that they could have used. Why did they pick those particular models? Why not divide the 24 models into 3 groups of 8 and see what difference it makes? How much uncertainty is introduced here? We don’t know … but it may be substantial.

STEP 5: Here we have the question of the uncertainties in the optimal fingerprinting. These uncertainties are said to have been established by Monte Carlo procedures … which makes me nervous. The generation of proper data for a Monte Carlo analysis is a very subtle and sophisticated art. As a result, the unsupported claim of a Monte Carlo analysis doesn’t mean much to me without a careful analysis of their “random” proxy data.

More importantly, the data does not appear to be suitable for “optimal fingerprinting” by their own criteria.

End result of the five steps?

While they have calculated the uncertainty of their final result and shown it in their graphs, they have not included most of the uncertainties I listed above. As a result, they have greatly underestimated the real uncertainty, and their results are highly questionable on that issue alone.

OVERALL CONCLUSIONS

1. They have neglected the uncertainties from:

•  the bad individual records in the original data

•  the homogenization of the original data

•  the averaging into gridcells

•  the incorrect assumption of increasing correlation with decreasing distance

•  the use of a 3 parameter fitted different probability function for each gridcell

•  the use of a PI average on top of a weighted raw data average

•  the use of non-Gaussian data for an “optimal fingerprint” analysis

•  the conversion of the model results to the HADEX grid

•  the selection of the models

As a result, we do not know if their findings are significant or not … but given the number of sources of uncertainty and the fact that their results were marginal to begin with, I would say no way. In any case, until those questions are addressed, the paper should not have been published, and the results cannot be relied upon.

2.  There are a number of major issues with the paper:

•  Someone needs to do some serious quality control on the data.

•  The use of the HADEX RX1day dataset should be suspended until the data is fixed.

•  The HADEX RX1day dataset also should not be used until gridcell averages can be properly recalculated without distance-weighting.

•  The use of a subset of models which are selected without any ex-ante criteria damages the credibility of the analysis

•  If a probability-based index is going to be used, it should be used on the raw data rather than on averaged data. Using it on grid-cell averages of raw data introduces spurious uncertainties.

•  If a probability-based index is going to be used, it needs to be applied uniformly across all gridcells rather than using different distributions a gridcell by gridcell basis.

•  No analysis is given to justify the use of “optimal fingerprinting” with non-Gaussian data.

3. Out of the 731 US stations with rainfall data, including Alaska, Hawaii and Puerto Rico, 91% showed no significant change in the extreme rainfall events, either up or down.

4. Of the 340 mainland US stations with 40 years or more of records, 92% showed no significant change in extreme rainfall in either direction.

As a result, I maintain that their results are contrary to the station records, that they have used inappropriate methods, and that they have greatly underestimated the total uncertainties of their results. Thus the conclusions of their paper are not supported by their arguments and methods, and are contradicted by the lack of any visible trend in the overwhelming majority of the station datasets. To date, they have not established their case.

My best regards to all, please use your indoor voices in discussions …

w.

[UPDATE] I’ve put the widely-cited paper by Allen and Tett about “optimal fingerprinting” online here.