Here I compare the L1 outputs (CEDs + the burning emissions from RCMIP) with the default emissions passed into Hector.

# Packages used here 
library(dplyr)
library(ggplot2)
library(hector)
library(kableExtra)
library(knitr)
library(tidyr)

# Using relative paths! 
BASE_DIR <- here::here()

# Set the default plot optics
theme_set(theme_bw())

# Import CEDs data! 
file.path(BASE_DIR, "data", "L1", "L1.incomplete_ceds_hector.csv") %>%  
    read.csv() %>% 
    mutate(source = "CEDs +") -> 
    ceds_emiss

# Run Hector emissions to compare with CEDS 
scns <- c("ssp119", "ssp245", "ssp585")
vars <- unique(ceds_emiss$variable)
yrs <- min(ceds_emiss$year):(max(ceds_emiss$year) + 20)

lapply(scns, function(scn){
    ini_file <- list.files(system.file(package = "hector", "input"), 
                           pattern = scn, full.names = TRUE)
    hc <- newcore(ini_file, name = scn)
    run(hc)
    fetchvars(core = hc, dates = yrs, vars = vars) %>%  
        mutate(source = "hector default") ->
        out

    return(out)
}) %>% 
    do.call(what = "rbind") -> 
    hector_emiss

hector_emiss %>% 
    filter(scenario == "ssp245") %>% 
    filter(year <= 2014) %>% 
    mutate(scenario = "historical") -> 
    hector_historical

all_emiss <- bind_rows(ceds_emiss, hector_emiss, hector_historical)

Relative Differences

all_emiss %>%  
    filter(scenario == "historical") %>% 
    mutate(source = if_else(source == "hector default", "H", "C")) %>% 
    select(year, value, variable, source) %>% 
    spread(source, value) %>% 
    na.omit %>% 
    mutate(diff = (C-H)/H) -> 
    differences

What is the relative difference between the default Hector and CEDS emissions? This will help us get a sense of what emissions have been converted properly and that the CEDS & RCMIP emissions are aggregated properly.

differences %>% 
    group_by(variable) %>% 
    summarise(min = min(diff), 
              mean = mean(diff), 
              max = max(diff), 
              sd = sd(diff)) %>% 
  kable(digits = 4, format = "html", caption = "Summary Stats for Relative Difference by Emission") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE,
    position = "center"
  )

Summary Stats for Relative Difference by Emission
variable	min	mean	max	sd
BC_emissions	-0.2047	0.0094	0.1141	0.0492
CH4_emissions	-0.0403	0.0174	0.0811	0.0308
CO_emissions	-0.1125	-0.0121	0.0182	0.0230
N2O_emissions	-0.2702	-0.1920	-0.0827	0.0531
NH3_emissions	-0.0797	-0.0525	-0.0020	0.0191
NMVOC_emissions	-0.1330	-0.0154	0.0244	0.0498
NOX_emissions	-0.5149	-0.3387	-0.0504	0.1873
OC_emissions	-0.1838	-0.0103	0.0250	0.0372
SO2_emissions	-0.1338	0.0066	0.1119	0.0516
ffi_emissions	-0.1058	-0.0239	0.0987	0.0568

ggplot(differences) + 
    geom_boxplot(aes(variable, diff)) + 
    labs(y = "(CEDS - Hector) / Hector", x = NULL) + 
      theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

BC Emissions

var <- EMISSIONS_BC()

all_emiss %>% 
    filter(variable == var) %>% 
    filter(year <= 2030) %>% 
    ggplot(aes(year, value, color = source)) + 
    geom_line() +
    facet_wrap("scenario") + 
    labs(title = var, y = getunits(var))

CH4 Emissions

CEDs+ = only CEDS, there were not additional CH4 burning emissions

var <- EMISSIONS_CH4()

all_emiss %>% 
    filter(variable == var) %>% 
    filter(year %in% 1950:2022) %>% 
    ggplot(aes(year, value, color = source)) + 
    geom_line() +
    facet_wrap("scenario", scales = "free") + 
    labs(title = var, y = getunits(var))

CO Emissions

var <- EMISSIONS_CO()

all_emiss %>% 
    filter(variable == var) %>% 
    filter(year <= 2022) %>% 
    ggplot(aes(year, value, color = source)) + 
    geom_line() +
    facet_wrap("scenario", scales = "free") + 
    labs(title = var, y = getunits(var))

CO2 Emissions (ffi only)

CEDs+ is actually only CEDs

var <- FFI_EMISSIONS()

all_emiss %>% 
    filter(variable == var) %>% 
    filter(year <= 2022) %>% 
    ggplot(aes(year, value, color = source)) + 
    geom_line() +
    facet_wrap("scenario", scales = "free") + 
    labs(title = var, y = getunits(var))

N2O Emissions

CEDs + is only CEDs no emissions from RCMIP

var <- EMISSIONS_N2O()

all_emiss %>% 
    filter(variable == var) %>% 
    filter(year %in% 1900:2030) %>% 
    ggplot(aes(year, value, color = source)) + 
    geom_line() +
    facet_wrap("scenario", scales = "free") + 
    labs(title = var, y = getunits(var))

NH3 Emissions

var <- EMISSIONS_NH3()

all_emiss %>% 
    filter(variable == var) %>% 
    filter(year <= 2022) %>% 
    ggplot(aes(year, value, color = source)) + 
    geom_line() +
    facet_wrap("scenario", scales = "free") + 
    labs(title = var, y = getunits(var))

NMVOC Emissions

var <- EMISSIONS_NMVOC()

all_emiss %>% 
    filter(variable == var) %>% 
    filter(year <= 2022) %>% 
    ggplot(aes(year, value, color = source)) + 
    geom_line() +
    facet_wrap("scenario", scales = "free") + 
    labs(title = var, y = getunits(var))

OC Emissions

var <- EMISSIONS_OC()

all_emiss %>% 
    filter(variable == var) %>% 
    filter(year <= 2022) %>% 
    ggplot(aes(year, value, color = source)) + 
    geom_line() +
    facet_wrap("scenario", scales = "free") + 
    labs(title = var, y = getunits(var))

NOx Emissions

var <- EMISSIONS_NOX()

all_emiss %>% 
    filter(variable == var) %>% 
    filter(year <= 2022) %>% 
    ggplot(aes(year, value, color = source)) + 
    geom_line() +
    facet_wrap("scenario", scales = "free") + 
    labs(title = var, y = getunits(var))

S2O Emissions

var <- EMISSIONS_SO2()

all_emiss %>% 
    filter(variable == var) %>% 
    filter(year <= 2022) %>% 
    ggplot(aes(year, value, color = source)) + 
    geom_line() +
    facet_wrap("scenario", scales = "free") + 
    labs(title = var, y = getunits(var))

Benchmarking L1 Outputs

Kalyn Dorheim

September 11, 2024