The purpose of this note book is to document how we came to decide the protocol for the Hector v3 calibration protocol. Because there have been a number of attempts at and piles of code & figures trying to address this notebook will only contain the critical results that led to our decision.

# Load the required packages. 
library(assertthat)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union
library(FME)
Loading required package: deSolve
Loading required package: rootSolve
Loading required package: coda
library(GGally)
Loading required package: ggplot2
Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2
library(ggplot2)
library(magrittr)
library(tidyr)

Attaching package: ‘tidyr’

The following object is masked from ‘package:magrittr’:

    extract
# While not the most stable way of loading Hector right now we need to make sure that the v3 
# branch is the one that is going to be used here!
HECTOR_BASE <- "/Users/dorh012/projects/Hector-Versions/luc/hector"
devtools::load_all(HECTOR_BASE)
ℹ Loading hector
# Vis setting
theme_set(theme_bw())

npp_flux0 is fixed from litterature!

The decision to fix npp_flux0 at a set value from the literature comes from the fact that while NPP is sensitive to changes in this parameter value other Hector variables are not that sensitive to it.

How does changing npp_flux0 affect Hector output?

# Set up 3 Hector cores each with one low, med, and high NPP0 this is to test 
# interactions between NPP0 and Hector variables. 
# Ito et al. 2011 - doi: 10.1111/j.1365-2486.2011.02450.x
# reports a range of 56.2±14.3 so this will be the range of the values we use in this. 
ini <- file.path(HECTOR_BASE, "inst", "input", "hector_ssp245.ini")
core_245_low <- newcore(ini, name = "ssp245 low")
setvar(core = core_245_low, dates = NA, var = NPP_FLUX0(), values = 41.9, unit = "Pg C/yr")
reset(core_245_low)
Hector core:    ssp245 low
Start date: 1745
End date:   2300
Current date:   1745
Input file: /Users/dorh012/Documents/Hector-Versions/luc/hector/inst/input/hector_ssp245.ini
core_245_mid <- newcore(ini, name = "ssp245 mid")
setvar(core = core_245_mid, dates = NA, var = NPP_FLUX0(), values = 56.2, unit = "Pg C/yr")
reset(core_245_mid)
Hector core:    ssp245 mid
Start date: 1745
End date:   2300
Current date:   1745
Input file: /Users/dorh012/Documents/Hector-Versions/luc/hector/inst/input/hector_ssp245.ini
core_245_high <- newcore(ini, name = "ssp245 high")
setvar(core = core_245_high, dates = NA, var = NPP_FLUX0(), values = 70.5, unit = "Pg C/yr")
reset(core_245_high)
Hector core:    ssp245 high
Start date: 1745
End date:   2300
Current date:   1745
Input file: /Users/dorh012/Documents/Hector-Versions/luc/hector/inst/input/hector_ssp245.ini
# Run Hector! 
run(core_245_low)
Hector core:    ssp245 low
Start date: 1745
End date:   2300
Current date:   2300
Input file: /Users/dorh012/Documents/Hector-Versions/luc/hector/inst/input/hector_ssp245.ini
run(core_245_mid)
Hector core:    ssp245 mid
Start date: 1745
End date:   2300
Current date:   2300
Input file: /Users/dorh012/Documents/Hector-Versions/luc/hector/inst/input/hector_ssp245.ini
run(core_245_high)
Hector core:    ssp245 high
Start date: 1745
End date:   2300
Current date:   2300
Input file: /Users/dorh012/Documents/Hector-Versions/luc/hector/inst/input/hector_ssp245.ini
# Extract and save the results 
dates <- 1850:2300
vars <- c(NPP(), NBP(), GLOBAL_TAS(), OCEAN_UPTAKE(), CONCENTRATIONS_CO2())
npp0_lmh <- rbind(fetchvars(core_245_low, dates, vars), 
                  fetchvars(core_245_mid, dates, vars),
                  fetchvars(core_245_high, dates, vars))
npp0_lmh %>% 
  filter(year <= 2020) %>% 
  ggplot() + 
  geom_line(aes(year, value, color = scenario)) + 
  facet_wrap("variable", scales = "free") + 
  labs(title = "Historcal Hector output at low, medium, and high NPP0", x = NULL, y = NULL)

npp0_lmh %>% 
  filter(year >= 2020) %>% 
  ggplot() + 
  geom_line(aes(year, value, color = scenario)) + 
  facet_wrap("variable", scales = "free") + 
  labs(title = "Future Hector output at low, medium, and high NPP0", x = NULL, y = NULL)

npp0_lmh %>% 
 # filter(year >= 2020) %>% 
  ggplot() + 
  geom_line(aes(year, value, color = scenario)) + 
  facet_wrap("variable", scales = "free") + 
  labs(title = "Hector output at low, medium, and high NPP0", x = NULL, y = NULL)

Therefore we decided to leave NPP as a value fixed from the literature.

ECS is fixed from the litterature!

Our thinking here is that for default Hector we should use the best scientific consensus aka the value from the IPCC AR6 report. When using Hector as an emulator for specific ESMs models this is key tune-able parameter that would need to be adjusted for Hector to act as an emulator.

Alpha & Volscl are fixed

Once again Hector models aerosol radiative forcing and uses the prescribed volcanic radiative forcing that contributes fully to Hector’s global temperature calculation. This is more of an adjustable parameter to account for uncertainty in volcanic and aerosol radiative forcing. We explored the idea of fitting alpha and volscl however but allowing those parameters to be free did not improve the fits. Therefore we set those to be fixee values.

LS0tCnRpdGxlOiAiRGV2bG9waW5nIHRoZSBIZWN0b3IgdjMgZGVmYXVsdCBjYWxpYnJhdGlvbiBwcm90b2NvbCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKVGhlIHB1cnBvc2Ugb2YgdGhpcyBub3RlIGJvb2sgaXMgdG8gZG9jdW1lbnQgaG93IHdlIGNhbWUgdG8gZGVjaWRlIHRoZSBwcm90b2NvbCBmb3IgdGhlIEhlY3RvciB2MyBjYWxpYnJhdGlvbiBwcm90b2NvbC4gQmVjYXVzZSB0aGVyZSBoYXZlIGJlZW4gYSBudW1iZXIgb2YgYXR0ZW1wdHMgYXQgYW5kIHBpbGVzIG9mIGNvZGUgJiBmaWd1cmVzIHRyeWluZyB0byBhZGRyZXNzIHRoaXMgbm90ZWJvb2sgd2lsbCBvbmx5IGNvbnRhaW4gdGhlIGNyaXRpY2FsIHJlc3VsdHMgdGhhdCBsZWQgdG8gb3VyIGRlY2lzaW9uLiAKCgpgYGB7ciwgd2FybmluZz1GQUxTRSwgbWVzc2FnZT1GQUxTRX0KIyBMb2FkIHRoZSByZXF1aXJlZCBwYWNrYWdlcy4gCmxpYnJhcnkoYXNzZXJ0dGhhdCkKbGlicmFyeShkcGx5cikKbGlicmFyeShGTUUpCmxpYnJhcnkoR0dhbGx5KQpsaWJyYXJ5KGdncGxvdDIpCmxpYnJhcnkobWFncml0dHIpCmxpYnJhcnkodGlkeXIpCgojIFdoaWxlIG5vdCB0aGUgbW9zdCBzdGFibGUgd2F5IG9mIGxvYWRpbmcgSGVjdG9yIHJpZ2h0IG5vdyB3ZSBuZWVkIHRvIG1ha2Ugc3VyZSB0aGF0IHRoZSB2MyAKIyBicmFuY2ggaXMgdGhlIG9uZSB0aGF0IGlzIGdvaW5nIHRvIGJlIHVzZWQgaGVyZSEKSEVDVE9SX0JBU0UgPC0gIi9Vc2Vycy9kb3JoMDEyL3Byb2plY3RzL0hlY3Rvci1WZXJzaW9ucy9sdWMvaGVjdG9yIgpkZXZ0b29sczo6bG9hZF9hbGwoSEVDVE9SX0JBU0UpCgojIFZpcyBzZXR0aW5nCnRoZW1lX3NldCh0aGVtZV9idygpKQpgYGAKCgojIG5wcF9mbHV4MCBpcyBmaXhlZCBmcm9tIGxpdHRlcmF0dXJlIQoKVGhlIGRlY2lzaW9uIHRvIGZpeCBgbnBwX2ZsdXgwYCBhdCBhIHNldCB2YWx1ZSBmcm9tIHRoZSBsaXRlcmF0dXJlIGNvbWVzIGZyb20gdGhlIGZhY3QgdGhhdCB3aGlsZSBOUFAgaXMgc2Vuc2l0aXZlIHRvIGNoYW5nZXMgaW4gdGhpcyBwYXJhbWV0ZXIgdmFsdWUgb3RoZXIgSGVjdG9yIHZhcmlhYmxlcyBhcmUgbm90IHRoYXQgc2Vuc2l0aXZlIHRvIGl0LiAKCgpIb3cgZG9lcyBjaGFuZ2luZyBucHBfZmx1eDAgYWZmZWN0IEhlY3RvciBvdXRwdXQ/CgpgYGB7ciwgbWVzc2FnZT1GQUxTRSwgd2FyaW5nPUZBTFNFfQojIFNldCB1cCAzIEhlY3RvciBjb3JlcyBlYWNoIHdpdGggb25lIGxvdywgbWVkLCBhbmQgaGlnaCBOUFAwIHRoaXMgaXMgdG8gdGVzdCAKIyBpbnRlcmFjdGlvbnMgYmV0d2VlbiBOUFAwIGFuZCBIZWN0b3IgdmFyaWFibGVzLiAKIyBJdG8gZXQgYWwuIDIwMTEgLSBkb2k6IDEwLjExMTEvai4xMzY1LTI0ODYuMjAxMS4wMjQ1MC54CiMgcmVwb3J0cyBhIHJhbmdlIG9mIDU2LjLCsTE0LjMgc28gdGhpcyB3aWxsIGJlIHRoZSByYW5nZSBvZiB0aGUgdmFsdWVzIHdlIHVzZSBpbiB0aGlzLiAKaW5pIDwtIGZpbGUucGF0aChIRUNUT1JfQkFTRSwgImluc3QiLCAiaW5wdXQiLCAiaGVjdG9yX3NzcDI0NS5pbmkiKQpjb3JlXzI0NV9sb3cgPC0gbmV3Y29yZShpbmksIG5hbWUgPSAic3NwMjQ1IGxvdyIpCnNldHZhcihjb3JlID0gY29yZV8yNDVfbG93LCBkYXRlcyA9IE5BLCB2YXIgPSBOUFBfRkxVWDAoKSwgdmFsdWVzID0gNDEuOSwgdW5pdCA9ICJQZyBDL3lyIikKcmVzZXQoY29yZV8yNDVfbG93KQoKY29yZV8yNDVfbWlkIDwtIG5ld2NvcmUoaW5pLCBuYW1lID0gInNzcDI0NSBtaWQiKQpzZXR2YXIoY29yZSA9IGNvcmVfMjQ1X21pZCwgZGF0ZXMgPSBOQSwgdmFyID0gTlBQX0ZMVVgwKCksIHZhbHVlcyA9IDU2LjIsIHVuaXQgPSAiUGcgQy95ciIpCnJlc2V0KGNvcmVfMjQ1X21pZCkKCmNvcmVfMjQ1X2hpZ2ggPC0gbmV3Y29yZShpbmksIG5hbWUgPSAic3NwMjQ1IGhpZ2giKQpzZXR2YXIoY29yZSA9IGNvcmVfMjQ1X2hpZ2gsIGRhdGVzID0gTkEsIHZhciA9IE5QUF9GTFVYMCgpLCB2YWx1ZXMgPSA3MC41LCB1bml0ID0gIlBnIEMveXIiKQpyZXNldChjb3JlXzI0NV9oaWdoKQoKIyBSdW4gSGVjdG9yISAKcnVuKGNvcmVfMjQ1X2xvdykKcnVuKGNvcmVfMjQ1X21pZCkKcnVuKGNvcmVfMjQ1X2hpZ2gpCgojIEV4dHJhY3QgYW5kIHNhdmUgdGhlIHJlc3VsdHMgCmRhdGVzIDwtIDE4NTA6MjMwMAp2YXJzIDwtIGMoTlBQKCksIE5CUCgpLCBHTE9CQUxfVEFTKCksIE9DRUFOX1VQVEFLRSgpLCBDT05DRU5UUkFUSU9OU19DTzIoKSkKbnBwMF9sbWggPC0gcmJpbmQoZmV0Y2h2YXJzKGNvcmVfMjQ1X2xvdywgZGF0ZXMsIHZhcnMpLCAKICAgICAgICAgICAgICAgICAgZmV0Y2h2YXJzKGNvcmVfMjQ1X21pZCwgZGF0ZXMsIHZhcnMpLAogICAgICAgICAgICAgICAgICBmZXRjaHZhcnMoY29yZV8yNDVfaGlnaCwgZGF0ZXMsIHZhcnMpKQpgYGAKCgpgYGB7cn0KbnBwMF9sbWggJT4lIAogIGZpbHRlcih5ZWFyIDw9IDIwMjApICU+JSAKICBnZ3Bsb3QoKSArIAogIGdlb21fbGluZShhZXMoeWVhciwgdmFsdWUsIGNvbG9yID0gc2NlbmFyaW8pKSArIAogIGZhY2V0X3dyYXAoInZhcmlhYmxlIiwgc2NhbGVzID0gImZyZWUiKSArIAogIGxhYnModGl0bGUgPSAiSGlzdG9yY2FsIEhlY3RvciBvdXRwdXQgYXQgbG93LCBtZWRpdW0sIGFuZCBoaWdoIE5QUDAiLCB4ID0gTlVMTCwgeSA9IE5VTEwpCmBgYAoKYGBge3J9Cm5wcDBfbG1oICU+JSAKICBmaWx0ZXIoeWVhciA+PSAyMDIwKSAlPiUgCiAgZ2dwbG90KCkgKyAKICBnZW9tX2xpbmUoYWVzKHllYXIsIHZhbHVlLCBjb2xvciA9IHNjZW5hcmlvKSkgKyAKICBmYWNldF93cmFwKCJ2YXJpYWJsZSIsIHNjYWxlcyA9ICJmcmVlIikgKyAKICBsYWJzKHRpdGxlID0gIkZ1dHVyZSBIZWN0b3Igb3V0cHV0IGF0IGxvdywgbWVkaXVtLCBhbmQgaGlnaCBOUFAwIiwgeCA9IE5VTEwsIHkgPSBOVUxMKQpgYGAKYGBge3J9Cm5wcDBfbG1oICU+JSAKICMgZmlsdGVyKHllYXIgPj0gMjAyMCkgJT4lIAogIGdncGxvdCgpICsgCiAgZ2VvbV9saW5lKGFlcyh5ZWFyLCB2YWx1ZSwgY29sb3IgPSBzY2VuYXJpbykpICsgCiAgZmFjZXRfd3JhcCgidmFyaWFibGUiLCBzY2FsZXMgPSAiZnJlZSIpICsgCiAgbGFicyh0aXRsZSA9ICJIZWN0b3Igb3V0cHV0IGF0IGxvdywgbWVkaXVtLCBhbmQgaGlnaCBOUFAwIiwgeCA9IE5VTEwsIHkgPSBOVUxMKQpgYGAKClRoZXJlZm9yZSB3ZSBkZWNpZGVkIHRvIGxlYXZlIE5QUCBhcyBhIHZhbHVlIGZpeGVkIGZyb20gdGhlIGxpdGVyYXR1cmUuIAoKCiMgRUNTIGlzIGZpeGVkIGZyb20gdGhlIGxpdHRlcmF0dXJlISAKCk91ciB0aGlua2luZyBoZXJlIGlzIHRoYXQgZm9yIGRlZmF1bHQgSGVjdG9yIHdlIHNob3VsZCB1c2UgdGhlIGJlc3Qgc2NpZW50aWZpYyBjb25zZW5zdXMgYWthIHRoZSB2YWx1ZSBmcm9tIHRoZSBJUENDIEFSNiByZXBvcnQuIFdoZW4gdXNpbmcgSGVjdG9yIGFzIGFuIGVtdWxhdG9yIGZvciBzcGVjaWZpYyBFU01zIG1vZGVscyB0aGlzIGlzIGtleSB0dW5lLWFibGUgcGFyYW1ldGVyIHRoYXQgd291bGQgbmVlZCB0byBiZSBhZGp1c3RlZCBmb3IgSGVjdG9yIHRvIGFjdCBhcyBhbiBlbXVsYXRvci4gCgojIEFscGhhICYgVm9sc2NsIGFyZSBmaXhlZCAKCk9uY2UgYWdhaW4gSGVjdG9yIG1vZGVscyBhZXJvc29sIHJhZGlhdGl2ZSBmb3JjaW5nIGFuZCB1c2VzIHRoZSBwcmVzY3JpYmVkIHZvbGNhbmljIHJhZGlhdGl2ZSBmb3JjaW5nIHRoYXQgY29udHJpYnV0ZXMgZnVsbHkgdG8gSGVjdG9yJ3MgZ2xvYmFsIHRlbXBlcmF0dXJlIGNhbGN1bGF0aW9uLiBUaGlzIGlzIG1vcmUgb2YgYW4gYWRqdXN0YWJsZSBwYXJhbWV0ZXIgdG8gYWNjb3VudCBmb3IgdW5jZXJ0YWludHkgaW4gdm9sY2FuaWMgYW5kIGFlcm9zb2wgcmFkaWF0aXZlIGZvcmNpbmcuIFdlIGV4cGxvcmVkIHRoZSBpZGVhIG9mIGZpdHRpbmcgYWxwaGEgYW5kIHZvbHNjbCBob3dldmVyIGJ1dCBhbGxvd2luZyB0aG9zZSBwYXJhbWV0ZXJzIHRvIGJlIGZyZWUgZGlkIG5vdCBpbXByb3ZlIHRoZSBmaXRzLiBUaGVyZWZvcmUgd2Ugc2V0IHRob3NlIHRvIGJlIGZpeGVlIHZhbHVlcy4gCgoKCgoK