This post should be the first in a series covering how we might derive an "empirical" estimate for climate sensitivity from ocean heat content (OHC) and surface temperature changes. I’ve touched on this topic a few times in previous posts, but my goal for this one is to have it be a more thorough examination.
The basic equation I want to use here is one we’ve seen quite a few times before:
ΔN = ΔF + λ*ΔT
Here, N is the TOA radiative flux, F is the forcing, T is the surface temperature, and λ is the radiative response term, also referred to as "effective sensitivity" in Soden and Held (2006) and "climate feedback parameter" in Forster and Gregory (2006), both of which can be confusing as they mean slightly different things than "sensitivity" and "feedback" in current climate lexicon. For more information you can see my page here, although I think parts of that may out-of-date enough to not totally reflect my evolving views. For instance, I see more and more evidence that this radiative response to inter-annual fluctuations is a poor predictor of the radiative response on the climate scale. This is why I hope to use differences in longer periods — e.g. 30 years — to determine a more relevant value for λ. Of course, we don’t have a single satellite that runs that long to compare TOA imbalance, but we can estimate it…from the OHC data.
One other thing worth mentioning here is that while theoretically it should be possible to determine the equilibrium climate sensitivity (temperature change with a doubling of CO2) by simply dividing the forcing (~3.7 W/m^2) by the radiative response, this assumes that λ is a constant for different timescales and forcing magnitudes, which is far from true in some models (whether this departs significantly from constant in the "real world" system is debatable). Winton et al. (2010) refer to this in terms of "Ocean Heat Uptake Efficacy", and Isaac Held has a discussion of it on his blog. . Paul_K also discussed this at The Blackboard. This is why the dividing the CO2 forcing by "effective sensitivity" in Soden and Held (2006) calculated from 100 years of the A1B scenario does not directly match the equilibrated temperature change from the idealized CO2x2 scenario. While the latter may be near impossible to accurately determine without thousands of years of well-constrained observations, the former is arguably much more useful anyhow, so I’ll set my sights on that.
2. Forcing Uncertainty
The start of the NOAA OHC data is 1955, so estimating the TOA imbalance from this data means we’re limited to 1955 and on (actually a bit later, but we’ll discuss that in the next section). So, what kind of forcing uncertainty are we looking at over this period? For this, we can first take a look at the GISS forcings, derived in the method described from Hansen et al. (2005), which are also very similar to the AR4 time evolution of forcings modeled using MIROC+SPRINTARS. Similarly, I digitized the GFDL CM2.1 forcing evolution from Held et al. (2010), and we get something similar. Moreover, since the emissions histories are all very similar, it seems the differences in these forcing histories can largely be determined based simply on an "aerosol efficacy factor" for anthropogenic aerosols. Here is a reconstruction of the GFDL CM2.1 forcing history from the GISS forcings, but using an "efficacy" of 0.75 for anthropegic aerosols:
As you can see, it is a pretty solid match. Thus, in order to determine the potential forcing histories used for “observations” in this experiment, I take the bounds for present-day aerosol estimates and compare that to the GISS forcing, and then use this factor for efficacy. The AR4 estimate for direct aerosol forcing is -0.5 +/- 0.4 W/m^2, and for indirect it is a most likely value of -0.7 W/m^2, with 5% to 95% range of -0.3 to -1.8 W/m^2. . From 1955, this produces the following uncertainty bounds for forcings if we use an aerosol forcing from -0.4 to -2.3 W/m^2:
With Monte Carlo, we can use a Gaussian distribution for the direct aerosol forcing and a triangle distribution for the indirect effect to get the following distribution for combined aerosol forcing difference (1955 – present day):
3. Inferred TOA Imbalance
As the measurements are very sparse down to 2000m prior to 2005, Levitus et al. (2012) provides only pentadal averages for OHC in the 0-2000m depths over that time period. We can calculate the approximate 5 year average TOA imbalance using the difference in 5 year OHC averages. For example, we estimate the 1957.5-1962.5 average TOA imbalance by subtracting the 1955-1959 average from the 1960-1964 average (the conversion from differenced Joules of 5 year averages to global average flux is ~0.124). However, we also need to note that not all of the extra energy is stored in the 0-2000m of the ocean…some goes into the atmosphere, into the cryosphere, or into the deeper ocean, so we’ll need to multiply this by some factor. We can use output from GCMs to estimate how well this method works:
First, for GISS-E2-R, a regression suggests only about 70% of the heat goes into this 0-2000m layer, but the reconstruction is quite good:
For GFDL CM2.1, the amount of heat is stored in 0-2000m ocean is about 85% of the imbalance, which seems more realistic. Levitus et al. (2012) estimates that 90% of the heat has gone into the ocean.
For a comparison of the 5 year running averages of OHC 0-2000m between the CMIP5 GFDL CM2.1 runs and the GISS-E2-R runs and the NOAA observations, I downloaded a bunch of ocean temperature data by layer from two different Earth System Grid Federation nodes: here and here. PLEASE NOTE: that this is NOT raw output data from the models, and it took a good amount of time on my part to download and process the data from the freely available kernel ocean temperatures into global heat content, so if you are using the data in the form I make available here (you can see the links in the scripts) I would request you acknowledge here as the source.
Additionally, NOAA provides 1 year averages for OHC 0-2000m from 2005.5-2011.5. By calculating the current imbalance using annual differences based on a regression dOHC/dT, we can estimate our imbalance up to the most recent full year.
4. Estimating Radiative Response Term using Monte Carlo
There are a number of uncertainties present, and to see their net effect I use Monte Carlo with 1000 iterations for each start year, from 1957 to 1985. The prior distribution for the forcing change was described in section 3. For estimating the heat going into the 0-2000m layer, I use a Gaussian distribution with standard deviation equal to the published standard error (from Levitus et al. 2012) for each year, with the OHC value for that year obviously centered on the most likely value presented in that data. This is then diff’d to determine the change in OHC, and converted into TOA imbalance by sampling from a uniform distribution that assumes 75%-85% of any heat imbalance is stored in that 0-2000m layer. Except for the “current” imbalance from 2005-2011, where the values are better constrained, I use 10-year averages (rather than 5) to limit the uncertainty for each start year.
Ultimately, the “radiative response” is determined from ((N2-N1)-(F2-F1))/(T2-T1), where N2 is the 2005-2011 Imbalance, N1 is the ten year average imbalance from an earlier start year period (starting at the StartYear), F2-F1 is the difference in forcing between the most recent 8 years and the average forcing over the earlier 10 year period, and T2-T1 is the difference in surface temperature between the most recent 8 years and the average surface temperature over the earlier 10 year period (the temperature used here is an average between the GISTemp and NCDC temperature datasets).
For those curious, the mean value I get for the TOA imbalance from 2005.5-2011.5 is ~0.57 W/m^2, which is pretty consistent with other estimates.
The figure above shows the 2.5%-97.5% uncertainty (and median) for the observations, versus the span of 5 runs (and their median) for GISS E2-R and GFDL CM2.1. A few things jump out: first, the uncertainty for GISS E2-R is extremely small even compared to GFDL CM2.1, which we could probably attribute to an underestimate of internal variability. Second, that even though the error bars overlap for a number of periods, it would appear that these models underestimate the radiative response, suggesting that they likely overestimate climate sensitivity.
Now, for kicks, if I were to ignore the issue of λ not being a constant (and assuming it was the same for all forcings), I could flip this and get the following pdf for climate sensitivity.
As you can see, such an approach yields “most likely” estimate of around 1.8 K that is outside of IPCC likely range. However, this method in itself fails to constrain the extremely high sensitivities due to the scenarios where negative aerosols have almost completely offset GHG warming, leaving a potential of > 10 K at around 1 in 200.
Furthermore, the observant reader may have noticed that in the graphs of the two GCMs above, this method has actually overestimated the sensitivity. The GISS ER CO2 doubling forcing is 4.06 W/m^2, but the median radiative response is only -1.11 W/m^2/K, yielding a sensitivity of 3.6 K that is greater than it’s known sensitivity of 2.7 K. The median radiative response for GFDL CM2.1 using this method is merely –1.02 W/m^2/K, and with its 3.5 W/m^2 forcing this yields a sensitivity of 3.4 K that is almost its actual equilibrium sensitivity. However, the latter is likely just luck, because the GFDL CM2.1 has an extremely high uptake efficiency factor that should cause an underestimate of the sensitivity in any short-term estimate like this. One need only compare to the “effective sensitivities” in Soden and Held (2006) to see that this method underestimates the radiative response (-1.64 W/m^2/K and –1.37 W/m^2/K for GISS and GFDL respectively), and I’m not sure why exactly; my best guess at this point would have to do with the response to the high volcanic activity during this time period.
I will compare the results we get here to those of other CMIP5 models in (hopefully) my next post in order to see how effective this method might be for determining century scale sensitivity.
When running my script again for the most recent post, I noticed I had used the mean of the 2004-2008 forcings rather than 2005-2011 forcings for the recent period. The change skews the PDF slightly to the left, but otherwise does not result in a huge change. The script has been updated in the link below. Here is the updated look:
Please acknowledge any code or data used from this post. Thanks!
Forster and Gregory (2006): http://journals.ametsoc.org/doi/pdf/10.1175/jcli3611.1
Soden and Held (2006): http://www.gfdl.noaa.gov/bibliography/related_files/bjs0601.pdf
Winton et al. (2010): http://www.met.igp.gob.pe/publicaciones/2010/efficacy.pdf
Hansen et al. (2005): http://pubs.giss.nasa.gov/docs/2005/2005_Hansen_etal_2.pdf
Held et al. (2010): http://journals.ametsoc.org/doi/abs/10.1175/2009JCLI3466.1
Levitus et al. (2012): http://www.agu.org/pubs/crossref/pip/2012GL051106.shtml
Hansen et al. (2011) http://www.columbia.edu/~jeh1/mailings/2011/20110415_EnergyImbalancePaper.pdf