This post is more or less a thinking-out-loud exercise, but I thought I would inflict it on you, since there are not enough blog posts about climate change. It’s a comparison of point estimates of expected global mean surface temperature change, as given by the IPCC AR5, and ones I computed using Transient Climate Sensitivities (TCS).
This was inspired by a recent cartoon in xkcd and a conversation with Brian McGill (starting here) over at Dynamic Ecology. It’s just a cartoon, and the author seems to try hard to get science and stat topics right. But on first look, I thought it was not even “ballpark” correct, except perhaps under the highest radiative forcing (RF) scenario (RCP85), whereas some other people thought it was ballpark OK. Of course, one person’s concept of “ballpark” might not be another’s.
The projections in the AR5 and CMIP5 are based on a set of very complex Global Climate Models that are spatio-temporally specific. There’s no way one can come close to mimicking that detail, but it is reasonable to expect a decent agreement, for global mean temperature estimates, between these GCM outputs and simpler estimates made using just (1) the given CMIP5/AR5 RF values used, and (2) estimates of the TCS. And that is what I’ve done here. Note that I’m not addressing the relative likelihoods of the world experiencing one RCP scenario over another in the future–that’s a very important, but separate, question.
Note also that this is not an indictment of IPCC AR5 numbers. I criticize a lot of stuff here, but this is not one of those criticisms. I get lower point estimates than they do, but there a number of possible reasons for this, and my working assumption is that if somebody’s numbers are less defensible than the others, those are mine, not the IPCC’s. I really just want to get a handle on why the numbers are different.
Update: The script (somewhat rearranged from what’s below) and all of the input data used can be downloaded here.
First the comments and notes:
## Script to evaluate expected T in relation to atm [CO2] forcings under varying RCP radiative forcing scenarios ## # June 2014, Jim Bouldin # # All GHG mixing ratios and CO2 equivalent forcings, are taken from the RCP Database, Version 2.0.5 # These data are at: http://tntcat.iiasa.ac.at:8787/RcpDb/dsd?Action=htmlpage&page=download, in section titled "CMIP5 recommended data, RCP concentrations" # HOWEVER, the formatting of the .dat (text) files there is very strange (appears to be either column-formatted or fixed width, but neither format reads properly) and a huge pain to read properly if working in R; therefore: # Use instead the Excel-formatted data found at http://www.pik-potsdam.de/~mmalte/rcps/ The entire data set in one zip file is at http://www.pik-potsdam.de/~mmalte/rcps/data/ALLDATA_30May2010.zip. Strip all header info first. # The zip file also includes data on all forcings used in all runs, and other data. # These analyses are based on the concepts of equivalent CO2 forcings for part, or all, of the anthropogenic forcings, as thoroughly discussed in AR5 WG1, especially chapters 8 and 12 # Note that TCR (= Transient Climate Response) is defined strictly as the T change resulting from a 1% increase in CO2, (and thus a linear increase in forcing) over 70 years. This defn. requires generalization in order to apply to different RCP's. # The analyses below follow exactly the idea of TCS as a simple generalization of TCR, as described in the first paragraph of AR5 WG1 section 10.8.1 "Transient Climate Response" # See AR5 WG1, ch 10, Fig 10.20, and ch. 12 Box 12.2, for more detail. # I computed a distribution of TCR values using WG1 Figure 9.5 point estimates, combined with the eyeballed modal values shown in Fig 10.20 for studies occurring since AR4
The following gives the relevant code, using the GHG forcing data only, over the period 1950-2100. I picked those dates because those are +/- the ones used to define the cartoon’s claims.
rm(list=ls())
# set working directory with setwd("directory name")
# filenames for radiative forcings are exactly as given in the zip file referenced above ("ALLDATA_30May2010.zip"), except that .csv replaces .xls
# Assign atm. [CO2] and [CO2eq] (CO2 equivalent) values (ppm) at reference time points:
f26 = read.csv("RCP26_MIDYEAR_RADFORCING.csv", header = T)
f45 = read.csv("RCP45_MIDYEAR_RADFORCING.csv", header = T)
f60 = read.csv("RCP60_MIDYEAR_RADFORCING.csv", header = T)
f85 = read.csv("RCP85_MIDYEAR_RADFORCING.csv", header = T)
year.start = 1950; year.stop = 2100
yr.1 = which(f26[,1]==year.start); yr.n = which(f26[,1]==year.stop)
options(digits=3)
FORCING = rbind(f26[c(yr.1,yr.n),], f45[c(yr.1,yr.n),], f60[c(yr.1,yr.n),], f85[c(yr.1,yr.n),])
FORCING$RCP = c(rep("rcp26",2),rep("rcp45",2),rep("rcp60",2),rep("rcp85",2))
n = ncol(FORCING);n
FORCING = FORCING[,c(n, 1:(n-1))]; FORCING = FORCING[c(seq(1,7,2),seq(2,8,2)),1:12]
# Read the GCM forcing parameters. These are in AR5 WG1 ch 9, Table 9.5. Copied exactly as given, first 7 columns only.
GCM.pars = read.csv("CMIP5 Model Forcing Parameters.csv")[,1:7]
medianTCR.CMIP5 = median(GCM.pars[,5]); meanTCR.CMIP5 = mean(GCM.pars[,5])
TCR.postAR4 = c(2,1.67,1.6,1.6,1.5,1.25,1.25); TCR.all = c(since.AR4,GCM.pars[,5])
medianTCR.all = median(GCM.pars[,5]); meanTCR.all = mean(GCM.pars[,5])
medianTCR.CMIP5; medianTCR.ALL
meanTCR.CMIP5; meanTCR.ALL
RCP = FORCING$RCP[1:4]; Period = paste(year.start,"-",year.stop,sep="")
GHG_RF = data.frame(RCP,FORCING$GHG_RF[5:8] - FORCING$GHG_RF[1])
dT.med = GHG_RF[,2]*(medianTCR.ALL/3.7)
dT.mean = GHG_RF[,2]*(meanTCR.ALL/3.7)
> data.frame(Period,RCP,dT.med, dT.mean)
Period RCP dT.med dT.mean
1 1950-2100 rcp26 0.945 0.971
2 1950-2100 rcp45 1.632 1.676
3 1950-2100 rcp60 2.219 2.278
4 1950-2100 rcp85 3.428 3.520
That last table gives the results. RCP85 is only giving about a 3.5 degree increase from GHG forcing, from 1950 to 2100. So that is indeed about 1.0 below what the xkcd cartoon claimed, assuming that cartoon was in fact assuming an RCP85 scenario. In my discussion with Brian I thought it was about 2.9, so that estimate was too low.
If I run the same exact same analysis using total anthropogenic forcing, which includes tropospheric ozone, aerosols and perhaps some other minor forcings I get:
# Total anthro. forcing: replace lines 27-30 (above) with the following:
RCP = FORCING$RCP[1:4]; Period = paste(year.start,"-",year.stop,sep="")
ANTHRO_RF = data.frame(RCP,FORCING$TOTAL_ANTHRO_RF[5:8] - FORCING$TOTAL_ANTHRO_RF[1])
dT.med = ANTHRO_RF[,2]*(medianTCR.all/3.7)
dT.mean = ANTHRO_RF[,2]*(meanTCR.all/3.7)
> data.frame(Period,RCP,dT.med, dT.mean)
Period RCP dT.med dT.mean
1 1950-2100 rcp26 0.853 0.876
2 1950-2100 rcp45 1.611 1.654
3 1950-2100 rcp60 2.180 2.238
4 1950-2100 rcp85 3.471 3.564
…again about 3.5 degrees for RCP85. That is just under the midpoint of the AR5 estimates, but the AR5 only estimates starting from 1986-2005, not 1950, and ends at 2081-2100, not 2100 itself. So my estimates are low compared to theirs. I'm not sure what the explanation is for it–there are a number of possibilities. If I use the midpoints of the AR5 time periods instead (1995, 2090), with these forcings as input:
RCP YEAR TOTAL_INCLVOLCANIC_RF VOLCANIC_ANNUAL_RF SOLAR_RF TOTAL_ANTHRO_RF GHG_RF KYOTOGHG_RF CO2CH4N2O_RF CO2_RF CH4_RF N2O_RF
231 rcp26 1995 1.64 0.0587 0.026 1.56 2.32 2.01 2.00 1.40 0.470 0.132
2311 rcp45 1995 1.64 0.0587 0.026 1.56 2.32 2.01 2.00 1.40 0.470 0.132
2312 rcp60 1995 1.64 0.0587 0.026 1.56 2.32 2.01 2.00 1.40 0.470 0.132
2313 rcp85 1995 1.64 0.0587 0.026 1.56 2.32 2.01 2.00 1.40 0.470 0.132
326 rcp26 2090 2.67 0.0000 0.174 2.50 3.05 2.95 2.82 2.30 0.284 0.231
3261 rcp45 2090 4.23 0.0000 0.174 4.06 4.43 4.34 4.24 3.51 0.427 0.306
3262 rcp60 2090 5.29 0.0000 0.174 5.12 5.52 5.42 5.34 4.45 0.505 0.386
3263 rcp85 2090 7.68 0.0000 0.174 7.51 7.78 7.69 7.48 5.97 1.055 0.454
I get these T anomalies as a result:
Period RCP dT.med dT.mean AR5.mid
1 1995-2090 rcp26 0.433 0.444 1.0
2 1995-2090 rcp45 1.145 1.176 1.8
3 1995-2090 rcp60 1.632 1.676 2.2
4 1995-2090 rcp85 2.726 2.799 3.7
If I use only WMGHG forcings instead of total anthro, I get:
Period RCP dT.med dT.mean AR5.mid
1 1995-2090 rcp26 0.334 0.343 1.0
2 1995-2090 rcp45 0.967 0.993 1.8
3 1995-2090 rcp60 1.465 1.504 2.2
4 1995-2090 rcp85 2.499 2.565 3.7
These values are quite a bit less than the midpoints of the AR5 estimates, for each scenario. I’m not sure what the cause is. Please contribute any ideas you may have.
*An idea of the estimates of transient response (TCR) is given by AR5 WG1 ch. 12 Box 12.2, here:
NOTE: I rearranged the code part of the script above in the text file in the zip file mentioned below, for clarity.
Ecologically Oriented 
Thanks for this post. It’s interesting. There are a couple hangups on data input which slow things down though. There is no RCP26_MIDYEAR_RADFORCING. I believe that should be RCP3PD_MIDYEAR_RADFORCING. There’s also no “CMIP5 Model Forcing Parameters” file. A person can make one via the IPCC table you mention, but it’s inconveniently formatted so that takes a bit of work. And as you pointed out, you have to modify the spreadsheets before reading them in.
None of that is hard to deal with, but it does take a bit of time. You could avoid that by saving your workspace and sharing that. That will save all of your variables (and custom functions) so they can be loaded with a single line. All you have to do is type:
save.image(“filename.RData”)
And upload the result. People can then get the same worksapce as you by downloaded it then replacing save.image with load. That’d let them quickly get on the same page as you while the data you’ve provided links for allows them to verify the accuracy of your workspace.
You might consider giving it a try. Either way, I’ll comment again when I’ve had more time to look through this (I need to clean up after dinner now).
Aha! I always completely forget about that function, very useful. Thanks Brandon.
You are correct on the original filenames. I changed RCP3 to RCP26, and created the Model Forcing Parameters file.
Here’s the link to a zip file containing the code (somewhat rearranged from what’s above), plus all necessary input forcing and TCR files:
https://drive.google.com/file/d/0B-_niykXh7REcVhubzBJSm8ycEk/edit?usp=sharing
No problem. I think I realized you could do that when someone made the same suggestion to me once. I’ve learned a lot about R from random conversations on blogs.
Anyway, I managed to get the data set up, and I got the same results as you. Then I looked at Chapter 9 a bit more closely, and I think I spotted an issue. The table you used, Table 9.5 only shows entries for AOGCMs. Table 9.1 has more models listed because it includes ESMs, and it indicates the CMIP5 ensemble used both types. It could be the AOGCMs had lower sensitivities than the ESMs. If so, not including the ESMs would result in lower estimates. I’m not sure offhand. I didn’t locate similar results for them yet.
I remember doing that when I first started with R, but never since for some reason.
OK, good spot on the ESM issue, I will check it. However, since ESMs differ from non-ESM’s mainly in their inclusion of long-term feedbacks (i.e. via the carbon cycle), I can’t see that it should make a difference. But who knows.
btw, you’re the first person to actually check anything since my tree ring posts last year. Which is fairly discouraging frankly.
You wouldn’t think the exclusion/inclusion of ESMs would cause much difference, but it’s tricky to know just how the IPCC got the results it got. Another issue I noticed is Chapter 12 indicates (e.g. Figure 12.5) a different number of models were available for each RCP. That shouldn’t cause the results for all RCPs to be higher than you got, but it’s another complexity when examining the results.
Another thing I noticed is Chapter 12 says they only used one model from each ensemble. Offhand, I’m not sure what that indicates (there’s a lot of text in the chapters, and I haven’t spent too much time on this). For example, I’m not sure if models like ACCESS 1.0 and ACCESS 1.3 were considered part of the same ensemble or not. I’m sure answers to questions like that can be found, but it bugs me they are conveniently available.
As for people checking your work, I feel your pain. It seems few people will check anything technical. Most people just like being spoon fed answers. I think your best bet is to try to find a handful of people who are interested in examining things in detail and keep them around. Don’t ask me how to do that though. I haven’t the slightest.
Thanks for those thoughts Brandon–there is indeed a lot of text (and graphics) to wade through. The discussion on feedbacks and forcing is scattered over 3 different chapters.
Anyway, this is getting kind of interesting; more tomorrow.
I submitted a comment earlier which hasn’t shown up. It doesn’t show up as in moderation either, so I’m guessing it flagged as spam. Could you check and see?
Also, that version is definitely easier to work with. Thanks for uploading it.
Sure enough, that’s where your 2 comments were. No idea what’s happening here now.
And no problem on the upload, good idea on your part.
My last comment went through with no problem. I’m guessing all the acronyms in my one comment made WordPress thing I was a spammer.
Now I’m going to check my spam folder. I can’t remember the last time I did.
In checking a few more things I’m finding nothing that increases the values I get here. One principal concern was that I was using the TAR/AR4 definition of RF, which is radiative forcing allowing only for the stratospheric T to equilibrate, which occurs in a few months time. But the AR5 uses “effective” RF (ERF), which allows the entire atmosphere, and the land surface, to come to equlibrium, (but no feedbacks therefrom, like water vapor). Table 8.6 shows the two metrics being equal from 1750 to present and text discussions also indicate that the two are very similar.
The computed values will increase by about 10% if you change the forcing attributed to CO2 doubling. In your R script, you use 3.7 Wm-2, which is the “canonical” value, as well as the mean in the first column of Table 9.5. However, Forster et al. Table 1 evaluates the mean over 23 GCMs as 3.4 Wm-2, which is also the mean of the second column of Table 9.5
I knew I could count on you for an insightful comment Harold!
Indeed, I was just noticing that this afternoon. The thing is, TCR values are computed for CO2 doubling based on the traditional RF of 3.71 W/m^2. So I’d have to correct the TCR values also, which would cancel out that 10% effect.
I’m really to the point where I just cannot explain this. These results are very different from the AR5 midpoints–they are near the bottom of the 90% CI’s in fact.
Pingback: Global temperature change computations using transient sensitivity, part two | Ecologically Orientated