BC IRF from the Sand et al Impulse Response

Purpose

The purpose of this document is to outline the process I used to get the BC IRF.

Format Data

# Import the required pacakges
library(dplyr)
library(tidyr)
library(ggplot2)
library(readxl)
# Define the directories
BASE_DIR  <- "/Users/dorh012/Documents/2019/cmip-arctic-data/SLCFimpulse"
INPUT_DIR <- file.path(BASE_DIR, 'data-raw')
# Import the Sand et al BC impulse resonse function and 
path_csv <- file.path(INPUT_DIR, 'TREFHT_fig1_toSteveSmith.xlsx')
raw_data <- read_excel(path_csv)

-
/
                                                                                                                                                 

The Sand et. al impulse response is the monthly temperature response but we want to work with the annual temperature response so I averaged the monthly data to get the average annual temperature response.

# Figure out how many year of data there are in the data set. 
n_yrs <- nrow(raw_data) / 12
# Then add the column of the years starting from 1 to the number of years 
# of data in monthly data set. Note that becasue the data is currently 
# monthly each year will need to be repreated 12 times in a row. 
data.frame(year = rep(1:n_yrs, each = 12),
           value = raw_data$`SAT BC`) %>%
  # Remove 1 from the year so that the temperature response is relative to the year of the step.
  mutate(year = year - 1) %>%
  # Find the annual average.
  group_by(year) %>%
  summarise(value = mean(value)) %>%
  ungroup ->
  data

Here is the annual temperature response to the BC step experiment. There is a lot of noise in the data so we will need to fit an exponential equation to the time sereis.

ggplot(data = data, aes(x = year, y = value)) + 
  geom_line() + 
  labs(y = 'deg C', 
       x = 'year', 
       title = 'Annual Temperature Response to BC Step') -> 
  temp_step_plot

The function we are fitting to the step temperature response is the following.

\[amplitude *(1 - {e}^{\frac{-year}{\tau}})\] where we need to solve for the amplitude and \(\tau\).

# Fit the exponetial results.
fit       <- nls(data$value ~ a *(1 - exp(-data$year / tau)), data = data, start = list(a = 2, tau = 1))
fit_rslts <- predict(fit)

How does the fit compare ti the annual data?

# plot the fit.
plot(data$year,data$value, xlab = 'Year', ylab = 'deg C', main = 'Annual Temperature Response to BC Step with Fit')
lines(data$year, fit_rslts, col="red")

What are our values for \(amplitude\) and \(\tau\)?

coefficients(summary(fit)) 

    Estimate Std. Error    t value      Pr(>|t|)
a   1.176718 0.01028081 114.457737 2.384155e-146
tau 2.079153 0.30534521   6.809188  2.291458e-10

amplitude <- coefficients(summary(fit))[1,1]
tau       <- coefficients(summary(fit))[2,1]

Does the amplitude match what we would expect the amplitude to look like? How close it is to the mean temp of the step after 50 years? (By this point in time the temperature has reached equilibrium).

# Subset the data so that the it only contains the temperature response after the first 50 years.
data %>% 
  filter(year > 50) %>% 
  pull(value) %>% 
  mean() -> 
  expected_amp
# Our expected value based on the equilibirum temp 
expected_amp 

[1] 1.179832

# The amplitude we got from the nonlinear fit
amplitude

[1] 1.176718

Hey they are pretty close to one another! I think that is pretty darn good!

Now that we have values for \(amplitude\) and \(\tau\) we can take a look at the deriative of the step fit which is equivalent to the temperature impulse response.

The derivative of the step fit looks like this.

\[ \frac{amplitude}{\tau}{e}^{\frac{-year}{\tau}}\]

OKay so now that we have the temperature impulse response we are going to have to normalize the temperature repsonse by RF to get the BC temperature response to a RF pulse because that is what we want the to used in HIRM.

So the units we want here are \(\frac{degC * m^2}{W}\)

If we start with \(deg C\) then we need to divide by the \(\frac\)

# Run Hector for some emissions pathwaway, we will use the BC emissions and RF output to determine Hecotr's RF/Tg BC relationship. 
core <- newcore(system.file('input/hector_rcp26.ini', package = 'hector')) 
run(core)
hector_output <- fetchvars(core, 1850:2100, vars = c(EMISSIONS_BC(), RF_BC()))

hector_output %>% 
  select(year, variable, value) %>% 
  spread(variable, value) %>% 
  mutate(value = FBC/BC_emissions) %>% 
  pull(value) %>% 
  mean -> 
  RF_Emiss

However the RF_Emiss is the radaitve foring response per TG of BC emissions and according to an email exchange

Use Sand et al.. That is temperature change for a 133 Tg step change in BC emissions. You’ll want to mathematically invert that to the change for 133 Tg impulse.

Therefore we need to multiply the RF_Emiss so that the RF_Emiss we use to normalize the BC temperature response is the appropriate size.

OKay so when we compare the BC and the generic IRF against one another they are similar the BC Sand et al IRF has less of a warming effect and declines faster. I guess that is the dynamic that we are interested in. But let’s see how it holds up when we use it in HIRM.

# First set up HIRM to reproduce Hector temperature. 
matrix       <- SLCFimpulse::Hector_ConfigMatrix
matrix[ , ]  <- 1
rcp45_rows   <- which(grepl(x = row.names(matrix), pattern = 'rcp45'))
config_rcp45 <- matrix[rcp45_rows, , drop = FALSE]
rlst_rcp45   <- core_IRM(config_matrix = config_rcp45, rf_list = Hector_RF, irf_list = Hector_IRF, match_agents = FALSE)
rcp45_temp   <- mutate(rlst_rcp45$total_temp, model = 'HIRM = Hector')
rcp45_conT   <- mutate(rlst_rcp45$temp_contributions, model = 'HIRM = Hector')

BC_IRF_list <- BC_RF_IRF
BC_IRF_list$source <- 'Sand et al'
BC_IRF_list <- split(BC_IRF_list, BC_IRF_list$source)

BC_RF_IRF$agent <- NA
BC_RF_IRF$name <- NA
input <- BC_RF_IRF
  
  Hector_IRF_new          <- Hector_IRF
  Hector_IRF_new[['BC']] <-  input
  rlst_BC   <- core_IRM(config_matrix = config_BC, rf_list = Hector_RF, irf_list = Hector_IRF_new, match_agents = FALSE)