require(LIM) require(tidyverse) require(GGally) require(gridExtra) require(grid) source("WoodStoich_functions.R") # Load in functions that we need. ## Model differences across model structures ----- prelim_gen <- Read("fullmodel_v4.lim") listofflows =list(grep('BASE', prelim_gen$flows$name), unique(c(grep('A', prelim_gen$flows$name), grep('S', prelim_gen$flows$name))), unique(c(grep('F', prelim_gen$flows$name), grep('B', prelim_gen$flows$name), grep('X', prelim_gen$flows$name), grep('Y', prelim_gen$flows$name)), grep('Z', prelim_gen$flows$name)), unique(c(grep('A', prelim_gen$flows$name), grep('S', prelim_gen$flows$name), grep('F', prelim_gen$flows$name), grep('B', prelim_gen$flows$name), grep('X', prelim_gen$flows$name), grep('Y', prelim_gen$flows$name), grep('Z', prelim_gen$flows$name))), unique(c(grep('M', prelim_gen$flows$name), grep('J', prelim_gen$flows$name), grep('G', prelim_gen$flows$name), grep('T', prelim_gen$flows$name))), unique(c(grep('X', prelim_gen$flows$name), grep('Y', prelim_gen$flows$name), grep('Z', prelim_gen$flows$name))), unique(c(grep('R', prelim_gen$flows$name), grep('W', prelim_gen$flows$name), grep('Q', prelim_gen$flows$name), grep('K', prelim_gen$flows$name), grep('U', prelim_gen$flows$name))), unique(c(grep('H', prelim_gen$flows$name), grep('K', prelim_gen$flows$name), grep('G', prelim_gen$flows$name), grep('T', prelim_gen$flows$name)))) modnames = c("Base", "No animals", "No micobes", "Only abiotic processes", "No MAOM", "No Symbiotes", "No Recalcitrance", "No physical protection") if(F){ debugonce(runscenario3) a = runscenario3(PRELIM_GEN = prelim_gen,conSTOICH = 1,spSTOICH = 1,spUE = 0, LISTOFFLOWS = listofflows, wf=1, MODNAMES = modnames) b = runscenario3(PRELIM_GEN = prelim_gen,conSTOICH = 4,spSTOICH = 4,spUE = 4, LISTOFFLOWS = listofflows, wf=1, MODNAMES = modnames) a %>% as_tibble() %>% select(ID, diff2) %>% rename(C4 = diff2) %>% left_join( b %>% as_tibble() %>% select(ID, diff2) %>% rename(C1 = diff2) ) %>% mutate(diff = abs(C1)-abs(C4)) %>% filter(abs(diff) > 1) %>% ggplot(aes(x=ID, y=diff)) + geom_point() } # Test one round of runscenario2 if(F){ dd = runscenario3(wf=2) dd %>% separate(ID, into=c("ID", "Element"),sep=-1) %>% # filter(parsimonious >0) %>% ggplot(aes(x=ID, y=parsimonious)) + geom_point(size=1) + geom_errorbar(aes(ymin=min, ymax=max)) + theme_classic() + facet_grid(Element~., scales="free") rm(dd) } # Check a single scenario designmat = data.frame(constoich = rep(c(1,2,3,4,4,4,4,4,4,4,4,4), length(modnames)), spstoich = rep(c(1,1,1,1,2,3,4,1,1,1,1,4), length(modnames)), spue = rep(c(0,0,0,0,0,0,0,1,2,3,4,4), length(modnames)), wf = rep(c(1:length(modnames)), each=12)) if(F){ designmat = data.frame(constoich = c(1,2,3,4,4,4,4,4,4,4,4,4), spstoich = c(1,1,1,1,2,3,4,1,1,1,1,4), spue = c(0,0,0,0,0,0,0,1,2,3,4,4), wf = rep(1, each=12)) } NN = dim(designmat)[1] listofresults <- vector(mode = "list", length = NN) for(ii in 1:NN){ listofresults[[ii]] = runscenario3(PRELIM_GEN = prelim_gen, conSTOICH = designmat$constoich[[ii]], spSTOICH = designmat$spstoich[[ii]], spUE = designmat$spue[[ii]], LISTOFFLOWS = listofflows, wf = designmat$wf[[ii]], MODNAMES = modnames) print(paste("Done", ii, "of",NN)) } output = do.call("rbind", listofresults) write.csv(output, "Data/output_18Aug2019.csv", row.names = F) # output = read.csv("Data/output_18Aug2019.csv") ## Model differences across STOICH and UE ---- modnames = c("Base", "No symbiont constraints", "No saprotrophs constraints", "No microbial constraints", "No animal constraints", "No biotic constraints", "No MAOM constraints", "No abiotic constraints") listofpools = list("ALLTOGETHER", "Zc|Yc|Xc|Zue|Yue|Xue", "Bc|Fc|Bue|Fue", "Zc|Yc|Xc|Bc|Fc|Zue|Yue|Xue|Bue|Fue", "Ac|Sc|Aue|Sue", "Zc|Yc|Xc|Bc|Fc|Ac|Sc|Zue|Yue|Xue|Bue|Fue|Aue|Sue", "Mc|Jc", "Mc|Jc|LDc") designmat = data.frame(constoich = rep(c(1,2,3,4,4,4,4,4,4,4,4,4), length(modnames)), spstoich = rep(c(1,1,1,1,2,3,4,1,1,1,1,4), length(modnames)), spue = rep(c(0,0,0,0,0,0,0,1,2,3,4,4), length(modnames)), wf = rep(c(1:length(modnames)), each=12)) # runscenario4(wf=1, conSTOICH = 4, spSTOICH = 1, spUE = 2) NN = dim(designmat)[1] listofresults <- vector(mode = "list", length = NN) for(ii in 1:NN){ listofresults[[ii]] = runscenario4(PRELIM_GEN = prelim_gen, conSTOICH = designmat$constoich[[ii]], spSTOICH = designmat$spstoich[[ii]], spUE = designmat$spue[[ii]], LISTOFPOOLS = listofpools, wf = designmat$wf[[ii]], MODNAMES = modnames) print(paste("Done", ii, "of",NN)) } output2 = do.call("rbind", listofresults) write.csv(output2, "Data/output_diffconst_18Aug2019.csv", row.names = F) # output2 <- read.csv("Data/output_diffconst_18Aug2019.csv") ## Model differences across PARAMETERS ---- designmat = data.frame(constoich = rep(c(1,2,3,4,4,4,4,4,4,4,4,4),100), spstoich = rep(c(1,1,1,1,2,3,4,1,1,1,1,4),100), spue = rep(c(0,0,0,0,0,0,0,1,2,3,4,4),100)) # runscenario5(conSTOICH = 4, spSTOICH = 4, spUE = 4) NN = dim(designmat)[1] listofresults <- vector(mode = "list", length = NN) for(ii in 1:NN){ listofresults[[ii]] = runscenario5(PRELIM_GEN = prelim_gen, conSTOICH = designmat$constoich[[ii]], spSTOICH = designmat$spstoich[[ii]], spUE = designmat$spue[[ii]]) print(paste("Done", ii, "of",NN)) } output3 = do.call("rbind", listofresults) write.csv(output3, "Data/output_diffparams_18Aug2019.csv", row.names = F) # output3 = read.csv("Data/output_diffparams_18Aug2019.csv") ## ....Plot results: Prep ------ nameMODS = data.frame(Group = c("General", "Broad \n C:N", "Broad \n C:N:P", "Broad \n C:N:P:S", "Specific \n C:N", "Specific \n C:N:P", "Specific \n C:N:P:S", "CUE", "CUE & NUE", "CUE, NUE \n & PUE", "CUE, NUE, \n PUE & SUE", "Efficiency & \n Stoichiometry"), conSTOICH = c(1,2,3,4,4,4,4,4,4,4,4,4), spSTOICH = c(1,1,1,1,2,3,4,1,1,1,1,4), spUE = c(0,0,0,0,0,0,0,1,2,3,4,4)) nameELEMENTS = data.frame(Element = c("c","n","p","s"), Elementname = c("Carbon", "Nitrogen", "Phosphorus", "Sulfur")) toplot = output %>% as_tibble() %>% select(-parsimonious) %>% mutate(diff2 = abs(diff2)) %>% separate(ID, into=c("ID", "Element"),sep=-1) %>% group_by(Element, conSTOICH, spSTOICH, spUE, Modification) %>% summarize(diff = sum(diff2)) %>% left_join( nameMODS ) %>% left_join( nameELEMENTS ) %>% mutate(diff = diff/1000) #convert to kg toplot2 = output2 %>% as_tibble() %>% select(-parsimonious) %>% mutate(diff2 = abs(diff2)) %>% separate(ID, into=c("ID", "Element"),sep=-1) %>% group_by(Element, conSTOICH, spSTOICH, spUE, Modification) %>% summarize(diff = sum(diff2)) %>% left_join( nameMODS ) %>% left_join( nameELEMENTS ) toplot3 = output3 %>% as_tibble() %>% mutate(diff2 = abs(diff2)) %>% separate(ID, into=c("ID", "Element"),sep=-1) %>% group_by(Run,Element, conSTOICH, spSTOICH, spUE) %>% summarize(diff = sum(diff2)) %>% group_by(Element, conSTOICH, spSTOICH, spUE) %>% summarise(lower = quantile(diff, 0.25)/1000, upper = quantile(diff, 0.75)/1000) %>% left_join( nameMODS ) %>% left_join( nameELEMENTS ) toplot3_1 = output3 %>% as_tibble() %>% mutate(diff2 = abs(diff2)) %>% separate(ID, into=c("ID", "Element"),sep=-1) %>% group_by(Run,Element, conSTOICH, spSTOICH, spUE) %>% summarize(diff = sum(diff2)) %>% left_join( nameMODS ) %>% left_join( nameELEMENTS ) %>% filter(Element=="c") toplot1_1 = toplot %>% left_join( toplot %>% filter(Group=="General") %>% rename(diffbase = diff) %>% ungroup() %>% select(Element, Modification, diffbase) ) %>% rename(diffunstd=diff) %>% mutate(diff = diffunstd/ diffbase) %>% # because of rounding errors, models without MAOM and microbes can look # artifically higher than base, but actual numbers are the same (i.e. 4300 units +/- 10^-13). # I clean this error out of the data here. mutate(diff = ifelse(diff > 1, 1, diff)) ## ...........Run the plot ------ STCON0A = toplot %>% filter(Element=="c") %>% filter(Modification == "Base") %>% left_join(toplot3 %>% filter(Element == "c")) %>% filter(Group %in% c("General", "Broad \n C:N", "Broad \n C:N:P", "Broad \n C:N:P:S", "Specific \n C:N", "Specific \n C:N:P", "Specific \n C:N:P:S")) %>% ggplot(aes(x=Group, y=diff, group=1)) + geom_vline(xintercept = c(1.5, 4.5), col="grey", linetype ="dashed") + geom_point(size=3, col="red") + geom_errorbar(aes(ymin=lower, ymax = upper), col="red", width=0.2) + theme_classic() + ylab(parse(text = "Uncertainty~(kg[C]~m^-2~year^-1)")) + xlab("Scenario") + scale_y_continuous(breaks = scales::pretty_breaks(3), limits=c(0,300)) + scale_x_discrete(limits = c("General", "Broad \n C:N", "Broad \n C:N:P", "Broad \n C:N:P:S", "Specific \n C:N", "Specific \n C:N:P", "Specific \n C:N:P:S")) STCON0B = toplot %>% filter(Element=="c") %>% filter(Modification == "Base") %>% left_join(toplot3 %>% filter(Element == "c")) %>% filter(Group %in% c("Broad \n C:N:P:S", "CUE", "CUE & NUE", "CUE, NUE \n & PUE", "CUE, NUE, \n PUE & SUE")) %>% ggplot(aes(x=Group, y=diff, group=1)) + geom_vline(xintercept = c(1.5), col="grey", linetype ="dashed") + geom_point(size=3, col="red") + geom_errorbar(aes(ymin=lower, ymax = upper), col="red", width=0.2) + theme_classic() + ylab(parse(text = "Uncertainty~(kg[C]~m^-2~year^-1)")) + xlab("Scenario") + scale_y_continuous(breaks = scales::pretty_breaks(3), limits=c(0,300)) + scale_x_discrete(limits = c("Broad \n C:N:P:S", "CUE", "CUE & NUE", "CUE, NUE \n & PUE", "CUE, NUE, \n PUE & SUE")) STCON1_1 = toplot %>% filter(Element=="c" & Group == "General") %>% ggplot(aes(x=1, y=diff, group=Modification, color=Modification, linetype=Modification, alpha = Modification)) + geom_jitter(aes(shape=Modification),size=3) + theme_bw() + ylab(parse(text = "kg[C]~m^-2~year^-1")) + xlab("") + scale_y_continuous(breaks = scales::pretty_breaks(3)) + scale_x_discrete(limits = c("General")) + scale_shape_manual(values=c(19,11,8,9,17,18,10,12)) + scale_alpha_manual(values=c(1, rep(0.5, 8))) + theme(legend.position = "none") g2 = ggplotGrob(STCON1_1) STCON1A = toplot1_1 %>% filter(Element=="c") %>% ggplot(aes(x=Group, y=diff, group=Modification, color=Modification, linetype=Modification, alpha = Modification)) + geom_vline(xintercept = c(1.5, 4.5), col="grey", linetype ="dashed") + geom_jitter(aes(shape=Modification),size=3, width=0.2) + theme_bw() + ylab("Uncertainty (% of General)") + xlab("Scenario") + scale_x_discrete(limits = c("General", "Broad \n C:N", "Broad \n C:N:P", "Broad \n C:N:P:S", "Specific \n C:N", "Specific \n C:N:P", "Specific \n C:N:P:S")) + scale_shape_manual(values=c(19,11,8,9,17,18,10,12)) + scale_alpha_manual(values=c(1, rep(0.5, 8))) + # theme(legend.position = "top") + scale_y_continuous(breaks = scales::pretty_breaks(3), labels = scales::percent_format(), limits = c(0,1)) + annotation_custom(grob = g2, xmin=0.5, xmax=3, ymin=0.01, ymax=0.7) + theme(legend.position = "top") STCON1B = toplot1_1 %>% filter(Element=="c") %>% ggplot(aes(x=Group, y=diff, group=Modification, color=Modification, linetype=Modification, alpha = Modification)) + geom_vline(xintercept = c(1.5), col="grey", linetype ="dashed") + geom_jitter(aes(shape=Modification),size=3, width=0.2) + theme_bw() + ylab("Uncertainty (% of General)") + xlab("Scenario") + scale_x_discrete(limits = c("Broad \n C:N:P:S", "CUE", "CUE & NUE", "CUE, NUE \n & PUE", "CUE, NUE, \n PUE & SUE")) + scale_shape_manual(values=c(19,11,8,9,17,18,10,12)) + scale_alpha_manual(values=c(1, rep(0.5, 8))) + # theme(legend.position = "top") + scale_y_continuous(breaks = scales::pretty_breaks(3), labels = scales::percent_format(), limits = c(0,1)) + theme(legend.position = "top") STCON2A = toplot2 %>% filter(Element=="c") %>% mutate(diff = diff/1000) %>% # convert to kg ggplot(aes(x=Group, y=diff, group=Modification, color=Modification, linetype=Modification, alpha = Modification)) + geom_vline(xintercept = c(1.5, 4.5), col="grey", linetype ="dashed") + geom_jitter(aes(shape=Modification),size=3, width=0.2) + theme_bw() + ylab(parse(text = "Uncertainty~(kg[C]~m^-2~year^-1)")) + xlab("Scenario") + scale_x_discrete(limits = c("General", "Broad \n C:N", "Broad \n C:N:P", "Broad \n C:N:P:S", "Specific \n C:N", "Specific \n C:N:P", "Specific \n C:N:P:S")) + scale_shape_manual(values=c(19,11,8,9,17,18,10,12)) + scale_alpha_manual(values=c(1, rep(0.5, 8))) + theme(legend.position = "top") + scale_y_continuous(limits=c(0,300)) STCON2B = toplot2 %>% filter(Element=="c") %>% mutate(diff = diff/1000) %>% # convert to kg ggplot(aes(x=Group, y=diff, group=Modification, color=Modification, linetype=Modification, alpha = Modification)) + geom_vline(xintercept = c(1.5), col="grey", linetype ="dashed") + geom_jitter(aes(shape=Modification),size=3, width=0.2) + theme_bw() + ylab(parse(text = "Uncertainty~(kg[C]~m^-2~year^-1)")) + xlab("Scenario") + scale_x_discrete(limits = c("Broad \n C:N:P:S", "CUE", "CUE & NUE", "CUE, NUE \n & PUE", "CUE, NUE, \n PUE & SUE")) + scale_shape_manual(values=c(19,11,8,9,17,18,10,12)) + scale_alpha_manual(values=c(1, rep(0.5, 8))) + theme(legend.position = "top") + scale_y_continuous(limits=c(0,300)) # Plot the results for the main text: png(paste0("Plots/Figure2_","c","_",Sys.Date(),".png"), width=9, height=4, units="in", res=600) ggpubr::ggarrange(STCON0A, STCON0B, labels = c("(A)", "(B)"), common.legend = T,hjust = 0) dev.off() png(paste0("Plots/Figure3_","c","_",Sys.Date(),".png"), width=9, height=4, units="in", res=600) ggpubr::ggarrange(STCON1A, STCON1B, labels = c("(A)", "(B)"), nrow=1, ncol = 2, common.legend = T,hjust = 0) dev.off() png(paste0("Plots/Figure4_","c","_",Sys.Date(),".png"), width=9, height=4, units="in", res=600) ggpubr::ggarrange(STCON2A, STCON2B, labels = c("(A)", "(B)"), nrow=1, ncol = 2, common.legend = T,hjust = 0) dev.off() # Plot the supplemental N, P, and S plots source("suppplots.R") ## Model differences across ecosystems ----- prelim_gen <- Read("fullmodel_v4.lim") pname = prelim_gen$pars$name fundecodata = read.csv("Data/fundamental_biomeXecosystem_data.csv") BXE = paste0(fundecodata$Biome, "_", fundecodata$Ecosystem) fundecodata = cbind(BXE, fundecodata[,c(-1,-2)]) # Load SRDB dataset and extract litterfall parameters srdb0 = read.csv( "Data/srdb-data-V3.csv" ) # Calculate proportion of litter: srdb0 %>% mutate(proproot = (TotDet_flux - Litter_flux)/TotDet_flux) %>% group_by(Biome, Ecosystem_type) %>% summarise(proproot = median(proproot, na.rm=T)) %>% View() fundecodata = fundecodata %>% left_join( # add the litter flux data data.frame(BXE = paste0(srdb0$Biome, "_", srdb0$Ecosystem_type), TotDet_flux = srdb0$TotDet_flux) %>% filter(!is.na(TotDet_flux)) %>% group_by(BXE) %>% summarise(TotDet_flux = mean(TotDet_flux)) ) %>% mutate(PropLitter2 = TotDet_flux/iTOT) # Transpose fundecodata = fundecodata %>% gather(key = Parameter, value = value, 2:ncol(fundecodata)) %>% spread_(key = names(fundecodata)[1],value = 'value') %>% as_tibble() # SUE parameterization # Load in Qiao data for CUE informaiton, use biome level only Qiao = read.csv("Data/Qiao_etal_2019_Data.csv") # parameterized the upper bound of SUE with the higher SUE from # Frey et al. 2011 (Glucose and Phenol), and the lower bound with # the lower SUE from the same MS SUEFrey2011 = data.frame( Substrate = rep(c("Glucose", "Glutamic acid", "Oxalic acid", "Phenol"), each=3), Type = rep(c("l", "l", "r","r"), each=3), key = rep(c("lmax", "lmin", "rmin", "rmax"), each = 3), Temp = rep(c(5,15,25), 4), SUE = c(75.7, 71.5, 70.2, 66, 59.8, 46.2, 4.5, 2.8, 3.2, 41.5, 21.3, 18.6) ) %>% mutate(SUE = SUE/100) # Groups of SUE = l vs r ; max vs min lmax = lm(SUE~Temp, data=subset(SUEFrey2011, Substrate == "Glucose")) lmin = lm(SUE~Temp, data=subset(SUEFrey2011, Substrate == "Glutamic acid")) rmax = lm(log(SUE)~Temp, data= subset(SUEFrey2011, Substrate == "Phenol")) rmin = lm(SUE~Temp, data=subset(SUEFrey2011, Substrate == "Oxalic acid")) dd = fundecodata %>% gather(-Parameter, key=BXE, value=param) %>% spread(key = Parameter, value=param) %>% select(BXE, Temp) %>% mutate(BXE2 = BXE) %>% separate(BXE2, into =c("Biome", "Ecosystem"), sep="_") %>% left_join( Qiao %>% group_by(broad.climate) %>% summarize(max = 1-min(CUE), min = 1-max(CUE)) %>% rename(Biome = broad.climate) ) %>% select(-Biome, -Ecosystem) # predict SUE based on Frey et al. data # We don't use precipitation data because the trends are context-dependent # We only constrain carbon SUE b/c the data aren't there for other SUE dd2 = dd %>% mutate(lCUEmin = (1-predict(lmax, dd)), # These are inverted so it is (1-SUE) lCUEmax = (1-predict(lmin, dd)), rCUEmin = (1-exp(predict(rmax, dd))), rCUEmax = (1-predict(rmin, dd))) %>% mutate(lCUEmin = ifelse(is.na(min), lCUEmin,ifelse(lCUEmin > min, min,lCUEmin)), lCUEmax = ifelse(is.na(min), lCUEmax,ifelse(lCUEmax < max & lCUEmax !=1, max,lCUEmax)))%>% select(-max, -min) dd2[c(4:6),3:4] = dd2[c(1:3),3:4] # Use Arctic data in Boreal # ADDED THIS TO CHECK ROBUSTNESS OF RESULTS # dd2[c(16:20),3] = 0.16 # dd2[c(16:20),4] = 0.62 dd2 = dd2 %>% gather(-BXE, -Temp, key=Type, value=param) %>% full_join( expand.grid(Species = c("B", "F", "Y", "Z", "X"), Type = c("lCUEmax","lCUEmin", "rCUEmax", "rCUEmin")) ) %>% mutate(Parameter = paste0(Species, Type)) %>% select(-Species, -Type, -Temp) # Animal CUE baseline comes from Reyes and Tiedje 1973 and Holtkamp 2011 with 0.2 for recal and 0.9 for labile. dd3 = expand.grid(Species = c("A", "S"), Qual = c("l", "r"), Type = c("CUEmin", "CUEmax"), BXE = unique(dd$BXE)) %>% mutate(Type2 = paste0(Qual, Type)) %>% left_join(data.frame( Type2 = c("lCUEmin","lCUEmax","rCUEmin","rCUEmax"), param = c(0.1, 0.5, 0.7, 0.9) # actual values before inversion --> (0.5, 0.9, 0.1, 0.3) )) %>% mutate(Parameter = paste0(Species, Qual, Type)) %>% select(BXE, Parameter, param) # Other efficiencies not well defined, so just doing broad range # Don't need to invert these b.c of symmetry! # We dont have these data, but less important if we know stoichiometry dd4 = expand.grid(Species = c("B", "F", "Y", "Z", "X","A", "S"), Qual = c("l", "r"), Element = c("NUE", "PUE", "SUE"), Type = c("min", "max"), BXE = unique(dd$BXE)) %>% mutate(param = ifelse(Type == "max", 0.95, 0.05), Parameter = paste0(Species, Qual, Element, Type)) %>% select(BXE, Parameter, param) # Put it all together into one data frame for running model fundecodatalong = fundecodata %>% gather(-Parameter, key=BXE, value=param) %>% spread(key = Parameter, value=param) %>% mutate(iLc = iTOT*PropLitter*PropLabile, iRc = iTOT*PropLitter*(1-PropLabile), iDc = (1-PropLitter -PropSymbiont)*iTOT, iYc = iTOT*PropSymbiont*PropEcto_Y, iZc = iTOT*PropSymbiont*PropAM_Z, iXc = iTOT*PropSymbiont*Propfix_X ) %>% gather(-BXE, key = Parameter, value=param) %>% bind_rows( dd2 ) %>% bind_rows( dd3 ) %>% bind_rows( dd4 ) # Invert substrate use efficiencies to reflect the parameterization of the model pdata = data.frame(pname,val = prelim_gen$pars$val) %>% as_tibble() %>% mutate(pname = as.character(pname)) %>% left_join( fundecodatalong %>% spread(key= BXE, value=param) %>% rename(pname = Parameter) %>% mutate(pname = as.character(pname)) ) %>% as.data.frame() # 68.4+102.6+76+133 pdata[c(12:108),c(3:27)] = pdata[c(12:108),c(2)] temp111 = pdata[2:3,c(-1,-2)] colnames(temp111)[grepl("Grassland", colnames(temp111)) | grepl("Agriculture", colnames(temp111)) | grepl("Savanna", colnames(temp111))] temproot = rep(0.5, dim(temp111)[2]) templeaf = rep(0.5, dim(temp111)[2]) temproot[grepl("Grassland", colnames(temp111)) | grepl("Agriculture", colnames(temp111)) | grepl("Savanna", colnames(temp111))] = 0.7 templeaf[grepl("Grassland", colnames(temp111)) | grepl("Agriculture", colnames(temp111)) | grepl("Savanna", colnames(temp111))] = 0.3 pdata[2:3,c(-1,-2)] = temp111*templeaf pdata[4:5,c(-1,-2)] = temp111*temproot rm(temp111, temproot, templeaf) # Verify order is correct for parameters all(prelim_gen$pars$name == pdata$pname) # Generate the better data for SUE and STOICH redrate = 0.45 pdataCUE = pdata aquarter = (pdata[grepl("CUEmax", pdata$pname),-1] - pdata[grepl("CUEmin", pdata$pname),-1])*redrate pdataCUE[grepl("CUEmax", pdata$pname),-1] = pdataCUE[grepl("CUEmax", pdata$pname),-1] - aquarter pdataCUE[grepl("CUEmin", pdata$pname),-1] = pdataCUE[grepl("CUEmin", pdata$pname),-1] + aquarter pdataSUE = pdata aquarter = (pdata[grepl("UEmax", pdata$pname),-1] - pdata[grepl("UEmin", pdata$pname),-1])*redrate pdataSUE[grepl("UEmax", pdata$pname),-1] = pdataSUE[grepl("UEmax", pdata$pname),-1] - aquarter pdataSUE[grepl("UEmin", pdata$pname),-1] = pdataSUE[grepl("UEmin", pdata$pname),-1] + aquarter pdataSTOICH = pdata pickmax = grepl("cnmax", pdata$pname) | grepl("cpmax", pdata$pname) | grepl("csmax", pdata$pname) pickmin = grepl("cnmin", pdata$pname) | grepl("cpmin", pdata$pname) | grepl("csmin", pdata$pname) aquarter = (pdata[pickmax,-1] - pdata[pickmin,-1])*redrate pdataSTOICH[pickmax,-1] = pdataSTOICH[pickmax,-1] - aquarter pdataSTOICH[pickmin,-1] = pdataSTOICH[pickmin,-1] + aquarter ## .....Actually run the ecosystems ------ pdatalist = list(pdata, pdataSTOICH, pdataCUE, pdataSUE) listofscenarios <- vector(mode = "list", length = 3) for(kkk in 1:3){ NN = dim(pdata)[2] - 2 listofecosystems <- vector(mode = "list", length = NN) for(ii in 1:NN){ listofecosystems[[ii]] = runecosystem(PRELIM_GEN = prelim_gen, ecosystem = ii, PDATA = pdatalist[[kkk]]) print(paste("Done", ii, "of",NN, "for scenario", kkk)) } itern = do.call("rbind", listofecosystems) itern = cbind(itern, rep(kkk, dim(itern)[1])) listofscenarios[[kkk]] = itern print(paste("Done scenario", kkk, "of 3")) } output = do.call("rbind", listofscenarios) colnames(output)[7] = "Scenario" ecosystem_out = output %>% filter(Scenario == 1) ## ....Plot results ------ toplot = output %>% as_tibble() %>% select(-diff2) %>% separate(ID, into=c("ID", "Element"),sep=-1) %>% group_by(Element, Ecosystem, Scenario) %>% summarize(par = sum(parsimonious)/1000, min = sum(min)/1000, max = sum(max)/1000) %>% filter(Ecosystem != "Baseline") %>% left_join( data.frame(Element = c("c","n","p","s"), Elementname = c("Carbon", "Nitrogen", "Phosphorus", "Sulfur")) ) %>% left_join( data.frame(Scenario = c(1,2,3,4), Scenarioname = c("Broad stoichiometry", "Reduced stoichiometry \n uncertainty", "Reduced carbon use \n efficiency uncertainty", "Reduced substrate use \n efficiency uncertainty")) ) %>% separate(Ecosystem, into=c("Biome", "Ecosystem"), sep="_") %>% mutate(XAXIS = as.numeric(as.factor(Ecosystem)) + 0.25*(as.numeric(Scenarioname)-2)) eco1 = toplot %>% # filter(Scenario != 4) %>% filter(Biome %in% c("Temperate", "Tropical")) %>% filter(Ecosystem %in% c("Forest", "Grassland")) %>% ggplot(aes(x=Ecosystem, y=par, color=Scenarioname, shape=Scenarioname)) + facet_grid(Elementname~Biome,scales="free") + # geom_point(size=3, position=position_dodge(width=0.5)) + geom_errorbar(aes(ymin = min, ymax = max), width=0.2, position=position_dodge(width=0.5)) + ylab(parse(text="Total~Cycling~(kg~m^-2~year^-1)")) + xlab("Ecosystem") + scale_color_manual(name="Constraint", values=c("black", "orange", "blue", "cyan")) + scale_shape_manual(name="Constraint", values=c(19,18,17,16)) + theme(axis.text.x = element_text(size=15, angle=0, color="black"), axis.text.y = element_text(size=15, angle=0, color="black"), axis.ticks = element_line(), text=element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(fill = "white"), panel.border = element_rect(colour = "black", fill=NA, size=0.5), legend.justification = c(0, 1), legend.position = "top",#c(0.025, 0.95), legend.background = element_rect(colour = NA), legend.key = element_rect(colour = "white", fill = NA), legend.title = element_blank(), axis.ticks.length=unit(0.25, "cm") ) + scale_y_continuous(breaks = scales::pretty_breaks(3), limits = c(0, NA)) png(paste0("Plots/Ecosystem_uncertainty_",Sys.Date(),".png"), width=14, height=12, units="in", res=300) eco1 dev.off() png(paste0("Plots/SUPP_Ecosystem_uncertainty_",Sys.Date(),".png"), width=14, height=12, units="in", res=300) toplot %>% filter(Scenario != 4) %>% ggplot(aes(x=Ecosystem, y=par, color=Scenarioname, shape=Scenarioname)) + facet_wrap(Biome~Elementname,scales="free", ncol=4) + geom_errorbar(aes(ymin = min, ymax = max), width=0.2, position=position_dodge(width=0.5)) + theme_classic() + ylab(parse(text="Total~Cycling~(kg~m^-2~year^-1)")) + xlab("Ecosystem") + theme(legend.position = "top") + scale_color_manual(name="Constraint", values=c("black", "orange", "blue", "cyan")) + scale_shape_manual(name="Constraint", values=c(19,18,17,16)) dev.off() ## Compare with SRDB (version 3.0) results ------ #Load base data srdb0 = read.csv( "Data/srdb-data-V3.csv" ) table(srdb0$Ecosystem_type) table(srdb0$Biome) table(srdb0$Biome,srdb0$Ecosystem_type) srdb = read.csv( "Data/srdb-data-V3.csv" ) %>% rename(Ecosystem = Ecosystem_type) %>% select(Biome, Ecosystem, TotDet_flux, Rh_annual) %>% mutate(Ecosystem = paste0(Biome, "_", Ecosystem)) %>% select(-Biome) %>% filter(Ecosystem %in% unique(ecosystem_out$Ecosystem)) %>% mutate(BXE = Ecosystem) %>% separate(Ecosystem, into=c("Biome", "Ecosystem"), sep="_") ecosystem_out2 = ecosystem_out %>% filter(Ecosystem %in% unique(srdb$BXE)) %>% mutate(BXE = Ecosystem) %>% separate(Ecosystem, into=c("Biome", "Ecosystem"), sep="_") pp1 = ecosystem_out2 %>% filter(grepl("Ec", ID)) %>% filter(ID != "DEc") %>% group_by(BXE, Biome, Ecosystem) %>% summarise(min = sum(min), parsimonious = sum(parsimonious), max = sum(max)) %>% ggplot(aes(x=Ecosystem)) + geom_boxplot(data=srdb,aes(y=Rh_annual)) + geom_pointrange(aes(ymin = min, y = parsimonious, ymax = max), color="red", linetype=2) + ylab("Heterotrophic respiration (g C / m2)") + facet_wrap(.~Biome, scales="free") + scale_y_continuous(breaks = scales::pretty_breaks(3)) + theme_classic() pp2 = ecosystem_out2 %>% filter(grepl("Lc", ID)|grepl("Rc", ID)) %>% filter(!(grepl("RLc", ID))) %>% group_by(BXE, Biome, Ecosystem) %>% summarise(min = sum(min), parsimonious = sum(parsimonious), max = sum(max)) %>% ggplot(aes(x=Ecosystem)) + geom_boxplot(data=srdb, aes(y=TotDet_flux)) + geom_pointrange(aes(ymin = min, y = parsimonious, ymax = max), color="red", linetype=2) + ylab("Total Litter Flux (g C / m2)") + facet_wrap(.~Biome, scales="free") + scale_y_continuous(breaks = scales::pretty_breaks(3)) + theme_classic() png("Plots/srdb_vs_model.png", width=12, height= 8, units="in", res=300) ggpubr::ggarrange(pp1,pp2, ncol = 1, nrow=2, labels = c("(A)", "(B)")) dev.off() # Flow matrix ----- lim <- Setup("fullmodel_v4.lim") nnn = substring(lim$Unknowns, 1,2) flowmat = t(lim$Flowmatrix[c(44:65, 88),c(44:65, 88)]) flowmat[flowmat !=0] = nnn[flowmat[flowmat !=0]] # Add microbial and animal respiration flowmat["EXT", "BP"] = "BE" flowmat["EXT", "FP"] = "FE" flowmat["EXT", "YP"] = "YE" flowmat["EXT", "XP"] = "XE" flowmat["EXT", "ZP"] = "ZE" flowmat["EXT", "AP"] = "AE" flowmat["EXT", "SP"] = "SE" # Mycorrhizal C inputs flowmat["ZP", "EXT"] = "EZ" flowmat["YP", "EXT"] = "EY" # Nitrogen fixer C and N inputs flowmat["XP", "EXT"] = "EX" # Nitrogen deposition flowmat["IP", "EXT"] = "EI" colmat = flowmat colmat[colmat !=0] ="black" colmat["EXT", "BP"] = "grey" colmat["EXT", "FP"] = "grey" colmat["EXT", "AP"] = "grey" colmat["EXT", "SP"] = "grey" colmat["IP",][colmat["IP",] !=0] = "salmon" colmat[,"IP"][colmat[,"IP"] !=0] = "salmon" colmat["ZP",][colmat["ZP",] !=0] = "salmon" colmat["YP",][colmat["YP",] !=0] = "salmon" colmat["ZP", "EXT"] = "grey" colmat["YP", "EXT"] = "grey" # colmat["EXT", "YP"] = "grey" # colmat["EXT", "ZP"] = "grey" colmat["IP", "EXT"] = "red" colmat["XP", "EXT"] = "rosybrown" colmat["EXT", "XP"] = "rosybrown" colmat["IP", "XP"] = "orangered" colmat["XP", "IP"] = "orange" namevec = substring(colnames(flowmat),1,1) boxcol = c("brown", "brown", "brown", "brown", "purple", "purple", "orange", "orange", "orange", "orange", "black","black", "blue", "blue", "blue", "blue", "blue", "cyan", "cyan", "cyan", "cyan", "magenta", "black") boxfill = rep("white", length(namevec)) boxfill[namevec %in% c("H", "G", "K", "T")] = "grey80" boxshape = rep("ellipse", length(namevec)) boxshape[c(5,6,13:17)] = "rect" boxshape[23] = "hexa" boxshape[c(2,4,9,10,12, 20,21)] = "diamond" # Entire System png(paste0("Plots/system_diagram_",Sys.Date(),".png"), width=12, height=10, units="in", res=300) layout(matrix(c(1,1,1,2,3,4,5,6,7), 3, 3, byrow = TRUE), heights=c(2,1,1)) par(mar=c(1,1,1,1)) plotmat(flowmat, arr.lcol = colmat, name = namevec, shadow.size = 0, box.size = 0.02, box.prop = 1, box.lcol = boxcol, box.col = boxfill, arr.length = 0.2, cex=0.01, box.type = boxshape) text(x=0.2, y=0.94, labels="A. Entire Web", cex=1.2) legend("topright", legend = c("C, N, P & S", "C & N", "N, P & S", "C","N", "P & S"), lty=1, lwd = 3, col=c("black", "rosybrown", "salmon", "grey", "red", "orange"), title="Nutrients involved") namesel = c("B. Leaf litter \n Decomposition", "C. Root litter \n Decomposition", "D. Soil animals", "E. Heterotrophic \n microbes", "F. Symbioses", "G. Mineral-associated \n organic matter") subselmat = list(c("LP", "RP"), c("VP","WP"), c("AP", "SP"), c("BP", "FP"), c("ZP", "YP", "XP"), c("MP", "JP") ) for(ii in 1:length(namesel)){ fv = subselmat[[ii]] if(length(fv) >1){ fv2 = unique(c(fv, colnames(flowmat)[(apply(flowmat[fv,] =="0", 2, all)==F)], colnames(flowmat)[(apply(flowmat[,fv] =="0", 1, all)==F)])) }else{ fv2 = unique(c(fv, colnames(flowmat)[(flowmat[fv,] !="0")], colnames(flowmat)[(flowmat[,fv] !="0")])) } subsel = seq(1,23, 1)[colnames(flowmat) %in% fv2] plotmat(flowmat[subsel,subsel], arr.lcol = colmat[subsel,subsel], arr.col = colmat[subsel,subsel], name = namevec[subsel], shadow.size = 0, box.size = 0.03, box.prop = 1, box.lcol = boxcol[subsel], box.col = boxfill[subsel], arr.length = 0.2, cex=0.01, box.type = boxshape[subsel]) text(x=0.2, y=0.94, labels=namesel[ii], cex=1.2) } dev.off()