Experiment 1
20 November, 2020
Description
In experiment 1 we apply the different FoRTE girdling treatments to the ED forest stand. This stand starts from near bare ground starting in 1900. After a series of tests, we had to exclude pines from the ED configuration.
See GitHub issue for some more details.
Observations
- These results make more sense espcially in the above ground biomass, the reduction in biomass matches the treatment gradient
- NPP and NEP results are interesting the difference in Rh between baseline and treatment groups are relatively small
Questions
- By 2030 NPP hasn’t completely recovered should we extend the length of the run?
- The cyclical NPP is most likely related to how the met data was looped over, do we want to keep it or swtich to a climate record like we do with the future climate? Part of me thinks that it is silly to waste too much time on this, from a writing perspective it will probably take more words to explain the three different sources of cliamte data. Where as using a climate record would decrease interannual variability which could help make the responses clearer, it depends what our objective is I think.
- How do we want to calculate resilience?
Set Up
Import the results
# Note that this is the object created with the drake frame work, there should be some better way to trigger this.
exp_object <- readRDS(file.path(BASE_DIR, 'ED-outputs', "exp-1.rds"))
exp_object <- map(exp_object, function(x){ x[year < 2030, ] }) # Note the year 2030 is not a complete run Event information
<year> 2019 </year>
<doy> 141 </doy>Meteorological 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 <- 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')## Warning in melt(met_data, measure.vars = c("nbdsf", "nddsf", "vbdsf", "vddsf", :
## The melt generic in data.table has been passed a data.frame and will attempt
## to redirect to the relevant reshape2 method; please note that reshape2 is
## deprecated, and this redirection is now deprecated as well. To continue using
## melt methods from reshape2 while both libraries are attached, e.g. melt.list,
## you can prepend the namespace like reshape2::melt(met_data). In the next
## version, this warning will become an error.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 MONTHLY_BACKGROUND(plot = ggplot(data = long_met_data)) +
THEME +
geom_line(aes(date, value)) +
facet_wrap("variable", scales = 'free') +
labs(title = 'Monthly Time Series of ED Met Data Inputs', y = NULL) 
Because of the variability in the seasonal cycle is rather large and hides longer trem trends. Let’s take a look at the annual min, max and mean.
annual_data <- long_met_data[ , list('value' = mean(value), 'min' = min(value), 'max' = max(value)),
by = list(year, variable)]
ANNUAL_BACKGROUND(plot = ggplot(data = annual_data)) +
geom_ribbon(aes(year, ymin = min, ymax = max), fill = 'grey', alpha = 0.8) +
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) +
THEME
Zooming in what does the actual observational data look like?
ggplot(data = annual_data[year %in% 1979:2008, ]) +
geom_ribbon(aes(year, ymin = min, ymax = max), fill = 'grey', alpha = 0.8) +
geom_line(aes(year, min), color = 'grey') +
geom_line(aes(year, max), color = 'grey') +
geom_line(aes(year, value)) +
facet_wrap('variable', scales = 'free') +
THEME +
labs(y = NULL)
Baseline Scenario
ANNUAL_BACKGROUND(plot = {exp_object$NPP[scn == 'harvest 0', ] %>%
.[ , value := sum(value), by = c('year', 'variable', 'unit', 'scn')] %>%
.[ , list(year, variable, unit, scn, value)] %>%
unique() %>%
ggplot()}) +
geom_line(aes(year, value)) +
THEME +
labs(title ='NPP',
x = 'Year',
y = unique(exp_object$NPP[scn == 'harvest 0', ][['unit']]), caption = NULL) ->
NPP_plot
MONTHLY_BACKGROUND(plot = ggplot(data = exp_object$LAI[scn == 'harvest 0', ] %>%
.[ , value := sum(value),
by = c('datetime', 'variable', 'unit', 'scn')] %>%
.[ , list(datetime, variable, unit, scn, value)] %>%
unique())) +
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
ANNUAL_BACKGROUND(plot = ggplot(exp_object$NEP[scn == 'harvest 0', ] %>%
.[ , value := sum(value), by = c('year', 'variable', 'unit', 'scn')] %>%
.[ , list(year, variable, unit, scn, value)] %>%
unique())) +
geom_line(aes(year, value)) +
THEME +
labs(title ='NEP',
x = 'Year',
y = unique(exp_object$NEP[scn == 'harvest 0', ][['unit']]), caption = NULL) ->
NEP_plot
ANNUAL_BACKGROUND(plot = ggplot(exp_object$ABG[scn == 'harvest 0', ] %>%
.[ , value := sum(value), by = c('year', 'variable', 'unit', 'scn')] %>%
.[ , list(year, variable, unit, scn, value)] %>%
unique())) +
geom_line(aes(year, value)) +
THEME +
labs(title ='ABG',
x = 'Year',
y = unique(exp_object$ABG[scn == 'harvest 0', ][['unit']]), caption = NULL) ->
ABG_plot
plot_grid(ABG_plot, NPP_plot, NEP_plot, LAI_plot)
The base line scenario looks good to me it is intresting just how much smoother the output looks when we swtich to the climatological record instead of weather, do we want to only use the climatological values.
Experiment 1
Above Ground Biomass
So the differences in total above ground biomass between the different sencaiors is faily small and is really only visible when we zoom in on the disturbance and recovery. What is intresting is that the as we would expect the no harvest scenario has a larger biomass until September and then the AGB starts to even out, is pretty similar during winter but then with the next growing season the patches with the more extreme treatments have a larger biomass.
exp_object$ABG$datetime <- date(exp_object$ABG$datetime)
data <- exp_object$ABG[ , value := sum(value), by = c("datetime", "scn", "unit", "year")]
data %>%
ggplot() +
geom_vline(xintercept = date('2019-05-01'), color = 'grey', size = 2) +
geom_point(aes(datetime, value, color = scn, group = interaction(scn, pft_name))) +
THEME +
labs(x = 'Year', y = unique(exp_object$ABG$unit), title = 'ABG') 
NPP
Annual values, what is odd is that the harvest 0% has NPP is more negative than the other treatments? Is that because of reading in the 0 harvest? Also does the monthly NPP have any other patterns?
ANNUAL_BACKGROUND(ggplot(exp_object$NPP[year %in% 2000:2030, ])) +
geom_vline(xintercept = 2019, color = 'grey', size = 2) +
geom_line(aes(year, value, color = scn, group = interaction(scn, pft_name)), size = 1.5) +
THEME +
coord_cartesian(xlim = c(2000, 2025)) +
labs(x = 'Year', y = unique(exp_object$NPP$unit), title = 'Annual NPP')
So the annual is some what hard to see the the pattern with it but after a decade there the NPP has not returned to the annual values.
h0_object <- readRDS(file.path(BASE_DIR, 'ED-outputs', 'exp-1', 'harvest_0_above.rds'))
h45_object <- readRDS(file.path(BASE_DIR, 'ED-outputs', 'exp-1', 'harvest_45_above.rds'))
h65_object <- readRDS(file.path(BASE_DIR, 'ED-outputs', 'exp-1', 'harvest_65_above.rds'))
h85_object <- readRDS(file.path(BASE_DIR, 'ED-outputs', 'exp-1', 'harvest_85_above.rds'))
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$scn <- 'harvest_0'
h45$scn <- 'harvest_45'
h65$scn <- 'harvest_65'
h85$scn <- 'harvest_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", "scn", "value", 'unit')]
NPP$year <- lubridate::year(NPP$datetime)
NPP$datetime <- date(NPP$datetime)
NPP_monthly <- NPP NPP_monthly <- NPP_monthly[, value := sum(value), by = c("datetime", "scn", "unit", "year")]
NPP_monthly %>%
ggplot(aes(datetime, value, color = scn)) +
# Note because of the fact that the results are recorded monthly the disturbance date is May 5th.
geom_vline(xintercept = date('2019-05-01'), color = 'grey', size = 2) +
geom_line() +
THEME +
labs(title = 'Monthly NPP', y = unique(NPP_monthly$unit), x = 'Time') 
NPP_monthly %>%
.[year %in% 2018:2022] %>%
ggplot(aes(datetime, value, color = scn)) +
geom_vline(xintercept = date('2019-05-01'), color = 'grey', size = 1.5) +
geom_line(size = 1) +
THEME +
labs(title = 'Monthly NPP', y = unique(NPP_monthly$unit), x = 'Time') 
NEP
exp_object$NEP %>%
dplyr::filter(year >= 2000 & year <= 2021) %>%
ggplot() +
geom_vline(xintercept = 2019, color = 'grey', size = 1.5) +
geom_line(aes(year, value, color = scn), size = 1) +
THEME +
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$scn <- 'harvest_0'
h45$scn <- 'harvest_45'
h65$scn <- 'harvest_65'
h85$scn <- 'harvest_85'
NEP <- rbind(h0, h45, h65, h85)
NEP$unit <- unique(exp_object$NEP$unit)
NEP <- NEP[, c("datetime", "scn", "value" = "MMEAN_NEP_PY")]
names(NEP) <- c("datetime", "scn", "value")
NEP$year <- lubridate::year(NEP$datetime)
NEP$datetime <- date(NEP$datetime)NEP %>%
dplyr::filter(year >= 2019 & year <= 2021) %>%
ggplot() +
geom_vline(xintercept = date('2019-05-01'), color = 'grey', size = 2) +
geom_line(aes(datetime, value, color = scn), size = 1) +
THEME +
labs(x = 'Year', y = unique(exp_object$NEP$unit), title = 'Monthly NEP')
LAI
exp_object$LAI$year <- year(exp_object$LAI$datetime)
exp_object$LAI %>%
dplyr::filter(year >= 2019 & year <= 2025) %>%
ggplot() +
geom_vline(xintercept = date('2019-05-01'), color = 'grey', size = 2) +
geom_line(aes(datetime, value, color = scn)) +
THEME +
labs(x = 'Year', y = unique(exp_object$LAI$unit), title = 'Monthly LAI') +
facet_wrap('pft_name', ncol = 1)
Soil 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$scn <- 'harvest_0'
h45$scn <- 'harvest_45'
h65$scn <- 'harvest_65'
h85$scn <- 'harvest_85'
RH <- rbind(h0, h45, h65, h85)
units <- 'kg/m2/yr'
RH %>%
ggplot(aes(datetime, MMEAN_RH_PY, color = scn)) +
geom_line() +
THEME +
labs(title = 'Monthly Heterotrophic Respiration',
y = units,
x = 'Date Time')
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: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.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] purrr_0.3.4 lubridate_1.7.8 cowplot_1.0.0 data.table_1.12.8
## [5] magrittr_1.5 ggplot2_3.3.2
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.5 plyr_1.8.6 pillar_1.4.4 compiler_3.6.3
## [5] tools_3.6.3 digest_0.6.25 evaluate_0.14 lifecycle_0.2.0
## [9] tibble_3.0.1 gtable_0.3.0 pkgconfig_2.0.3 rlang_0.4.6
## [13] yaml_2.2.1 xfun_0.14 withr_2.2.0 dplyr_1.0.0
## [17] stringr_1.4.0 knitr_1.28 generics_0.0.2 vctrs_0.3.1
## [21] rprojroot_1.3-2 grid_3.6.3 tidyselect_1.1.0 glue_1.4.1
## [25] here_0.1 R6_2.4.1 rmarkdown_2.3.1 farver_2.0.3
## [29] reshape2_1.4.3 scales_1.1.1 backports_1.1.8 ellipsis_0.3.1
## [33] htmltools_0.5.0 colorspace_1.4-1 labeling_0.3 stringi_1.4.6
## [37] munsell_0.5.0 crayon_1.3.4