In this third post of this series, I’ll demonstrate various outputs of the method/approach described in parts one and two of this series. You may need to look at those to follow this. The key point is that the methods described there allow one to estimate the model-predicted, global mean surface temperature (GMST) at any point in time, if one knows the equilibrium climate sensitivity and the temperature trajectory (i.e. response and lag) of a pulsed increase in radiative forcing. Both of these are obtained from CMIP5 climate model output. This can be used to both look at future temperature predictions, and also to examine historic temperature changes in light of estimated RFs for various agents. This post concentrates on the former, the next one on the latter.

A few more points first though. The key one is that Good et al. (2011) demonstrated that a series of linearly scaled RF pulses (at 1% of the instant 4X CO2 experiment), when continued for 140 years, gave a *very* good estimate of that produced by the annual 1% CO2 increase experiments, for the nine models’ outputs available in 2011. Caldeira and Myhrvold (2013) confirm this, and note two and three parameter negative exponential models performed about equally well in this prediction (based on RSMSE criteria; actually these are four- and six-parameter fits–see the paper). Information criteria tell us we should go with the simpler model in such cases, as it will often be a better predictor for cases beyond the original sample range. Note that everything is fully deterministic here–I’m not including any random variation, i.e. systemic variability, measurement error, etc. Also, CO2 is the lone RF evaluated, and thus an important underlying assumption is that the actual RF change is known with high confidence.

First, what do the 20 CMIP5 AOGCM model temperature trajectories of the instant 4X CO2 experiment look like? For the CM2013 two-parameter model fits, they look like this over 1000 years (starting from 1850, which I take as +/- end of the pre-industrial period):

and like this over the first 50 years:

The differences between these trajectories are due entirely to the differing physics/chemistry embedded in each model. Importantly, in every case, the temperature rise in the first ten years varies between 2 and 4 degrees C. [Note: I show only the CM2013 two-parameter fits here, even though the three-parameter fits are reported to be slightly better, because when the former are scaled to the much smaller impulse forcings of the 1% per year experiment and added year upon year (i.e. the Good et al. approach described in the last post), the authors report that the results fit the CMIP5 1% annual increase expt. results better than do those from the three-parameter models. And it’s those more gradual RF increase situations that we’re really interested in.]

So, that’s what the models think it would look like, for the GMST, if the earth instantly went from 280 to 4X (1160) ppm CO2. But our objective is to see what they predict for realistic scenarios, and we can readily conceive of several such, using e.g. the four RCP scenarios, and/or historical trajectories. First, that 1% per year CO2 increase scenario, for 140 years, and then stable thereafter. Graphed again over 1000 years:

There is a slight acceleration near the beginning, followed by a nearly linear increase, and at t = 140 (1990, the year of the last 1% impulse), an abrupt transition to a decelerating increase towards an asymptote, as determined by the ECS for the model.

Secondly, and more realistically still, we can mimic historic, and/or expected future RF increases; there are numerous possibilities here. One often finds comparisons referencing the middle of the 20th century as a baseline, including in the AR5. Since 1942, atm. CO2 has increased, on average, at about a third of a percentage point per year [= (400/310)^(1/72) = 0.0035]. I picked 1942 because it’s at the mid-point of a relative stasis in atm. CO2 that extended from about 1934 to 1950. That rate, for 140 years starting from 1850, followed by stasis, gives these trajectories:

The curves are very similar in overall shape to those for the 1% increase, just of lesser magnitude.

[**EDIT**: I’ve discovered math mistakes for the next three graphs but don’t have time to correct the code and graphs right now. The estimates are *too high*. These will be corrected in the next day or so when I get time; ignore for now.]

**[Fixed]**

Although 0.35% per year is a good mean estimate for CO2 annual increase since 1942, the rate since then (or since 1850 for that matter), has not in fact been constant. It has instead been accelerating, and over the last couple decades or so has been over 0.5%. This must in turn lead to a slightly accelerating CO2 RF, assuming the standard relationship between CO2 and RF. A reasonable estimate for this accelerations, starting from 1850 (delta CO2 assumed to be zero then, rising to 0.55% for ~ the last decade), would be: 1.0055/164 = 0.000032 percent per year. If that atm. CO2 increase rate then went to 0 immediately (i.e. now), and was stable until t = 1000–completely unrealistic–the T trajectories would look like this:

Which, from 1850 to 2014 only would look like this:

Projecting that acceleration all the way to year 2100 gives a 0.84% annual rate of increase at that time (with a mean rate of about 0.4% over the entire period). Drastic CO2 mitigation actions taken at that time would thus give trajectories like this over 1000 years (note that this scenario does not exactly match any of the AR5/CMIP5 RCP scenarios, but it’s closest to RCP 8.5):

Lastly, suppose a scenario in which the atm. CO2 increase and acceleration since 1850 is exactly mirrored by a decrease/deceleration over the next 164 years, thereby representing a plan in which the world starts CO2 mitigation by simply limiting the atmospheric *acceleration* to zero (it continues to grow, but does not accelerate), followed over the 164 years by a steadily decreasing growth rate. Such a scenarios is at least realistic in the sense that the world seems more likely to attempt stopping the heavy bleeding first and then proceed from there: triage. Over 1000 years, those trajectories would now look like this:

And over only the period of symmetry (the next 164 years), like this:

Note that the abrupt T transition points from the previous graphs, in which the atmospheric increase has halted suddenly, have now disappeared as expected, and the temperatures at t = 1000 are accordingly higher.

The point of these exercises is to explore what the implications for the temporal (i.e. “transient”) trajectories of GMST are, as a function of the physics embedded in the AOGCMs used in the CMIP5 experiments, as summarized by the equations given by Caldeira and Myhrvold (2013), and presumed to apply, when scaled, to much lower pulses of atm. CO2 than the 4X CO2 experiment (as per Good et al., 2010, 2011). Keep in mind that these results apply to CO2 forcings only, and thus do not even consider all the well mixed GHGs, let alone other agents not considered such (e.g. aerosols, tropospheric ozone, CH4). In the next post I want to look at what the implications from the models are when considering data on historical temperature changes and various RF agents.

The R code, and CM2013 model fit parameters, are given below. It was easier to copy the code and adjust the values and graph titles for each scenario than to try to write a loop, so I did.

### Computations of expected transient temperature change (GMST) using Caldeira and Myhrvold's (2013) equations and various assumptions about future atm. CO2 dynamics ### ## Prelims: make file, "CM2013.4X.fit.pars.csv", containing CM2013 parameters for fits to 4X CO2 RF (their Tables S4 and S5) (pars = read.csv("CM2013.4X.pars.csv")) # temperatures and forcings, from Table S1 (fitpars = read.csv("CM2013.4X.fit.pars.csv")) # model fits, from Tables S4 and S5 fitpars[is.na(fitpars)]=0 # For the two parameter models, third Theta and Tau params. = 0 fitpars$n.pars = rep(c(4,6),each=20) ## 1. Compute yearly T for nyears from CM2013 2 and 3 parameter fits to CMIP5 4X instant CO2: t1 = 1850; nyears = 1000 results1 = matrix(NA, nrow=nyears, ncol=40) for (i in 1:40){ t = 1:nyears results1[,i] = 1 - ((fitpars$Theta.0[i] * exp(-t/fitpars$Tau.0[i])) + (fitpars$Theta.1[i] * exp(-t/fitpars$Tau.1[i])) + (fitpars$Theta.2[i] * exp(-t/fitpars$Tau.2[i]))) } (ECS = rep(pars[,5],2)) results1 = t(t(results1) * ECS) head(results1) ## 2. Plot same, over full nyears and over first x years only: colrs = round(seq(20,500,length.out=19)); colors()[colrs] # Strictly to upset Ed Hawkins, but note, not rainbow colors here. # a. Two par. models windows() plot(t1:(t1+nyears-1), results1[,1], ylim = c(0,max(ECS)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013, 2 parameter fits", pch=20) for (j in 2:20) points(t1:(t1+nyears-1), results1[,j], col = colors()[colrs[j]], pch=20) # b. Three par. models windows() plot(t1:(t1+nyears-1), results1[,21], ylim = c(0,max(ECS)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013, 3 parameter fits", pch=20) for (j in 22:40) points(t1:(t1+nyears-1), results1[,j], col = colors()[colrs[j-20]], pch=20) nyears = 50 # a. Two par. models windows() plot(t1:(t1+nyears-1), results1[1:nyears,1], ylim = c(0,max(ECS)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013, 2 parameter fits", pch=20) for (j in 2:20) points(t1:(t1+nyears-1), results1[1:nyears,j], col = colors()[colrs[j]], pch=20) for (j in 2:20) lines(t1:(t1+nyears-1), results1[1:nyears,j], col = colors()[colrs[j]], lwd=1) # b. Three par. models windows() plot(t1:(t1+nyears-1), results1[1:nyears,21], ylim = c(0,max(ECS)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013, 3 parameter fits", pch=20) for (j in 22:40) points(t1:(t1+nyears-1), results1[1:nyears,j], col = colors()[colrs[j-20]], pch=20) for (j in 22:40) lines(t1:(t1+nyears-1), results1[1:nyears,j], col = colors()[colrs[j-20]], lwd=1) ## 3. Compute various transient sensitivities (TCS): nyears=1000 results2 = matrix(NA, nrow=nyears, ncol=40) # a. 1% per year CO2 increase, to 4X: rates = 1.01; dt = 140 for (j in 1:40){ results3 = matrix(0, nrow=nyears, ncol=dt) # temporary matrix, not saved for (i in 0:(dt-1)) results3[(i+1):nyears, i+1] = results1[1:(nyears-i), j] results3 = t(t(results3) * log(rates)/log(4)) results2[,j] = apply(results3,1,sum) } # Graph two par. models only: windows() plot(t1:(t1+nyears-1), results2[,1], ylim = c(0,max(ECS)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013 2 par fit, 1% CO2 to 4X", pch=20) for (j in 2:20) points(t1:(t1+nyears-1), results2[,j], col = colors()[colrs[j]], pch=20) nyears = 140 windows() plot(t1:(t1+nyears-1), results2[1:nyears,1], ylim = c(0,max(ECS)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013 2 par fit, 1% CO2 to 4X", pch=20) for (j in 2:20) points(t1:(t1+nyears-1), results2[1:nyears,j], col = colors()[colrs[j]], pch=20) for (j in 2:20) lines(t1:(t1+nyears-1), results2[1:nyears,j], col = colors()[colrs[j]], lwd=1) # b. 0.4% per year CO2 increase, for 140 years: nyears=1000 results2 = matrix(NA, nrow=nyears, ncol=40) rates = 1.004; dt = 140 for (j in 1:40){ results3 = matrix(0, nrow=nyears, ncol=dt) # temporary matrix, not saved for (i in 0:(dt-1)) results3[(i+1):nyears, i+1] = results1[1:(nyears-i), j] results3 = t(t(results3) * log(rates)/log(4)) results2[,j] = apply(results3,1,sum) } # Two par. models windows() plot(t1:(t1+nyears-1), results2[,1], ylim = c(0,max(ECS)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013 2 par fit, 0.4% annual CO2 increase for 140 years", pch=20) for (j in 2:20) points(t1:(t1+nyears-1), results2[,j], col = colors()[colrs[j]], pch=20) nyears = 140 windows() plot(t1:(t1+nyears-1), results2[1:nyears,1], ylim = c(0,max(ECS)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013 2 par fit, 0.4% annual CO2 increase for 140 years", pch=20) for (j in 2:20) points(t1:(t1+nyears-1), results2[1:nyears,j], col = colors()[colrs[j]], pch=20) for (j in 2:20) lines(t1:(t1+nyears-1), results2[1:nyears,j], col = colors()[colrs[j]], lwd=1) # c1. Accelerating RF increases to 2014, via increasing rate of annual CO2 increase: nyears=1000 results2 = matrix(NA, nrow=nyears, ncol=40) dt = 164; rates = 1 + ((1:dt) * (21/379/10/dt)) for (j in 1:40){ results3 = matrix(0, nrow=nyears, ncol=dt) # temporary matrix, not saved for (i in 0:(dt-1)) results3[(i+1):nyears, i+1] = results1[1:(nyears-i), j] results3 = t(t(results3) * log(rates)/log(4)) results2[,j] = apply(results3,1,sum) } # Two par. models windows() plot(t1:(t1+nyears-1), results2[,1], ylim = c(0,max(results2)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013 2 par fit, accelerating CO2 increase to present", pch=20) for (j in 2:20) points(t1:(t1+nyears-1), results2[,j], col = colors()[colrs[j]], pch=20) nyears = 164 windows() plot(t1:(t1+nyears-1), results2[1:nyears,1], ylim = c(0,max(results2)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013 2 par fit, accelerating CO2 increase to present", pch=20) for (j in 2:20) points(t1:(t1+nyears-1), results2[1:nyears,j], col = colors()[colrs[j]], pch=20) for (j in 2:20) lines(t1:(t1+nyears-1), results2[1:nyears,j], col = colors()[colrs[j]], lwd=1) # c2. Accelerating RF increases, to 2100: nyears=1000 results2 = matrix(NA, nrow=nyears, ncol=40) dt = 250; rates = 1 + ((1:dt) * (21/379/10/dt)) for (j in 1:40){ results3 = matrix(0, nrow=nyears, ncol=dt) # temporary matrix, not saved for (i in 0:(dt-1)) results3[(i+1):nyears, i+1] = results1[1:(nyears-i), j] results3 = t(t(results3) * log(rates)/log(4)) results2[,j] = apply(results3,1,sum) } # Two par. models windows() plot(t1:(t1+nyears-1), results2[,1], ylim = c(0,max(results2)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013 2 par fit, accelerating CO2 increase to 2100", pch=20) for (j in 2:20) points(t1:(t1+nyears-1), results2[,j], col = colors()[colrs[j]], pch=20) # d. Accelerating and then decelerating RF, symmetrical, centered on 2014: nyears=1000 results2 = matrix(NA, nrow=nyears, ncol=40) dt1 = 164; incr.rates = 1 + ((1:dt1) * (21/379/10/164)) dt2 = dt1; decr.rates = rev(incr.rates) rates = c(incr.rates,decr.rates) dt = length(rates) for (j in 1:40){ results3 = matrix(0, nrow=nyears, ncol=dt) # temporary matrix, not saved for (i in 0:(dt-1)) results3[(i+1):nyears, i+1] = results1[1:(nyears-i), j] results3 = t(t(results3) * log(rates)/log(4)) results2[,j] = apply(results3,1,sum) } # Two par. models windows() plot(t1:(t1+nyears-1), results2[,1], ylim = c(0,max(results2)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013 2 par fit, symmetrical CO2 rise/fall", pch=20) for (j in 2:20) points(t1:(t1+nyears-1), results2[,j], col = colors()[colrs[j]], pch=20) nyears = dt windows() plot(t1:(t1+nyears-1), results2[1:nyears,1], ylim = c(0,max(results2)), xlab = "Year", ylab = "dT (degr. C)", main = "CM2013 2 par fit, symmetrical CO2 rise/fall", pch=20) for (j in 2:20) points(t1:(t1+nyears-1), results2[1:nyears,j], col = colors()[colrs[j]], pch=20) for (j in 2:20) lines(t1:(t1+nyears-1), results2[1:nyears,j], col = colors()[colrs[j]], lwd=1)

# The parameters for the 2 and 3 parameter models, as given by CM2013 (note that it's actually 4 and 6 parameter models): # = fitpars in code above Model Theta.0 Theta.1 Theta.2 Tau.0 Tau.1 Tau.2 n.pars 1 BCC-CSM1.1 0.566 0.434 0.000 3.150 130.10 0.0 2 2 BCC-CSM1.1(m) 0.587 0.413 0.000 3.000 109.20 0.0 2 3 CanESM2 0.577 0.423 0.000 3.520 160.80 0.0 2 4 CSIRO-Mk3.6.0 0.386 0.614 0.000 3.380 191.90 0.0 2 5 FGOALS-g2 0.482 0.518 0.000 4.540 245.90 0.0 2 6 FGOALS-s2 0.525 0.475 0.000 3.530 287.60 0.0 2 7 GFDL-CM3 0.446 0.554 0.000 3.710 174.00 0.0 2 8 GFDL-ESM2G 0.560 0.440 0.000 2.390 291.90 0.0 2 9 GFDL-ESM2M 0.542 0.458 0.000 2.850 272.20 0.0 2 10 INM-CM4 0.668 0.332 0.000 3.500 462.90 0.0 2 11 IPSL-CM5A-LR 0.569 0.431 0.000 7.380 363.00 0.0 2 12 IPSL-CM5A-MR 0.529 0.471 0.000 5.360 237.40 0.0 2 13 IPSL-CM5B-LR 0.567 0.433 0.000 3.050 149.70 0.0 2 14 MIROC5 0.631 0.369 0.000 2.590 232.10 0.0 2 15 MIROC-ESM 0.541 0.459 0.000 4.920 290.10 0.0 2 16 MPI-ESM-LR 0.571 0.429 0.000 3.050 153.80 0.0 2 17 MPI-ESM-MR 0.590 0.410 0.000 3.070 151.20 0.0 2 18 MPI-ESM-P 0.588 0.412 0.000 2.660 141.80 0.0 2 19 MRI-CGCM3 0.602 0.398 0.000 3.670 132.10 0.0 2 20 NorESM1-M 0.499 0.501 0.000 3.120 201.60 0.0 2 21 BCC-CSM1.1 0.235 0.352 0.412 0.691 6.19 140.7 3 22 BCC-CSM1.1(m) 0.303 0.334 0.363 0.596 8.63 130.2 3 23 CanESM2 0.458 0.245 0.298 2.130 26.98 322.3 3 24 CSIRO-Mk3.6.0 0.197 0.212 0.591 0.801 8.83 208.0 3 25 FGOALS-g2 0.333 0.227 0.440 1.611 27.56 322.2 3 26 FGOALS-s2 0.079 0.453 0.468 0.194 4.71 297.7 3 27 GFDL-CM3 0.181 0.284 0.535 0.745 7.22 185.7 3 28 GFDL-ESM2G 0.130 0.432 0.438 0.295 3.17 294.6 3 29 GFDL-ESM2M 0.160 0.385 0.455 0.372 4.13 275.6 3 30 INM-CM4 0.197 0.481 0.322 0.322 5.31 543.1 3 31 IPSL-CM5A-LR 0.216 0.394 0.390 0.000 16.23 461.3 3 32 IPSL-CM5A-MR 0.185 0.379 0.436 0.442 10.27 294.5 3 33 IPSL-CM5B-LR 0.292 0.316 0.393 0.482 8.81 176.3 3 34 MIROC5 0.259 0.384 0.356 0.639 4.71 253.7 3 35 MIROC-ESM 0.204 0.364 0.432 0.690 9.34 348.8 3 36 MPI-ESM-LR 0.278 0.315 0.407 0.904 6.70 169.0 3 37 MPI-ESM-MR 0.230 0.380 0.390 0.406 5.83 165.0 3 38 MPI-ESM-P 0.302 0.317 0.380 0.577 7.07 162.1 3 39 MRI-CGCM3 0.305 0.356 0.339 0.679 10.55 170.0 3 40 NorESM1-M 0.223 0.297 0.480 0.515 6.88 221.6 3