Hector3_GCB

Project: Hector v3 manuscript (calibration investigation)

Objective: take a look at the Global Carbon Budget and compare it to Hector inputs/outputs. While we don’t think there will be 100% match between everything we would ideally want Hector in the same ball park range as their results or understand how and why Hector & GCB diverge from one another.

Main Take Away: My main take away from this exercise is that while the inputs for LUC for the historical do different this difference is not large enough to be the cause of the discrepancy between the observed \(CO_2\) concentrations and Hector output during the historical period.

Set Up

# Load required packages, note that Hector is not included in this list because need the v3.
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(ggplot2)
library(tidyr)
library(readxl)

# 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/Documents/Hector-Versions/v3/hector"
devtools::load_all(HECTOR_BASE)

## ℹ Loading hector

# Vis setting
theme_set(theme_bw())

Import Data

# 1. Import Data -----------------------------------------------------------------------------------------
# Load of the historical data 
GCB_data <- read_xlsx(here::here("misc_data", "Global_Carbon_Budget_2021v1.0.xlsx"), 
                      sheet = "Global Carbon Budget") 
GCB_data_historical <- read_xlsx(here::here("misc_data", "Global_Carbon_Budget_2021v1.0.xlsx"), 
                                 sheet = "Historical Budget", skip = 15) 
GCB_land_sink <- read_xlsx(here::here("misc_data", "Global_Carbon_Budget_2021v1.0.xlsx"), 
                           sheet = "Terrestrial Sink", skip = 24) 

## New names:
## * `` -> ...3
## * `` -> ...21

# Format the data so easier to work with. 
names(GCB_data) <- gsub(pattern = " |-", replacement = "_", x = names(GCB_data))
names(GCB_data_historical) <- gsub(pattern = " |-", replacement = "_", x = names(GCB_data_historical))
names(GCB_land_sink) <- gsub(pattern = " |-", replacement = "_", x = names(GCB_land_sink))


# Historical hector results. 
dates <- 1750:2020
vars <- c(FFI_EMISSIONS(), LUC_EMISSIONS(), GLOBAL_TAS(), CONCENTRATIONS_CO2(), NBP(), NPP())

# Emission driven run. 
ini <- file.path(HECTOR_BASE, "inst", "input", "hector_ssp245.ini")
core <- newcore(ini)
run(core)

## Hector core: Unnamed Hector core
## Start date:  1745
## End date:    2300
## Current date:    2300
## Input file:  /Users/dorh012/Documents/Hector-Versions/v3/hector/inst/input/hector_ssp245.ini

out1 <- fetchvars(core, dates, vars)
out1$scenario <- "hist emission"

# Concentration driven runs. 
out2 <- read.csv(file.path(HECTOR_BASE, "inst", "input", "tables", "ssp434_emiss-constraints_rf.csv"), comment.char = ";")
setvar(core, out2$Date, var = CO2_CONSTRAIN(), values = out2$CO2_constrain, unit = getunits(CO2_CONSTRAIN()))
reset(core)

## Hector core: Unnamed Hector core
## Start date:  1745
## End date:    2300
## Current date:    1745
## Input file:  /Users/dorh012/Documents/Hector-Versions/v3/hector/inst/input/hector_ssp245.ini

run(core)

## Hector core: Unnamed Hector core
## Start date:  1745
## End date:    2300
## Current date:    2300
## Input file:  /Users/dorh012/Documents/Hector-Versions/v3/hector/inst/input/hector_ssp245.ini

out2 <- fetchvars(core, dates, vars)
out2$scenario <- "hist conc"

hector_data <- rbind(out1, out2)

Compare Inputs

LUC Emissions

# 2. Compare LUC emissions ---------------------------------------------------------------------------------------
hector_data %>% 
  filter(variable == LUC_EMISSIONS() & scenario == "hist emission") %>% 
  filter(year %in% GCB_data_historical$Year) -> 
  hector_luc

ggplot() + 
  geom_line(data = GCB_data_historical, aes(Year, land_use_change_emissions, color = "GCB")) + 
  geom_line(data = hector_luc, aes(year, value, color = "hector")) + 
  labs(title = "LUC emissions")

## Warning: Removed 100 row(s) containing missing values (geom_path).

FFI Emissions

hector_data %>% 
  filter(variable == FFI_EMISSIONS() & scenario == "hist emission") %>% 
  filter(year %in% GCB_data_historical$Year) -> 
  hector_luc

GCB_data_historical %>% 
  mutate(cement_carbonation_sink = if_else(is.na(cement_carbonation_sink), 0, cement_carbonation_sink)) %>% 
  mutate(ffi = fossil_emissions_excluding_carbonation + cement_carbonation_sink) -> 
  GCB_to_plot

ggplot() + 
  geom_line(data = GCB_to_plot, aes(Year, ffi, color = "GCB")) + 
  geom_line(data = hector_luc, aes(year, value, color = "hector"), linetype = 2) + 
  labs(title = "GCB Fossil Fuel Emissions")

GCB Driven Hector

GCB_data_historical %>% 
  filter(!is.na(land_use_change_emissions)) -> 
  new_inputs

ini <- file.path(HECTOR_BASE, "inst", "input", "hector_ssp245.ini")
core <- newcore(ini)
setvar(core, dates = new_inputs$Year, values = new_inputs$land_use_change_emissions, 
       var = LUC_EMISSIONS(), unit = "Pg C/yr")
reset(core)

## Hector core: Unnamed Hector core
## Start date:  1745
## End date:    2300
## Current date:    1745
## Input file:  /Users/dorh012/Documents/Hector-Versions/v3/hector/inst/input/hector_ssp245.ini

run(core)

## Hector core: Unnamed Hector core
## Start date:  1745
## End date:    2300
## Current date:    2300
## Input file:  /Users/dorh012/Documents/Hector-Versions/v3/hector/inst/input/hector_ssp245.ini

out3 <- fetchvars(core, dates, vars)
out3$scenario <- "GCB driven"

hector_data <- rbind(out3, hector_data)

hector_data %>% 
  filter(variable == LUC_EMISSIONS() & year <= 2020) %>% 
  ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
  geom_line() + 
  labs(title = "LUC Emissions", 
       y = NULL, x = NULL, 
       caption = "Comparing Hector intputs with GCB inputs this fig shows that the Hector cores were set up correctly.")

hector_data %>% 
  filter(variable == CONCENTRATIONS_CO2() & year <= 2020) %>% 
  ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
  geom_line()

hector_data %>% 
  filter(variable == NBP() & year <= 2020) %>% 
  ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
  geom_line() + 
  labs(title = "Hector NBP", 
       x = NULL, y = NULL,
       caption = "Comparison of Hector NBP output from 3 different Hector runs, two emission driven 
       runs with LUC emissionns from different sources and then CO2 concentration driven Hector which would 
       be our approximation of historical NBP")

hector_data %>% 
  filter(variable == GLOBAL_TAS() & year <= 2020) %>% 
  ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
  geom_line() + 
  labs(title = "Hector Global Tas", 
       x = NULL, y = NULL, 
       caption = "Comparison of Hector temp output from 3 different Hector runs, two emission driven 
       runs with LUC emissionns from different sources and then CO2 concentration driven Hector which would 
       be our approximation of historical tas")

hector_data %>% 
  filter(variable == CONCENTRATIONS_CO2()) %>% 
  split(f = .$scenario) %>% 
  lapply(function(x){
    value <- c(diff(x[["value"]]), NA)
    x[["value"]] <- value
    return(x)
  }) %>%
  bind_rows -> 
  hector_growth

# The units for the ppmv CO2 -- could it be a difference in units? 

co2_obs <- read.csv(here::here("data", "co2_annmean_mlo.csv"), comment.char = "#")
co2_obs[["value"]] <- c(diff(co2_obs$mean), NA)


ggplot() + 
  geom_line(data = hector_growth, aes(year, value, color = scenario)) + 
  geom_line(data = GCB_data_historical, aes(Year, atmospheric_growth, color = "GCB")) + 
  geom_line(data = co2_obs, aes(year, value, color = "obs")) + 
  labs(title = "Amtospheric Growth", 
       x = NULL, y = NULL, 
       caption = "Comparison of Hector chagne in [CO2] output from 3 different Hector runs & reported by GCB. 
       What is interesting is the all three Hector runs are within one another which it good to know and tells
       us that the relative input time series for the emission driven runs are good in that they fall within 
       the rate of chagne in the observational record.")

## Warning: Removed 3 row(s) containing missing values (geom_path).

## Warning: Removed 1 row(s) containing missing values (geom_path).