Symmetry beta and q10

Purpose

During the concentration driven experiments we saw that \(S\) and \(\kappa\) were able to trade off with one another, making it impossible for the optimization routine to identify unique solutions of \(S\) and \(\kappa\) for the different model fits. We are also concerned that something similar is occurring with carbon cycle parameters in the emission driven calibration. Which is why we opted to apply a penalty to constrain the carbon cycle parameters. Here are the “symmetry tests”, where we fix one of the symmetrical parameters and solve for the remaining Hector parameters for the concentration and emission driven run.

Set Up

library(dplyr)
library(tidyr)
library(ggplot2)
library(knitr)
BASE_DIR <- "/Users/dorh012/Documents/2019/hectorcal/analysis/paper_1"

# Import the data from the carbon cycle symetry results and add info about the pentaly. 
beta_q10 <- read.csv(list.files(file.path(BASE_DIR, 'output', 'emiss_beta_q10'), '.csv', full.name = TRUE), 
                     stringsAsFactors = FALSE)
beta_q10_penalty <- read.csv(list.files(file.path(BASE_DIR, 'output', 'emiss_beta_q10_penalty2'), '.csv',
                                        full.name = TRUE), stringsAsFactors = FALSE)
beta_q10$penalty         <- 'no penalty'
beta_q10_penalty$penalty <- 'penalty'
beta_q10_df <- bind_rows(beta_q10, beta_q10_penalty)
# Import the data for the cliamte system symetry test runs. 
S_kappa <- read.csv(list.files(file.path(BASE_DIR, 'output', 'conc_S_kappa_temp'), '.csv', full.name = TRUE), 
                     stringsAsFactors = FALSE)
S_kappa_HF <- read.csv(list.files(file.path(BASE_DIR, 'output', 'conc_S_kappa_tempHF'), '.csv',
                                        full.name = TRUE), stringsAsFactors = FALSE)
S_kappa_df <- bind_rows(S_kappa, S_kappa_HF)

Concentration Driven Tests

For the concentration driven symmetry tests \(\kappa\) was fixed at values ranging from 0 to 5, the remaining climate variables were optimized by our calibration protocol that uses temp-only or a combination of temp and heat flux. The first plot compares the \(\kappa\) values vs the minimized temperature MSE. The second plot shows how the relationship between \(\kappa\) and \(S\) change when the ocean heat flux constraint it added.

ggplot(data = S_kappa_df, 
       aes(kappa, min, color = model))+ 
    geom_line() + 
    geom_point() + 
    facet_wrap('comp_data')  + 
    labs(x = expression(kappa), 
         y = 'minmized temp MSE')

The relationship between \(\kappa\) and temp MSE changes from being constant in the temperature-only run to having a defined minimum in the temperature-heat flux set up.

ggplot(data = S_kappa_df, 
       aes(kappa, S, color = model))+ 
    geom_line() + 
    geom_point() + 
    facet_wrap('comp_data') + 
        labs(x = expression(kappa), 
         y = 'S')

The correlation between \(\kappa\) and \(S\) reverses when we add the ocean heat flux constraint. This is similar to what we are seeing in the mcmc runs which is dope.

The symmetry about \(S\) and \(\kappa\) is fairly obvious in this experiment. I would hope for something similar in the set up for the emission driven runs with vs without the penalty.

Shout out to BBL

So I talked to BBL about observational values of \(\beta\) and \(Q_{10}\). And this is what he passed along.

Q10 global scale papers:
* Mahecha et al. 2010, found convergence around 1.5. http://dx.doi.org/10.1126/science.1189587
* Bond-Lamberty and Thomson 2010, global soil respiration Q10 of 1.4 http://dx.doi.org/10.1038/nature08930
* Johnston and Sibly 2018, latitudinal gradient of 2.3 to 2.7. http://dx.doi.org/10.1038/s41559-018-0648-6
* The E3SM team has found that a global Q10 of 1.5 leads to best performance vis-a-vis ILAMB benchmarks and observational records.


The CO2 fertilization literature is more scattered (or maybe I just don’t know it as well). Tree rings, satellites, greenhouse studies, FACE studies–they seem to all give wildly divergent estimates of possible beta strength.

When I spent some time looking at the \(\beta\) literature what I could find didn’t report values, if it did it was rarely a global value and since Hector’s beta is unit less it was hard to get a sense of what values would be comparable to one another.

New constraint, the following is pulled from BBL’s nature paper.

These data suggest a moderate response of global RS to temperature: a Q10 (rate of change of RS with an increase in temperature of 10 °C) of 1.5. This value matches, within confidence limits, global Q10 values for RS (2.1 ± 0.7 and 1.9 ± 0.4) constrained by the observed interannual variability in atmospheric CO2 using the UK Met Office Hadley Centre coupled global model11. 

Emission Driven Tests

Now we suspect that there is trade off occurring between \(\beta\) and \(Q_{10}\). So now we fix \(\beta\) values and solve for the remaining climate and carbon cycle parameters. Not all of the sampled values of \(\beta\) converged, it took a lot more iterations to get the emission driven runs to converge. In the penalty set up \(Q_{10}\) was set constrained using the log mesa function where a = 1.5, b = 2.8, dig = 0.10 from BBL’s paper. This range is very different from the range we found from the DOE PI meeting MCMC.

ggplot(data = beta_q10_df, 
       aes(beta, min_value, color = model))+ 
    geom_line() + 
    geom_point() + 
    facet_wrap('penalty')  + 
    labs(x = expression(beta), 
         y = 'MSE')

As we were expecting without the penalty the MSE is constant across the fixed \(\beta\) values however it does look like CanESM2 has a bit of a trend in it. I don’t really have a sense of what sort of slope in the parameter vs MSE line would mean that optim would be able to solve for a unique solution. My gut says that the gradient of the line is not enough for optim to identify unique solutions but what do you guys think? And that even with the penalty applied there is some un-identifiablity with \(\beta\) less than 0.25.

Pay attention to the scales on these!

ggplot(data = beta_q10_df, 
       aes(beta, q10_rh, color = model))+ 
    geom_line() + 
    geom_point() + 
    facet_wrap('penalty', scales = 'free')  + 
    labs(x = expression(beta), 
         y = expression(Q[10]))

Based off of these two plots I am inclined to say that the penalty function does help resolve some of the the Identifiability issues we are seeing with \(\beta\) and \(Q_{10}\) but I suspect that optim is still going to be used to really low values of \(\beta\).

Conclusions

Apply the penalty function to the calibration set up somewhat works but I am not sure if it works entirely.