Correcting Aerosol RF

Objective: Take a look at how the correction to the forcing scalings affect Hector output.

The resulting changes are to the following 5 aerosol constants:

SI Eq 12 -ρ𝐵𝐶 (new) = ρ𝐵𝐶(old) * ( (0.37+0.07) / (0.28 + 0.07) )**

SI Eq 13 through SI Eq 15 (for all the reflecting forcings) multiply by ( 0.65/0.63 )

SI Eq 16 -ρaci (new) = ρaci (old) * ( 0.88/ 0.81 )

The result will be a about a 5% larger net negative forcing, and a larger positive ARI from BC.

library(dplyr)
library(ggplot2)
library(hector)
library(knitr)

theme_set(theme_bw())

Load and format the data

# The Hector comp files were generated by the raw-data/generate-comp-data.R script
# before an after implementing the new aerosol correction scalar values, note that
# it was not until...
og_data <- read.csv("~/Documents/Hector/tests/testthat/compdata/hector_comp-dev-branch.csv")
new_data <- read.csv("~/Documents/Hector/tests/testthat/compdata/hector_comp-new-aerosol.csv")


og_data$version <- "V3"
new_data$version <- "Updated Aerosols"

# Save a long formatted data frame for easy plotting
long_data <- rbind(og_data, new_data)


# Save a copy of a wide data frame to compute the differences. 
og_data %>%
    select(-version) %>%
    rename(v3_value = value) %>%
    inner_join(new_data) %>%
    rename(aero_value = value) %>%
    mutate(SE = (v3_value - aero_value) ^2,
           percent = 100 * abs(abs(v3_value - aero_value)/v3_value) ) ->
    wide_data

Joining with `by = join_by(scenario, year, variable, units, commit)`

wide_data %>%
    group_by(variable) %>%
    summarise(MSE = mean(SE), percent = mean(percent, na.rm = TRUE)) %>% 
    arrange(desc(percent)) -> 
    summary_stats

kable(summary_stats)

variable	MSE	percent
RF_CO2	0.0001510	26.7411263
RF_BC	0.0010348	25.7142857
sst	0.0000203	21.8943594
global_tas	0.0000381	13.7688273
RF_tot	0.0001193	12.0917674
RF_aci	0.0010591	8.6419753
RF_SO2	0.0000287	3.1746032
RF_NH3	0.0000075	3.1746032
RF_OC	0.0000018	3.1746032
heatflux	0.0000406	2.6692331
CO2_concentration	0.6836167	0.1841800
atmos_co2	3.1015006	0.1841800
HL_pH	0.0000007	0.0092686
ocean_c	14.3051887	0.0087844

Global Temperature

long_data %>%
    filter(variable == GLOBAL_TAS()) %>%
    ggplot(aes(year, value, color = version)) +
    geom_line() + 
    labs(y = "deg C")

summary_stats %>%
    filter(variable == GLOBAL_TAS())

NA

Aerosol RFs

vars <- c(RF_BC(), RF_OC(), RF_SO2(), RF_NH3(), RF_ACI())

long_data %>%
    filter(variable %in% vars) %>%
    ggplot(aes(year, value, color = version)) +
    geom_line() + 
    labs(y = "W") + 
    facet_wrap("variable", scales = "free")

wide_data %>% 
    filter(variable %in% vars) %>% 
    ggplot(aes(year, percent, color = variable)) + 
    geom_line()

CO2 Variables

vars <- c("CO2_concentration", "RF_CO2", "ocean_c", "HL_pH")

long_data %>%
    filter(variable %in% vars) %>%
    ggplot(aes(year, value, color = version)) +
    geom_line() + 
    facet_wrap("variable", scales = "free")

vars <- c("CO2_concentration", "RF_CO2", "ocean_c", "HL_pH")

wide_data %>%
    filter(variable %in% vars) %>%
    filter(year > 1900) %>% 
    ggplot(aes(year, percent, color = version)) +
    geom_line() + 
    facet_wrap("variable", scales = "free")

Question for Steve

What was your vision for how these changes were going to be added to Hector. Does it mean adding a new corrective parameter per aerosol equation? Or would it be as simple as updating the \(\rho\) scalers with the new ones you’ve come up with that are the product of two scalars?