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.

Table of contents

Calculation of carbon pools

Description of input data

Preprocessing

Preprocessing steps not covered here:

  • AGB Map
    • convert to resolution of WorldClim 5 Minutes (i.e. 0.08333333 deg, 1920x4320); values in Mg Biomass / ha
  • Ecological zones
    • rasterized to WorldClim 5 Minutes
    • expand into NA-regions to avoid missing data
# 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

Calculation of total living biomass

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

  • 0.5: Initial Conversion Biomass -> Carbon (Note: will be updated below)
  • RootExpansionFactor: Values from Robinson (2007), LivingC = AbovegroundC * RootExpansionFactor
The biome-specific factors to caculate total living biomass carbon from aboveground biomass.
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)")

Calculation of deadwood pools

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:

  • Alaska was assigned the ratio of Canada instead of the value for US
  • Russia was split in an European and Asian part, Kalinigrad was assigned the average value of Europe
  • the value for China was changed, since the value in Pan et al is dubiously low. We used the value from FRA 2015 (for the year 2007), 15.319%.
  • Biome mean values were used for countries with missing data

Moreover, according to Martin et al (2021) deadwood C fractions deviate from living biomass. The deadwood C fraction per biome were (Table 1):

  • Boreal: 48.84
  • Subtropical / Mediterrean: 46.24 (not used due to small sample size)
  • Temperate : 49.29
  • Tropical: 47.16

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)

Carbon pools summary

Total Living Carbon (including roots, PgC)

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 %

Deadwood (>10cm, PgC)

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 %

Preparations for upscaling

Including deadwood with smaller DBH

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")

Splitting of deadwood into broadleaves / conifers

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)

Correcting for deadwood diameter

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.

Correcting for slower decay in standing deadwood

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:

  • Tropics: decay in standing deadwood has the same speed as downed deadwood. No correction is applied.
  • Temperate forests: fraction of standing dead wood between 20% and 40%, we assume here 30%
  • Boreal forests: fraction of standing dead wood between 20% and 30%, we assume here 25%

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")

Global climate data

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[])

Upscaling

load data and define functions

# 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

# 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)

Summarise per biome

# 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]))

Results

Total deadwood decomposition and insect effect

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

Figures

# 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()

Figure of total loss / insect affect along a latitudinal gradient

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()