EPA Hector Task: Updating RF to reflect AR6
03 September, 2021
1 Objective
Explore the effects associated with changing RF associated Hector paramters from AR5 to AR6.
library(hector) # Note this is a locally modified branch of hector that allows for rho to be read in via R
library(dplyr)
library(tidyr)
library(ggplot2)2 Halocarbons
Updating the halocarbons is the most straight forward task. We can read in the new life times (tau) and radiative forcing efficienies (rho) via the ini file. AR6 adds an additional step, adjusting CFC RF by some value to account for tropospheric interactions. Here I do that adjustment after running Hector. There are two different ways this could be addressed in Hector
2.1 Unresolved
- incorperating trop adjustments
- add an RF_adjs factor for all of the halocarbons, for most of the speicies this this will be set equal to 0 (probably easier to implent & would allow those paramters to be updated & tested in uncertainty analyses).
- add the the trop adjustment factor to only the specific halocarbons
- leave out the trop adjustments to the CFCs
- HFC4310_halocarbon is a halocarbon read into Hector but it is not listed in AR6 so we cannot update the parameters, where did this halo carbon come from? Do we know if it has been grouped into something else as part of AR6?
- are the trop adjustments paramters we would want to fix? or would want to be able to vary as part of a senstivity analysis?
- Pre industiral CH4 & CO2 paramters are defined in the SARF functions for CH4 & N2O, does that mean that we want to fix those values in Hector? If the pre-inndustiral concentrations change do the other paramters need to be adjusted?
2.2 Results
# Make a list of all the halocarbonns modeled by Hecotr
vars <- c(RF_HFC32(), RF_CF4(), RF_C2F6(), RF_HFC23(), 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())
# Certain RF values from AR6 must be adjusted, define the function that does the tropsoheric
# applies the appropriate adjustment.
ar6_adjustments <- function(data){
# for the values the assessed adjustment to CFCs is therefore 12 % ± 13%
data %>%
filter(grepl(pattern = "CFC", variable)) ->
CFC_data
# adjust the CFC data by the trop adjustment
CFC_data$value <- (CFC_data$value * 0.12) + CFC_data$value
data %>%
filter(!grepl(pattern = "CFC", variable)) %>%
bind_rows(CFC_data) ->
out
return(out)
}
# Run the Hector 2.5 with the original parameters
ini <- system.file("input/hector_rcp45.ini", package = "hector")
core45 <- newcore(ini)
run(core45)## Hector core: unnamed hector core
## Start date: 1745
## End date: 2300
## Current date: 2300
## Input file: /Users/dorh012/Library/R/3.6/library/hector/input/hector_rcp45.inirslt_rcp45 <- fetchvars(core45, 1850:2100 ,vars)
rslt_rcp45$run <- 'hector_2.5'
# Run RCP 4.5 but us the ini file with the updated parameters related to
# halocarbon forcing (tau and rho).
ini <- system.file("input/hector_rcp45_new.ini", package = "hector")
core_ar6 <- newcore(ini)
run(core_ar6)## Hector core: unnamed hector core
## Start date: 1745
## End date: 2300
## Current date: 2300
## Input file: /Users/dorh012/Library/R/3.6/library/hector/input/hector_rcp45_new.inirslt_ar6 <- ar6_adjustments(fetchvars(core_ar6, 1850:2100 ,vars)) # Make sure to adjust forcing to account for the trop adjustments
rslt_ar6$run <- 'AR6'Compare the results from Hector 2.5 and Hector AR6
rslt_rcp45 %>%
bind_rows(rslt_ar6) %>%
ggplot(aes(year, value, color = run, linetype = run)) +
geom_line() +
facet_wrap('variable', scales = "free") +
labs(title = "Comparison of Hector 2.5 params & AR6 params",
subtitle = "RCP 4.5",
y = "W m-2") +
theme_bw()
rslt_rcp45 %>%
bind_rows(rslt_ar6) %>%
select(year, variable, value, run) %>%
spread(run, value) %>%
mutate(mse = (AR6 - hector_2.5)^2) %>%
ggplot(aes(year, mse, color = variable)) +
geom_line() +
theme_bw() +
labs(title = 'Squared Error between Hector 2.5 & AR6 RF',
y = 'SE')
Looking at the MSE the RF for the two CFCs had the largest change, this all together not thatsurprising to me considering that these were the two variables where the trop adjustment factor came into play. Even then the max MSE was a value of 7x10^4 W/m^2. Which seems like a small number.
rslt_rcp45 %>%
bind_rows(rslt_ar6) %>%
select(year, variable, value, run) %>%
spread(run, value) %>%
mutate(se = (AR6 - hector_2.5)^2) %>%
group_by(variable) %>%
summarise(MSE = mean(se)) %>%
ungroup ->
MSE_table## `summarise()` ungrouping output (override with `.groups` argument)MSE_table %>%
arrange(desc(MSE)) %>%
knitr::kable()| variable | MSE |
|---|---|
| FCFC12 | 0.0001142 |
| FCFC11 | 0.0000564 |
| FHCFC22 | 0.0000205 |
| FCCl4 | 0.0000082 |
| FCFC113 | 0.0000031 |
| FHFC143a | 0.0000013 |
| FCF4 | 0.0000008 |
| FHFC134a | 0.0000005 |
| FHCFC141b | 0.0000002 |
| FCFC115 | 0.0000002 |
| FHFC245fa | 0.0000002 |
| FSF6 | 0.0000002 |
| FHFC32 | 0.0000001 |
| FCH3Cl | 0.0000001 |
| FCFC114 | 0.0000000 |
| FHFC125 | 0.0000000 |
| FCH3CCl3 | 0.0000000 |
| FHCFC142b | 0.0000000 |
| FHFC23 | 0.0000000 |
| Fhalon1301 | 0.0000000 |
| FHFC227ea | 0.0000000 |
| Fhalon2402 | 0.0000000 |
| FC2F6 | 0.0000000 |
| FCH3Br | 0.0000000 |
| Fhalon1211 | 0.0000000 |
| FHFC4310 | 0.0000000 |
For context how does this compare to the total over RF?
total_rf <- fetchvars(core45, 1850:2100, RF_TOTAL())mean(total_rf$value)## [1] 1.619367Compared to the average value of RF this error is a somewhat small fraction.
max(MSE_table$MSE)/mean(total_rf$value)## [1] 7.053875e-05rslt_rcp45 %>%
bind_rows(rslt_ar6) %>%
select(year, variable, value, run) %>%
spread(run, value) %>%
mutate(diff = (AR6 - hector_2.5)) %>%
group_by(year) %>%
summarise(diff = sum(diff)) %>%
ungroup ->
RF_change## `summarise()` ungrouping output (override with `.groups` argument)fetchvars(core45, 1850:2100, RF_TOTAL()) %>%
left_join(RF_change) %>%
mutate(frac = diff/value) %>%
ggplot(aes(year, frac)) +
geom_line() +
theme_bw() +
labs(y = "Fraction of change in halo RF vs. total RF")## Joining, by = "year"
For the most part the change in halocarbon RF represnts a small fraction of Hector’s total RF. There are some instances of rather larges peaks however those peaks correspond to when Hector’s total RF approaches 0 so might be an artifact of dividing by a small number.
3 GHGs
Updating the RF for the GHGs is going to be slighlty more complicated. In AR6 they have introdcuted new functional forms for CO2, CH4, and N2O RF. To see how much of an affect this has on the RF times series here I am going to use the Hector output for concentrations to calculate RF using the H25 (hector 25) and AR6 approach.
Run Hecctor 2.5 & get some data from it.
ini <- system.file("input/hector_rcp45.ini", package = "hector")
core <- newcore(ini)
run(core)## Hector core: unnamed hector core
## Start date: 1745
## End date: 2300
## Current date: 2300
## Input file: /Users/dorh012/Library/R/3.6/library/hector/input/hector_rcp45.inico2_conc <- fetchvars(core, 1850:2100, vars = ATMOSPHERIC_CO2())
n2o_conc <- fetchvars(core, 1850:2100, vars = ATMOSPHERIC_N2O())
ch4_conc <- fetchvars(core, 1850:2100, vars = ATMOSPHERIC_CH4())
co2_rf <- fetchvars(core, 1850:2100, vars = RF_CO2())
n2o_rf <- fetchvars(core, 1850:2100, vars = RF_N2O())
ch4_rf <- fetchvars(core, 1850:2100, vars = RF_CH4())
hector_rf <- rbind(co2_rf, n2o_rf, ch4_rf)
hector_rf$run <- 'hector_25'3.1 CO2
Set up the functions to get RF CO2 using the two different approaches.
# define the funtion that calcualtes CO2 RF in Hector forcings cpp file.
AR5_co2 <- function(C){
assertthat::assert_that(unique(C['variable']) == "Ca")
# Extract the paramters from Hector
C0 <- fetchvars(core, NA, PREINDUSTRIAL_CO2())[['value']]
rho <- fetchvars(core, NA, RHO_CO2())[['value']]
forcings <- rho * log(C$value/C0)
out <- data.frame(year = C['year'],
value = forcings,
variable = RF_CO2(),
run = "H25")
return(out)
}# define the funtions used to get the CO2 RF using the AR6 approach using
# equations from Table 7.SM.1 and then adjusting for the ERF
AR6_co2 <- function(co2_conc, n2o_conc){
# Define the coefficents from AR6 equations
a1 = -2.478e-7
b1 = 7.590e-4
c1 = -2.149e-3
d1 = 5.248
C0 = 277.1
C <- co2_conc$value
N <- n2o_conc$value
C_alpha_max = C0 - ((b1)/(2*a1))
N2O_alpha = c1 * sqrt(N)
# Define the function that calculates alpha prime, it is a piecewise function
# so there is a in statement.
alpha_prime_func <- function(single_c){
if(single_c > C_alpha_max){
return(d1 - ((b1^2)/ 2 * a1))
} else if (C0 < single_c & single_c < C_alpha_max){
return(d1 + a1 * (single_c - C0)^2 + b1 * (single_c - C0))
} else if (single_c < C0){
return(d1)
}
}
alpha_prime_values <- unlist(lapply(C, alpha_prime_func))
sarf_co2 <- (alpha_prime_values + N2O_alpha) * log(C / C0)
# Adjust SARF using the trop adjustments
value <- sarf_co2 + 0.05 * sarf_co2
out <- data.frame(year = co2_conc$year,
value = value,
variable = RF_CO2(),
run = "AR6")
return(out)
}Plot the results.
bind_rows( AR5_co2(co2_conc),
AR6_co2(co2_conc, n2o_conc)) ->
co2_rf
co2_rf %>%
ggplot(aes(year, value, color = run)) +
geom_line() +
theme_bw() +
labs(title = "RF CO2: Hector 2.5 vs AR 6")
What is the difference between the two?
co2_rf %>%
spread(run, value) %>%
mutate(SE = (H25 - AR6) ^2) ->
co2_SE
co2_SE %>%
ggplot(aes(year, SE)) +
geom_line() +
theme_bw() +
labs(title = "Squared Error bewteen AR6 & Hector 2.5 CO2 RF",
y = "Squared Error")
How does this change compare the total rf?
co2_rf %>%
spread(run, value) %>%
mutate(error = (AR6 - H25)) ->
co2_error
mean(co2_error$error)/mean(total_rf$value)## [1] 0.02176857This change makes up may be 2% of total RF.
3.2 CH4
# Get the CH4 RF from the hector output
ch4_rf$run <- "H25"Define the function(s) used to get the CH4 RF using the AR6 approach.
# Get the AR6 CH4 values
AR6_CH4_func <- function(ch4_conc, n2o_conc){
# get the sarf value based on the equaitons from Table 7.SM.1
AR6_SARF_CH4 <- function(M, N){
a3 <- -8.960e-5
b3 <- -1.246e-4
d3 <- 0.04519
M0 <- 731.41
return((a3 * sqrt(M) + b3 * sqrt(N) + d3) * (sqrt(M) - sqrt(M0)))
}
arf_value <- AR6_SARF_CH4(ch4_conc$value, n2o_conc$value)
# Adjust the CH4 RF for ERF
value <- arf_value * (1-.14)
out <- data.frame(year = ch4_conc$year,
value = value,
variable = RF_CH4(),
run = "AR6")
return(out)
}
CH4_ar6 <- AR6_CH4_func(ch4_conc, n2o_conc)bind_rows(CH4_ar6, ch4_rf) %>%
ggplot(aes(year, value, color = run)) +
geom_line() +
theme_bw() +
labs(title = "RF CH4: Hector 2.5 vs AR 6")
What is the MSE between this output?
bind_rows(CH4_ar6, ch4_rf) %>%
select(year, variable, run, value) %>%
spread(run, value) %>%
mutate(SE = (AR6 - H25) ^2) %>%
pull(SE) %>%
mean()## [1] 0.00068854273.2.1 Probelm
- the pre industrial CH4 between HEctor v 2.5 & AR6 are different from one another!
AR6 CH4 pre = 731.41 Hector v 2.5 pre = 653
3.3 N2O
n2o_rf$run <- "H25"AR6_N2O_func <- function(ch4_conc, n2o_conc, co2_conc){
AR6_SARF_N20 <- function(C, N, M){
a2 = -3.419e-4
b2 = 2.545e-4
c2 = -2.435e-4
d2 = 0.1217
N0 = 273.8
(a2 * sqrt(C) + b2 * sqrt(N) + c2 * sqrt(M) + d2) * (sqrt(N) - sqrt(N0))
}
rf <- AR6_SARF_N20(M = ch4_conc$value, N = n2o_conc$value, C = co2_conc$value)
values <- rf + (0.07) * rf
out <- data.frame(year = n2o_conc$year,
value = values,
variable = RF_N2O(),
run = "AR6")
return(out)
}
ar6_n2o <- AR6_N2O_func(ch4_conc, n2o_conc, co2_conc)bind_rows(ar6_n2o,
n2o_rf) %>%
ggplot(aes(year, value, color = run)) +
geom_line() +
theme_bw() +
labs(title = "RF N2O: Hector 2.5 vs AR6")
What is the MSE between this output?
bind_rows(ar6_n2o,
n2o_rf) %>%
select(year, variable, run, value) %>%
spread(run, value) %>%
mutate(SE = (AR6 - H25) ^2) %>%
pull(SE) %>%
mean()## [1] 3.125535e-05