Wrapping Up Hector EPA RF task
01 October, 2021
1 Objective
library(dplyr)
library(tidyr)
library(ggplot2)
library(data.table)
library(hector)
theme_set(theme_bw())
DIR <- "/Users/dorh012/projects/2021/EPA_Hector" Load the AR6 RF results.
format_data <- function(path){
x <- gsub(pattern = ".csv", replacement = "", x = basename(path))
xx <- unlist(strsplit(x, split = "_"))
if (length(xx) == 3 ){
out <- as.data.frame(matrix(xx, nrow = 1, dimnames = list(NULL, c("output", "scn", "time"))))
}
if (length(xx) == 4 ){
out <- as.data.frame(matrix(xx, nrow = 1, dimnames = list(NULL, c("output", "scn", "time", "percent"))))
}
return(out)
}
get_cmip_data <- function(path){
dat <- read.csv(path, stringsAsFactors = FALSE)
info <- format_data(path)
dat %>%
melt(id.vars = c( "year"),
variable.name = "variable", value.name = "value") %>%
cbind(info) ->
out
return(out)
}
AR6_dir <- "/Users/dorh012/Documents/2021/Chapter-7/data_output"
files <- list.files(file.path(AR6_dir,"SSPs"), pattern = "ERF", full.names = TRUE)
files %>%
lapply(get_cmip_data) %>%
bind_rows() %>%
select(year, variable, value, scenario = scn, percent) %>%
as.data.table() %>%
filter(is.na(percent)) ->
ar6_data_og
ar6_names <- c("co2","ch4","n2o","o3","h2o_stratospheric","aerosol.cloud_interactions",
"volcanic","total", "HFC.134a","HFC.23","HFC.32","HFC.125","HFC.143a",
"HFC.227ea", "HFC.245fa", "SF6","CF4","CFC.11", "CFC.113",
"CFC.114", "CFC.115", "land_use")
hector_names <- c(RF_CO2(), RF_CH4(), RF_N2O(), RF_O3_TROP(), RF_H2O_STRAT(), RF_ACI(),
RF_VOL() ,RF_TOTAL(), RF_HFC134A(),RF_HFC23(), RF_HFC32(), RF_HFC125(),
RF_HFC143A(), RF_HFC227EA(), RF_HFC245FA(), RF_SF6(), RF_CF4(),
RF_CFC11(), RF_CFC113(), RF_CFC114(), RF_CFC115(), RF_T_ALBEDO())
ar6_data_og %>%
filter(variable %in% ar6_names) %>%
left_join(data.frame(variable = ar6_names,
hvar = hector_names)) %>%
select(year, variable = hvar, value, scenario) ->
ar6_results vars <- c(RF_CF4(), RF_C2F6(), RF_HFC23(), RF_HFC32(), RF_HFC4310(), RF_HFC125(), RF_HFC134A(),
RF_HFC143A(), RF_HFC227EA(), RF_HFC245FA(), RF_SF6(), RF_CFC11(), RF_CFC12(),
RF_CFC113(), RF_CFC114() , RF_CFC115(), RF_CCL4(), RF_CH3CCL3(), RF_HCFC22(),
RF_HCFC141B(), RF_HCFC142B(), RF_HALON1211(), RF_HALON1301(), RF_HALON2402(),
RF_CH3CL(), RF_CH3BR())Run Hector with different scenarios driven with different conditions.
run_hector_ssp <- function(ini){
vars <- c(RF_CO2(), RF_CH4(), RF_N2O(), RF_O3_TROP(), RF_H2O_STRAT(), RF_ACI(),
RF_VOL(), RF_TOTAL(), RF_HFC134A(), RF_HFC23(),
RF_HFC32(), RF_HFC125(), RF_HFC143A(), RF_HFC227EA(), RF_HFC245FA(),
RF_SF6(), RF_CF4(), RF_CFC11(),
RF_CFC113(), RF_CFC114(), RF_CFC115(), RF_BC(), RF_SO2(),
RF_OC(), RF_NH3(), ATMOSPHERIC_N2O(), GLOBAL_TEMP(), RF_T_ALBEDO(),
RF_CF4(), RF_C2F6(), RF_HFC23(), RF_HFC32(), RF_HFC4310(), RF_HFC125(), RF_HFC134A(),
RF_HFC143A(), RF_HFC227EA(), RF_HFC245FA(), RF_SF6(), RF_CFC11(), RF_CFC12(),
RF_CFC113(), RF_CFC114() , RF_CFC115(), RF_CCL4(), RF_CH3CCL3(), RF_HCFC22(),
RF_HCFC141B(), RF_HCFC142B(), RF_HALON1211(), RF_HALON1301(), RF_HALON2402(),
RF_CH3CL(), RF_CH3BR(), ATMOSPHERIC_CH4())
name <- gsub(pattern = "hector_|.ini", replacement = "", x = basename(ini))
core <- newcore(ini, name = name)
run(core, runtodate = 2100)
out <- fetchvars(core, 1745:2100, vars)
return(out)
}
list.files("/Users/dorh012/projects/2021/EPA_Hector/input", pattern = "ini", full.names = TRUE) %>%
lapply(run_hector_ssp) %>%
bind_rows() ->
hector_out_og
hector_out_og %>%
mutate(type = ifelse(grepl(pattern = "-tot", x = scenario), 'rf const', 'emiss')) %>%
mutate(type = ifelse(grepl(pattern = "constCO2", x = scenario), 'conc CO2', type)) %>%
mutate(type = ifelse(grepl(pattern = "constN2O", x = scenario), 'conc N2O', type)) %>%
mutate(type = ifelse(grepl(pattern = "constCH4", x = scenario), 'conc CH4', type)) %>%
mutate(type = ifelse(grepl(pattern = "constGHG", x = scenario),
'conc CO2, CH4, N2O', type)) %>%
mutate(type = ifelse(grepl(pattern = "missing", x = scenario), 'conc const + extra RF', type)) %>%
mutate(scenario = ifelse(grepl(pattern = "_constCH4|_constN2O|_constCO2|-tot|_constGHG|_missing",
x = scenario),
gsub(pattern = "_constCH4|_constN2O|_constCO2|-tot|_constGHG|_missing",
x = scenario, replacement = ""), scenario)) ->
hector_outFormat the data into a single data frame.
ar6_results$source <- "AR6"
ar6_results$type <- "AR6"
hector_out$source <- "Hv3"
comp_data <- bind_rows(ar6_results, hector_out)
scns <- c("ssp119", "ssp126", "ssp245", "ssp460", "ssp534-over")2 AR6 vs Hector Comparison
So it is possible to run Hector with prescribed RF. When we read in AR6 total RF as we expect Hector’s RF matches our expected value. Which is great and means that we have some estimate of global mean temperature.
types <- c("AR6", "rf const")
comp_data %>%
filter(variable %in% c(RF_TOTAL())) %>%
filter(year <= 2100) %>%
filter(scenario %in% scns) %>%
filter(type %in% types) %>%
mutate(source = paste0(source, " ", type)) %>%
ggplot(aes(year, value, color = source, linetype = source)) +
geom_line() +
facet_wrap("scenario", scales = "free") +
labs(y = "RF Wm-2", title = "Comparing Total RF")
Now there are several configurations we can run Hector in
- emission driven - as noted inn N2O, CO2, & CH4 sections there will be discrepancies due to issues with emissions –> concentrations step
- concentration driven - Hector is driven with N2O, CO2, CH4 concentrations to avoid the issues with the emission –> concentrations step
When driven with concentrations Hector’s total RF is pretty close to the AR6 total!
types <- c("AR6", "emiss", "conc CO2, CH4, N2O")
comp_data %>%
filter(variable == RF_TOTAL()) %>%
filter(year <= 2100) %>%
filter(scenario %in% scns) %>%
filter(type %in% types) %>%
mutate(source = paste0(source, " ", type)) %>%
ggplot(aes(year, value, color = source, linetype = source)) +
geom_line() +
facet_wrap("scenario", scales = "free") +
labs(y = "RF Wm-2", title = "Comparing Total RF")
What is the MSE between the AR6 & concentration driven Hector runs?
types <- c("AR6", "conc CO2, CH4, N2O")
comp_data %>%
filter(variable == RF_TOTAL()) %>%
filter(year <= 2100) %>%
filter(scenario %in% scns) %>%
filter(type %in% types) %>%
select(year, value, scenario, type) %>%
spread(type, value) %>%
na.omit() %>%
mutate(dif = (`conc CO2, CH4, N2O` - AR6)^2) %>%
group_by(scenario) %>%
summarise(value = mean(dif, na.rm = TRUE)) %>%
knitr::kable(digits = 3, format = "markdown")## `summarise()` ungrouping output (override with `.groups` argument)| scenario | value |
|---|---|
| ssp119 | 0.103 |
| ssp126 | 0.103 |
| ssp245 | 0.102 |
| ssp460 | 0.105 |
| ssp534-over | 0.105 |
How different is our estimated AR6 temperature from forcing driven Hector with the concentration driven Hector?
types <- c("rf const", "conc CO2, CH4, N2O")
comp_data %>%
filter(variable == GLOBAL_TEMP()) %>%
filter(year <= 2100) %>%
filter(scenario %in% scns) %>%
filter(type %in% types) %>%
mutate(type = if_else(type == "rf const", "estimated AR6", type)) %>%
ggplot(aes(year, value, color = type)) +
geom_line() +
facet_wrap("scenario", scales = "free") +
labs(y = "deg C", title = "Global Mean Temp")
types <- c("rf const", "conc CO2, CH4, N2O")
comp_data %>%
filter(variable == GLOBAL_TEMP()) %>%
filter(year <= 2100) %>%
filter(scenario %in% scns) %>%
filter(type %in% types) %>%
mutate(type = if_else(type == "rf const", "estimated AR6", type)) %>%
select(year, scenario, value, type) %>%
spread(type, value) %>%
mutate(error = (`conc CO2, CH4, N2O` - `estimated AR6`)^2) %>%
group_by(scenario) %>%
summarise(MSE = mean(error)) %>%
knitr::kable(digits = 3, format = "markdown")## `summarise()` ungrouping output (override with `.groups` argument)| scenario | MSE |
|---|---|
| ssp119 | 0.036 |
| ssp126 | 0.036 |
| ssp245 | 0.035 |
| ssp460 | 0.036 |
| ssp534-over | 0.036 |
The average MSE between our estimated AR6 temp & best Hector result are <= 0.01.
3 Larger Problems
3.1 Cloud Interactions
var <- RF_ACI()
comp_data %>%
filter(scenario %in% scns) %>%
filter(variable %in% var) %>%
filter(type %in% c("AR6", "emiss")) %>%
filter(year <= 2100) ->
cloud_rf
cloud_rf %>%
ggplot(aes(year, value, color = source, linetype = source)) +
geom_line() +
facet_wrap("scenario", scales = "free") +
labs(title = "Aerosol Cloud Interactions", y = "") 
The indirect cloud RF is our largest discrepancy!
cloud_rf %>%
filter(type %in% c("emiss", "AR6")) %>%
select(year, scenario, type, value) %>%
spread(type, value) %>%
na.omit %>%
mutate(dif = emiss - AR6) ->
diff_bw_cloudsHow much does this difference matter? In terms of total forcing?
diff_bw_clouds %>%
ggplot(aes(year, dif, color = scenario)) +
geom_line() +
labs(y = "Hector - AR6 Cloud Aerosol RF")
comp_data %>%
filter(variable == RF_TOTAL()) %>%
filter(type == "AR6") %>%
select(year, total = value, scenario) %>%
left_join(diff_bw_clouds, by = c("scenario", "year")) %>%
mutate(fraction = dif/total) %>%
filter(year <= 2100) %>%
ggplot(aes(year, fraction, color = scenario)) +
geom_line() +
coord_cartesian(ylim = c(-1, 1)) +
labs(title = "Difference in cloud RF/ Total RF",
y = "Frac")## Warning: Removed 2106 row(s) containing missing values (geom_path).
3.2 Not included in Hector
AR6 includes some RF that are not modeled by Hector such as sun, bc on snow or contrails.
ar6_data_og %>%
filter(variable %in% c("contrails", "bc_on_snow", "solar")) ->
missing_rf
scns## [1] "ssp119" "ssp126" "ssp245" "ssp460" "ssp534-over"missing_rf %>%
dplyr::filter(scenario %in% scns) %>%
filter(year <= 2100) %>%
group_by(year, scenario) %>%
summarise(value = sum(value)) ->
excluded_from_hector## `summarise()` regrouping output by 'year' (override with `.groups` argument)excluded_from_hector %>%
ggplot(aes(year, value, color = scenario)) +
geom_line() +
labs(title = "Forcings Excluded from Hector")
How much does this matter when compared with total RF?
ar6_data_og %>%
filter(variable == "total") %>%
select(year, total = value, scenario) ->
total_rf
excluded_from_hector %>%
left_join(total_rf, on = c("year", "scenario")) %>%
mutate(fraction = value / total) %>%
ggplot(aes(year, fraction, color = scenario)) +
geom_line() +
coord_cartesian(ylim = c(-1, 1)) +
geom_hline(yintercept = .3) +
geom_hline(yintercept = -.3) +
labs(caption = "horizontal +/- 0.3",
y = "missing from Hector / total RF")## Joining, by = c("year", "scenario")
It looks like it is a fairly small fraction of total RF, although it does add more inter annual variability. If we really want we include it we can pretty easily add a miscellaneous external rf.
4 Some mismatch (minor)
4.1 Other GHGs aka halocarbons
The AR6 has a category that they call “other GHGs” which I took to mean as the halocarbons. There wasn’t a fundamental change in halocarbon modeling with AR6 it continued to look like
RF = effective radiative forcing * change in concentrations
AR6 provided new values for the effective radiative forcing and life time so those are the changes we have made here.
Hector’s models 26 different hydrocarbons, that is the Hector values shown here, note that is higher than the expected value from AR6.
hvar <- c("FHFC134a", "FHFC23", "FHFC32", "FHFC125", "FHFC143a", "FHFC227ea",
"FHFC245fa", "FSF6", "FCF4", "FCFC11", "FCFC113", "FCFC114", "FCFC115",
"FC2F6", "FHFC4310", "FCFC12", "FCCl4", "FCH3CCl3", "FHCFC22",
"FHCFC141b", "FHCFC142b", "Fhalon1211", "Fhalon1301", "Fhalon2402",
"FCH3Cl", "FCH3Br")
hector_out %>%
dplyr::filter(variable %in% hvar) %>%
group_by(year, scenario, source, type) %>%
summarise(value = sum(value)) %>%
ungroup() ->
hector_other_ghgs## `summarise()` regrouping output by 'year', 'scenario', 'source' (override with `.groups` argument)ar6_data_og %>%
filter(scenario %in% scns & variable == "other_wmghg") %>%
mutate(source = "AR6") %>%
bind_rows(hector_other_ghgs) %>%
filter(scenario %in% scns) %>%
filter(year <= 2100) ->
halocarbon_compare
halocarbon_compare %>%
ggplot(aes(year, value, color = source)) +
geom_line() +
facet_wrap("scenario", scales = "free") +
labs(y = "", title = "other ghg")
How different are these values?
halocarbon_compare %>%
select(year, value, scenario, source) %>%
distinct() %>%
spread(source, value) %>%
na.omit() %>%
mutate(dif = Hv3 - AR6) ->
halocarbon_diff halocarbon_diff %>%
ggplot(aes(year, dif, color = scenario)) +
geom_line() +
labs(title = "Difference in halocarbon RF",
y = "Hv3 - AR6")
halocarbon_diff %>%
group_by(scenario) %>%
summarise(max = max(dif),
mean = mean(dif),
sd = sd(dif)) %>%
knitr::kable(digits = 3, markdown = "markdown")## `summarise()` ungrouping output (override with `.groups` argument)| scenario | max | mean | sd |
|---|---|---|---|
| ssp119 | 0.117 | 0.032 | 0.041 |
| ssp126 | 0.126 | 0.034 | 0.045 |
| ssp245 | 0.166 | 0.048 | 0.064 |
| ssp460 | 0.195 | 0.050 | 0.069 |
| ssp534-over | 0.206 | 0.048 | 0.066 |
How much does <.2 RF represent of the total annual RF value?
halocarbon_diff %>%
select(year, halocarbon = dif, scenario) %>%
left_join(hector_out %>%
filter(variable == RF_TOTAL())) %>%
filter(source == "Hv3" & type == "emiss") %>%
na.omit() %>%
mutate(frac = halocarbon / value,
frac = replace_na(frac, replace = 0)) %>%
ggplot(aes(year, frac, color = scenario)) +
geom_line() +
facet_wrap("scenario") +
coord_cartesian(ylim=c(0,.5))## Joining, by = c("year", "scenario")
It looks like it is a relatively small fraction of total RF.
4.2 O3
So I did not make any changes to Hector’s troposphere O3 RF modeling, while the AR6 did provide an update to O3 because the AR6 O3 RF includes the effects from ozone depleting halocarbons & would add complexity to how we model halocarbons. Since it is not as straightforward would like to consider all our options.
var <- RF_O3_TROP()
comp_data %>%
filter(scenario %in% scns) %>%
filter(variable %in% var) %>%
filter(type %in% c("AR6", "emiss")) %>%
filter(year <= 2100) %>%
ggplot(aes(year, value, color = source, linetype = source)) +
geom_line() +
facet_wrap("scenario", scales = "free") +
labs(title = var, y = "") 
4.3 N2O
Like CO2 & CH4 AR6 changed the N2O RF calculation, since this calculation depends on CO2, CH4, and N2O which we know these concentrations do have some discrepancies all three had to be constrained at once. When we do so the N2O RF better represents the expected RF from AR6. Emission driven Hector N2O RF is more different than the AR6 RF and that is due to differences in concentrations which is caused by the natural N2O emissions.
var <- RF_N2O()
comp_data %>%
filter(scenario %in% scns) %>%
filter(variable %in% var) %>%
filter(type %in% c("AR6", "emiss", "conc CO2, CH4, N2O")) %>%
filter(year <= 2100) %>%
ggplot(aes(year, value, color = type, linetype = type)) +
geom_line(size = 1) +
facet_wrap("scenario", scales = "free" ) +
labs(title = var, y = "") 
N2O Concentrations.
var <- ATMOSPHERIC_N2O()
comp_data %>%
filter(scenario %in% scns) %>%
filter(variable %in% var) %>%
filter(type %in% c("AR6", "emiss", "conc CO2, CH4, N2O")) %>%
filter(year <= 2100) %>%
ggplot(aes(year, value, color = type, linetype = type)) +
geom_line(size = 1) +
facet_wrap("scenario", scales = "free" ) +
labs(title = var, y = "") 
Look at how sensitive N2O RF is to the natural emissions, changing the N2O natural emissions can really change the shape of the N2O RF & N2O concentrations. However this has a somewhat small effect on the temperature output.
vars <- c(ATMOSPHERIC_N2O(), RF_N2O(), NAT_EMISSIONS_N2O(), GLOBAL_TEMP())
yrs <- 1745:2100
ini <- file.path(DIR, "input", "hector_ssp245.ini")
core <- newcore(ini, name = "default N2O emiss")
run(core) ## Hector core: default N2O emiss
## Start date: 1745
## End date: 2300
## Current date: 2300
## Input file: /Users/dorh012/Documents/2021/EPA_Hector/input/hector_ssp245.iniout1 <- fetchvars(core, yrs, vars)
default_n2o_emiss <- fetchvars(core, 1745:2100, NAT_EMISSIONS_N2O())
vals <- 9.9 * default_n2o_emiss$value / default_n2o_emiss$value
setvar(core, default_n2o_emiss$year, values = vals, var = NAT_EMISSIONS_N2O(), unit = "Tg N")
reset(core)## Hector core: default N2O emiss
## Start date: 1745
## End date: 2300
## Current date: 1745
## Input file: /Users/dorh012/Documents/2021/EPA_Hector/input/hector_ssp245.inirun(core)## Hector core: default N2O emiss
## Start date: 1745
## End date: 2300
## Current date: 2300
## Input file: /Users/dorh012/Documents/2021/EPA_Hector/input/hector_ssp245.iniout2 <- fetchvars(core, yrs, vars)
out2$scenario <- "9.9 N2O emiss"comp_data %>%
filter(variable == RF_N2O(), source == "AR6") %>%
filter(scenario %in% "ssp245") %>%
mutate(scenario = "ar6 ssp245") ->
ar6_data_n2o
out1 %>%
bind_rows(out2) %>%
bind_rows(ar6_data_n2o) %>%
filter(year <= 2100) %>%
ggplot(aes(year, value, color = scenario, linetype = scenario)) +
geom_line(size = 1) +
facet_wrap("variable", scales = "free") +
labs(y = NULL)
5 Some (minor) mismatch
5.1 CH4
AR6 completely changed the functions used to calculate CO2 RF, here is Hector driven with emissions and with CH4 concentrations constrained. When CH4 concentrations are directly read in Hector CH4 RF matches the expected AR6 RF. I suspect that the discrepancy is related to the fact that the CH4 concentrations (see below).
comp_data %>%
filter(scenario %in% scns) %>%
filter(variable %in% RF_CH4()) %>%
filter(type %in% c("AR6", "emiss", "conc CH4")) %>%
filter(year <= 2100) %>%
ggplot(aes(year, value, color = type, linetype = type)) +
geom_line(size = 1) +
facet_wrap("scenario", scales = "free" ) +
labs(title = RF_CH4(), y = "") 
var <- ATMOSPHERIC_CH4()
comp_data %>%
filter(scenario %in% scns) %>%
filter(variable %in% var) %>%
filter(type %in% c("AR6", "emiss", "conc CH4")) %>%
filter(year <= 2100) %>%
ggplot(aes(year, value, color = type, linetype = type)) +
geom_line(size = 1) +
facet_wrap("scenario", scales = "free" ) +
labs(title = var, y = "") 
How much of an affect does the change in CH4 concentrations have on global mean temperature? For the most part the difference between the global mean temp output are less than 4%.
var <- GLOBAL_TEMP()
comp_data %>%
filter(scenario %in% scns) %>%
filter(variable %in% var) %>%
filter(type %in% c("emiss", "conc CH4")) %>%
ggplot() +
geom_line((aes(year, value, color = type, linetype = type)), size = 1) +
facet_wrap("scenario", scales = "free") +
labs(title = var,
y = "")
var <- GLOBAL_TEMP()
comp_data %>%
filter(scenario %in% scns) %>%
filter(variable %in% var) %>%
filter(type %in% c("emiss", "conc CH4", "rf const")) %>%
select(year, value, scenario, type) %>%
spread(type, value) %>%
mutate(diff = (`conc CH4` - `emiss`)/emiss) %>%
ggplot() +
geom_line((aes(year, diff)), size = 1) +
geom_hline(yintercept = 4, color = "pink") +
geom_hline(yintercept = -4, color = "pink") +
facet_wrap("scenario") +
coord_cartesian(ylim = c(-5,10)) +
labs(title = "Percent Diff in Global Mean Temp",
y = "%",
caption = "horizontal lines +/-4%")
5.2 H2O
AR6 did not update the functions used to model stratospheric H2O, in Hector 2.5 we modeled stratospheric H2O as a fraction of the “pure” CH4 RF as they did in Joos et al 2001. However because with the updates to AR6 CH4 RF that is no longer an option. Instead we model stratospheric H2O RF proportional to the change in CH4 concentrations to be consistent with AR6.
var <- RF_H2O_STRAT()
comp_data %>%
filter(scenario %in% scns) %>%
filter(variable %in% var) %>%
filter(type %in% c("AR6", "conc CH4", "emiss")) %>%
filter(year <= 2100) %>%
ggplot(aes(year, value, color = type, linetype = type)) +
geom_line(size = 1) +
facet_wrap("scenario", scales = "free") +
labs(title = var, y = "") 
6 Looks Good
6.1 CO2
AR6 completely changed the functions used to calculate CO2 RF, here is Hector driven with emissions and with CO2 concentrations constrained. Without the carbon-cycle temperature feedbacks turned on (CO2 concentrations constrained) Hector v3 is able to replicated the AR6 CO2 RF.
var <- RF_CO2()
comp_data %>%
filter(scenario %in% scns) %>%
filter(variable %in% var) %>%
filter(type %in% c("AR6", "emiss")) %>%
filter(year <= 2100) %>%
ggplot(aes(year, value, color = type, linetype = type)) +
geom_line(size = 1) +
facet_wrap("scenario", scales = "free") +
labs(title = var, y = "") 
6.2 Volcanic
var <- RF_VOL()
comp_data %>%
filter(scenario %in% scns) %>%
filter(variable %in% var) %>%
filter(type %in% c("AR6", "emiss")) %>%
filter(year <= 2100) %>%
ggplot(aes(year, value, color = source, linetype = source)) +
geom_line(size = 1) +
facet_wrap("scenario", scales = "free") +
labs(title = var, y = "") 
6.3 Land use
var <- RF_T_ALBEDO()
comp_data %>%
filter(scenario %in% scns) %>%
filter(variable %in% var) %>%
filter(type %in% c("AR6", "emiss")) %>%
filter(year <= 2100) %>%
ggplot(aes(year, value, color = source, linetype = source)) +
geom_line(size = 1) +
facet_wrap("scenario", scales = "free") +
labs(title = var, y = "") 
6.4 Direct Aerosol
hector_out %>%
filter(variable %in% c(RF_BC(), RF_SO2(), RF_OC(), RF_NH3())) %>%
group_by(year, scenario, type) %>%
summarise(value = sum(value)) %>%
ungroup %>%
mutate(variable = "aerosol.radiation_interactions") %>%
mutate(source = "H3") ->
hector_dir## `summarise()` regrouping output by 'year', 'scenario' (override with `.groups` argument)ar6_data_og %>%
filter(variable == "aerosol.radiation_interactions") %>%
mutate(source = "AR6", type = "AR6") %>%
bind_rows(hector_dir) %>%
filter(scenario %in% scns) %>%
filter(type %in% c("AR6", "emiss")) %>%
filter(year <= 2100) %>%
ggplot(aes(year, value, color = source, linetype = source)) +
geom_line() +
facet_wrap("scenario", scales = "free") +
labs(title = "Direct Aerosol", y = "")
7 Unresovled problems
- The additional forcing aka (contrails, black carbon on snow, solar) do we want the ability to easily read in miscellaneous rf?
- What do we want to do about O3RF?
- Update the CH4, N2O, SO2 natural emissions, What is going on with the SO2 natural emissions?
- How are we feeling about the aerosol cloud interactions?