Disturbance Treatments

Description

On 2020-11-12 BBL and I looked over the experiment 1 results but have some concerns about the treatments. Here we let’s look at different harvest event set ups

• the original treatment only harvest the storage C
• uniform harvest acccross all four C pools
• apply harvest to the storage C and read in 0s for all three C pools

Set Up

Results Discussed on 2020-11-16

AGB

The above ground biomass from the xml files I used in experiment 1 a few weeks a go. Our big concern was that the different treatment groups all decreased in biomass by baout 90%

object$ABG[grepl('original|baseline', scn), ] %>% 
  .[, value := sum(value), by = list(datetime, variable, description, unit, scn)] %>% 
  .[, list(datetime, value, scn, variable, unit)] %>%  
  ggplot() + 
  geom_line(aes(datetime, value, color = scn)) + 
  THEME + 
  labs(title = 'ABG for the Harvest XML files used for 2020-11-11 results', y = "kgC/m2")

Now let’s apply the harvests uniformly to all of the tree carbon storage pools (above ground biomass, below ground biomass, folliage, and storage)

object$ABG[grepl('uniform', scn), ] %>% 
  .[, value := sum(value), by = list(datetime, variable, description, unit, scn)] %>% 
  .[, list(datetime, value, scn, variable, unit)] %>%  
  ggplot() + 
  geom_line(aes(datetime, value, color = scn)) + 
  THEME + 
  labs(title = 'ABG for the Uniform Harvest', y = "kgC/m2")

This is more along the lines of what we would expect to see!

Now let us apply the harvest to only the storage C pool.

data <- object$ABG
data <- cbind(data, 'year' = year(data$datetime))
data[grepl(pattern = 'other|baseline', x = scn), ] %>% 
  .[, value := sum(value), by = list(datetime, variable, description, unit, scn, year)] %>% 
  .[, list(datetime, value, scn, variable, unit, year)] %>%  
  unique %>% 
  .[year %in% 2019:2022, ] %>% 
  ggplot() + 
  geom_line(aes(datetime, value, color = scn)) + 
  THEME + 
  labs(title = 'ABG C Storage Harvest Treatment Only', y = "kgC/m2")

While yes the pattern is along the lines of what we would expect but the differences are pretty small, which begs the question is the difference meandingful or trivial?

cbind(object$ABG, year = date(object$ABG$datetime)) %>% 
  .[grepl(pattern = 'other|baseline', x = scn), ] %>% 
  .[, list(datetime, value, scn, variable, unit, year)] %>%  
  .[, value := sum(value), by = list(datetime, variable, unit, scn, year)] %>% 
  unique %>% 
  .[year %in% 2019:2022, ] %>% 
  dcast(year + datetime + variable + unit ~ scn, value.war = 'value') %>% 
  melt(id.vars = c("year", "datetime", "variable","unit","baseline"), variable.name = 'treatment') %>% 
  .[ , diff := baseline - value] %>%  
    ggplot() + 
  geom_line(aes(datetime, value, color = treatment)) + 
  THEME + 
  labs(title = 'Difference in ABG C Storage Harvest Treatment Only', y = "Baseline - Treatment (kgC/m2)")

Or perhaps AGB is not the right variable to look at, since we are not havesting the AGB it is not surprising that the treatments have a minimal impact on this output variable.

NPP

Start by checking out the annual values.

object$NPP %>% 
  .[grepl(pattern = 'other|baseline', x = scn), ] %>% 
  .[ , scn := gsub(pattern = ' other 0 ', replacement = '', x = scn)] %>% 
  .[ , treat := 'C only'] -> 
  annual_NPP_other
  
object$NPP %>% 
  .[grepl(pattern = 'uniform|baseline', x = scn), ] %>% 
  .[ , scn := gsub(pattern = ' uniform', replacement = '', x = scn)] %>% 
  .[ , treat := 'uniform'] -> 
  annual_NPP_uniform
rbind(annual_NPP_other, annual_NPP_uniform) %>%  
  .[, value := sum(value), by = list(year, scn, year, description, unit, treat)] %>%
  unique %>% 
  ggplot(aes(year, value, color = scn)) + 
  geom_line() + 
  THEME + 
  facet_wrap('treat')

Import & format the uniform disturbance treatments.

h0 <- readRDS(file.path(ED_OUTPUT_DIR, 'disturbence-treatments', 'harvest_0_uniform.rds'))
h45 <- readRDS(file.path(ED_OUTPUT_DIR, 'disturbence-treatments', 'harvest_45_uniform.rds'))
h65 <- readRDS(file.path(ED_OUTPUT_DIR, 'disturbence-treatments', 'harvest_65_uniform.rds'))
h85 <- readRDS(file.path(ED_OUTPUT_DIR, 'disturbence-treatments', 'harvest_85_uniform.rds'))
h0 <- as.data.table(h0$df_cohort[[1]])[  ,list(datetime, MMEAN_NPP_CO, PFT, DBH, NPLANT)]
h45 <- as.data.table(h45$df_cohort[[1]])[  ,list(datetime, MMEAN_NPP_CO, PFT, DBH, NPLANT)]
h65 <- as.data.table(h65$df_cohort[[1]])[  ,list(datetime, MMEAN_NPP_CO, PFT, DBH, NPLANT)]
h85 <- as.data.table(h85$df_cohort[[1]])[  ,list(datetime, MMEAN_NPP_CO, PFT, DBH, NPLANT)]
h0$scn  <- 'baseline'
h45$scn <- 'harvest_45'
h65$scn <- 'harvest_65'
h85$scn <- 'harvest_85'
h0$treat <- 'Uniform Harvest'
h45$treat <- 'Uniform Harvest'
h65$treat <- 'Uniform Harvest'
h85$treat <- 'Uniform Harvest'
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 <- NPP[, c("datetime", "PFT", "scn", "value", "treat")]
NPP$year <- lubridate::year(NPP$datetime)
NPP$datetime <- date(NPP$datetime)
NPP_uniform <- NPP 

Import & format the storage specific treatments

h45 <- readRDS(file.path(ED_OUTPUT_DIR, 'disturbence-treatments', 'harvest_45_other0.rds'))
h65 <- readRDS(file.path(ED_OUTPUT_DIR, 'disturbence-treatments', 'harvest_65_other0.rds'))
h85 <- readRDS(file.path(ED_OUTPUT_DIR, 'disturbence-treatments', 'harvest_85_other0.rds'))
h45 <- as.data.table(h45$df_cohort[[1]])[  ,list(datetime, MMEAN_NPP_CO, PFT, DBH, NPLANT)]
h65 <- as.data.table(h65$df_cohort[[1]])[  ,list(datetime, MMEAN_NPP_CO, PFT, DBH, NPLANT)]
h85 <- as.data.table(h85$df_cohort[[1]])[  ,list(datetime, MMEAN_NPP_CO, PFT, DBH, NPLANT)]
h0$scn <- 'baseline'
h45$scn <- 'harvest_45'
h65$scn <- 'harvest_65'
h85$scn <- 'harvest_85'
h0$treat <- 'C Storage Harvest'
h45$treat <- 'C Storage Harvest'
h65$treat <- 'C Storage Harvest'
h85$treat <- 'C Storage Harvest'
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 <- NPP[, c("datetime", "PFT", "scn", "value", "treat")]
NPP$year <- lubridate::year(NPP$datetime)
NPP$datetime <- date(NPP$datetime)
NPP_other0 <- NPP 

rbind(NPP_uniform, NPP_other0) %>% 
  .[year %in% 2019:2022] %>% 
  .[, value := sum(value), by = list(datetime, scn, year, treat)] %>%
  unique -> 
  NPP_to_plot 
NPP_to_plot %>%  
  ggplot(aes(datetime, value, color = scn)) + 
  geom_line() +
  facet_wrap('treat') + 
  THEME + 
  labs(y = 'NPP kgC/m2/yr per month')

NPP_to_plot %>%  
  ggplot(aes(datetime, value, color = treat, linetype = treat)) + 
  geom_line() +
  facet_wrap('scn') + 
  THEME + 
  labs(y = 'NPP kgC/m2/yr per month')

Conclusion

• Turns out that reading the version of ED’s fortran code that if a carbon disturbance is not specific then ED assigns the disturbance of 1 (complete disturbance) which is explains the large decline in AGB we saw in the results from last week.
• Now we have runs that are more realistic where we apply the same disturbance treatment to all of the different plant carbon pools and runs that only apply the disturbance to the C storage as per Jackie’s suggestion. However the difference in AGB and carbon fluexes (implied by NPP) are smaller between the C storage treatments than I was expecting. It is unclear if this is a problem or not. Looking at the Pest ED paper applything phlome disturbances of varying degrees had a substantial impact on carbon fluxes (NEE) which is what makes the ED results seem a tad suspect

Ideas for next steps

• BBL which one seems better/more realistic? My concern is that the C only disturbances are too small but perhaps that is realistic.
• Run the complete runs with the two different harvests, does it matter how mature the forest is?
• Followup with Jackie and Alexey about the results

January Disturbance

So it looks like starting in before the growing season does not have the intended impact.

Above Ground Biomass Disturbance

object$ABG[grepl('above', scn), ] %>% 
  .[, value := sum(value), by = list(datetime, variable, description, unit, scn)] %>% 
  .[, list(datetime, value, scn, variable, unit)] %>%  
  ggplot() + 
  geom_line(aes(datetime, value, color = scn, linetype = scn)) + 
  THEME + 
  labs(title = 'ABG Above Ground Biomass', y = "kgC/m2")

---

sessionInfo()

R version 3.6.3 (2020-02-29)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ed4forte_0.0.0.9000 tidyselect_1.1.0    purrr_0.3.4         cowplot_1.0.0       magrittr_1.5       
 [6] ggplot2_3.3.2       assertthat_0.2.1    lubridate_1.7.8     data.table_1.12.8   tibble_3.0.1       
[11] drake_7.9.0        

loaded via a namespace (and not attached):
 [1] storr_1.2.1      xfun_0.14        reshape2_1.4.3   lattice_0.20-41  testthat_2.3.2   colorspace_1.4-1
 [7] vctrs_0.3.1      generics_0.0.2   htmltools_0.5.0  yaml_2.2.1       base64enc_0.1-3  utf8_1.1.4      
[13] rlang_0.4.6      pillar_1.4.4     txtq_0.2.0       glue_1.4.1       withr_2.2.0      sp_1.4-4        
[19] plyr_1.8.6       lifecycle_0.2.0  stringr_1.4.0    munsell_0.5.0    gtable_0.3.0     evaluate_0.14   
[25] labeling_0.3     knitr_1.28       parallel_3.6.3   fansi_0.4.1      Rcpp_1.0.5       readr_1.3.1     
[31] scales_1.1.1     backports_1.1.8  filelock_1.0.2   desc_1.2.0       pkgload_1.1.0    jsonlite_1.6.1  
[37] farver_2.0.3     hms_0.5.3        digest_0.6.25    stringi_1.4.6    dplyr_1.0.0      grid_3.6.3      
[43] rprojroot_1.3-2  here_0.1         rgdal_1.4-8      cli_2.0.2        tools_3.6.3      base64url_1.4   
[49] tidyr_1.1.0      crayon_1.3.4     pkgconfig_2.0.3  ellipsis_0.3.1   rmarkdown_2.3.1  rstudioapi_0.11 
[55] R6_2.4.1         igraph_1.2.4.2   compiler_3.6.3