10. Breakpoint adjustment and its effect on the scaling behaviour of anomalies

In my previous post (see here) I demonstrated that the noise or random temporal fluctuations in the temperature anomaly from weather station data were not consistent with white noise, but were instead characteristic of a one-dimensional time-dependent fractal with a 1/√ω power spectrum (where ω is the frequency of the oscillations in the noise signal). I then showed that these temperature fluctuations exhibit self-similarity with a fractal dimension of 0.25 (approximately). I also showed that one consequence of this is that the expected long term temperature fluctuations (i.e. those with timescales greater than 100 years) would be at least six times greater in magnitude (almost 1.0 °C) than most climate scientists probably realize. This analysis was done on the original raw unadulterated data (or at least what I assume is the original raw unadulterated data). In this post I shall apply the same technique to data that has been processed and adjusted by Berkeley Earth and look at the difference.

One of the most frequent claims made against climate science by sceptics is that the temperature data is manipulated. It is of course a claim that those climate scientists deny, but ever since the Climategate email scandal there has been a bad smell. It was, after all, that bad smell that prompted Richard Muller to set up Berkeley Earth in order to provide some independent verification of the results that the other groups (NOAA, GISS, HADCRUT) were producing.

That said, the claim that some climate scientists manipulate the data has never been completely rejected by climate scientists themselves. In fact they justify much of it on the grounds that the data is unreliable and needs to be “cleaned up”. Their rationale for these statistical interventions is this: some temperatures records are very long and therefore it is inevitable that something will go wrong with the data in that time. People will change, instruments will change (e.g. changes from Stevenson screens to electronic systems), procedures will change (e.g. time of data collection), and locations will change, so error checking needs to be put in place to compensate for any potential errors that may arise. This is particularly important since the primary purpose of most climate research is to determine the temperature change over time. This in turn requires the scientist to measure two temperatures, under identical conditions, at the same location, but a long time apart. It is therefore no good having very accurate measurements now if what you are comparing them against is uncertain.

In addition, data is often missing or incorrectly entered (e.g. in Fahrenheit rather than Celsius). The result of all this is that the climate scientists believe the data needs to be subject to a strict quality control during the analysis process, and two of the “data cleaning” techniques that they rely on to do this are homogenization and breakpoint adjustment.

Homogenization is a primarily used for two outcomes. The first is to use extrapolation or interpolation to fill in missing data points in the record. The second is to compare records from neighbouring stations and construct a regional average that can be used as a more reliable reference dataset.

Breakpoint adjustment is a data manipulation technique that first uses homogenization to construct an idealized reference for the region and then compares the actual station data to the reference. It then looks for points in the dataset of the target station where the difference between the target and the reference changes over a short time interval by an amount that is much larger than is normally seen. At this point a cut or break in placed in the data and data each side of the cut is shifted up or down. An example (taken from Berkeley Earth) is shown in Fig. 10.1 below.

Fig. 10.1: Anomaly data for Waiouru Airstrip.

In the above data for Waiouru Airstrip (Berkeley Earth ID - 172950) a breakpoint was detected by Berkeley Earth between October and November 1986, as shown below in Fig. 10.2. The red lines in Fig. 10.2 indicate the amount by which the data on the right (post October 1986) is deemed to be too low in value relative to the data before November 1986 when compared with a regional average.

Fig. 10.2: Difference in station anomaly and regional anomaly for Waiouru Airstrip.

This regional average appears to be constructed from the surrounding station data weighted by the square of their correlation coefficient with the target station (Waiouru Airstrip). There is only one problem with this, though; there is only one other station within 100 km with data spanning 1986 and that is at Ohakea 87 km away with only 234 months of data (see below). So how reliable is the weighting process? That is a major issue, particularly in the Southern Hemisphere where station density is much lower than in the USA or Europe.

 Fig. 10.3: Anomaly data for Ohakea.

Nevertheless, it can be seen that the breakpoint correction in Fig. 10.2 is an impressive 1.2 °C. However, this adjustment is not applied evenly. First, it equates to a different amount for each month (January - December) as shown in Fig. 10.4 (see below).

 Fig. 10.4: Total breakpoint correction for Waiouru Airstrip.

This effectively means that while the breakpoint is set at the same point for each of the twelve monthly datasets that make up the whole, the adjustment is different for each. Why? I’m still not sure about this and can see no obvious rationale. But surely if the cause is the same same (e.g. a station move), then the impact on all monthly readings should be the same.

Fig. 10.5: Actual breakpoint correction for Waiouru Airstrip.

The second feature of note is that adjustments are made to data on both sides of the breakpoint, but by different amounts (see Fig. 10.5 above). Those months before November 1985 are adjusted down, while those after are adjusted up. These adjustments are not equal in magnitude, but are scaled in proportion to the number of data points on the opposite side of the breakpoint so that the mean reading remains unchanged but the difference in adjustments equals the total adjustment. For example, the total adjustment for January is 1.21 °C, but there are 12 January data readings before November 1986, and 25 after. So the January readings before November 1986 are reduced by 0.82 °C ( = 1.21x25÷37), with the January readings after November 1986 being increased by 0.39 °C ( = 1.21x12÷37). The mean adjustment is therefore zero but the difference between the two is 0.39 - (-0.82) = 1.21 °C.

 Fig. 10.6: Berkeley Earth adjusted anomaly for Waiouru Airstrip.

The impact of the adjustment on the anomaly data in shown above in Fig. 10.6, the outcome of which is to change the gradient of the best-fit trend line from a negative cooling of -3.24 °C per century (the green line in Fig. 10.1) to a positive warming of +0.82 °C per century (the red line in Fig. 10.6).

Variants on breakpoint adjustment are used by most of the major climate science groups that analyse the global temperature record. In the case of NOAA they are called changepoints and are identified by pair-wise comparisons of neighbouring stations. The effect of these adjustments is to reduce long temperature records to a series of short ones. In the case of NOAA, the records are reduced to a typical length of 15-20 years; for Berkeley Earth it is 13.5 years. You can see this in the New Zealand data I have analysed in the last three posts. Of the ten long records (1200+ months of data) and seventeen medium ones (400+ months of data), all had breakpoints inserted by the Berkeley Earth group, while of the thirty short records only nine did. In the long records alone, there are a total of 46 breakpoints. That is about one breakpoint every 28.6 years and does not include the record gaps and documented station moves that also effectively serve as breakpoints. So the obvious question is, are they justified?

Fig. 10.7: Scaling behaviour of the stemperature anomaly noise for Auckland with and without breakpoint adjustment.

The assumption is that breakpoint adjustments correct the data and return it to its true values and behaviour that would be the case before human error was introduced. Yet if we look at the scaling behaviour I described in the last post, this suggests that the impact may be more negative than positive.

The data in Fig. 10.7 above compares the scaling behaviour of the temperature anomaly for Auckland-Albert Park (Berkeley Earth ID - 157062) before and after the breakpoint adjustments were added. It can clearly be seen that the original data is less scattered than the adjusted data. The gradient for the original data is -0.248 while it is -0.266 for the adjusted. But it is the scatter of the data that is most striking. If we calculate the percentage error (or residual) by comparing the r.m.s. value of the difference between the y-values and the best fit line with the r.m.s. value of the y-values themselves, the results are 4.8% for the original data and 8.1% for the adjusted data.

Fig. 10.8: Scaling behaviour of the temperature anomaly noise for Wellington with and without breakpoint adjustment.

A similar pattern is seen for Wellington-Kelburn (Berkeley Earth ID - 18625) above. The gradient of the best fit line is -0.235 for the original data and -0.239 for the adjusted data. While the gradients are almost identical the residuals are again different, with the adjusted data having a residual (3.3%) that is more than double that of the original data (1.6%).

Fig. 10.9: Scaling behaviour of the temperature anomaly noise for Christchurch with and without breakpoint adjustment.

Finally, if we look at the data above (Fig. 10.9) from the station at Christchurch AP-Harewood (Berkeley Earth ID - 157045), similar differences are seen with the gradient of the best fit to the original data being -0.275 while it is -0.382 for the adjusted data. Here, though, the adjusted data has a residual (20.3%) that is more than seven times that of the original data (2.6%), an astonishing disparity.

Fig. 10.10: Scaling behaviour of the mean temperature anomaly noise for all ten long stations in New Zealand.

It is, though, not just individual stations where this disparity between the scaling behaviour of the original and adjusted data is seen. It is seen in the regional trends as well. So in Fig. 10.10 the trend data from Fig. 7.6 in Post 7 is analysed. This data represents the mean anomaly of all ten long stations in New Zealand. Here the scaling follows a best fit line with a slope of -0.219 and a residual of 1.4%, while the breakpoint adjusted data for all ten long stations (see Fig. 7.8 in Post 7) has a scaling behaviour with a best fit line of slope -0.198 and a residual of 3.2% (see Fig. 10.11 below). So again the residual for the adjusted data is more than double that of the original data.

Fig. 10.11: Scaling behaviour of the mean temperature anomaly noise for all ten long stations in New Zealand after breakpoint adjustment.

What I think this shows is that breakpoints are unsubtle in their impact. The claims that they are surgical in their precision and benign with regard to the general temperature trend are dubious at the very least. At best they may be seen as blunt statistical interventions that correct some of the worst measurement errors, and at worst as being a greater corrupting influence on the data than the errors that they were intended to eliminate.

9. Fooled by randomness

Is global warming real? That is probably a justifiable question given what I revealed in the last post about breakpoint alignment. But what I am going to demonstrate here and over the next two or three posts should also make you question everything you think you know about climate change. The first topic I am going to explore is a concept that most physicists and mathematicians are all too familiar with, but which appears to be totally off the radar of climate scientists: chaos theory and fractal geometry.

Fig. 9.1:  Record 1.

First a test. Look at the dataset above (Fig. 9.1) and the one below (Fig. 9.2). Can you tell which one is a real set of temperature data and which one is fake?

Fig. 9.2:  Record 2.

Okay, so actually it was a trick question because they are both real sets of data. In fact they are both from the same set of station data, and they are partially from the same time period as well, but there is clearly a difference. The difference is that the data in Fig. 9.1 above is only a small part of the actual temperature record but the data from Fig. 9.2 is from the entire record. The data in Fig. 9.1 is taken from the Christchurch station (Berkeley Earth ID - 157045) and is monthly data for the period 1974 - 1987. The data in Fig. 9.2 is from the same record but for the time interval 1864 - 2013: it has also been smoothed with a 12 month moving average. Yet they look the very similar in terms of the frequency and height of their fluctuations - why? Well, what you are seeing here is an example of self-similarity or fractal behaviour. The temperature record for Christchurch is a one-dimensional fractal, and so for that matter is every other temperature record.

Self-similarity is common in nature. You see it everywhere from fern leaves and cauliflowers to clouds and snowflakes. It is observed when you magnify some objects and look at them in greater detail, only to find, to your surprise, that the detail looks just like a smaller version of the original object. This is known as self-similarity: the object looks like itself but in microcosm. It is also an example of scaling behaviour. There is usually a fixed size ratio between the original and the smaller copies from which it is made.

In order to make the smoothed data in Fig. 9.2 look similar to the original data in Fig. 9.1 two scaling adjustments were made. First the time scale on the horizontal axis in Fig. 9.2 was shrunk by a factor of twelve. This is to compensate for the smoothing process which effectively combines twelve points into one. The second was to scale up the temperature axis in Fig. 9.2 by a factor 12 ^0.275. The reason for the power of 0.275 will become apparent shortly, but it is important as it has profound implications for the noise level we see in temperature records over long time periods (i.e. centuries).

To demonstrate the scaling behaviour of the temperature record we shall do the following. First we smooth the data with a moving average of length say two points and then calculate the standard deviation of the smoothed data. Then we repeat this for the original data, but with a different number of data points in the moving average and again calculate the standard deviation of the new smoothed data. After doing this for six or seven different moving averages we plot a graph of the logarithm of the standard deviation versus log(N) where N is the number of points used each time for the moving average. The result is shown below in Fig. 9.3.

Fig. 9.3:  Plot of the standard deviation of the smoothed  anomaly data against the smoothing interval N for temperature data from Christchurch (1864-2013).

The important feature of the graph in Fig. 9.3 is that the data lies on an almost perfect straight line of slope -0.275 (remember that number)? I have to confess that even I was shocked by how good the fitting was when I first saw it, particularly given how imperfect temperature data is supposed to be. What this graph is illustrating is that as we smooth the data by a factor N, the noise level is reducing by a factor N^-0.275. But is this reproducible for other data? Well the answer appears to be, yes.

Fig. 9.4:  Plot of the standard deviation of the smoothed  anomaly data against the smoothing interval N for temperature data from Auckland (1853-2013).

The graph above (Fig. 9.4) shows the same scaling behaviour for the station at Auckland (Berkeley Earth ID = 157062) while the one below (Fig. 9.5) illustrates it for the station at Wellington (Berkeley Earth ID = 18625). The gradients of the best fit lines (i.e. the power law index in each case) are -0.248 and -0.235 respectively. This suggests that the real value is probably about -0.25.

Fig. 9.5:  Plot of the standard deviation of the smoothed  anomaly data against the smoothing interval N for temperature data from Wellington (1863-2005).

But it is the implications of this that are profound. Because the data is such a perfect fit in all three cases, we can extrapolate to longer smoothing operations such as one hundred years. That corresponds to a scaling term of 1200 (because it is equal to 1200 months and thus is 1200 greater in period than the original data) and a noise reduction of 1200 ^0.25 = 5.89. In other words, the noise level on the underlying one hundred year moving average is expected to be about six times less than for the monthly data. This sounds like a lot but the monthly data for Christchurch has a noise range of up to 5 °C (see Fig. 9.6 below), so this implies that the noise range on a 100 year trend will still be almost 1 °C. Now if that doesn’t grab your attention, I have to wonder what will? Because it implies that the anthropogenic global warming (AGW) that climate scientists think they are measuring is probably all just low frequency noise resulting from the random fluctuations of a chaotic non-linear system.

Fig. 9.6:  The temperature anomaly data from Christchurch (1864-2013) plus a 5-year smoothing average.

What we are seeing here is a manifestation of the butterfly effect which, put simply, says that there is no immediate causal link between some current phenomena such as the temperature fluctuations we see today and current global events. This is because the fluctuations are actually the result of dynamic effects that played out long ago but which are only now becoming visible.

Fig. 9.7:  Typical mean station temperatures for each decade over time.

To illustrate the potential of this scaling behaviour further we can use it to make other predictions. Because the temperature record exhibits self-similarity on all timescales, it must do so for long timescales as well, such as centuries. So we can predict what the average temperature over hundreds of years might look like (qualitatively but not precisely) just by taking the monthly data in Fig. 9.6, expanding the time axis by a factor of 120 and shrinking the amplitude of the fluctuations by a factor of 120 ^0.25 = 3.310. The result is shown in Fig. 9.7 above. Because of the scaling by a factor of 120, each monthly data point in Fig. 9.6 becomes a decade in Fig. 9.7. The data in Fig. 9.7 thus indicates that the average temperature for each decade can typically fluctuate by about ±0.5 °C or more over the course of time.

Fig. 9.8:  Typical mean station temperatures over 100 years over time.

Then, if we smooth the data in Fig. 9.7, we can determine the typical fluctuations over even longer timescales. So, smoothing with a ten point moving average will yield the changes in mean temperature for 100 year intervals as shown in the graph above (Fig. 9.8). This again shows large fluctuations (up to 0.5 °C) over large time intervals. But what we are really interested in from a practical viewpoint is the range of possible fluctuations over 100 years as this corresponds to the timeframe most quoted by climate scientists.

To examine this we can subtract from the value at current time t the equivalent value from one hundred years previous, i.e. ∆T = T(t) - T(t-100).

Fig. 9.9:  Typical change in the 100-year mean temperature for a time difference of 100 years.

So, as an example we may wish to look at the change in mean temperature from different epochs, say from one century to the next. Well the data in Fig. 9.9 shows just that. Each data point represents difference between the mean temperature over a hundred years at that point in time with the same value for a hundred years previous. Despite the large averaging periods we still see significant temperature changes of ± 0.25 °C or more. However, if we compare decades in different centuries it is even more dramatic.

For example, Fig. 9.10 below predicts the range of changes in the average decadal temperatures from one century to the next, in other words, the difference between the 10-year mean temperature at a given time t and the equivalent decadal mean for a time one hundred years previous. What Fig. 9.10 indicates is that there is a high probability that the mean temperature in the 1990s could be 0.5 °C higher or lower that the mean temperature in the 1890s, and this is just as a consequence of low frequency noise.

Fig. 9.10:  Typical change in mean decadal temperature for a time difference of 100 years.

So why have climate scientists not realized all this? Maybe it's because their cadre comprise more geography graduates and marine biologists than people with PhDs in quantum physics. But perhaps it is also due to the unique behaviour of the noise power spectrum.

If the noise in the temperature record behaved like white noise it would have a power spectrum that is independent of frequency, ω. If we define P(ω) to be the total power in the noise below a frequency, ω, then the power spectrum is the differential of P(ω). For white noise this is expected to be constant across all frequencies up to a cutoff frequency ω_o.

(9.1)

This in turn means that P(ω) has the following linear form up to the cutoff frequency ω_o.

P(ω) = aω

(9.2)

where a is a constant. The cutoff frequency is the maximum frequency in the Fourier spectrum of the data and is set by the inverse of the temporal spacing of the data points. If the data points are closer together then the cutoff frequency will be higher. Graphically P(ω) looks like the plot shown below in Fig. 9.11, a continuous horizontal line up to the cutoff frequency ω_o.

Fig. 9.11: The frequency dependent power function P(ω) for white noise.

The effect of smoothing with a moving average of N points is to effectively reduce the cutoff frequency by a factor of N because you are merging N points into one. And because the noise power is proportional to the noise intensity, which is proportional to the square of the noise amplitude, this means that the noise amplitude (as well as the standard deviation of the noise) will reduce by a factor equal to √N when you smooth by a factor of N.

For a 100-year smoothing the scaling factor compared to a monthly average is 1200, and so the noise will therefore reduce by a factor of 1200 ^0.5 = 34.64 . That means the temperature fluctuations will be typically less than 0.1 °C. This is probably why climate scientists believe that the long term noise will always be smoothed or averaged out, and therefore why any features that remain in the temperature trend must be "real". The problem is, this does not appear to be true.

Instead the standard deviation varies as N ^-0.25. So the intensity of the noise varies as N ^-0.5 and P(ω) will increase as √N. It therefore follows that the power spectrum is not independent of frequency as is the case for white noise, but instead varies with frequency as

(9.3)

and P(ω) will look like the curve shown in Fig. 9.12 below.

Fig. 9.12:  The frequency dependent power function P(ω) for temperature data.

The net result is that the random fluctuations in temperature seen over timescales of 100 years or more are up to six times greater in magnitude than most climate scientists probably think they will be. So the clear conclusions is this: most of what you see in the smoothed and averaged temperature data is noise not systemic change (i.e. warming). Except, unfortunately, most people tend to see what they want to see.

5. Combining temperature records into local trends

To see how temperature records can be added or combined it is first necessary to identify the separate components that go to comprise each record.

If we define the temperature record at a position r_i (where i just denotes a label for that location or weather station) on the Earth's surface at time t to be T_i(r_i,t), then each temperature record T_i(r_i,t) can be thought of as the sum of four distinct components as shown below.

 (5.1)

The first of the terms to the right of the equality is the one we are interested in: G(t). This is the global warming trend for the whole planet. It is in effect two terms in one. The first is the benchmark term G(0), which is the global temperature before the start of our measurement, T_i(r_i,t). The second is the change in G(t) over time which we can denote as ∆G(t). So

 (5.2)

The other terms after G(t) in Eq. (5.1) are (in order), the local warming trend (L_i), the seasonal temperature change (S_i), and the local random variations (W_i). These last three terms all vary with time (t) and the location (r_i). The term L_i(r_i,t), like G(t), is actually the sum of two terms, one that denotes the value of L_i(r_i,t) at the start of the temperature record, L_i(r_i,0), and a second that represents the change of L_i(r_i,t) with time, ∆L_i(r_i,t).

  L_i(r_i,t)  =  L_i(r_i,0)  +   ∆L_i(r_i,t).

(5.3)

The term W_i(r_i,t) is just the local weather. It is therefore random and should average to zero, either over time via a temporal moving average on the dataset, T_i(r_i,t), or over position if we combine records from different stations across the globe.

The term S_i(r_i,t) is related to the monthly reference temperature (MRT) that was discussed in the previous post, which we can denote as M_i(r_i,t). So we know how to calculate this, by adding a sufficient amount of data from the same month but different years, in the same record and over multiple years, and finding the mean. However this value will inevitably contain the benchmark values of G(t) and L_i(r_i,t), namely G(0) and L_i(r_i,0) respectively.

(5.4)

Given that T_i(r_i,t) is our data and ∆G(t) or G(t) is what we are trying to measure, the task then is to be able to identify and remove the other three terms L_i , S_i and W_i. We can do this as follows.

The term L_i(r_i,t) is the local warming trend. This is how much the local temperature data is changing with time in a manner that is different from G(t). This could be due to differences in latitude, altitude or local geography (such as whether the station is inland or close to the sea, or surrounded by mountains). However, what is also true of L_i(r_i,t) is that, in absence of any global warming G(t), it should average to a fixed value. But that fixed value would just in effect be a second global warming term that should by rights be part of G(t). In which case it would be logical for the sum of all the local warming terms to equate to zero.

(5.5)

But as already mentioned, the same should be true for the weather component as well:

(5.6)

These two constraints allow us to find the global warming term G(t). To do so we first need to find the temperature anomaly for station or location i.

(5.7)

This, as explained in the last post, is just the difference between the station temperature, T_i(r_i,t), and the MRT. If we substitute for T_i(r_i,t) and M_i(r_i,t) in Eq. 5.7 from the expressions in Eq. 5.1 and 5.4 respectively, together with Eq. 5.2 and Eq. 5.3, after cancellation of terms we get the expression

(5.8)

Now when we sum the anomalies we get

(5.9)

But we know that the summation over W_i(r_i,t) will tend to zero as the number of stations increases, or if a smoothing process is utilised, as will the summation over ∆L_i(r_i,t). So the net result is that if we sum over all possible stations, then

(5.10)

where N is the total number of stations in the summation. Hence the global warming trend is just the mean of the temperature anomalies.

(5.11)

or

(5.12)

If we just sum over all local stations within a region R we will get a similar result, but one that includes the regional trend ∆L_R(t).

(5.13)

As the number of regions in the summation increases, the ∆L_R(t) term will tend to zero and the result will tend to that in Eq. 5.11. It is important to note, however, that when combining anomalies, those anomalies should be derived using MRTs calculated over the same time period. Typically this is chosen to be from 1961 to 1990, mainly because this period contains data from the largest number of temperature records. A second consideration, though, when choosing a suitable period is to select the one where the temperature is most stable. However, given that most temperature records are fairly recent and recent temperature anomalies tend to show the largest rate of warming, it follows that these two criteria are often mutually exclusive. 

What the above analysis demonstrates is that the global warming trend, either globally or locally, is just the sum of the temperature anomalies. None of this requires the climate scientist to determine the local climatic warming L_i(r_i,t). That is just an additional layer of complexity via which the data may possibly be corrupted or even cynically distorted.

The major caveat that one should apply to the above analysis is that it assumes that all stations are equally important in the summation in Eq. 5.13 or Eq. 5.11. This is generally not true as each station effectively represents the area between itself and its nearest neighbours. Hence stations that are more isolated are, therefore, responsible for a larger surface area and should command a stronger weighting.

To account for this we can give each station a weighting coefficient a_i that reflects the proportion of the land area (A_i) that surrounds station i, in comparison to the total land area A_tot. Thus if

(5.14)

then

(5.15)

It is also worth noting the because the areas A_i sum to A_tot, it follows that the coefficients a_i must sum to unity.

So are these weighting coefficients an important consideration? Well most of the major climate science groups that have derived curves for the global warming (NOAA, NASA-GISS, Berkeley Earth and UK Met. Office Hadley Centre/UEA CRU) seem to calculate the area around each station to great precision and use this in the weighting. Personally I think this is generally a waste of effort. I would argue that in most countries the weather stations appear to be fairly uniformly distributed and any variation is relatively small compared to the much larger measurement uncertainties seen in the temperature data. However, there are certainly enormous differences in station densities between countries and also regions (as noted previously), so it is mainly when combining results from different countries and regions that this becomes a major issue, in my opinion.

The main point, however, is this: To get the global warming trend, you just need to average the anomalies. If you do this for a country or a region, the trend will be unique for that country or region. But this will also show you which countries and regions have the greatest warming.