This document covers the technical steps of the the upscaling to global scale in:
Seibold et al (2021): Deadwood decomposition mediated by insects influences global carbon fluxes.
Aboveground biomass data: GlobBiom dataset (http://globbiomass.org/), aggregated to 5min resolution
Ecological biome map: FAO http://www.fao.org/docrep/017/ap861e/ap861e00.pdf
Ratios of belowground (root) biomass to total living biomass: Robinson 2007 (https://royalsocietypublishing.org/doi/pdf/10.1098/rspb.2007.1012, Table 3)
Deadwood to total living biomass fractions (Pan et al 2011, http://www.sciencemag.org/cgi/doi/10.1126/science.1201609)
Preprocessing steps not covered here:
# load data
library(raster)
library(ggplot2)
library(gridExtra)
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.6.3
library(mgcv)
agb5m <- raster("data/globbiom_agb5min.tif")
ecozones.expanded <- raster("data/FAO_ecozones_expanded.tif") # expanded into NAs to avoid missing data
We extend the above-ground-biomass from GlobBiom with biome-specifc root biomass expansion factors, and convert total living biomass to carbon.
The expansion factors are based on Robinson (2007): https://royalsocietypublishing.org/doi/pdf/10.1098/rspb.2007.1012
these factors were also used by Liu et al 2015 (NCC): http://web.science.unsw.edu.au/~jasone/publications/liuetal2015.pdf
Living Biomass Carbon (MgC/ha) = MgBiomass/ha * 0.5 * RootExpansionFactor
biome | Above(PgC) | Below(PgC) | Total(PgC) | BGFactor | BGFrac(ofTotal) |
---|---|---|---|---|---|
tropical forest | 266 | 74 | 340 | 1.278 | 21.8% |
temperate | 109 | 30 | 139 | 1.275 | 21.6% |
Boreal | 42 | 15 | 57 | 1.357 | 26.3% |
Total | 417 | 119 | 536 | 1.285 | 22.2% |
# reclassify the original FAO Dataset, groups:
# Boreal: 1, Temperate: 5, Subtropical: 2, Tropical: 3, Polar: 9, Water: 0
ecozones.expanded[][ecozones.expanded[] %in% c(4,6,7,8)] <- 5
# table(ecozones.expanded[])
# root exapansion map
ecozone.multiplier <- ecozones.expanded
ecozone.multiplier[] <- 1 # default
# Robinson 2007, "previous best estimate" (Saugier 2001)
ecozone.multiplier[] [ecozones.expanded[] %in% c(1,9)] <- 1.357142857 * 0.5 # Boreal forests
ecozone.multiplier[] [ecozones.expanded[] %in% c(2,3) ] <- 1.278195489 * 0.5 # subtropical / tropical forests
ecozone.multiplier[] [ecozones.expanded[] == 5] <- 1.275229358 * 0.5 # Temperate
plot(ecozone.multiplier, main="Expansion Factor (AGB->TotalLivingC)")
The amount of deadwood in each cell is derived from estimates in Pan et al (2011). Pan et al. provide values for dead wood and total living biomass for countries/regions for the year 2007 (Table S3). The ratio of deadwood and total living biomass was calculated, and used to create a map (based on a shape file for the countries of the world). This map was rasterized, and some additional changes were made:
Moreover, according to Martin et al (2021) deadwood C fractions deviate from living biomass. The deadwood C fraction per biome were (Table 1):
Martin, A. R., Domke, G. M., Doraisami, M., & Thomas, S. C. (2021). Carbon fractions in the world’s dead wood. Nature Communications, 12(1), 889. https://doi.org/10.1038/s41467-021-21149-9
dwfrac.fix <- raster("data/Pan_deadwood_frac_modified.tif")
plot(dwfrac.fix, main="Deadwood fraction from Pan et al (modified)")
ocean.mask <- dwfrac.fix
ocean.mask[] <- ifelse( is.na(ocean.mask[]), 1, 0)
# Do the actual calculation
# correction factor for C fraction in deadwood (Martin et al 2021)
deadwoodC.corr <- dwfrac.fix
deadwoodC.corr[] <- 1
deadwoodC.corr[] [ecozones.expanded[] %in% c(1,9)] <- 0.4884 / 0.5 # Boreal forests
deadwoodC.corr[] [ecozones.expanded[] %in% c(2,3) ] <- 0.4716 / 0.5 # subtropical / tropical forests
deadwoodC.corr[] [ecozones.expanded[] == 5] <- 0.4929 / 0.5 # Temperate
# total Living carbon (MgC/ha)
agb5m_totalC <- agb5m * ecozone.multiplier
agb5m_totalC[agb5m_totalC[] == 0] <- NA
#writeRaster(agb5m_totalC, "total_biomassC_tha.tif", overwrite=T)
plot(agb5m_totalC, main="Living C (MgC/ha)")
dw10cmC <- agb5m_totalC * dwfrac.fix * deadwoodC.corr
plot(deadwoodC.corr)
plot(dw10cmC, main="Deadwood (>10cm) C Mg/ha")
agb.area <- area(agb5m) # area for each cell (km2)
# calc total (note: the area of each cell is taken into account!)
dw10cmC.t <- dw10cmC * agb.area * 100 / 10^9 # PgC
zDW10cm_Pg <- zonal(dw10cmC.t, ecozones.expanded, fun="sum")
totalDWPg <- sum(dw10cmC.t[], na.rm=T)
# total living carbon
agb5m_totalC.t <- agb5m_totalC * agb.area * 100 / 10^9 # PgC
zTotalLivingC_Pg <- zonal(agb5m_totalC.t, ecozones.expanded, fun="sum")
totalLivingCPg <- sum(agb5m_totalC.t[], na.rm=T)
Biome | Pan et al | Our results |
---|---|---|
Boreal | 53.9 | 60.2340642 |
Temperate | 46.6 | 53.5038084 |
Tropical | 261.1 | 272.7034687 |
World | 362.6 | 386.4802662 |
Difference to Pan et al: 6.5858428 %
Biome | Pan et al | Our results |
---|---|---|
Boreal | 16.1 | 17.7828725 |
Temperate | 3.3 | 9.3141611 |
Tropical | 53.6 | 49.3826089 |
World | 73 | 76.4801064 |
Difference to Pan et al: 4.767269 %
Calculation based on:
Christensen, M., Hahn, K., Mountford, E.P., Ódor, P., Standovár, T., Rozenbergar, D., Diaci, J., Wijdeven, S., Meyer, P., Winter, S., Vrska, T., 2005. Dead wood in European beech (Fagus sylvatica) forest reserves. For. Ecol. Manage. 210, 267–282. https://doi.org/10.1016/j.foreco.2005.02.032
We derived the expansion factor from observed deadwood distributions in Christensen et al. (2005). We found a factor of 1.21104 to convert from “deadwood >10cm” to “deadwood>2cm”, i.e. to include also smaller diameters.
dw2cmC <- dw10cmC * 1.21104
plot( dw2cmC, main="deadwood (>=2cm) tC/ha")
The approach for calculating the share of conifers/broadleaved is based on this global landcover product: Kobayashi, T., Tateishi, R., Alsaaideh, B., Sharma, R.C., Wakaizumi, T., Miyamoto, D., Bai, X., Long, B.D., Gegentana, G., Maitiniyazi, A. (2017). Production of Global Land Cover Data - GLCNMO2013. Journal of Geography and Geology, Vol. 9, No. 3, 1-15, 2017, http://dx.doi.org/10.5539/jgg.v9n3p1
https://globalmaps.github.io/glcnmo.html
The map resolution is 1/240 degree (i.e. 20x the resolution of 5min), and pixels are in one of 20 land cover classes. The map was reclassified to “Broadleaved”, “Conifer”, and “Mixed forest”, and aggregated to 5min cells. The result is a map of “conifer-share” and “broadleave-share” (summing up to 1); Mixed forests were treated as 50% conifers.
Furthermore, these maps (located in data/GLCNMO_blshare5min.tif, data/GLCNMO_conshare5min.tif) were expanded in order not to loose data points (cells that have values in the deadwood map, but which are not covered by the landcover product). See the section “Pre-processing of maps for the relative share of broadleaves and conifers” in auxiliary_calculations.R.
# load pre-processed broadleave/conifer share
bl.share.fix <- raster("data/GLCNMO_blshare5min_expanded.tif")
con.share.fix <- raster("data/GLCNMO_conshare5min_expanded.tif")
# calculate deadwood by clade
dw2cm.con.C <- dw2cmC * con.share.fix
dw2cm.bl.C <- dw2cmC * bl.share.fix
# check if the process of splitting up is not loosing something:
# plot(dw2cm.bl.C + dw2cm.con.C - dw2cmC, zlim=c(-0.0001, 0.0001), main="BL + Con = Total DW") # yeah
dw2cm.total.C <- dw2cm.bl.C+dw2cm.con.C
# this is also used as intermediate product, e.g. for the auxiliary climate calculations
writeRaster(dw2cm.total.C, "products/total_deadwood_2cm_tha.tif", overwrite=T)
writeRaster(dw2cm.bl.C, "products/total_bl_2cm_tha.tif", overwrite=T)
writeRaster(dw2cm.con.C, "products/total_con_2cm_tha.tif", overwrite=T)
Coarse deadwood decays slower than fine woody debris; our measurements from the sample plots are FWD only. We therefore developed a correction factor to decrease the decay rates for CWD.
To that end, we extracted 11 studies from the literature that measured both FWD and CWD decay over a wide ecological gradient and calculated the proportion of FWD-decay-rate / CWD-decay-rate. The values were mean= 0.54,stddev: 0.10, median: 0.526315789
We applied the median value as the CWD correction factor for the carbon pool.
Standing deadwood decays in some regions slower than downed wood. To account for this in our estimate, we assume (i) relationship between decay rates for standing and downed wood per biome, and (ii) estimate the fraction of standing deadwood on total deadwood. The assumptions are:
We used empirical data on mass loss in standing and downed logs from the temperate und boreal biome (Harmon et al 2011) to derive a factor that corrects for the slower decay in standing dead wood. Given the high variability in the data, we assumed a 50% slow down for decay in standing dead wood, representing a middle reduction of annual mass loss.
The factor \(SD_{corr}\) is applied as: \[ k_{mod} = f_{standing} \cdot SD_{corr} \cdot k_{predicted} + (1-f_{standing}) \cdot k_{predicted} \] \[ k_{mod} = k_{predicted} \cdot ( f_{standing} \cdot SD_{corr} + (1 - f_{standing}) ) \]
SD.corr <- 0.5 # the decay rate k in standing dead wood is 50% of decay in downed wood
# map with correction factors for decay rate k
sd.multiplier <- ecozones.expanded
sd.multiplier[] <- 1 # default
sd.multiplier[] [ecozones.expanded[] %in% c(1,9)] <- 0.25*SD.corr + 0.75 # Boreal forests
sd.multiplier[] [ecozones.expanded[] %in% c(2,3) ] <- 1 # subtropical / tropical forests
sd.multiplier[] [ecozones.expanded[] == 5] <- 0.3*SD.corr + 0.7 # Temperate
plot(sd.multiplier, main="decay rate reduction for standing deadwood")
To avoid extrapolation beyond the climate conditions found at our sample plots, we limited the climate data to an area in climate space given by a convex hull around our plot points. Specifically, we used for each point outside of the climatic range given by the plots the nearest point along the hull. This was done separately for gymnosperms and angiosperms (since the spatial coverage differs).
See “auxiliary_calculations.R” for details on how the modified grids were prodcued.
# load modified climate data for gymno/angio
angio.fix.temp <- raster("data/mean_temp_angio.tif")
angio.fix.map <- raster("data/map_angio.tif")
gymno.fix.temp <- raster("data/mean_temp_gymno.tif")
gymno.fix.map <- raster("data/map_gymno.tif")
cdat.angio <- data.frame(Temp = angio.fix.temp[], Prec=angio.fix.map[])
cdat.gymno <- data.frame(Temp = gymno.fix.temp[], Prec=gymno.fix.map[])
# load the "ret_n" dataframe, i.e. the plots for measurements
# the data file is produced by the Seibold_et_al_2021.R script.
load("data/ret_natural.rda")
ret_n$y_u <- with(ret_n, pmin(rate_n_uncaged, .99))
ret_n$y_c <- with(ret_n, pmin(rate_n_closed, .99))
### The gam model for upscaling:
### Temperatur in Celcius, Precipitation in cm
m_u <- gam(y_u ~ s(Temp, Prec, by = clade),
data = ret_n, family = betar)
m_c <- gam(y_c ~ s(Temp, Prec, by = clade),
data = ret_n, family = betar)
# predictModel: apply the model on a global deadwood map
# model: the statistical model (m_u, m_c)
# cl: clade ("angio", "gymno")
# clim.dat: a data frame consisting of MAT/Prec for each cell
# dwmap: deadwood map to use (DW in Mg/ha)
# cwd.decay.corr: correction factor for the decay for coarse deadwood (see explanation above)
# Returns: a map with the predicted value per cell
predictModel <- function(model, cl, clim.dat, dwmap, cwd.decay.corr = 0.526315789) {
# limit predictions to non-NA cells
predict.mask <- !is.na(dwmap[])
pred.dat <- 1 - predict.gam( model,
newdata = cbind(clim.dat, clade=cl)[predict.mask,],
type="response") # 1 - predict: loss rate
# apply the reduced decay factor (to correct for slower decay of CWD)
pred.dat <- pred.dat * cwd.decay.corr
# create a map from the predicted data
pred.map <- rep(NA, length(dwmap[]))
pred.map[predict.mask] <- pred.dat
pred.map <- setValues( dwmap, pred.map)
# apply the correction factor for standing dead wood
pred.map <- pred.map * sd.multiplier
}
# run predictions for both gymno- and angiosperms and two models (=on/off)
# the function calculates the "effect" = difference between the two given models
# returns a stack with prediction values for angio/gymno and effect.
predictAllModels <- function(model1, model2) {
# make predictions for both angio and gymnosperm
angio1 <- predictModel(model1, 'angio', cdat.angio, dw2cm.bl.C)
angio2 <- predictModel(model2, 'angio', cdat.angio, dw2cm.bl.C)
gymno1 <- predictModel(model1, 'gymno', cdat.gymno, dw2cm.con.C)
gymno2 <- predictModel(model2, 'gymno', cdat.gymno, dw2cm.con.C)
# the effect is the difference of both predictions
angio.effect <- angio1 - angio2
gymno.effect <- gymno1 - gymno2
# return as a raster stack
res <- stack(angio1, angio2, angio.effect, gymno1, gymno2, gymno.effect)
names(res) <- c("angio1", "angio2", "angioeffect", "gymno1", "gymno2", "gymnoeffect")
res
}
# total effect / worldmaps
# the input rasters are always in MgC/ha. We convert
# to absolute numbers by considering the actual area of each cell
# r.angio/r.gymno: maps for angio/gymnosperm
# returns: vector with two rasters
makeWorldMap <- function(r.angio, r.gymno) {
angio_t_ha <- dw2cm.bl.C * r.angio
gymno_t_ha <- dw2cm.con.C * r.gymno
total_t_ha <- angio_t_ha + gymno_t_ha
total_pg <- total_t_ha * agb.area * 100 / 10^9 # PgC
c("t_ha" = total_t_ha, "pg" = total_pg)
}
# Run the model (takes a while)
s_uncaged_closed <- predictAllModels( m_u, m_c )
# total decay ("m_u" model)
map.total <- makeWorldMap( s_uncaged_closed[[1]], s_uncaged_closed[[4]] )
writeRaster(map.total$t_ha, "products/total_flux_t_ha.tif", overwrite=T)
# total decay (insect effect, "m_c"-"m_u")
map.uc <- makeWorldMap( s_uncaged_closed[[3]], s_uncaged_closed[[6]] )
writeRaster(map.uc$t_ha, "products/total_uc_t_ha.tif", overwrite=T)
# total values by biome / world
rzone.total <- zonal(map.total$pg, ecozones.expanded, fun="sum")
rzone.insect <- zonal(map.uc$pg, ecozones.expanded, fun="sum")
effect.ecozone <- data.frame( zone=rep(c("Boreal", "Temperate", "Tropical"),each=3),
type=rep(c("Total", "Other", "Insect"), 3),
value=c(rzone.total[1,2]+rzone.total[5,2], rzone.total[1,2]+rzone.total[5,2]-rzone.insect[1,2]-rzone.insect[5,2], rzone.insect[1,2]+rzone.insect[5,2],
rzone.total[4,2], rzone.total[4,2]-rzone.insect[4,2], rzone.insect[4,2],
rzone.total[2,2]+rzone.total[3,2], rzone.total[2,2]+rzone.total[3,2]-rzone.insect[2,2]-rzone.insect[3,2], rzone.insect[2,2]+rzone.insect[3,2]))
Biome | Total flux (Pg) | Insect effect (Pg) |
---|---|---|
Boreal | 0.4440975 | 0.003222 |
Temperate | 0.2798908 | -0.0094478 |
Tropical | 10.2183523 | 3.1767946 |
World | 10.9424202 | 3.1705909 |
# from sensitivity analysis:
#dwPg loss.wood.Pg loss.insect.Pg loss.wood.boreal loss.wood.tropical loss.wood.temperate insect.bor insect.tr insect.temp
#14.09411 3.143694 0.8811415 0.5093355 2.805989 0.1924267 0.03125649 0.897586 0.008361141
# In order to produce the graph with error bars, we need to break the flow, and take already now
# values from the sensitivty analysis!
effect.ecozone$sd <- c(0.5093355, # total boreal
0,
0.03125649, # insect oreal
0.1924267, # total temperate
0,
0.008361141, # insect temperate
2.805989, # total tropical
0,
0.897586) # insect tropical
#pdf("products/figure3_part.pdf", width=4, height=4)
ggplot(effect.ecozone %>% filter(type!="Other"), aes(x=zone, y=value, fill=type)) +
geom_col(colour="black", position = "dodge") + geom_errorbar(aes(ymin=value-sd, ymax=value+sd), width=0.2, position=position_dodge(0.9), lwd=1) +
ylab(expression(Pg~C~a^{-1})) + theme_bw() + scale_fill_brewer(palette="BrBG") +
theme(panel.grid.minor = element_blank()) +
theme(legend.position = c(0.25, 0.8), legend.background = element_rect(fill="white", colour="darkgrey") ) +
xlab("Biome") + labs(fill="Annual deadwood\n carbon release")
#dev.off()
zones <- init(map.total$pg, v='y')
z_tha <- zonal(map.total$t_ha, zones, 'mean')
zinsect_tha <- zonal(map.uc$t_ha, zones, 'mean')
#pdf("products/tha_over_lat.pdf", width=2, height=4)
ggplot(as.data.frame(z_tha), aes(x=zone, y=mean)) + geom_area(lwd=2, fill="lightgray") +
geom_area(data=as.data.frame(zinsect_tha), fill="darkred", lwd=2) +
coord_flip() + theme_bw() + xlab("Latitude") + ylab("C released (Mg ha-1)")
## Warning: Removed 29 rows containing missing values (position_stack).
## Warning: Removed 29 rows containing missing values (position_stack).
#dev.off()