Reflections on Identifiability of Climate Sensitivity and Diffusivity

library('hector')
library('ggplot2')
library('ggthemes')

What we want to know is, how individually identifiable are the climate sensitivity \(S\) and the ocean thermal diffusivity \(\kappa\), given just temperature data over the 21st century. To do this, we will look at a range of \(\kappa\) values, and for each one we will find the \(S\) value that best reproduces the 21st century temperatures for the baseline run. Then we will look at how different those runs are over the 21st century, and whether they diverge with longer runs.

Here is our function to return output data for a single \(\kappa\) value.

compyears <- 2001:2100
kopt <- function(hcore, kappa, temp2100)
{
    setvar(hcore, NA, DIFFUSIVITY(), kappa, 'cm2/s')
    optfun <- function(S) {
        hfail <- function(e) {1e6*(S-2)}  # return a large value on parameter out of range
        tryCatch({
            setvar(hcore, NA, ECS(), S, 'degC')
            reset(hcore)
            run(hcore, max(compyears))
            sum((fetchvars(hcore, compyears, GLOBAL_TEMP())$value - temp2100)^2)
        }, error=hfail) 
    }
    op <- optimize(optfun, c(0.1, 50))
    # if(op$objective > 1e-4) {
    #     warning('Optimization failed. kappa= ', kappa, ' val= ', op$objective)
    # }
    ## Have the location of the optimum, now get the data
    setvar(hcore, NA, ECS(), op$minimum, 'degC')
    runtodate <- function(date) {
        reset(hcore)
        run(hcore,date)
        df <- fetchvars(hcore, 2000:date, GLOBAL_TEMP())
        df$kappa <- kappa
        df$S <- op$minimum
        df
    }
    ## Some models won't be able to run to 2300, so back those down as necessary
    rslt <- NULL
    date <- 2301
    hfail2 <- function(e) {NULL}
    while(is.null(rslt)) {
        date <- date-10
        rslt <- tryCatch({
            runtodate(date)
        }, error=hfail2)
    }
    rslt
}

To run this we’ll need a Hector core. My hunch is that the high emissions scenario will give us the best discriminatory power, so we’ll use that. We also need the 21st century temperatures.

inifile <- system.file('input/hector_rcp85.ini', package='hector')
hc <- newcore(inifile)
invisible(run(hc, max(compyears)))
temps <- fetchvars(hc, compyears, GLOBAL_TEMP())$value
kvals <- c(0.01, 2.3, 10.0, 20.0)
rslts <- lapply(kvals,
                function(k) {
                    r <- kopt(hc, k, temps)
                    message('k= ', k, ' S= ', r$S[1])
                    r
                })

k= 0.01 S= 2.00909721511787
k= 2.3 S= 2.99999511064658
k= 10 S= 6.83593956125711
k= 20 S= 49.9999397962838

Now, let’s take a look at these runs. We’ll first plot only out to 2100, then we’ll plot everything.

pltdata <- dplyr::bind_rows(rslts)
pltdata$kappa <- factor(pltdata$kappa, ordered=TRUE)
pltdata2100 <- dplyr::filter(pltdata, year<=2100)
ggplot(data=pltdata2100, aes(x=year, y=value, color=kappa)) + geom_line(size=1.1) +
    ylab('Temperature') + theme_solarized_2(light=FALSE) + scale_color_solarized()

ggplot(data=pltdata, aes(x=year, y=value, color=kappa)) + geom_line(size=1.1) +
    ylab('Temperature') + theme_solarized_2(light=FALSE) + scale_color_solarized()

As expected, \(S\) and \(\kappa\) can to some extent trade off against one another, but the tradeoff isn’t perfect. They look wildly different for the longer runs, of course. What matters is how this partial tradeoff compares to our ability to fit the ESMs. We know the ESMs have noise that we won’t ever be able to fit. If the magnitude of that noise is comparable to or larger than the difference (in Hector temperatures) between the largest and smallest \(\kappa\) runs, then all the \(\kappa\), \(S\) pairings effectively look the same, and the parameters are not going to be identifiable. Therefore, we’ll calculate the mean squared error relative to the \(\kappa = 2.3\) run for each of the other runs.

sapply(rslts, 
       function(r){
           r2100 <- dplyr::filter(r, year<=2100)
           cmp2100 <- dplyr::filter(rslts[[2]], year<=2100)
           mean(
               (r2100$value - cmp2100$value)^2
           )
       }) -> mse
names(mse) <- kvals
mse

       0.01         2.3          10          20 
0.009771107 0.000000000 0.006692442 0.022865701 

The MSE for the total range of \(\kappa\) values is only about 0.033. That’s pretty small. I bet the noise in the ESMs is already higher than that. We can estimate that noise by calculating the MSE between the ESM data and a smoothed version. We’ll use loess smoothing to make it easy on ourselves.

esmnoise <- function(mod) {
    dd <- dplyr::filter(cmip_individual, model==mod, experiment=='rcp85')
    ll <- loess(data=dd, value~year, span=0.5)
    mean((predict(ll) - dd$value)^2)
}
mods <- unique(dplyr::filter(cmip_individual, experiment=='rcp85')$model)
n <- sapply(mods, esmnoise)
n

     ACCESS1-0      ACCESS1-3        CanESM2          CCSM4      CESM1-BGC     CESM1-CAM5    CESM1-WACCM 
  7.218283e-03   5.839079e-03   1.302654e-02   7.927348e-03   5.569663e-03   1.141015e-02   1.114608e-02 
     CMCC-CESM        CMCC-CM       CMCC-CMS       CNRM-CM5  CSIRO-Mk3-6-0       EC-EARTH      FGOALS-g2 
  3.857231e-02   8.656464e-03   1.597073e-02   5.739298e-03   1.798573e-02   1.019948e-02   4.136632e-03 
      GFDL-CM3     GFDL-ESM2G   GISS-E2-H-CC      GISS-E2-H   GISS-E2-R-CC      GISS-E2-R         inmcm4 
  9.435532e-03   9.271329e+04   2.819891e-03   1.150993e-01   2.735334e-03   1.454114e-01   2.602672e-03 
  IPSL-CM5A-MR   IPSL-CM5B-LR         MIROC5 MIROC-ESM-CHEM      MIROC-ESM     MPI-ESM-LR     MPI-ESM-MR 
  7.553371e-03   1.092388e-02   3.148565e-02   5.838751e-03   6.169950e-03   2.025255e-02   1.304763e-02 
     MRI-CGCM3     NorESM1-ME 
  5.150667e-03   6.506770e-03 

summary(n)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
    0.00     0.01     0.01  3090.00     0.02 92710.00 

Something clearly went wrong for GFDL-ESM2G. Other than that the results look interesting, though not good for identifiability. The smallest MSE is 0.003 for inmcm4. If the true value of \(\kappa\) were near the Hector canonical value of 2.3, we would have at least a fighting chance of ruling out \(\kappa \approx 0\) or \(\kappa \ge 10\), assuming that high-frequency variability is the only source of error (i.e, no fitting error and no compromise between experiments).
For 3/4 of the models the MSE is above 0.006, which is comparable to the difference between \(\kappa \approx 2.3\) and \(\kappa \approx 10\). The range of allowable \(\kappa\) values covers essentially all physically plausible values.

By way of example, here’s CESM1-CAM5 overlaid on the \(\kappa\) plots.

cesm_all <- dplyr::filter(cmip_individual, model=='CESM1-CAM5', experiment=='rcp85')
ggplot(data=pltdata2100, aes(x=year, y=value, color=kappa)) + geom_line(size=1.1) +
    ylab('Temperature') + theme_solarized_2(light=FALSE) + scale_color_solarized() +
    geom_point(data=cesm_all, color='lightgrey')

Well, ok, not quite the illustration I’d hoped because these model runs weren’t optimized to that CESM run, but hopefully it’s clear that when you include all the ensemble members, the band easily spans a width comparable to the width of the band of \(\kappa\) values.

What next?

One thing that might help would be to smooth the ESMs so that we have Hector trying to fit their smooth trend, which is all it can really hope to do anyhow. I’m curious to see if doing the fits in principal component space might accomplish this in a less ad-hoc way than using something like loess (more about the principal components angle tomorrow). Of course, the discrepancies between experiments and ensemble members may still lead to non-identifiability between \(S\) and \(\kappa\). Even if it doesn’t, the fact that one \(\kappa\) is favored doesn’t mean it’s the right value, so I’m not sure how much stock we should put in fits done with smoothed results. However, I do think we want to know what happens when we try those other methods so that we can plan our next move, as well as so we can start to figure out what our findings are going to be for this paper.