The Mystery of CNRM 3.0

I’ve “completed” downloading the data I need for what I described in my last post. However, I cut some corners by deleting/not downloading files for the first few decades of the 20c3m runs (whenever possible). I still have to do another huge purge on my external hard drive to make room if I ever want to have the entire data set!

But here’s what I’ve done so far. I noted that Steve M. found some weird behavior in the CNRM 3.0 20c3m run. Before I wrote any code I saw that the 20c3m run was broken up into three files covering 1860-1899, 1900-1959 and 1960-1999. I immediately thought the second file may have been mixed up with some other model run since that jump and drop in anomaly occurs over that time interval. The figure below shows this feature for the tropical troposphere.

I took the three netCDF files and loaded them into Matlab and computed area-weighted averages for the first 12 pressure levels. Steve M. reports here the T2LT weights by altitude via Christy. I multiplied the 12 time series by their respective weights, summed them up, took anomalies relative to 1980-1999 and here’s what I found:

Keep in mind that this is for the globe as a whole, not just the tropics. My replication doesn’t show the exact same behavior as the Santer version.  To find what exactly is causing this jump, let’s look at time series at the first six pressure levels.

The discontinuity begins with the second level. The jump continues up to the last level. From this, I can’t see how the data drops back off in 1960. I smell a fish.

CMIP Data or What I’ve been up to.

I’ve spent most of the day downloading atmospheric temperature and skin surface temperature data from the CMIP archive. At last count I’m up to 36 GB of data and I’m pretty close to being done. My purpose is to generate synthetic MSU temperatures like those used by Santer that are readily available on the PCMDI website. Steve M. recently posted on a strange anomaly (pun) with the CNRM3.0 model showing a maximum value of over 13 °C:

I don’t know whether the error rests with Santer’s collation, with the archiving at PCMDI or with the underlying model.

I’ll find out soon enough. First, I need to obtain the weighting function which I’m in the process of doing right now. Second, I need to pick a model whose data (from the Santer release) doesn’t have any obvious errors so I can compare the results I get as a check on my calculations. Then I’ll see what happens when I process the CNRM 3.0 data. More to come later.

HadCrut vs. HadCM3 – The Battle Concludes

Following up on my last post documenting tensions between HadCRUT and HadCM3, I’ll quickly recap.

Computational Steps

1. Load HadCRUT and HadCM3 data on the inteval 1980-present. Compute ensemble for the 20C3M portion of HadCM3.

2. Since HadCRUT and HadCM3 had different grid spacing, I smoothed the HadCM3 data onto a grid map matching HadCRUT’s grid spacing.

3. Compute anomalies for HadCRUT and HadCM3 relative to the period 1980-1999.

4. Calculate HadCRUT minus HadCM3.

5. Compute the Nychka-corrected OLS trend for each grid point as well as the relevent statistics that are generated in the process (AR(1) coefficient, standard error, margin of error, results of testing the null hypothesis and p-values.) In order for the regression to proceed forward, at least 30 real-valued points were required.

The Results

Figure 1 shows a map of anomalies. I adjusted the color scale because a few extreme values forced the scale up to about 2°C/year effectively suppressing nearly all of the variability in the data.

Figure 1

We can calculate the trend line’s 5-95% confidence interval by computing area-weighted averages of the upper and lower limits on the OLS results. Doing so we get [-0.2577,  0.1579] °C/year. I used this method rather than computing it assuming a normal distribution because such an assumption would be false as is obvious from figure 2.

Figure 2

Now we turn to figure 3 for the AR(1) coefficient for each grid box.

Figure 3

The areas showing a high persistance in the residuals should be expected to be most likely to fail to reject the null hypothesis. Figure 4 shows the results of the hypothesis tests. The black boxes are rejection regions, the white boxes are fail-to-reject regions and the grey boxes signify missing data. A significance level of 5% was used. Interestingly, the areas with high AR(1) values (most of the ocean in the tropics) didn’t reject as much as I thought they would since their effect would be to strongly widen the confidence intervals.

Figure 4

Figure 5 shows the p-values.

Figure 5

It appears that HadCM3 isn’t terribly bad. I say “appears” because with autocorrelated noise comes large type I errors. But when we correct for it, we exchange it for large type II errors. We may not be rejecting the null hypothesis for precisely that reason. That being said, considering the fact that the regions with high AR(1) coefficients did not strongly correlate (at least visually) with regions where we fail-to-reject, it may not be as simple as attributing these results to type II errors.

We’ll see how this pans out for other models. Pending any critical comments on this analysis, I’ll proceed with reproducing these calculations for other models.

UPDATE(5/5/09): If I understood Lucia correctly, I’ve added a graph of the area-weighted difference and its trend. The results were the exact opposite, though we have much less information to deal with because, well, it’s a global average.

HadCrut vs. HadCM3- Let the battle begin! Part 1

In the last two posts I gave a small preview with some pretty graphs to show what I consider a more appropriate though more time-consuming way to compare model data to observations. For this post, I’ll compare HadCrut and their in-house model HadCM3. I obtained the HadCrut and HadCM3 data from here and here, respectively. The code I used for these calculations is not the same that I used for my last two posts. Because I’m dealing with so much data and the calculations were very time consuming, I had to re-write everything to optimize the code’s performance. The difference is like night and day (5 minutes vs. 30 minutes)!

I wrote a script to extract the HadCrut gridded data. The data conatins NaN to flag missing values. I used those missing values to determine the extent of the spatial coverage over time. The next step is to process the HadCM3 data. I used three files. Two corresponded to two 20C3M runs spanning Jan. 1860 – Dec. 1999. The third file was for the SRES A1B scenario. I averaged the two 20C3M runs and appended the SRES A1B data. To validate my code, I compared the output from Climate Explorer  to my output. In the recent post on replicating Climate Explorer , I found the a difference of about 0.002 °C between my and Climate Explorer’s output. Figure 1 shows the difference between the numbers reported on CE and what I calculated in Matlab. This difference is much larger than I had found previously.

Figure 1

After failing to find any obvious errors in my code, I checked Climate Explorer and it took me all of one minute to find the problem. If you download the SRES A1B data , Climate Explorer takes the first 20C3M run and appends SRES A1B. It does not compute an average for 20C3M. When computing the annal cycle to calculate anomalies based on this data, from 1860-1999, it partially anomalizes the data. But from 2000 on, a sinusoidal shape appears. This is the left over annual cycle that was accounted for in my calculation, but neglected by Climate Explorer. Anyone using data from Climate Explorer must first get the 20c3m data and have Climate Explore compute the ensemble average. Then manually append the SRES A1B data. To double check that my code really was performing properly, I checked the 20C3M ensemble on CE to what I had already calculated (1860-1999) and found the difference to be exceedingly small with a magnitude of at most 0.0015 °C as shown in figure 2.

Figure 2

Now let’s turn to the HadCRUT data. I calculated a simple area-weighted average and at first look, it appeared to match what you can download from the Met Office. I took the difference to check and found significant discrepencies as shown in figure 3.

Figure 3

I initially suspected the source of the error may have been in the way I calculated the area map but quickly discounted that because the error would roughly be constant. I thought the problem may lay in the fact that this data is in anomalies not absolute temperature and some additional processing was done. I attempted to reconstruct the absolute temperature data by adding the annual cycle but found riduculous results. Figure 4 shows the botched results on the left and the correct results on the right. The correct results were found by downloading the area-weighted global average from the Met Office and adding the seasonal cycle.

Figure 4

So what caused the calculated temperatures to be so high? Computationally, I’m taking the gridded data at each time step and multiplying it by a matrix containing fractions of the Earth’s surface area to weight the grid points. If we had complete coverage, then we’d just have to sum up the matrix in both dimensions and we’re done. If there is missing data, then we have to divide the sum by the fraction of the Earth’s surface that is represented by the data. This number will always be smaller than 1 so the calculated sum will increase. The only way the data can avoid being scaled up so much is if there was total coverage! So I took the matrix with HadCRUT anomalies, replaced the NaNs with zeros and proceeded to add the annual cycle. Figure 5 shows the difference between the correct results shown on the right in figure 4 and the results found by replacing the missing values.

(EDIT: 5/11/09 , I made the previous sentence in bold because the figure below it to which it is referring to is mis-titled, rather than re-run the program and generate a new graph.)

Figure 5

I think we can safely say that the mystery has (edit: almost) been solved. Figure 6 shows a map of the temperature data during the first month of the base period.

Figure 6

The only way the entire map can be populated with data is if they interpolated from an incomplete set. I think they do this to maximize the number of grid points with real data.  If they didn’t then the maximum coverage HadCrut as whole could have would be what it was in 1961, about 79%.

In order to make a more level playing field, when calculating anomalies for the HadCM3 data, I will use the unmasked data to get the seasonal cycle and then remove this cycle from the masked data.

I await your feedback.

A Look at What’s to Come

As we saw in my previous post, I was able to replicate Climate Explorer’s ouput for model data to a sufficiently high degree of accuracy. Now that I know the routine works, I can move onto the next step. This post and the ones to follow were inspired when I downloaded HadCrut’s anomalies in their netCDF format. I plotted the data and was amazed to see just how much of the Earth was missing! Figure 1 shows the percent of the Earth’s surface accounted for by HadCrut over the course of the time series.

Figure 1. Note the dips in coverage between 1914-1918 and 1938-1945 coinciding with the World Wars.

It dawned on me that comparing HadCrut to climate model simulations may be a bit unfair. Because of polar amplification, it should warm faster in the Arctic and a bit faster in Antarctica than the rest of the globe. Since the climate models have total coverage, the effects of polar amplification are not lost. Since the surface temperature record is missing a significant portion of the Arctic and practically all of Antarctica, the effects of amplification may be lost. To illustrate the difference in coverage, figure 2 shows HadCrut and the HadCM3 ensemble for the month of December 1999.  I chose this date because I have only so far unpacked the 20c3m data and this is the most recent date I can compare the two. I will piece it together with the SRES A1B data later.

Figure 2

To compare like-with-like, I used isfinite in Matlab to produce a layered matrix of logical values (ones and zeros) for the data at each time step in HadCrut. Ones indicating finite values and zeros indicating non-finite values (NaN). I will use this matrix to mask out the values in HadCM3 so both observational and model data have similar coverage. HadCM3 has higher resolution than HadCrut (2.5° x 3.75° versus 5° x 5°). I need to smooth the HadCM3 data on the same lat/lon vectors so I can do the matrix multiplication to simulate the same missing coverage. I’m still working on that.

After solving that problem, I can then begin calculating global averages for the simulation data and proceed with comparing it to HadCrut. A month ago or so, I downloaded atmospheric temperature data and was going to compute the tropospheric temperature to compare to the MSU data. I intend to pick up on that with the same idea of comparing the two only after their difference in coverage is accounted for, but that’s much further down the road.

UPADTE (4/24 9:01 pm): I solved the smoothing problem and produced the following graph to show the difference between HadCrut and HadCM3.

It’s fairly evident that the model for this particular point in time produces a pretty good match throughout the oceans and on land in the tropics but diverges over land outside of the tropics. The model errors in the Northern Hemisphere are very cold and the errors in the South are very warm. Interesting.

UPDATE (4/24 10:38 pm): This’ll be the last update. I’ve masked the missing values for HadCM3 and produced the following graph.

Now we can compare like with like.

Replicating Climate Explorer

While writing code for analysis for an upcoming post, I was advised by Lucia to check whether or not my calculation of global mean surface temperature derived from raw model data matches that which can be found on Climate Explorer. I had done a quick comparison a while back and noticed the differences to be extremely small. I want to take this moment to show that my Matlab routine does work properly. I downloaded the two runs for HadCM3 running the 20c3m scenario from Climate Explorer and the PCMDI website. I noticed run 1 had a missing value (-999.9000) in Dec. 1999 so I replaced it with NaN. Figure 1 shows the difference between each run as reported by Climate Explorer and my routine. Figure 1.

The differences are at most are about 0.002 °C. The difference is due to that fact that I’m simplifying (out of concerns for memory consumption) the procedure for calculating the mean temperature. The routine generates an area map based on the variables ‘lat_bnds’ and ‘lon_bnds’ contained in the netcdf file. These two variables represent the upper and lower bounds of lat/lon for each grid point. I then extract the actual temperature data into a matrix. I multiply the area map by the temperature data and sum in both dimensions to get the area-weighted mean. The way Climate Explorer computes the area-weighted average was explained to me by Geert Jan van Oldenborgh:

[Climate Explorer] first computes the climatology for each grid point, subtracts that from the full values to obtain anomalies.  Next the anomalies are area-averaged, and the result added to the climatology.

If time, patience and motivation are in sufficient supply tomorrow, I might write a script to calculate the mean temperature using this method. I am concerned however that this may require more memory than I can spare. Matlab has a tendency to close when large memory demands are imposed on it (at least on my computer.)

I am considering publishing my scripts if there are any interested Matlab users who’d like to do their own analysis. I’ve written many scripts that greatly simply things like downloading data sets like GISS, HADCRUT, etc. It really beats copy and pasting from their websites!

UPDATE (4/24 2:30am): I implemented the procedure Geert Jan van Oldenborgh recommended and to my surprise, it didn’t cause too many problems in Matlab as far as memory was concerned. Sure, I had to pull the power cord out of my computer because Matlab froze everything up when I implemented a poorly thought out for-loop for averaging multiple runs from HadCM3. But hey, you can’t have it all. The figure below shows the results for this alternate method.

It’s an improvement over the method I had used previously and it’s not terribly hard to implement. It’s just a matter of writing the necessary routines to make doing the calculations a breeze.

Present-Day Control Runs

At Lucia’s suggestion, I’m took a look at some of the IPCC AR4 control runs from the model data archive. The particular run we’re looking at today is called “PDcntrl“. The linked website give the following description.

The forcing agents of the “PDspup” and “PDcntrl” experiments, greenhouse gases, sulfate aerosol direct effects and solar forcing, are fixed at the ‘present-day’ (CO2 = 348 ppmv, CH4 = 1650 ppbv, N2O = 306 ppbv, solar constant = 1367.0 Wm?2, sulfate aerosol from natural and anthropogenic source).

I found only four models with PDcntrl runs. Those models are ECHO G, MRI CGCM 2.3.2a, NCAR CSCM 3 and NCAR PCM 1. NCAR CSCM 3 had two runs. All of the models do not simulate over the same interval of time as you’ll see in the following figures. NCAR CSCM 3 runs start at year 0 which I assume is 1801. MRI CGCM had a run spanning 1801-1900 in one file and another run from 1901-1950. However, the second file also started from 1801. I simply appended the second file onto the first to get the complete run. I inquired at the website about this discrepency and was told that it didn’t really matter. The reason being that

The year numbers are arbitrary because the climate-forcing factors (like atmospheric CO2) don’t change, so the control run years not linked to any actual calendar years.

Figure one shows all five models with HadCrut for comparison.

Figure 1.

Just eyeballing it, CGCM comes the closest to at least being very close to HadCrut. Since the climate forcing is held constant, we should expect the trend in the data to be either very small or indistinguishable from zero. Since I haven’t figured out how to insert tables that don’t suck, I’ll just copy and paste the output from my program:

The CCSM1 shows a trend of 0.0094 ± 0.0022 °C/century.
The CCSM2 shows a trend of -0.019 ± 0.0023 °C/century.
The CGCM shows a trend of 0.039 ± 0.005 °C/century.
The ECHOG shows a trend of 0.069 ± 0.0018 °C/century.
The PCM shows a trend of 0.016 ± 0.0025 °C/century.

All the trends are statistically significant, but are practically equal to zero! As a side note, the trend is calculated with OLS, but the confidence intervals are calculated in a different-than-usual way. The short answer is: I calculated the statistically significant lags in the residuals and directly removed the autocorrelation structure and recalcualted the standard error. I’ll give a more thorough explanaition in my next post.

Only two of the models had enough data to directly compare to HadCrut (i.e. same baseperiod).  Figure 2 shows the two runs from CCSM and the one run from PCM.

Figure 2.

Figure 3 shows the anomalies for the remaining two models.

Figure 3.

The next thing I looked at were Fourier transforms of the deseasonalized data. Since it’s been deseasonalized, it shouldn’t show any significant component at certain frequencies (like 1 cycle/year). But to my surprise, I found the opposite as shown in these two figures.

Figure 4.

Figure 5

CCSM run 1 shows the signal has some significant components at frequencies 1,2,4 and 5 cycles/year. CCSM run 2 and PCM show a strong peak at 1 cycle/year. ECHO G shows a fairly uniform distribution across the spectrum. CGCM shows a nice peak at 3 cycles/year. All the models, except ECHO G, show significant contribution at low frequencies. This indicates there is red noise in the signal. Comments, criticisms, suggestions?

A Slightly New Approach to Comparing Model/Empirical Data

In my previous comparisons of model simulations to the surface temperature records (STR), the STR’s uncertainty was left out for lack of data. I had mentioned in my first post the intention to try to incorporate uncertainties into my analysis, but the issue was left on the back burner. So far I have only been able to obtain error estimate data for HadCrut (description and data). I contacted James Hansen to get similar data for the GISS dataset and he referred me to Hansen et al. 2006. This paper doesn’t give error estimates for monthly data, but it does for annual data. It states (p. 14288),

Estimated 2(sigma) error (95% confidence) in comparing nearby years of global temperature (Fig. 1A), such as 1998 and 2005, decreases from 0.1°C at the beginning of the 20th century to 0.05°C in recent decades (4). Error sources include incomplete station coverage, quantified by sampling a model- generated data set with realistic variability at actual station locations (7), and partly subjective estimates of data quality problems (8).

This figure of 0.05 °C is the same value I got from “eyeing” the last error bar in the following graph (fig. A2). This gives an annual standard deviation of 0.025 °C.

I’ve put in a request with NCDC to see if they have any error estimates comparable to the one I found at HadCrut. Turning to RSS/UAH, we find in Christy et al. 2000 (PDF), p. 1163

How accurate are the annual anomalies and trends of version D? We will show below that the 95% confidence interval (CI) for annual anomalies T2 · D and T2LT · D is about ± 0.10 K and that the CI of the trend is ± 0.06 K/decade. It is important to understand we are only describing one type of error in these ranges: measurement error.

This implies the standard deviation for the annual anomalies should be 0.10 K/2 or 0.05 K. Unfortunately, this error estimate isn’t as encompassing in sources of error as the data available from HadCrut. Figure 1 shows HadCrut monthly data from Jan. 2001-Dec. 2008 and the 5-95% confidence interval.

Figure 1.

The purpose of obtaining the error estimates is to calculate the standard deviation for use in a Monte Carlo approach in obtaining a probability distribution for the trendline. The standard deviation was found by subtracting the lower 5% combined estimate from the upper 95% combined estimate and dividing by 2*1.96. Looking at the standard deviation, we find a strange cyclical pattern as shown in figure 2. You can see the series from 2001-2008 here to better visualize the seasonal pattern.

Figure 2.

Figure 3 shows the results of a Fourier transform on the standard deviation series. The signal is composed mostly at the frequencies 1,2 and 3 cycles/year. Any ideas why the standard deviation shows this cyclical behavior?

UPDATE 2/07: I’ve found conflicting results using two different implementations of the Fourier transform in Matlab. It should be obvious enough that the signal has a cycle of about 1/year and  3 or 4/year. I first used this implementation and then this one and got different results.  The second one appears to give the right answer, but the plotted results look lousy to say the least.

UPDATE 2/07: I found a Fourier transform implementation that works. I’ve updated the figure.

Figure 3.

Some Preliminary Results

Using Matlab, I ran 10,000 simulations using the standard deviations shown in figure 2 to generate randomized versions of HadCrut using normally distributed random numbers. Assuming the data represent an AR(1) process, I ran a script to calculate the trendline using Cochrane-Orcutt.  I took a look at few of the randomized datasets and found tha there are statistically significant correlations at lags greater than 1. It may be necessary to write another script to check all significant lags up to some point and run Cochrane-Orcutt to get more accurate results. The histogram in figure 4 shows the relative frequency of the calculated trendlines.

Figure 4.

The temperature anomaly trend is found to be -1.4 ± 0.73 °C/century. This is a quick and dirty post intended to introduce a methodology that I’m still working on. I’m sure there are plenty of issues to iron out before taking it all out.

Comments? Criticisms?

I should have written this earlier

I haven’t posted anything for quite a while now because I was busy packing to go back to school. I moved back yesterday and I’m settling in and trying to come up with some ideas for the next post. For the past few days, I’ve been studying simplified climate models. No good climate blog worth its salt can go for long without posting on even a simple energy balance model! So that’ll probably be the topic of my next post. I found a simplified approach to calculating the absorption and emmission of radiation by CO2 and H20 in my old heat transfer book. Who knew such things could be so useful after you’ve finished the class?

Model Data and Falsification

Lucia recently posted another in a series of hypothesis tests comparing the central tendency of the IPCC projections for warming to the surface temperature record. It was her website that inspired me to start my own blog. What really motivated me was my thought that it wasn’t ideal to get a trend line from seven years of surface temperature data and compare it to a trend line with no confidence interval based on 30 years of model data. I took up the issue in my third post but I found that my use of ttest2 was flawed. Well, I modified the function to perform a two-sample t-test with unequal variances. See the Wikipedia entry for a brief rundown.

I took data from 15 models running scenario A1B and calculated two temperature trends based on Jan-2001/July-2008 and Jan-2001/2031. The figure below shows the distribution for 7 and 30-year trends. The surface temperature record (STR) that’s plotted is a merge of GISS/HADCRUT/NCDC. The data for the 7-year and 30-year trends were baselined on their respective 7-year interval and 30-year intervals. The method used to calculate the temperature trends was Cochrane-Orcutt.

As I suspected, the temperature trends are more concentrated in the oft-quoted 2 °C/century region when we trend line out 30 years. Just eyballing the first graph we should expect that a few models should fail to reject in a hypothesis test. When I was writing the code to generate the graphs, I hadn’t even written anything for the hypothesis tests. I didn’t see the point in doing them for the 30-year trends since it is obvious from the graphs they would all fail! But I did them anyway. The first table shows the calculated temperature trends (95% CI). The second table summarizes the hypothesis test results. The null hypothesis is that the temperature trend in the models equals the temperature trend in the surface temperature record. The alternate hypothesis is that they are not equal. I tested with a 5% significance level.

Sorry if the second table looks dithered (damn WordPress!) but I’m sure you can read it. FGOALS, ECHO G and ECHAM 5 performed excellently. FGOALS, aka. the Alternating Current Model had a VERY large spread! This model’s results makes me glad many models are used because its mean is about 0 °C/century with half the confidence interval on either side! So is it warming or cooling? Echo G is the best fit to the STR. ECHAM 5 shows behavior similar to FGOALS. If those two models were politicians and you asked them where they stand on the whole global warming debate, they’d straddle both sides of the issue! (That was a joke.) For the most part, the p-values corresponding to rejecting the null hypothesis were extremely small. But in some cases where reject was the result, it would have failed to reject if we used 1% as our significance level. Comparing to GISS, GISS EH and BCM 2.0 rejected at 2.51% and 1.13%. Compared to HadCrut, INM CM 3.0 rejects at 2.55%. Compared to NCDC, BCM 2.0 and GISS EH reject at 1.43% and 4.89%. Compared to Merged, CM 3.0 passes at 5.25%. Finally, compared with Merged+RSS, CM 3.0 rejects at 3.65%.

As I suspected, all the models rejected at 5% significance with p-values on the order of 10^-12 or worse. If the raw numbers weren’t enough for you, I generated the following plot to show the relative distributions of the STRs used in the hypothesis tests.

Like most people would suspect, GISS is on the higher end. Interestingly, NCDC is a bit higher. I’ll leave it at that.

Update (8/25/08): I’ve added UAH to the mix.