Hector: Initalized Ocean C

1 Set Up

library(data.table)
library(dplyr)
library(tidyr)
library(ggplot2)
library(magrittr)

BASE_DIR <- here::here()

theme_set(theme_bw())

1.1 Objective

Take a look at how different the Hector output is when different initial C pools are used.

out <- as.data.table(rbind(read.csv(file.path(BASE_DIR, "hector_new.csv")),
                           read.csv(file.path(BASE_DIR, "hector_old.csv"))))

2 How did the inital C values change?

df <- data.table(box = c("HL", "LL", "IO", "DO"), 
           old = c(140, 770, 8400, 26000), 
           new = c(135, 765, 9080, 28019))
df

Surface

data.table(old = sum(df[df$box %in% c("HL", "LL")]$old), 
           new = sum(df[df$box %in% c("HL", "LL")]$new))

Deep / Interior

data.table(old = sum(df[df$box %in% c("IO", "DO")]$old), 
           new = sum(df[df$box %in% c("IO", "DO")]$new))

Total Ocean C changes by about 7%.

data.table(old = sum(df$old), 
           new = sum(df$new))

3 Ocean Output

4 pH

out  %>%
    filter(grepl(x = tolower(variable), pattern = "ph")) ->
    ph_data

ph_data %>%
    ggplot(aes(year, value, color = pre, linetype = scenario)) +
    geom_line() +
    facet_wrap("variable", scales = "free", ncol = 1)

ph_data %>%
    distinct() %>%
    dcast(scenario + year + units + variable ~ pre, value.var = "value", fun.aggregate = min) %>%
    mutate(percent =  100 *(`new` - `old`) / old) %>%
    group_by(variable) %>%
    summarise(min = min(percent),
              max = max(percent),
              mean = mean(percent))

5 CO3

out  %>%
    filter(grepl(x = tolower(variable), pattern = "co3")) ->
    co3_data

co3_data %>%
    ggplot(aes(year, value, color = pre, linetype = scenario)) +
    geom_line() +
    facet_wrap("variable", scales = "free", ncol = 1)

co3_data %>%
    distinct() %>%
    dcast(scenario + year + units + variable ~ pre, value.var = "value") %>%
    mutate(percent =  100 *(`new` - `old`) / `old`) %>%
    group_by(variable) %>%
    summarise(min = min(percent),
              max = max(percent),
              mean = mean(percent))

6 DIC

out  %>%
    filter(grepl(x = tolower(variable), pattern = "dic")) ->
    dic_data

dic_data %>%
    ggplot(aes(year, value, color = pre, linetype = scenario)) +
    geom_line() +
    facet_wrap("variable", scales = "free", ncol = 1)

dic_data %>%
    distinct() %>%
    dcast(scenario + year + units + variable ~ pre, value.var = "value") %>%
    mutate(percent =  100 *(`new` - `old`) / `old`) %>%
    group_by(variable) %>%
    summarise(min = min(percent),
              max = max(percent),
              mean = mean(percent))

7 Carbon

out  %>%
    filter(grepl(x = tolower(variable), pattern = "carbon_")) ->
    carbon_out

carbon_out %>%
    ggplot(aes(year, value, color = pre, linetype = scenario)) +
    geom_line() +
    facet_wrap("variable", scales = "free", ncol = 1)

carbon_out %>%
    distinct() %>%
    dcast(scenario + year + units + variable ~ pre, value.var = "value") %>%
    mutate(percent =  100 *(`new` - `old`) / `old`) %>%
    group_by(variable) %>%
    summarise(min = min(percent),
              max = max(percent),
              mean = mean(percent))

8 SST

out  %>%
    filter(variable == "Tgav_ocean_ST") ->
    tgav_out

tgav_out %>%
    ggplot(aes(year, value, color = pre, linetype = scenario)) +
    geom_line() +
    facet_wrap("variable", scales = "free", ncol = 1)

tgav_out %>%
    distinct() %>%
    dcast(scenario + year + units + variable ~ pre, value.var = "value") %>%
    mutate(percent =  100 *(`new` - `old`) / `old`) %>%
    group_by(variable) %>%
    summarise(min = min(percent),
              max = max(percent),
              mean = mean(percent))

9 The Big 4

big_four_vars <- c("Tgav", "Ftot", "Ca", "FCO2")
out %>%
    filter(variable %in% big_four_vars) -> 
    big_four_out
    
big_four_out %>% 
    ggplot(aes(year, value, color = pre, linetype = scenario)) +
    geom_line() +
    facet_wrap("variable", scales = "free")

9.1 Percent Difference in the Big Four

big_four_out %>%
    distinct() %>%
    dcast(scenario + year + units + variable ~ pre, value.var = "value") %>%
    mutate(percent =  100 *(`new` - `old`) / `old`) %>%
    group_by(variable) %>%
    summarise(min = min(percent),
              max = max(percent),
              mean = mean(percent))

10 Log files

10.1 Hector using the old values

Tue Feb 15 16:30:46 2022:DEBUG:stashCValues: annualflux_sum=0.0032725 Pg C
Tue Feb 15 16:30:46 2022:DEBUG:log_state: ----- State of LL box -----
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    carbon = 783.646 -> 783.65
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    T=20.9 degC, surfacebox=1, active_chemistry=1
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    CarbonAdditions = 69.2904 (8.84205%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    CarbonSubtractions = 68.2904 (8.71444%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    FPgC = -0.9965 (-0.127162%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    ao_flux = 0 (0%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    oa_flux = 0.9965 (0.127162%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    K0 = 0.0324716 mol/l/atm Tr = 0.0337689 gC/m/month/atm
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    Surface DIC = 2076.28 umol/kg
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    Connection #0 to HL, k=0.0742532, window=1
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    Connection #1 to intermediate, k=0.0128912, window=1
Tue Feb 15 16:30:46 2022:DEBUG:log_state: ----- State of HL box -----
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    carbon = 140.525 -> 140.525
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    T=1.6 degC, surfacebox=1, active_chemistry=1
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    CarbonAdditions = 98.3685 (70.0008%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    CarbonSubtractions = 99.3685 (70.7124%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    FPgC = 0.999773 (0.711456%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    ao_flux = 0.999773 (0.711456%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    oa_flux = 0 (0%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    K0 = 0.0609349 mol/l/atm Tr = 0.0368937 gC/m/month/atm
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    Surface DIC = 2109.82 umol/kg
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    Connection #0 to deep, k=0.707124, window=1
Tue Feb 15 16:30:46 2022:DEBUG:log_state: ----- State of intermediate box -----
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    carbon = 8418.94 -> 8418.94
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    T=18 degC, surfacebox=0, active_chemistry=0
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    CarbonAdditions = 273.472 (3.2483%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    CarbonSubtractions = 273.472 (3.24829%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    K0 = 0 (undefined) Tr = 0 (undefined)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    Connection #0 to LL, k=0.0082303, window=1
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    Connection #1 to HL, k=0.0047726, window=1
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    Connection #2 to deep, k=0.01948, window=1
Tue Feb 15 16:30:46 2022:DEBUG:log_state: ----- State of deep box -----
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    carbon = 25991.8 -> 25991.8
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    T=18 degC, surfacebox=0, active_chemistry=0
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    CarbonAdditions = 263.369 (1.01328%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    CarbonSubtractions = 263.37 (1.01328%)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    K0 = 0 (undefined) Tr = 0 (undefined)
Tue Feb 15 16:30:46 2022:DEBUG:log_state:    Connection #0 to intermediate, k=0.0101328, window=1

10.2 Hector using the new values

Tue Feb 15 16:16:24 2022:DEBUG:log_state: ----- State of LL box -----
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    carbon = 844.184 -> 844.186
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    T=20.9 degC, surfacebox=1, active_chemistry=1
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    CarbonAdditions = 74.5659 (8.8329%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    CarbonSubtractions = 73.5659 (8.71444%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    FPgC = -0.997694 (-0.118184%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    ao_flux = 0 (0%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    oa_flux = 0.997694 (0.118184%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    K0 = 0.0324716 mol/l/atm Tr = 0.0337689 gC/m/month/atm
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    Surface DIC = 2236.67 umol/kg
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    Connection #0 to HL, k=0.0742532, window=1
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    Connection #1 to intermediate, k=0.0128912, window=1
Tue Feb 15 16:16:24 2022:DEBUG:log_state: ----- State of HL box -----
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    carbon = 151.208 -> 151.208
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    T=1.6 degC, surfacebox=1, active_chemistry=1
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    CarbonAdditions = 105.923 (70.0511%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    CarbonSubtractions = 106.923 (70.7124%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    FPgC = 1.00002 (0.661357%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    ao_flux = 1.00002 (0.661357%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    oa_flux = 0 (0%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    K0 = 0.0609349 mol/l/atm Tr = 0.0368937 gC/m/month/atm
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    Surface DIC = 2270.22 umol/kg
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    Connection #0 to deep, k=0.707124, window=1
Tue Feb 15 16:16:24 2022:DEBUG:log_state: ----- State of intermediate box -----
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    carbon = 9059.93 -> 9059.93
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    T=18 degC, surfacebox=0, active_chemistry=0
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    CarbonAdditions = 294.293 (3.24829%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    CarbonSubtractions = 294.293 (3.24829%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    K0 = 0 (undefined) Tr = 0 (undefined)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    Connection #0 to LL, k=0.0082303, window=1
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    Connection #1 to HL, k=0.0047726, window=1
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    Connection #2 to deep, k=0.01948, window=1
Tue Feb 15 16:16:24 2022:DEBUG:log_state: ----- State of deep box -----
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    carbon = 27969.5 -> 27969.5
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    T=18 degC, surfacebox=0, active_chemistry=0
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    CarbonAdditions = 283.41 (1.01328%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    CarbonSubtractions = 283.411 (1.01328%)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    K0 = 0 (undefined) Tr = 0 (undefined)
Tue Feb 15 16:16:24 2022:DEBUG:log_state:    Connection #0 to intermediate, k=0.0101328, window=1
Tue Feb 15 16:16:24 2022:DEBUG:record_state: Recording component state at t= 1746

11 Comparison with obs

read.csv(file.path(BASE_DIR, "obs", "HOTS_PH.csv"), skip = 3) %>% 
    rename(variable = vtag) %>% 
    select(year = Year, value, variable) %>% 
    mutate(scenario = "HOTS") -> 
    obs1 

read.csv(file.path(BASE_DIR, "obs", "ensemble_mean_model.csv")) %>% 
    select(year, value, variable) %>% 
    mutate(scenario = "CMIP") -> 
    obs2

obs_out <- bind_rows(obs1, obs2) %>% 
        mutate(variable = tolower(variable)) 

out %>% 
    filter(variable %in% c("pH_HL", "pH_LL", "pH")) %>% 
    filter(pre == "new") %>% 
    mutate(variable = tolower(variable)) -> 
    hector_ph

out %>% 
    filter(variable %in% c("pH_HL", "pH_LL", "pH")) %>% 
    filter(pre == "old") %>% 
    mutate(variable = tolower(variable)) -> 
    hector_ph_old

ggplot() + 
    geom_point(data = hector_ph, aes(year, value, color = 'new hector'), size = 1) + 
    geom_point(data = hector_ph_old, aes(year, value, color = 'old hector'), size = 1) + 
    geom_point(data = obs_out, aes(year, value, color = scenario), size = 1) + 
    facet_wrap("variable") + 
    coord_cartesian(xlim = c(1850, 2000), ylim = c(7.45, 8.25)) + 
    labs(y = "ph total value")

Warning: Removed 4 rows containing missing values (geom_point).

Well this is less than ideal. Hmmmmmm…