#' Mass of living roots
#'
#' This function returns the mass of the living roots between two layers. Given a mass per area of roots and a maximum rooting depth, this calculates the total mass between two layers with a specified distribution. Currently only a linear distribution is implemented. \cr
#' The linear algorithm is as follows: \itemize{
#' \item #mass_per_depth = slope * depth + intercept
#' \item slope <- -2 * totalRootMass / (rootDepthMax^2)
#' \item intercept <- 2 * totalRootMass / rootDepthMax
#' \item #mass = integral(mass_per_depth, depth)
#' \item rootMass <- intercept * (layerBottom-layerTop) + slope/2 * (layerBottom ^2-layerTop^2)
#' }
#'
#' @param layerBottom an array of depths from soil top to the bottom of the layer, generally in cm
#' @param layerTop an array of depths from soil top to the top of the layer, generally in cm
#' @param totalRootMassPerArea an integer that is the total mass per area of the roots, generally in g cm-3
#' @param soilLength unit length of soil area, generally 1
#' @param soilWidth unit width of soil area, generally 1.
#' @param shape flag of root shape, only \code{linear} is implimented.
#' @param expDecayRatePerMaxDepth a numeric, value for exponential decay rate per maximum root depth
#' @param ...
#' @param rootDepthMax an integer that is the maximum root depth, generally in cm
#'
#' @return an array of the associated live root mass for each layer specified.
#' @export
#'
#' @examples
#' massLiveRoots(layerBottom = 1:10, layerTop = 0:9,
#' totalRootMassPerArea = 0.3, rootDepthMax=5,
#' soilLength=1, soilWidth=1, shape='linear')
massLiveRoots <- function(layerBottom, layerTop,
totalRootMassPerArea,
rootDepthMax,
soilLength=1, soilWidth=1,
shape,
expDecayRatePerMaxDepth = log(0.05),
...){
totalRootMass <- soilLength*soilWidth*totalRootMassPerArea
if (totalRootMass == 0) {
rootMass <- rep(0, length(layerBottom))
return(rootMass)
} else {
layerBottom[is.na(layerBottom)] <- 0
layerTop[is.na(layerTop)] <- 0
if(any(layerBottom < layerTop)){
#print(data.frame(layerBottom, layerTop))
stop('Bad layer definition.')
}
##reset the layers that are beyond the root inputs to have 0 depths
layerBottom[layerBottom > rootDepthMax] <- rootDepthMax
layerTop[layerTop > rootDepthMax] <- rootDepthMax
if(shape == 'linear'){
#mass_per_depth(depth) = slope * depth + intercept; mass(depth) = slope/2*depth^2+intercept*depth
#Given
# mass_per_depth(maxDepth) = 0 => intercept = -slope * maxDepth
# mass(maxDepth) = TotalMass =>
# TotalMass = slope/2*maxDepth^2 + intercept*maxDepth
# => slope = 2*TotalMass/maxDepth^2
# intercept = -2 * TotalMass/maxDepth
#
slope <- -2 * totalRootMass / (rootDepthMax^2)
intercept <- 2 * totalRootMass / rootDepthMax
#mass = integral(mass_per_depth, depth)
rootMass <- intercept*(layerBottom-layerTop) + slope/2*(layerBottom ^2-layerTop^2)
}else{
if(shape == 'exponential'){
# root mass with depth x: r(x) = a * exp(b * x) - m
# Let R(x) be the antiderivative: R(x) = a / b * exp(b * x) - m * x + C0
# Assume
# r( x = rootDepthMax) => m = a * exp(b * rootDepthMax)
# Total root biomass is integral of r from 0 to rootDepthMax
# TotalMass = a / b * exp(b*rootDepthMax) - m * rootDepthMax - a / b
# = a / b * exp(b*rootDepthMax) - a * exp(b * rootDepthMax) * rootDepthMax - a / b
# a = TotalMass * [1 / b * exp(b * rootDepthMax) - exp(b * rootDepthMax) * rootDepthMax - 1 / b]^-1
#
# mass between two layers is the integral of r from x1 to x2
# mass = (a / b * exp(b * x2) - m * x2) - (a / b * exp(b * x1) - m * x1)
# mass = a / b * (exp(b * x2)-exp(b * x1)) + m *(x1 - x2)
#layerTop <- 0; layerBottom = 1; expDecayRate_perRootDepthMax <- log(0.05); totalRootMass <- 0.2; rootDepthMax <- 30
b <- expDecayRatePerMaxDepth / rootDepthMax #convert from total profile to per cm
a <- totalRootMass * (1 / b * exp(b * rootDepthMax) -
exp(b * rootDepthMax) * rootDepthMax - 1 / b)^-1
m <- a * exp(b * rootDepthMax)
rootMass <- a / b * (exp(b * layerBottom)-exp(b * layerTop)) + m*(layerTop - layerBottom)
}else{
stop(paste('Unknown shape specified:', shape))
}
}
return(rootMass)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.