2020/12/08 FoRTE ED Disturbance Meeting 1

Objective

Ben and I would like to talk over some of our results from ED, start discussing and thinking about a paper that uses the model and field results, we have started outlining a manuscript here.

Met Inputs

• dlwrf — downward long wave radiation (W m-2)
• nbdsf — near infrared beam downward solar radiation (W m-2)
• nddsf — near IR diffuse downward solar radiation (W m-2)
• vbdsf — visible beam downward solar radiation (W m-2)
• vddsf — visible diffuse downward solar radiation (W m-2)
• prate — precipitation rate (kgH2O m-2 s-1)
• pres — atmospheric pressure (Pa)
• hgt — geopotential height (m)
• ugrd — zonal wind (m s-1)
• vgrd — meridional wind (m s-1)
• sh — specific humidity (kgH2O kgAir-1)
• tmp — air temperature (K)

met_data <- as.data.table(read.csv(file.path(BASE_DIR, 'A.inputs', 'monthly_mean_met.csv')))
met_data$date <- ymd(paste(met_data$year, met_data$month, '01', sep = '/'))
long_met_data <- melt(met_data, measure.vars = c("nbdsf", "nddsf", "vbdsf", "vddsf",
                                                 "prate", "dlwrf", "pres", "hgt", "ugrd", 
                                                 "vgrd", "sh", "tmp"),
                      variable.name = 'variable', value.name = 'value')
long_met_data <- as.data.table(long_met_data[, c("date", "variable", "value", "year")])
var_info <- data.table(variable = c('dlwrf', 'nbdsf', 'nddsf', 'vbdsf',
                                    'vddsf', 'prate', 'pres', 'hgt', 
                                    'ugrd', 'vgrd', 'sh', 'tmp'), 
                       description = c('downward long wave radiation', 
                                       'near infrared beam downward solar radiation', 
                                       'near IR diffuse downward solar radiation', 
                                       'visible beam downward solar radiation',
                                       'visible diffuse downward solar radiation', 
                                       'precipitation rate', 
                                       'atmospheric pressure', 
                                       'geopotential height', 
                                       'zonal wind',
                                       'meridional wind', 
                                       'specific humidity',
                                       'air temperature'),  
                       units = c('W m-2', 'W m-2', 'W m-2', 'W m-2',
                                 'W m-2', 'kgH2O m-2 s-1', 'Pa', 'm',
                                 'm s-1', 'm s-', 'kgH2O kgAir-1', 'K'))
long_met_data <- long_met_data[var_info, on = 'variable']
long_met_data <- long_met_data[variable != 'hgt', ] # exclude the geo spatial height on cause it doesn't 

annual_data <- long_met_data[ , list('value' = mean(value), 'min' = min(value), 'max' = max(value)),
                              by = list(year, variable)]
ADD_BACKGROUND(plot = ggplot(data = annual_data)) +
  geom_ribbon(aes(year, ymin = min, ymax = max), fill = 'grey') + 
  geom_line(aes(year, min), color = 'grey') + 
  geom_line(aes(year, max), color = 'grey') + 
  geom_line(aes(year, value)) +
  facet_wrap('variable', scales = 'free') + 
  labs(y = NULL, title = 'Annual met data with annual max and min')

Control Results

The control, baseline, or 0 harvest ED results.

# Process the data to plot, the annual values. 
exp_object$NPP[scn == 'harvest 0', ] %>%  
    .[  , value := sum(value), by = c('year', 'variable', 'unit', 'scn')] %>% 
    .[ , list(year, variable, unit, scn, value)] %>% 
  unique() -> 
  NPP_to_plot
exp_object$LAI[scn == 'harvest 0', ] %>%
  .[  , value := sum(value), by = c('datetime', 'variable', 'unit', 'scn')] %>% 
  .[ , list(datetime, variable, unit, scn, value)] %>%  
  unique() -> 
  LAI_to_plot
exp_object$NEP[scn == 'harvest 0', ] %>%  
  .[  , value := sum(value), by = c('year', 'variable', 'unit', 'scn')] %>% 
  .[ , list(year, variable, unit, scn, value)] %>%  
  unique() -> 
  NEP_to_plot
exp_object$ABG[scn == 'harvest 0', ] %>%  
  .[  , value := sum(value), by = c('year', 'variable', 'unit', 'scn')] %>% 
  .[ , list(year, variable, unit, scn, value)] %>%  
  unique() -> 
  ABG_to_plot
exp_object$GPP[scn == "harvest 0", ] %>%  
  .[  , value := sum(value), by = c('year', 'variable', 'unit', 'scn')] %>% 
  .[ , list(year, variable, unit, scn, value)] %>%  
  unique() -> 
  GPP_to_plot
ADD_BACKGROUND(plot = ggplot(data = NPP_to_plot)) + 
  geom_line(aes(year, value), size = 1) +
  THEME +
  labs(title ='NPP',
       x = 'Year', 
       y = unique(exp_object$NPP[scn == 'harvest 0', ][['unit']]), caption = NULL) -> 
  NPP_plot
ADD_BACKGROUND(plot = ggplot(data = LAI_to_plot)) + 
  geom_point(aes(datetime, value)) + 
  THEME + 
  labs(title ='LAI',
       x = 'Year', 
       y = unique(exp_object$LAI[scn == 'harvest 0', ][['unit']]), caption = NULL, size = 0.01) -> 
  LAI_plot
ADD_BACKGROUND(ggplot(NEP_to_plot)) + 
  geom_line(aes(year, value), size = 1) +
  THEME + 
  labs(title ='NEP',
       x = 'Year', 
       y = unique(exp_object$NEP[scn == 'harvest 0', ][['unit']]), caption = NULL) -> 
  NEP_plot
ADD_BACKGROUND(ggplot(ABG_to_plot)) + 
  geom_line(aes(year, value), size = 1) + 
  THEME + 
  labs(title ='ABG',
       x = 'Year', 
       y = unique(exp_object$ABG[scn == 'harvest 0', ][['unit']]), caption = NULL) -> 
  ABG_plot
ADD_BACKGROUND(ggplot(GPP_to_plot)) + 
  geom_point(aes(year, value), size = 1) + 
  THEME + 
  labs(title = 'GPP', 
       x = 'Year', 
       y = 'kg/m2/yr', caption = NULL) -> 
  GPP_plot
plot_grid(ABG_plot, NPP_plot, NEP_plot, LAI_plot, GPP_plot)

Experiment 1

Near bare ground runs starting in 1900 until 2030. The 0, 45, 65, 85 disturbance severity treatments are applied in late May 2018.

Above Ground Biomass

exp_object$ABG$datetime <- date(exp_object$ABG$datetime)
data <- exp_object$ABG[  , value := sum(value), 
                         by = c("datetime", "scn", "unit", "year", "severity")] 
ggplot(data = data) + 
  ADD_DISTURBANCE_LINE_MONTH + 
  geom_line(aes(datetime, value, color = severity, 
                group = interaction(scn, pft_name)),
            size = 1.5) + 
  THEME + 
  ADD_COLORS + 
  labs(x = 'Year', y = unique(exp_object$ABG$unit), title = 'ABG') 

GPP

GPP <- exp_object$GPP
GPP <- dplyr::bind_rows(add_disturbance_severity(list(GPP)))
GPP$datetime <- date(GPP$datetime)
GPP$PFT <- as.character(GPP$PFT)
GPP <- GPP[year %in% 2000:2030, ] 
# Add the cohort labels to the data
split(GPP, interaction(GPP$datetime, GPP$severity, GPP$PFT), drop = TRUE) %>% 
  lapply(function(x){
    
    x$cohort <- LETTERS[1:nrow(x)]
    return(x)
    
  }) %>%  
  rlist::list.rbind() -> 
  GPP_cohort
to_plot <- GPP_cohort[year %in% 2018:2020, ] 
to_plot$value <- to_plot$value / 12
to_plot %>%  
  ggplot(aes(datetime, value, group = interaction(cohort, severity), color = severity)) + 
  ADD_DISTURBANCE_LINE_MONTH +
  geom_line(size = 2) + 
  facet_grid(. ~ pft_name) + 
  THEME + 
  ADD_COLORS + 
  labs(title = 'GPP per Cohort', 
       y = 'kgC/m2/month', 
       x = 'Date')

NPP

ggplot(exp_object$NPP[year %in% 2010:2030, ]) +
  ADD_DISTURBANCE_LINE_ANNUAL +
  geom_line(aes(year, value, color = severity, 
                group = interaction(scn, pft_name)), size = 1.5) + 
  THEME + 
  coord_cartesian(xlim = c(2010, 2025)) + 
  labs(x = 'Year', y = unique(exp_object$NPP$unit), title = 'Annual NPP') + 
  ADD_COLORS

h0 <- as.data.table(h0_object$df_cohort[[1]])[  ,list(datetime, MMEAN_NPP_CO, PFT, DBH, NPLANT)]
h45 <- as.data.table(h45_object$df_cohort[[1]])[  ,list(datetime, MMEAN_NPP_CO, PFT, DBH, NPLANT)]
h65 <- as.data.table(h65_object$df_cohort[[1]])[  ,list(datetime, MMEAN_NPP_CO, PFT, DBH, NPLANT)]
h85 <- as.data.table(h85_object$df_cohort[[1]])[  ,list(datetime, MMEAN_NPP_CO, PFT, DBH, NPLANT)]
h0$severity  <- '0'
h45$severity <- '45'
h65$severity <- '65'
h85$severity <- '85'
NPP <- rbind(h0, h45, h65, h85)
# Convert the NPP units from per plant to unit area
NPP$value <- NPP$MMEAN_NPP_CO * NPP$NPLANT
NPP$unit  <- unique(exp_object$NPP$unit)
NPP <- NPP[, c("datetime", "PFT", "severity", "value", 'unit')]
NPP$year <- lubridate::year(NPP$datetime)
NPP$datetime <- date(NPP$datetime)
monthly_NPP <- NPP
NPP_monthly <- monthly_NPP[, value := sum(value), by = c("datetime", "severity", "unit", "year")]
NPP_monthly$value <- NPP_monthly$value / 12
NPP_monthly %>% 
  .[year %in% 2018:2020] %>%  
  ggplot(aes(datetime, value, color = severity)) + 
  ADD_DISTURBANCE_LINE_MONTH +
  ADD_COLORS +
  geom_line(size = 1.5) + 
  THEME + 
  labs(title = 'Monthly NPP', y = 'kgC/m2/month', x = 'Time') 

NEP

exp_object$NEP %>% 
  dplyr::filter(year >= 2010 & year <= 2021) %>% 
  ggplot() + 
  geom_vline(xintercept = 2019, color = 'grey', size = 1.5) + 
  geom_line(aes(year, value, color = severity), size = 1.5) + 
  THEME + 
  ADD_COLORS +
  labs(x = 'Year', y = unique(exp_object$NEP$unit), title = 'Annual NEP')

h0 <- as.data.table(h0_object$df_scalar[[1]])[ ,list(datetime, MMEAN_NEP_PY)]
h45 <- as.data.table(h45_object$df_scalar[[1]])[ ,list(datetime, MMEAN_NEP_PY)]
h65 <- as.data.table(h65_object$df_scalar[[1]])[ ,list(datetime, MMEAN_NEP_PY)]
h85 <- as.data.table(h85_object$df_scalar[[1]])[ ,list(datetime, MMEAN_NEP_PY)]
h0$severity  <- '0'
h45$severity <- '45'
h65$severity <- '65'
h85$severity <- '85'
NEP <- rbind(h0, h45, h65, h85)
NEP$unit  <- unique(exp_object$NEP$unit)
NEP <- NEP[, c("datetime", "severity", "value" = "MMEAN_NEP_PY")]
names(NEP) <- c("datetime", "severity", "value")
NEP$year <- lubridate::year(NEP$datetime)
NEP$datetime <- date(NEP$datetime)
NEP %>% 
  dplyr::filter(year >= 2018 & year <= 2020) %>% 
  dplyr::mutate(value = value / 12, 
                units = 'kg/m2/month') %>% 
  ggplot() + 
  ADD_DISTURBANCE_LINE_MONTH + 
  geom_line(aes(datetime, value, color = severity), size = 1.5) + 
  THEME + 
  ADD_COLORS +
  labs(x = 'Year', y ='kg/m2/month', title = 'Monthly NEP')

Heterotrophic Respiration

h0 <- as.data.table(h0_object$df_scalar[[1]])[  ,list(datetime, MMEAN_RH_PY)]
h45 <- as.data.table(h45_object$df_scalar[[1]])[  ,list(datetime, MMEAN_RH_PY)]
h65 <- as.data.table(h65_object$df_scalar[[1]])[  ,list(datetime, MMEAN_RH_PY)]
h85 <- as.data.table(h85_object$df_scalar[[1]])[  ,list(datetime, MMEAN_RH_PY)]
h0$severity  <- '0'
h45$severity <- '45'
h65$severity <- '65'
h85$severity <- '85'
RH    <- rbind(h0, h45, h65, h85)
RH$datetime <- as_date(RH$datetime )
units <- 'kg/m2/month'
RH$MMEAN_RH_PY <- RH$MMEAN_RH_PY / 12
RH[datetime >= ymd('2018-06-01') & datetime <= ymd('2021-01-01'), ] %>% 
ggplot(aes(datetime, MMEAN_RH_PY, color = severity)) + 
  ADD_DISTURBANCE_LINE_MONTH +
  geom_line(size = 1.5) + 
  THEME + 
  labs(title = 'Monthly Heterotrophic Respiration', 
       y = units,
       x = 'Date Time') + 
  ADD_COLORS

LAI

exp_object$LAI$year <- year(exp_object$LAI$datetime)
exp_object$LAI %>% 
  dplyr::filter(year >= 2018 & year <= 2020) %>% 
  dplyr::filter(pft_name == 'Temperate broadleaf, early successional') %>% 
  ggplot() + 
  ADD_DISTURBANCE_LINE_MONTH + 
  geom_line(aes(datetime, value, color = severity), size = 1.5) + 
  ADD_COLORS +
  THEME + 
  labs(x = 'Year', y = unique(exp_object$LAI$unit), title = 'Monthly LAI') + 
  facet_wrap('pft_name', ncol = 1)

Carbon Cycling Resistance

Carbon Cycling Resistance

\[r = \ln \left( {\frac{{F_{{\text{dist}}} }}{{F_{{\text{cont}}} }}} \right)\] Equation 1 from Gough 2020.

exp_object$NPP[year %in% 2018:2030, ] %>%  
  .[, list(pft_name, year, variable, unit, value, severity)] %>% 
  # Make the table wide by disturbance severity treatment. 
  dcast(data = ., pft_name + year + variable + unit ~ severity) %>%
  # Now gather all of the distubrance treatments into a single column so that the 
  # control values are recorded in their own column. 
  melt(data = ., id.vars = c("pft_name", "year", "variable", "unit", "0"), 
       measure.vars = c('45', '65', '85'), variable.name = 'severity', 
       value.name = 'treatment_value') %>%  
  # Rename the columns 
  .[ , list(pft_name, year, variable, unit, control = `0`, severity, dist = treatment_value)] %>% 
  # Now calculate r 
  .[ , r := log(dist/control)] -> 
  NPP_r
  
NPP_r %>%   
ggplot(aes(year, r, color = severity)) + 
  geom_line(size = 1.5) + 
  ADD_COLORS + 
  THEME + 
  labs(title = 'NPP Resistance Values', 
       y = expression('r'[npp]), 
       x = 'Year')