Compare RCMIP & AR6 historical

1 Set Up

library(data.table)
library(dplyr)
library(ggplot2)
library(hector)
library(tidyr)

BASE_DIR <- here::here()

theme_set(theme_bw())

1.1 Objective

Here we want to compare the concentrations from RCMIP & the AR6 driven Hector runs. Does using the historical AR6 concentrations improve how Hector RF compares with with the the AR6 RF?

2 Run Hector

Run historical with RCMIP inputs.

GHG_cmip5_ini <- "/Users/dorh012/Documents/2021/2021-hector-ar6rf/input/hector_ssp585-GHG.ini"
core <- newcore(GHG_cmip5_ini)  
run(core, runtodate = 2015)  

Hector core:    Unnamed Hector core
Start date: 1745
End date:   2300
Current date:   2015
Input file: /Users/dorh012/Documents/2021/2021-hector-ar6rf/input/hector_ssp585-GHG.ini

ghg_rcmip_out <- fetchvars(core, 1850:2015, vars = c(GLOBAL_TEMP(), RF_CO2(), RF_N2O(), RF_CH4(), RF_TOTAL(),
                                    ATMOSPHERIC_CO2(), ATMOSPHERIC_N2O(), ATMOSPHERIC_CH4()))  
ghg_rcmip_out$scenario <- "conc hector rcmip"
  

emiss_cmip5_ini <- "/Users/dorh012/Documents/2021/2021-hector-ar6rf/input/hector_ssp585.ini"
core <- newcore(emiss_cmip5_ini)  
run(core, runtodate = 2015)  

Hector core:    Unnamed Hector core
Start date: 1745
End date:   2300
Current date:   2015
Input file: /Users/dorh012/Documents/2021/2021-hector-ar6rf/input/hector_ssp585.ini

emiss_rcmip_out <- fetchvars(core, 1850:2015, vars = c(GLOBAL_TEMP(), RF_CO2(), RF_N2O(), RF_CH4(), RF_TOTAL(), 
                                    ATMOSPHERIC_CO2(), ATMOSPHERIC_N2O(), ATMOSPHERIC_CH4()))  
emiss_rcmip_out$scenario <- "emiss hector rcmip"

Ru historical with AR6 inputs.

emiss_ini <- "/Users/dorh012/Documents/2021/2021-hector-ar6rf/input/hector_ssp585.ini"
core <- newcore(emiss_ini)  
ar6_conc <- read.csv(here::here("ar6", "ar6-table_A3.1_historical_abundances.csv"), skip = 1)
setvar(core, ar6_conc$YYYY, var = CO2_CONSTRAIN(), values = ar6_conc$CO2, unit =  getunits(CO2_CONSTRAIN()))
setvar(core, ar6_conc$YYYY, var = CH4_CONSTRAIN(), values = ar6_conc$CH4, unit =  getunits(CH4_CONSTRAIN()))
setvar(core, ar6_conc$YYYY, var = N2O_CONSTRAIN(), values = ar6_conc$N2O, unit =  getunits(N2O_CONSTRAIN()))
reset(core)

Hector core:    Unnamed Hector core
Start date: 1745
End date:   2300
Current date:   1745
Input file: /Users/dorh012/Documents/2021/2021-hector-ar6rf/input/hector_ssp585.ini

run(core, runtodate = 2015)  

Hector core:    Unnamed Hector core
Start date: 1745
End date:   2300
Current date:   2015
Input file: /Users/dorh012/Documents/2021/2021-hector-ar6rf/input/hector_ssp585.ini

ar6_out <- fetchvars(core, 1850:2015, vars = c(GLOBAL_TEMP(), RF_CO2(), RF_N2O(), RF_CH4(), RF_TOTAL(),
                                    ATMOSPHERIC_CO2(), ATMOSPHERIC_N2O(), ATMOSPHERIC_CH4()))  
ar6_out$scenario <- "ar6 hector"

2.1 Import the effective radiative forcing

ERF_data <- read.csv(here::here("SSPs", "ERF_ssp434_1750-2500.csv"))
ERF_data %>% 
  select(year, FN2O = n2o, FCO2 = co2, FCH4 = ch4, Ftot = total) %>% 
  as.data.table() %>% 
  melt(id.vars = "year") %>% 
  mutate(scenario = "ar6 results") %>% 
  filter(year %in% ar6_out$year) -> 
  ar6_rf

hector_out <- bind_rows(ghg_rcmip_out, emiss_rcmip_out, ar6_out, ar6_rf)

3 Compare Results

3.1 CO2

How different is the CO2 RF form one another? What does that mean for the CO2 RF?

hector_out %>%  
  filter(grepl(pattern = "CO2|Ca", x = variable)) %>% 
  ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
  geom_line(size = 1.5, alpha = 0.5) + 
  facet_wrap("variable", scales = "free") + 
  labs(y = NULL, x = NULL) 

3.2 N2O

hector_out %>%  
  filter(grepl(pattern = "N2O", x = variable)) %>% 
  ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
  geom_line(size = 1.5, alpha = 0.5) + 
  facet_wrap("variable", scales = "free") + 
  labs(y = NULL, x = NULL) 

3.3 CH4

hector_out %>%  
  filter(grepl(pattern = "CH4", x = variable)) %>% 
  ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
  geom_line(size = 1.5, alpha = 0.5) + 
  facet_wrap("variable", scales = "free") + 
  labs(y = NULL, x = NULL) 

3.4 Global Meann Temp

hector_out %>%  
  filter(variable %in% c(GLOBAL_TEMP(), RF_TOTAL())) %>% 
  ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
  geom_line(size = 1.5, alpha = 0.5) + 
  facet_wrap("variable", scales = "free") + 
  labs(y = NULL, x = NULL) 

What is the difference between the ar6 hector and the conc hector rcmip temperature? Does this make a difference in the long term?

hector_out %>%  
  filter(variable %in% c(GLOBAL_TEMP(), RF_TOTAL())) %>% 
  select(-units) %>% 
  as.data.table %>% 
  dcast(year + variable  ~ scenario) %>%  
  mutate(diff = 100 * (`conc hector rcmip` - `ar6 hector`) / `ar6 hector`) %>%  
  ggplot(aes(year, diff)) + 
  geom_hline(yintercept = -5, size = 1, color = "red", linetype = 2, alpha = 0.5) + 
  geom_hline(yintercept = 5, size = 1, color = "red", linetype = 2, alpha = 0.5) + 
  geom_line(size = 1.5, alpha = 0.5) + 
  facet_wrap("variable", scales = "free") + 
  labs(y = NULL, x = NULL, title = "% Difference between conc hector rcmip & ar6 hector", 
       caption = "red lines at +/- 5%") + 
  coord_cartesian(ylim = c(-50, 50))

4 Take