Objective

During the emissions-driven run calibration \(\beta\) and \(Q_{10}\) test we saw that that \(\beta\) and \(Q_{10}\) were able to compensate for one another such that they were able to trade off with one another to get the same answer. Even when we tried penalize unreasonable values of \(Q_{10}\), we still saw \(Q_{10}\) and \(\beta\) trade off with one another. What happens to NPP with these parametrizations? Could this variable be used to properly constrain Hector’s carbon cycle?

Run Hector

# Install packages
library(dplyr)
library(tidyr)
library(ggplot2)
library(knitr)
library(hectorcal)
library(cowplot)
devtools::load_all('~/Documents/hector')
# Set up directory
BASE_DIR <- "/Users/dorh012/Documents/2019/hectorcal/analysis/paper_1"
# Import results from the beta and q10 symmetry tests.
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)
# Add info about the penalty
beta_q10$penalty         <- 'no penalty'
beta_q10_penalty$penalty <- 'penalty'
beta_q10_df <- bind_rows(beta_q10, beta_q10_penalty)
# Set up the Hector cores
tibble::tibble( file = c(system.file('input/hector_rcp85.ini', package = 'hector'),
                         system.file('input/hector_rcp85.ini', package = 'hector')),
                experiment = c('esmHistorical', 'esmrcp85')) %>%
    select(ini_file = file, core_name = experiment, experiment) ->
    ini_files_tib
core_list <- mapply(newcore, inifile = ini_files_tib$ini_file, name = ini_files_tib$core_name )
# Run Hector
split(beta_q10_df, interaction(beta_q10_df$model, beta_q10_df$penalty, beta_q10_df$beta, drop = TRUE)) %>%
lapply(function(input){
    input <- data.frame(input)
    param_names <- c(ECS(), BETA(), DIFFUSIVITY(), VOLCANIC_SCALE(), AERO_SCALE(), PREINDUSTRIAL_CO2(), Q10_RH())
    params      <- select(input, param_names)
    info        <- select(input, -param_names)
    lapply(core_list, parameterize_core, params = params)
    lapply(core_list, run, runtodate = 2100)
    lapply(core_list, fetchvars, dates = 1850:2100, vars = c(GLOBAL_TEMP(), HEAT_FLUX(), 
                                                             ATMOSPHERIC_CO2(),  NPP())) %>%
        bind_rows() %>%
        cbind(info) %>%
        cbind(params)
    }) %>%
    bind_rows() %>%
    # Only save results for a single experiment for a single model fit.
    filter(!(variable == "esmHistorical" & year <= 2005)) %>%
    filter(model  == 'CanESM2') ->
    df
# Plot each of the variables
split(df, df$variable) %>%
    lapply(function(input){
        ylab  <- unique(input[['units']])
        title <- unique(input[['variable']])
        ggplot(data = input) +
            geom_line(aes(year, value, group = interaction(beta, q10_rh))) +
            labs(y = ylab, title = title) +
            facet_wrap('penalty')
    }) ->
    list

Results

list$Ca

Without the penalty the various combinations of \(\beta\) and \(Q_{10}\) have results that are almost identical. With the parameter penalty function atmospheric CO\(_2\) results are a bit different. Which is more or less what we were expecting.

list$heatflux

With and without the penalty the various combinations of \(\beta\) and \(Q_{10}\) have results that are almost identical. Which is not surprising.

list$Tgav

With and without the penalty the various combinations of \(\beta\) and \(Q_{10}\) have results that are almost identical. Which is also somewhat expected.

list$npp

With and without the penalty the various combinations of \(\beta\) and \(Q_{10}\) the npp results are drastically different from one another! Are we surprised, nope! This is exactly what we expected!

Conclusions

If it was available NPP could be used to constrain \(\beta\) and \(Q_{10}\) so that they are no longer able to compensate for one another.

LS0tCnRpdGxlOiAiV2hhdCBoYXBwZW5zIHRvIE5QUD8iCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCiMjIyBPYmplY3RpdmUgCgpEdXJpbmcgdGhlIGVtaXNzaW9ucy1kcml2ZW4gcnVuIGNhbGlicmF0aW9uICRcYmV0YSQgYW5kICRRX3sxMH0kIHRlc3Qgd2Ugc2F3IHRoYXQgdGhhdCAkXGJldGEkIGFuZCAkUV97MTB9JCB3ZXJlIGFibGUgdG8gY29tcGVuc2F0ZSBmb3Igb25lIGFub3RoZXIgc3VjaCB0aGF0IHRoZXkgd2VyZSBhYmxlIHRvIHRyYWRlIG9mZiB3aXRoIG9uZSBhbm90aGVyIHRvIGdldCB0aGUgc2FtZSBhbnN3ZXIuIApFdmVuIHdoZW4gd2UgdHJpZWQgcGVuYWxpemUgdW5yZWFzb25hYmxlIHZhbHVlcyBvZiAkUV97MTB9JCwgd2Ugc3RpbGwgc2F3ICRRX3sxMH0kIGFuZCAkXGJldGEkIHRyYWRlIG9mZiB3aXRoIG9uZSBhbm90aGVyLiBXaGF0IGhhcHBlbnMgdG8gTlBQIHdpdGggdGhlc2UgcGFyYW1ldHJpemF0aW9ucz8gQ291bGQgdGhpcyB2YXJpYWJsZSBiZSB1c2VkIHRvIHByb3Blcmx5IGNvbnN0cmFpbiBIZWN0b3LigJlzIGNhcmJvbiBjeWNsZT8gCgojIyMgUnVuIEhlY3RvcgoKYGBge3IsIHdhcm5pbmc9RkFMU0UsIG1lc3NhZ2U9RkFMU0V9CiMgSW5zdGFsbCBwYWNrYWdlcwpsaWJyYXJ5KGRwbHlyKQpsaWJyYXJ5KHRpZHlyKQpsaWJyYXJ5KGdncGxvdDIpCmxpYnJhcnkoa25pdHIpCmxpYnJhcnkoaGVjdG9yY2FsKQpsaWJyYXJ5KGNvd3Bsb3QpCmRldnRvb2xzOjpsb2FkX2FsbCgnfi9Eb2N1bWVudHMvaGVjdG9yJykKCiMgU2V0IHVwIGRpcmVjdG9yeQpCQVNFX0RJUiA8LSAiL1VzZXJzL2RvcmgwMTIvRG9jdW1lbnRzLzIwMTkvaGVjdG9yY2FsL2FuYWx5c2lzL3BhcGVyXzEiCgojIEltcG9ydCByZXN1bHRzIGZyb20gdGhlIGJldGEgYW5kIHExMCBzeW1tZXRyeSB0ZXN0cy4KYmV0YV9xMTAgPC0gcmVhZC5jc3YobGlzdC5maWxlcyhmaWxlLnBhdGgoQkFTRV9ESVIsICdvdXRwdXQnLCAnZW1pc3NfYmV0YV9xMTAnKSwgJy5jc3YnLCBmdWxsLm5hbWUgPSBUUlVFKSwKICAgICAgICAgICAgICAgICAgICAgc3RyaW5nc0FzRmFjdG9ycyA9IEZBTFNFKQpiZXRhX3ExMF9wZW5hbHR5IDwtIHJlYWQuY3N2KGxpc3QuZmlsZXMoZmlsZS5wYXRoKEJBU0VfRElSLCAnb3V0cHV0JywgJ2VtaXNzX2JldGFfcTEwX3BlbmFsdHkyJyksICcuY3N2JywKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZ1bGwubmFtZSA9IFRSVUUpLCBzdHJpbmdzQXNGYWN0b3JzID0gRkFMU0UpCgojIEFkZCBpbmZvIGFib3V0IHRoZSBwZW5hbHR5CmJldGFfcTEwJHBlbmFsdHkgICAgICAgICA8LSAnbm8gcGVuYWx0eScKYmV0YV9xMTBfcGVuYWx0eSRwZW5hbHR5IDwtICdwZW5hbHR5JwpiZXRhX3ExMF9kZiA8LSBiaW5kX3Jvd3MoYmV0YV9xMTAsIGJldGFfcTEwX3BlbmFsdHkpCgojIFNldCB1cCB0aGUgSGVjdG9yIGNvcmVzCnRpYmJsZTo6dGliYmxlKCBmaWxlID0gYyhzeXN0ZW0uZmlsZSgnaW5wdXQvaGVjdG9yX3JjcDg1LmluaScsIHBhY2thZ2UgPSAnaGVjdG9yJyksCiAgICAgICAgICAgICAgICAgICAgICAgICBzeXN0ZW0uZmlsZSgnaW5wdXQvaGVjdG9yX3JjcDg1LmluaScsIHBhY2thZ2UgPSAnaGVjdG9yJykpLAogICAgICAgICAgICAgICAgZXhwZXJpbWVudCA9IGMoJ2VzbUhpc3RvcmljYWwnLCAnZXNtcmNwODUnKSkgJT4lCiAgICBzZWxlY3QoaW5pX2ZpbGUgPSBmaWxlLCBjb3JlX25hbWUgPSBleHBlcmltZW50LCBleHBlcmltZW50KSAtPgogICAgaW5pX2ZpbGVzX3RpYgoKY29yZV9saXN0IDwtIG1hcHBseShuZXdjb3JlLCBpbmlmaWxlID0gaW5pX2ZpbGVzX3RpYiRpbmlfZmlsZSwgbmFtZSA9IGluaV9maWxlc190aWIkY29yZV9uYW1lICkKCiMgUnVuIEhlY3RvcgpzcGxpdChiZXRhX3ExMF9kZiwgaW50ZXJhY3Rpb24oYmV0YV9xMTBfZGYkbW9kZWwsIGJldGFfcTEwX2RmJHBlbmFsdHksIGJldGFfcTEwX2RmJGJldGEsIGRyb3AgPSBUUlVFKSkgJT4lCmxhcHBseShmdW5jdGlvbihpbnB1dCl7CgogICAgaW5wdXQgPC0gZGF0YS5mcmFtZShpbnB1dCkKICAgIHBhcmFtX25hbWVzIDwtIGMoRUNTKCksIEJFVEEoKSwgRElGRlVTSVZJVFkoKSwgVk9MQ0FOSUNfU0NBTEUoKSwgQUVST19TQ0FMRSgpLCBQUkVJTkRVU1RSSUFMX0NPMigpLCBRMTBfUkgoKSkKICAgIHBhcmFtcyAgICAgIDwtIHNlbGVjdChpbnB1dCwgcGFyYW1fbmFtZXMpCiAgICBpbmZvICAgICAgICA8LSBzZWxlY3QoaW5wdXQsIC1wYXJhbV9uYW1lcykKCiAgICBsYXBwbHkoY29yZV9saXN0LCBwYXJhbWV0ZXJpemVfY29yZSwgcGFyYW1zID0gcGFyYW1zKQogICAgbGFwcGx5KGNvcmVfbGlzdCwgcnVuLCBydW50b2RhdGUgPSAyMTAwKQogICAgbGFwcGx5KGNvcmVfbGlzdCwgZmV0Y2h2YXJzLCBkYXRlcyA9IDE4NTA6MjEwMCwgdmFycyA9IGMoR0xPQkFMX1RFTVAoKSwgSEVBVF9GTFVYKCksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgQVRNT1NQSEVSSUNfQ08yKCksICBOUFAoKSkpICU+JQogICAgICAgIGJpbmRfcm93cygpICU+JQogICAgICAgIGNiaW5kKGluZm8pICU+JQogICAgICAgIGNiaW5kKHBhcmFtcykKCiAgICB9KSAlPiUKICAgIGJpbmRfcm93cygpICU+JQogICAgIyBPbmx5IHNhdmUgcmVzdWx0cyBmb3IgYSBzaW5nbGUgZXhwZXJpbWVudCBmb3IgYSBzaW5nbGUgbW9kZWwgZml0LgogICAgZmlsdGVyKCEodmFyaWFibGUgPT0gImVzbUhpc3RvcmljYWwiICYgeWVhciA8PSAyMDA1KSkgJT4lCiAgICBmaWx0ZXIobW9kZWwgID09ICdDYW5FU00yJykgLT4KICAgIGRmCgojIFBsb3QgZWFjaCBvZiB0aGUgdmFyaWFibGVzCnNwbGl0KGRmLCBkZiR2YXJpYWJsZSkgJT4lCiAgICBsYXBwbHkoZnVuY3Rpb24oaW5wdXQpewoKICAgICAgICB5bGFiICA8LSB1bmlxdWUoaW5wdXRbWyd1bml0cyddXSkKICAgICAgICB0aXRsZSA8LSB1bmlxdWUoaW5wdXRbWyd2YXJpYWJsZSddXSkKCiAgICAgICAgZ2dwbG90KGRhdGEgPSBpbnB1dCkgKwogICAgICAgICAgICBnZW9tX2xpbmUoYWVzKHllYXIsIHZhbHVlLCBncm91cCA9IGludGVyYWN0aW9uKGJldGEsIHExMF9yaCkpKSArCiAgICAgICAgICAgIGxhYnMoeSA9IHlsYWIsIHRpdGxlID0gdGl0bGUpICsKICAgICAgICAgICAgZmFjZXRfd3JhcCgncGVuYWx0eScpCiAgICB9KSAtPgogICAgbGlzdApgYGAKCgojIyMgUmVzdWx0cyAKCmBgYHtyfQpsaXN0JENhCmBgYAoKV2l0aG91dCB0aGUgcGVuYWx0eSB0aGUgdmFyaW91cyBjb21iaW5hdGlvbnMgb2YgJFxiZXRhJCBhbmQgJFFfezEwfSQgaGF2ZSByZXN1bHRzIHRoYXQgYXJlIGFsbW9zdCBpZGVudGljYWwuIFdpdGggdGhlIHBhcmFtZXRlciBwZW5hbHR5IGZ1bmN0aW9uIGF0bW9zcGhlcmljIENPJF8yJCByZXN1bHRzIGFyZSBhIGJpdCBkaWZmZXJlbnQuIFdoaWNoIGlzIG1vcmUgb3IgbGVzcyB3aGF0IHdlIHdlcmUgZXhwZWN0aW5nLgoKCmBgYHtyfQpsaXN0JGhlYXRmbHV4CmBgYAoKV2l0aCBhbmQgd2l0aG91dCB0aGUgcGVuYWx0eSB0aGUgdmFyaW91cyBjb21iaW5hdGlvbnMgb2YgJFxiZXRhJCBhbmQgJFFfezEwfSQgaGF2ZSByZXN1bHRzIHRoYXQgYXJlIGFsbW9zdCBpZGVudGljYWwuIFdoaWNoIGlzIG5vdCBzdXJwcmlzaW5nLiAKCgpgYGB7cn0KbGlzdCRUZ2F2CmBgYAoKV2l0aCBhbmQgd2l0aG91dCB0aGUgcGVuYWx0eSB0aGUgdmFyaW91cyBjb21iaW5hdGlvbnMgb2YgJFxiZXRhJCBhbmQgJFFfezEwfSQgaGF2ZSByZXN1bHRzIHRoYXQgYXJlIGFsbW9zdCBpZGVudGljYWwuIFdoaWNoIGlzIGFsc28gc29tZXdoYXQgZXhwZWN0ZWQuCgoKYGBge3J9Cmxpc3QkbnBwCmBgYAoKCldpdGggYW5kIHdpdGhvdXQgdGhlIHBlbmFsdHkgdGhlIHZhcmlvdXMgY29tYmluYXRpb25zIG9mICRcYmV0YSQgYW5kICRRX3sxMH0kIHRoZSBucHAgcmVzdWx0cyBhcmUgZHJhc3RpY2FsbHkgIGRpZmZlcmVudCBmcm9tIG9uZSBhbm90aGVyISBBcmUgd2Ugc3VycHJpc2VkLCBub3BlISBUaGlzIGlzIGV4YWN0bHkgd2hhdCB3ZSBleHBlY3RlZCEgCgoKIyMjIENvbmNsdXNpb25zCgpJZiBpdCB3YXMgYXZhaWxhYmxlIE5QUCBjb3VsZCBiZSB1c2VkIHRvIGNvbnN0cmFpbiAkXGJldGEkIGFuZCAkUV97MTB9JCBzbyB0aGF0IHRoZXkgYXJlIG5vIGxvbmdlciBhYmxlIHRvIGNvbXBlbnNhdGUgZm9yIG9uZSBhbm90aGVyLgoK