#' @title calcSOM
#' @description calculates Soil Organic Matter Pool, accounting for the management history.
#' We assume carbon Stocks from LPJml natural vegetation as a starting point.
#' Here we use the upper 30cm soil layer (0-20cm of + 1/3 of 30-50 cm).
#' We then correct carbon pools by lost c-share depending on the climate region, using default
#' factors of IPCC Guidelines 2006 table 5.5.
#' We assume that this IPCC-corrected value is the target long-term equilibrium value for the soil stocks.
#' Because soil decline and build-up slowly, we assume that in every year, the carbon pools move 15% towards
#' this new equilibrium.
#' This assumption is in line with IPCC saying that the process will take 20 years: with our assumption,
#' after 5 years 44% of the carbon pool is gone, after 10 years 80% and after 20 years 96%.
#' We determine a carbon stock for cropland soils and non-cropland soils in every cell.
#' If the cropland area expands, the carbon stock of noncropland is proportionally assigned to the cropland pool
#' and vice versa. The outputs of the function are the soilc stocks for cropland and non-cropland.
#' Relevant for the release of N by SOM loss is also the change in carbon stocks per ha, as this relases or binds N.
#' This is done in delta cropland soilc.
#' @param climatetype Switch between different climate scenarios (default on "historical")
#' @param subtype "stock" (default) for absoulte values, "density" for per hectar values
#' @param cells "magpiecell" for 59199 cells or "lpjcell" for 67420 cells
#'
#' @return List of magpie object with results on country or cellular level,
#' weight on cellular level, unit and description.
#' @author Benjamin Leon Bodirsky, Kristine Karstens
#' @examples
#' \dontrun{
#' calcOutput("SOM2")
#' }
#'
calcSOM <- function(climatetype = "historical", subtype = "stock", cells = "lpjcell") {
years <- seq(1951, 2010, 1)
soilc <- calcOutput("LPJmL_new", version = "LPJmL4_for_MAgPIE_44ac93de",
climatetype = "GSWP3-W5E5:historical", subtype = "soilc_layer",
stage = "raw", aggregate = FALSE, years = years)
soilc <- setNames(soilc[, , 1] + 1 / 3 * soilc[, , 2], "soilc")
states <- readSource("LUH2v2", subtype = "states", convert = "onlycorrect")[, years, ]
crops <- c("c3ann", "c4ann", "c3per", "c4per", "c3nfx")
cropArea <- dimSums(states[, , crops], dim = 3)
noncropArea <- dimSums(states, dim = 3) - cropArea
rm(states)
cropshare <- toolFillYears(calcOutput("Croparea", sectoral = "kcr", physical = TRUE, cells = "lpjcell",
cellular = TRUE, irrigation = FALSE, aggregate = FALSE), years)
cropshare <- toolConditionalReplace(cropshare / dimSums(cropshare, dim = 3), "is.na()", 0)
carbshare <- calcOutput("SOCLossShare", aggregate = FALSE, subsystems = TRUE,
rate = "change", factor = "ipccReduced", cells = "lpjcell")
cshare <- dimSums(cropshare * carbshare, dim = 3)
cshare[cshare == 0] <- 1 # target for cropland in cells without cropland equal to nat veg just as backup.
# in principle possible to add begr/betr area based on LUH2v2
# crpbf_c3per: C3 perennial crops grown as biofuels
# crpbf_c4per: C4 perennial crops grown as biofuels
targetCcrop <- soilc * cshare * cropArea
targetCNoncrop <- soilc * noncropArea
transitions <- cropArea
stopifnot(length(years) >= 2)
transitions[, years[2:length(years)], ] <- (cropArea[, years[2:length(years)], ]
- setYears(cropArea[, years[seq_along(years) - 1], ],
years[2:length(years)]))
abandonnedland <- newland <- transitions
abandonnedland[abandonnedland > 0] <- 0
abandonnedland <- abandonnedland * (-1)
newland[newland < 0] <- 0
cropC <- cropCha <- deltaCcrop <- targetCcrop
noncropC <- noncropCha <- deltaCnoncrop <- targetCNoncrop
cropC[, 2:length(years), ] <- NA
noncropC[, 2:length(years), ] <- NA
cropCha[, , ] <- deltaCcrop[, , ] <- NA
noncropCha[, , ] <- deltaCnoncrop[, , ] <- NA
cropCha[, 1, ] <- cropC[, 1, ] / cropArea[, 1, ]
cropCha[is.nan(cropCha)] <- 0
noncropCha[, 1, ] <- noncropC[, 1, ] / noncropArea[, 1, ]
noncropCha[is.nan(noncropCha)] <- 0
for (yearX in (2:length(years))) {
cropC[, yearX, ] <- (setYears(cropC[, yearX - 1, ], NULL)
+ newland[, yearX, ] * setYears(noncropCha[, yearX - 1, ], NULL)
- abandonnedland[, yearX, ] * setYears(cropCha[, yearX - 1, ], NULL))
noncropC[, yearX, ] <- (setYears(noncropC[, yearX - 1, ], NULL)
- newland[, yearX, ] * setYears(noncropCha[, yearX - 1, ], NULL)
+ abandonnedland[, yearX, ] * setYears(cropCha[, yearX - 1, ], NULL))
# assumption on transition: 15% of the soil difference per year.
# 44% after 5 years, 20% after 10 years, 4% after 20 years
deltaCcrop[, yearX, ] <- (targetCcrop[, yearX, ] - cropC[, yearX, ]) * 0.15
deltaCnoncrop[, yearX, ] <- (targetCNoncrop[, yearX, ] - noncropC[, yearX, ]) * 0.15
# to avoid infs in division, a rounding is required
cropC[, yearX, ] <- round(cropC[, yearX, ] + deltaCcrop[, yearX, ], 10)
noncropC[, yearX, ] <- round(noncropC[, yearX, ] + deltaCnoncrop[, yearX, ], 10)
cropCha[, yearX, ] <- setYears(cropC[, yearX, ] / cropArea[, yearX, ], NULL)
cropCha[is.nan(cropCha)] <- 0
cropCha[abs(cropCha) == Inf] <- 0
noncropCha[, yearX, ] <- setYears(noncropC[, yearX, ] / noncropArea[, yearX, ], NULL)
noncropCha[is.nan(noncropCha)] <- 0
noncropCha[abs(noncropCha) == Inf] <- 0
}
# deltaC is not equivalent to the difference in carbon_cropland_soils over time, as the area changes
if (subtype == "stock") {
out <- mbind(setNames(cropC, "cropland.soilc"),
setNames(noncropC, "noncropland.soilc"),
setNames(deltaCcrop, "cropland.delta_soilc"),
setNames(deltaCnoncrop, "noncropland.delta_soilc"),
setNames(targetCcrop, "cropland.target_soilc"),
setNames(targetCNoncrop, "noncropland.target_soilc"))
unit <- "Mt C"
weight <- NULL
} else if (subtype == "density") {
deltaCCropHa <- toolNAreplace(deltaCcrop / cropArea)$x
deltaCNoncropHa <- toolNAreplace(deltaCnoncrop / noncropArea)$x
targetCCropHa <- toolNAreplace(targetCcrop / cropArea)$x
targetCNoncropHa <- toolNAreplace(targetCNoncrop / noncropArea)$x
out <- mbind(setNames(cropCha, "cropland.soilc"),
setNames(noncropCha, "noncropland.soilc"),
setNames(deltaCCropHa, "cropland.delta_soilc"),
setNames(deltaCNoncropHa, "noncropland.delta_soilc"),
setNames(targetCCropHa, "cropland.target_soilc"),
setNames(targetCNoncropHa, "noncropland.target_soilc"))
unit <- "t C per ha"
weight <- mbind(setNames(cropArea, "cropland"),
setNames(noncropArea, "noncropland"))[, -c(1:10), ]
} else {
stop(paste("Subtype", subtype, "does not exist yet."))
}
# delete first 20 years of spin-up
out <- out[, -c(1:10), ]
if (cells == "magpiecell") out <- toolCoord2Isocell(out)
return(list(
x = out,
weight = weight,
unit = unit,
description = paste("Carbon in cropland and non-cropland soils, as well as change",
"over time due to built-up or loss. Change is not equivalen to",
"the difference in carbon_cropland_soils over time, as the area changes."),
isocountries = FALSE))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.