1 Objective

As part of the EPA task to update Hector RF to be consistent with the AR6 RF I have encountered something funky going on with the \(CH_{4}\) and \(N_{2}O\) emissions to concentration step. Because of the relations between [CH4] & [N2O] and their RF one this has to be a joint investigation.

library(dplyr)
library(tidyr)
library(ggplot2)
library(data.table)
library(hector) # or devtools::load_all() depending on the R project space being used. 

theme_set(theme_minimal())

# Define the project home directory, note that because we are working with an devloping version of hector this script 
# is not always run from within the EPA_hector project but within the R Hector dev project so cannot use here::here
DIR <- "/Users/dorh012/projects/2021/EPA_Hector" # define the analysis directory. 

# Run hector for a ssp scenario emission driven and then again with concentrations. 
years <- 1750:2050
vars <- c(EMISSIONS_N2O(), EMISSIONS_CH4(), ATMOSPHERIC_N2O(),
          ATMOSPHERIC_CH4(), GLOBAL_TEMP(), RF_CH4(), RF_N2O(), RF_TOTAL())
   

# Read in the hector inputs for the corresponding run. 
file.path(DIR, "input", "emissions", "ssp245_emiss-constraints.csv") %>% 
    read.csv(stringsAsFactors = FALSE, comment.char = ";") %>% 
    select(year = Date, N2O = N2O_constrain, N2O_emissions, 
           CH4 = CH4_constrain, CH4_emissions) %>%  
    melt(id.vars = "year") %>% 
    filter(year %in% years) %>% 
    mutate(scenario = "Hector inputs") -> 
    hector_inputs


ini <- file.path(DIR, "input", "hector_ssp245.ini")
ssp245_core <- newcore(ini, name = "Hv3")
run(ssp245_core, max(years))
## Hector core: Hv3
## Start date:  1745
## End date:    2300
## Current date:    2050
## Input file:  /Users/dorh012/Documents/2021/EPA_Hector/input/hector_ssp245.ini
ssp245_hector_out <- fetchvars(ssp245_core, years, vars)

    
ini <- file.path(DIR, "input", "hector_ssp245.ini")
ssp245_core_ch4 <- newcore(ini, name = "Hv3 CH4 const")
atm_ch4_vals <- filter(hector_inputs, variable == ATMOSPHERIC_CH4())
setvar(ssp245_core_ch4, dates = years, values = atm_ch4_vals$value, var = CH4_CONSTRAIN(), unit = "ppbv CH4")
reset(ssp245_core_ch4)
## Hector core: Hv3 CH4 const
## Start date:  1745
## End date:    2300
## Current date:    1745
## Input file:  /Users/dorh012/Documents/2021/EPA_Hector/input/hector_ssp245.ini
run(ssp245_core_ch4, max(years))
## Hector core: Hv3 CH4 const
## Start date:  1745
## End date:    2300
## Current date:    2050
## Input file:  /Users/dorh012/Documents/2021/EPA_Hector/input/hector_ssp245.ini
ssp245_core_ch4out <- fetchvars(ssp245_core_ch4, years, vars)

ini <- file.path(DIR, "input", "hector_ssp245.ini")
ssp245_core_n2o <- newcore(ini, name = "Hv3 N2O const")
atm_n2o_vals <- filter(hector_inputs, variable == ATMOSPHERIC_N2O())
setvar(ssp245_core_n2o, dates = years, values = atm_n2o_vals$value, var = N2O_CONSTRAIN(), unit = "ppbv N2O")
reset(ssp245_core_n2o)
## Hector core: Hv3 N2O const
## Start date:  1745
## End date:    2300
## Current date:    1745
## Input file:  /Users/dorh012/Documents/2021/EPA_Hector/input/hector_ssp245.ini
run(ssp245_core_n2o, max(years))
## Hector core: Hv3 N2O const
## Start date:  1745
## End date:    2300
## Current date:    2050
## Input file:  /Users/dorh012/Documents/2021/EPA_Hector/input/hector_ssp245.ini
ssp245_core_n2oout <- fetchvars(ssp245_core_n2o, years, vars)

ini <- file.path(DIR, "input", "hector_ssp245.ini")
ssp245_core_n2o_ch4 <- newcore(ini, name = "Hv3 N2O CH4 const")
atm_n2o_vals <- filter(hector_inputs, variable == ATMOSPHERIC_N2O())
setvar(ssp245_core_n2o_ch4, dates = years, values = atm_n2o_vals$value, var = N2O_CONSTRAIN(), unit = "ppbv N2O")
atm_n2o_vals <- filter(hector_inputs, variable == ATMOSPHERIC_N2O())
setvar(ssp245_core_n2o_ch4, dates = years, values = atm_n2o_vals$value, var = N2O_CONSTRAIN(), unit = "ppbv N2O")
reset(ssp245_core_n2o_ch4)
## Hector core: Hv3 N2O CH4 const
## Start date:  1745
## End date:    2300
## Current date:    1745
## Input file:  /Users/dorh012/Documents/2021/EPA_Hector/input/hector_ssp245.ini
run(ssp245_core_n2o_ch4, max(years))
## Hector core: Hv3 N2O CH4 const
## Start date:  1745
## End date:    2300
## Current date:    2050
## Input file:  /Users/dorh012/Documents/2021/EPA_Hector/input/hector_ssp245.ini
ssp245_core_n2o_ch4out <- fetchvars(ssp245_core_n2o_ch4, years, vars)


# Read in the AR6 data 
AR6_dir <- "/Users/dorh012/Documents/2021/Chapter-7/data_output"
file <- file.path(AR6_dir, "SSPs", "ERF_ssp245_1750-2500.csv")
read.csv(file, stringsAsFactors = FALSE) %>% 
    select(year, FCH4 = ch4, FN2O = n2o)  %>%
    melt(id.vars = c( "year"), variable.name = "variable", value.name = "value") %>% 
    mutate(scenario = "IPCC") %>% 
    filter(year %in% years)-> 
    ipcc_rf_output


COLORS <- c("Hv3" = "#88CCEE", 
            "Hv3 CH4 const" = "#CC6677", 
            "Hector inputs" = "#999933", 
            "IPCC" = "#888888", 
            "Hv3 N2O const"="#44AA99", 
            "Hv3 N2O CH4 const"= "#661100")

2 The Probelm

2.1 CH4 RF

When we run Hector with the updated CH4 RF we get negative CH4 radiative forcing!

vars_to_plot <- RF_CH4()
ssp245_hector_out %>% 
    bind_rows(ssp245_core_ch4out, hector_inputs, ipcc_rf_output) %>% 
    filter(variable %in% vars_to_plot) %>% 
    ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
    geom_vline(xintercept = 2014, color = 'grey') + 
    geom_hline(yintercept = 0, color = 'black') + 
    geom_line(size = 1) + 
    facet_wrap('variable', scales = 'free') +
    labs(y = "radiative forcing W m-2", 
         caption = "gery vertical at 2015, balck horizontal at 0")  + 
    scale_color_manual(values = COLORS)

So when Hector CH4 concentrations are constrained the Hector output is nearly identical to the IPCC results, which is awesome, this suggests that Hector’s implementation of the AR6 CH4 RF equation is correct. However when Hector is driven with CH4 emissions the RF is cooler than expected and before 1800 there is even a period of time when CH4 RF is negative! Suggesting a cooling effect because CH4 concentrations decreased! Which leads be to suspect that something funky is going on during the emissions to concentrations step.

2.2 CH4 concentrations

Here is a comparison of the CH4 concentration the blue line what

vars_to_plot <- c(ATMOSPHERIC_CH4())
ssp245_hector_out %>% 
    bind_rows(hector_inputs) %>% 
    filter(variable %in% vars_to_plot) %>% 
    ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
    geom_vline(xintercept = 2014, color = 'grey') + 
    geom_line(size = 1) + 
    facet_wrap('variable', scales = 'free') + 
    labs(y = "ppbv CH4", title = "CH4 concentration") + 
    scale_color_manual(values = COLORS)

I didn’t expect the difference between the prescribed CH4 concentrations to be all that different from the CH4 concentrations calculated by emission driven Hector. Something is the CH4 concentrations to decrease…

2.3 N2O RF

When we run Hector with the updated CH4 RF we get negative CH4 radiative forcing!

vars_to_plot <- RF_N2O()
ssp245_hector_out %>% 
    bind_rows(ssp245_core_n2oout, hector_inputs, ipcc_rf_output, ssp245_core_n2o_ch4out) %>% 
    filter(variable %in% vars_to_plot) %>% 
    ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
    geom_vline(xintercept = 2014, color = 'grey') + 
    geom_hline(yintercept = 0, color = 'black') + 
    geom_line(size = 1) + 
    facet_wrap('variable', scales = 'free') +
    labs(y = "radiative forcing W m-2", 
         caption = "gery vertical at 2015, balck horizontal at 0")  + 
    scale_color_manual(values = COLORS)

So when Hector N2O, CH4 concentrations are constrained the Hector output is closer but not identical to the IPCC results. The behavior of the Hv3 RF is funky…

2.4 N2O concentrations

Here is a comparison of the N2O concentrations and the shape of the N2O concentrations

vars_to_plot <- c(ATMOSPHERIC_N2O())
ssp245_hector_out %>% 
    bind_rows(hector_inputs) %>% 
    filter(variable %in% vars_to_plot) %>% 
    ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
    geom_vline(xintercept = 2014, color = 'grey') + 
    geom_line(size = 1) + 
    facet_wrap('variable', scales = 'free') + 
    labs(y = "ppbv N2O", title = "N2O concentration") + 
    scale_color_manual(values = COLORS)

The differences between the N2O concentrations prescribed as an input or internally calculated internally be Hector are larger than I would have expected.

Possible causes

  • Could the RF be coded incorrectly?
  • Could the CH4 & N2O emissions inputs be formatted incorrectly?
  • Could it be related to the natural CH4 or N2O emissions?
  • Could something be wrong with the life time of CH4 or N2O?
  • Could it be related to the postindustrial values?

3 Trying to identify the problem

3.1 Could the RF be coded incorrectly?

We can rule out this possibility since when Hector v3 is driven with concentrations Hector v3 is able to reproduce IPCC RF because Hector can do this we can rule out an issue with how the new AR6 RF code is being implemented.

vars_to_plot <- RF_CH4()
ssp245_hector_out %>% 
    bind_rows(ssp245_core_ch4out, hector_inputs, ipcc_rf_output) %>% 
    filter(variable %in% vars_to_plot) %>% 
    ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
    geom_vline(xintercept = 2014, color = 'grey') + 
    geom_hline(yintercept = 0, color = 'black') + 
    geom_line(size = 1) + 
    facet_wrap('variable', scales = 'free') +
    labs(y = "radiative forcing W m-2", 
         caption = "gery vertical at 2015, balck horizontal at 0")  + 
    scale_color_manual(values = COLORS)

3.2 Could the CH4 & N2O emissions inputs be formatted incorrectly?

This is certainly a possibility we are using outputs from hector data for the ssp scenarios, however one of the hectordata checks that that we are able to reproduce the rcp outputs using the rcmip conversion, which is proof that the emissions cannot be off, we check by comparing the emissions inputs generated from hector data with emissions commited to github.

Let’s run RP 45 again anyways and check to make sure that the new results are reasonably close to the Hector 2.5 RF values. Which they are… although the N2O RF time series has a bit of a funky shape if you as me…

Notice that even with the rcp45 emissions we still see the negative FCH4 in H with the AR6 RF functions.

vars <- c(RF_CH4(), RF_N2O())
vars <- EMISSIONS_N2O()

ini <- file.path(DIR, "input", "hector_rcp45.ini") 
core <- newcore(ini, name = "Hv3")
run(core, runtodate = 2050)
## Hector core: Hv3
## Start date:  1745
## End date:    2300
## Current date:    2050
## Input file:  /Users/dorh012/Documents/2021/EPA_Hector/input/hector_rcp45.ini
hector_rcp45 <- fetchvars(core, years, vars)

read.csv(file.path(DIR, "Hector_25_rcp45.csv")) %>% 
    filter(year <= 2050) %>% 
    select(scenario, year, variable, value, units) %>% 
    mutate(scenario = 'old rcp45') -> 
    Hector_25_results

hector_rcp45 %>% 
    bind_rows(Hector_25_results) %>% 
    filter(variable %in% vars) %>% 
    ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
    geom_line() + 
    facet_wrap('variable', scales = "free")

3.3 Could it be realted to the natural CH4 or N2O emissions?

Unclear as to how one would answer this, but seeing that in above in the RCP Hv3 FCH4 time series there is are negative values after 1750 I even though the two different Hectors are receiving the same emissions inputs both anthropocentric and natural make me think that that the natural emissions are not the culprits here.

3.4 Could something be wrong with the life time of CH4 or N2O?

Looking at the CH4 component there is no one single life time for CH4 as described by in Hartin et al 2015

ch4 fig

So it seems unlikely that the cause of this the lifetimes, unless there have changes in the science in the atmospheric concentrations. Does Steve or Ben know?

What if it is related to life time of N2O?? So the default tau N2O is

TN2O0=132               ; initial lifetime of N2O,

What does is mean by the initial life time of N2O? How is that any different than the life time of N2O reported in the IPCCC report which has a value of 109.

years <- 1750:2100
vars <- c(RF_N2O(), RF_CH4(), RF_TOTAL(), GLOBAL_TEMP())

ini <- file.path(DIR, "input", "hector_rcp45.ini") 
core <- newcore(ini, name = "Hv3")
run(core, runtodate = 2050)
## Hector core: Hv3
## Start date:  1745
## End date:    2300
## Current date:    2050
## Input file:  /Users/dorh012/Documents/2021/EPA_Hector/input/hector_rcp45.ini
out1 <- fetchvars(core, years, vars)

ini <- file.path(DIR, "input", "hector_rcp45_n2o.ini") 
core <- newcore(ini, name = "Hv3 AR6 N2O tau")
run(core, runtodate = 2050)
## Hector core: Hv3 AR6 N2O tau
## Start date:  1745
## End date:    2300
## Current date:    2050
## Input file:  /Users/dorh012/Documents/2021/EPA_Hector/input/hector_rcp45_n2o.ini
out2 <- fetchvars(core, years, vars)

bind_rows(out1, out2) %>%  
    ggplot(aes(year, value, color = scenario, linetype = scenario)) + 
    geom_line() + 
    facet_wrap("variable", scales = "free")

When updating to use the AR6 N2O life time makes the N2O RF negative! And it is not clear what how the lifetimes are different from one another. Also not clear if this is related to the

4 My take aways

  • there seems to be something funky going on when we change the pre industrial value for CH4, although I don’t really have the slightest idea as to what is going on