sorcering | R Documentation |
SORCERING
can be used to model the fate of soil organic carbon (SOC) and soil organic nitrogen (SON) and to calculate N mineralisation
rates.
It provides a framework that numerically solves differential equations
of SOC models based on first-order kinetics.
An SOC model can be simply defined or a predefined existing SOC model can be chosen and then run to predict the temporal development of SOC. Beyond this, SORCERING
determines the fluxes of SON and N mineralisation / immobilisation.
Basic inputs are
(1) the model parameters of a given SOC model expressed as the C transfer matrix (including information on decomposition and transfer rates between model pools),
(2) either the initial distributions of C and N among model pools as a direct input or
time series of at least three C and N measurement points with which these initial distributions
can be calculated using linear regression,
and
(3) time series of C and N inputs and rate modifying environmental factors.
In case a predefined SOC model is used, instead of model parameters and time series of rate modifying factors,
model-specific environmental and stand data must be passed for the calculation of decomposition and transfer rates.
The fourth-order Runge-Kutta algorithm is used to numerically solve the system of differential equations.
\loadmathjax
sorcering( A = NULL,
tsteps = "monthly",
C0 = NULL,
N0 = NULL,
Cin = NULL,
Nin = NULL,
Cin_wood = NULL,
Nin_wood = NULL,
wood_diam = NULL,
xi = NULL,
env_in = NULL,
site = NULL,
theta = NULL,
theta_unc = NULL,
theta_n_unc = 1,
meas_data = NULL,
A_sl = NULL,
C0_sl = NULL,
N0_sl = NULL,
Cin_sl = NULL,
Nin_sl = NULL,
Cin_wood_sl = NULL,
Nin_wood_sl = NULL,
wood_diam_sl = NULL,
xi_sl = NULL,
env_in_sl = NULL,
site_sl = NULL,
sitelist = NULL,
meas_data_sl = NULL,
calcN = FALSE,
calcNbalance = FALSE,
calcN0 = FALSE,
calcC0 = FALSE,
calcCN_fast_init = FALSE,
CTool_input_raw = FALSE,
RothC_Cin4C0 = FALSE,
C0_fracts = NULL,
multisite = FALSE,
pooltypes = NULL,
CN_fast_init = 40,
CN_bio = 9,
CN_spin = NULL,
CN_fast_init_sl = NULL,
CN_bio_sl = NULL,
CN_spin_sl = NULL,
init_info = FALSE,
model = "",
spinup = FALSE,
t_spin = 2,
t_spin_sl = 2)
A |
square matrix. Transfer matrix typical for SOC modelling. Defines number of pools, decomposition and transfer rates. n\mjeqn\timesxn elements with n = number of pools. Diagonal values are decomposition rates [yr-1]. Off-diagonals represent the transfer between pools . Only used when |
tsteps |
character string indicating the type of simulation time steps. Valid options are |
C0 |
either vector with a length equal to the number of pools or scalar. If vector, initial soil organic carbon per pool [tC ha-1].
If scalar, initial total soil organic carbon [tC ha-1]. In the latter case, either |
N0 |
vector with a length equal to the number of pools. Contains initial soil organic nitrogen per pool [tN ha-1]. If |
Cin |
either matrix with a number of columns equal to the number of pools and a number of rows corresponding to simulation time steps (if |
Nin |
either matrix with a number of columns equal to the number of pools and a number of rows corresponding to simulation time steps (if |
Cin_wood |
list of lengths of different wood diameter classes. Each list element must be in |
Nin_wood |
list of lengths of different wood diameter classes. Each list element must be in |
wood_diam |
vector with wood diameter [cm]. The first element corresponds to the first list element of |
xi |
either matrix with a number of columns equal to the number of pools and a number of rows corresponding to simulation time steps (if |
env_in |
matrix with a model-specific number of columns and a number of rows corresponding to simulation time steps (if |
site |
vector of model-specific length. Contains site-specific information to calculate rate modifying factors (instead of passing them with |
theta |
either vector with model parameters for predefined models or matrix with rows of such parameters. If it is a matrix, each row is expected to represent a stochastic repetition that covers input uncertainties. If uncertainties are defined by another argument, e.g. |
theta_unc |
either number or vector of percentage values. If it is a vector, the same model-specific lengths as described for |
theta_n_unc |
number of stochastic repetitions when model parameters for predefined models should be determined from a random distribution. Only used when the number of stochastic repetitions is not defined by another argument (e.g. |
meas_data |
matrix with a number of rows equal to the number of measurement points. The first column defines the time of measurement, the metric of which is based on simulation time steps. The second row must contain values of measured soil organic carbon stock. The third row must contain values of measured soil organic nitrogen and is only used when |
A_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
C0_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
N0_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
Cin_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
Nin_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
Cin_wood_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
Nin_wood_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
wood_diam_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
xi_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
env_in_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
site_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
sitelist |
list with names of sites to simulate. Only used when |
meas_data_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
calcN |
logical indicating whether soil organic nitrogen should be modeled. |
calcNbalance |
logical indicating whether the balance of nitrogen cycling should be calculated. |
calcN0 |
logical indicating whether |
calcC0 |
logical indicating whether |
calcCN_fast_init |
logical indicating whether to calculate the initial CN ratio for fast pools (using |
CTool_input_raw |
logical defining of which type |
RothC_Cin4C0 |
logical defining whether the SOC input should be used for the calculation of initial SOC. If |
C0_fracts |
numerical vector of a length equal to the number of pools. Contains initial fractions of SOC in pools, the sum of which must be 1. Only used when |
multisite |
logical indicating whether multiple sites should be calculated with one program call. Then, |
pooltypes |
integer vector with a length equal to the number of pools. Contains information necessary for the calculation of |
CN_fast_init |
number that defines the initial CN ratio for fast pools ( |
CN_bio |
number that defines the initial CN ratio for slow pools ( |
CN_spin |
vector with a length equal to the number of pools. Defines the initial CN ratios for spin-up runs. For the case where the spinup starts from bare ground without soil organic components, the CN ratios of the pools must be defined. Since the CN ratios would then only be influenced by external inputs, the CN ratios for slow target pools without input would exceptionally be defined by the CN ratios of fast source pools due to a lack of alternatives. To prevent this, it is necessary to define initial CN ratios for spin-up runs in that case.
Only used when |
CN_fast_init_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
CN_bio_sl |
list with a length of number of sites to simulate. Each list element represents a site and must be in |
CN_spin_sl |
list of vectors that defines the initial CN ratios for spin-up runs. Each list element represents a site and must be in |
init_info |
logical indicating whether additional information about the calculation of initial carbon, initial nitrogen, and CN ratio should be printed out during the simulations. Only used when |
model |
character string specifying a predefined soil organic carbon model to use. Valid options are |
spinup |
logical indicating whether the simulations should run in spin-up mode. Then, from all time-depending input ( |
t_spin |
integer number of spin-up time steps. |
t_spin_sl |
list with a length of number of sites to simulate or integer number of spin-up time steps. If list, each list element represents a site and must be in |
SORCERING
is a general model framework to describe soil organic carbon (SOC) dynamics and soil organic nitrogen (SON) dynamics based on models of first-order kinetics.
It can be applied to any given SOC first-order kinetics model.
The approach has already been successfully tested to describe SOC dynamics of Yasso \insertCiteTuomi2009,Viskari2020,Viskari2022sorcering, RothC \insertCiteColeman1996sorcering and
C-Tool \insertCiteTaghizadeh-Toosi2014sorcering.
Moreover, it additionally offers the possibility of modelling N immobilisation and mineralisation by enhancing given SOC models by an additional N module.
SORCERING
was created using the C++ interface Rcpp
\insertCiteEddelbuettel2021sorcering and can handle multiple sites and multiple stochastic representations with just one function call. This makes SORCERING
a computationally efficient SOC and SON modelling tool.
In the following a description of each output value (see section 'Value') is given.
Detailed mathematical descriptions of the SOC and SON calculation, the optional extensions of the SORCERING
function and the predefined models used can be found in the extended R documentation at
browseVignettes("sorcering")
.
Value C
SORCERING
calculates SOC applying a given SOC model for every simulation time step defined by passing tsteps
and the number of rows of Cin
(or number of rows of matrix elements in Cin_sl
).
SOC models applied here are defined by a number of pools, each characterised by specific decomposition and turnover rates.
The model-specific decomposition kinetics and SOC fluxes among pools are described by a set of partial differential equations represented by the transfer matrix
\mjeqnAA (as passed with A
or provided by model
).
Each row and column of \mjeqnAA represent SOC pools. Off-diagonal elements of \mjeqnAA describe SOC fluxes and diagonal elements describe SOC decomposition.
The differential equations furthermore contain the boundary condition \mjeqnCin(t)Cin(t) (as passed with Cin
) and the model-specific generated rate modifying factor series \mjeqnxi(t)xi(t)
(as passed with xi
or calculated for a predefined model
).
The change of SOC concentration in time is then defined as:
dC(t) dt = Cin(t)+A_e(t) \cdot C(t)dC(t)/dt = Cin(t) + A_e(t) * C(t)
with \mjdeqnA_e(t) = A \cdot diag(xi(t))Ae(t) = A * diag(xi(t))
Initial conditions must be defined for every SOC pool by passing C0
or by using the capabilities of SORCERING
to calculate it.
A description of the numerical solution can be found in the extended pdf documentation at browseVignettes("sorcering")
.
For more information on the functioning and possibilities of solving first-order kinetics SOC models see \insertCiteSierra2012;textualsorcering.
Value N
As an extension to SOC modelling, SORCERING
allows the modelling of SON coupled to the modelling of SOC. Its implementation is based on the following simplifying assumptions:
(1) Nitrogen transfer and turnover rates are equal to carbon rates.
(2) There is no N limitation in the soil, i.e. mineral N is always available for N immobilisation processes.
(3) CN ratios of single pools are only affected by external inputs of N and C. The transfer of organic matter among pools does not affect CN ratios.
As for SOC, the development of SON depends on initial and boundary conditions.
As N decomposition is proportional to C decomposition, SON is calculated based on the results of the SOC calculations and
input conditions (for details see the extended pdf documentation at browseVignettes("sorcering")
).
Values Nloss, Nmin, Nmin.sink<1>, ..., Nmin.sink<n>
Along with modelling SON, further quantities are determined. Nitrogen losses are calculated as:
\mjdeqnNloss(t) = N(t-1) + Nin(t-1) - N(t) Nloss(t) = N(t-1) + Nin(t-1) - N(t)
In contrast, mineralisation rates contain information about sources and sinks of SON.
They are calculated based on the CN ratios in the pools and the turnover rates (for details see the extended pdf documentation at browseVignettes("sorcering")
).
Pool-specific N mineralisation \mjeqnNmin.sink\left\langle 1 \right\rangle, ..., Nmin.sink\left\langle n \right\rangleNmin.sink_(1, ..., n) and N mineralisation \mjeqnNminNmin are related as follows:
Nmin_j(t) = \sum_p=1^n Nmin.sink \left\langle j \right\rangle_p(t) Nmin_j(t) = sum over p=1..n of Nmin.sink<j>_p(t)
for each simulation time point \mjeqntt, each pool \mjeqnj = 1, ..., nj = 1, ..., n and each pool \mjeqnp = 1, ..., np = 1, ..., n and \mjeqnnn total pools. Or in other words, the row sum of \mjeqnNmin.sink\left\langle j \right\rangleNmin.sink<j> at one simulation time point equals the jth column of \mjeqnNminNmin at that time point.
As changes in SON must match the sums of all mineralisation paths, the sums over soil pools of Nloss
and Nmin
, respectively, must be approximately equal for all simulation time points:
matrix \sum_p=1^n Nloss_p(t) \approx \sum_p=1^n Nmin_p(t)\endmatrixsum over p=1..n of Nloss_p(t) ~ sum over p=1..n of Nmin_p(t) \mjtdeqn\sum_p=1^n Nloss_p(t) \approx \sum_p=1^n Nmin_p(t)
A verification of this relation is given by "Nbalance"
(see below).
Value Nbalance
The overall N change between two time steps is calculated as: \mjdeqn \Delta N (t) = \sum_p=1^n N_p(t-1) - \sum_p=1^n N_p(t) dN(t) = sum over p=1..n of N_p(t) - sum over p=1..n of N_p(t-1)
The total system N balance serves as a verification output. Both of the following equations should always give results close to zero: \mjdeqn N_bal1(t) = \sum_p=1^n Nin_p(t-1) + \Delta N (t) - \sum_p=1^n Nloss_p(t) \approx 0 N_bal1(t) = sum over p=1..n of Nin_p(t-1) - dN(t) - sum of p=1..n of Nloss_p(t) ~ 0
\mjdeqnN_bal2(t) = \sum_p=1^n Nin_p(t-1) + \Delta N (t) - \sum_p=1^n Nmin_p(t) \approx 0 N_bal2(t) = sum over p=1..n of Nin_p(t-1) - dN(t) - sum of p=1..n of Nmin_p(t) ~ 0
\mjeqn\Delta N (t)dN (t) is saved in the first column, \mjeqnN_bal1(t)N_(bal1)(t) in the second and \mjeqnN_bal2(t)N_(bal2)(t) in the third column of "Nbalance"
.
Model parameters
If a predefined model has been specified (model is not NULL
) the following standard parameters are used.
They can be changed using theta
within the program call.
RothC | |||||
k_dpm | 10 | Decomposition rate for DPM pool [yr-1] | |||
k_rpm | 0.3 | Decomposition rate for RPM pool [yr-1] | |||
k_bio | 0.66 | Decomposition rate for BIO pool [yr-1] | |||
k_hum | 0.02 | Decomposition rate for HUM pool [yr-1] | |||
k_iom | 0 | Decomposition rate for IOM pool [yr-1] | |||
R_W_max | 1 | Maximum rate modifying factor for soil moisture | |||
R_W_min | 0.2 | Minimum rate modifying factor for soil moisture | |||
C-Tool | |||||
C-Tool | C-Tool-org | ||||
k_fom_t | 1.44 | 1.44 | Decomposition rate for FOM pool (topsoil) [yr-1] | ||
k_hum_t | 0.0336 | 0.0336 | Decomposition rate for HUM pool (topsoil) [yr-1] | ||
k_rom_t | 0.000463 | 0 | Decomposition rate for ROM pool (topsoil) [yr-1] | ||
k_fom_s | 1.44 | 1.44 | Decomposition rate for FOM pool (subsoil) [yr-1] | ||
k_hum_s | 0.0336 | 0.0336 | Decomposition rate for HUM pool (subsoil) [yr-1] | ||
k_rom_s | 0.000463 | 0 | Decomposition rate for ROM pool (subsoil) [yr-1] | ||
tf | 0.03 | 0 | Fraction going to downward transport | ||
f_co2 | 0.628 | 0.628 | Fraction of CO2 released | ||
f_rom | 0.012 | 0 | Fraction of fresh organic matter going to ROM pool | ||
f_hum | 0 | 0.358 | Fraction of input going to HUM pool | ||
Yasso | |||||
Yasso07 | Yasso15 | Yasso20 | |||
kA | 0.66 | 0.49 | 0.51 | Base decomposition rate for pool A [yr-1] | |
kW | 4.3 | 4.9 | 5.19 | Base decomposition rate for pool W [yr-1] | |
kE | 0.35 | 0.25 | 0.13 | Base decomposition rate for pool E [yr-1] | |
kN | 0.22 | 0.095 | 0.1 | Base decomposition rate for pool N [yr-1] | |
kH | 0.0033 | 0.0013 | 0.0015 | Base decomposition rate for pool H [yr-1] | |
p1 | 0.32 | 0.44 | 0.5 | Transference fraction from pool A to pool W | |
p2 | 0.01 | 0.25 | 0 | Transference fraction from pool A to pool E | |
p3 | 0.93 | 0.92 | 1 | Transference fraction from pool A to pool N | |
p4 | 0.34 | 0.99 | 1 | Transference fraction from pool W to pool A | |
p5 | 0 | 0.084 | 0.99 | Transference fraction from pool W to pool E | |
p6 | 0 | 0.011 | 0 | Transference fraction from pool W to pool N | |
p7 | 0 | 0.00061 | 0 | Transference fraction from pool E to pool A | |
p8 | 0 | 0.00048 | 0 | Transference fraction from pool E to pool W | |
p9 | 0.01 | 0.066 | 0 | Transference fraction from pool E to pool N | |
p10 | 0 | 0.00077 | 0 | Transference fraction from pool N to pool A | |
p11 | 0 | 0.1 | 0.163 | Transference fraction from pool N to pool W | |
p12 | 0.02 | 0.65 | 0 | Transference fraction from pool N to pool E | |
pH | 0.04 | 0.0046 | 0.0015 | Transference fraction from AWEN pools to pool H | |
beta_1 | 0.076 | 0.091 | 0.158 | 1st-order temperature parameter for AWE pools [degrees C-1] | |
beta_2 | -0.00089 | -0.00021 | -0.002 | 2nd-order temperature parameter for AWE pools [degrees C-2] | |
beta_N1 | - | 0.049 | 0.17 | 1st-order temperature parameter for N pool [degrees C-1] | |
beta_N2 | - | -0.000079 | -0.005 | 2nd-order temperature parameter for N pool [degrees C-2] | |
beta_H1 | - | 0.035 | 0.067 | 1st-order temperature parameter for H pool [degrees C-1] | |
beta_H2 | - | -0.00021 | 0 | 2nd-order temperature parameter for H pool [degrees C-2] | |
gamma | -1.27 | -1.8 | -1.44 | Precipitation impact parameter for AWE pools [yr mm-1] | |
gamma_N | - | -1.2 | -2 | Precipitation impact parameter for N pool [yr mm-1] | |
gamma_H | - | -13 | -6.9 | Precipitation impact parameter for H pool [yr mm-1] | |
theta_1 | - | -0.44 | -2.55 | 1st-order impact parameter for wood size [cm-1] | |
theta_2 | - | 1.3 | 1.24 | 2nd-order impact parameter for wood size [cm-2] | |
r | - | 0.26 | 0.25 | Exponent parameter for wood size |
SORCERING
returns either a list of carbon and nitrogen output values or, when multisite = TRUE
, a list broken down by site with result lists for each site.
When modelling uncertainties (as can be defined by passing e.g. Cin
, Nin
, xi
or theta
),
the output is even extended to include another list dimension that covers these uncertainties.
The lowest output list-level contains the following components:
C |
matrix with a number of rows corresponding to simulation time steps (number of rows of |
N |
matrix with a number of rows corresponding to simulation time steps (number of rows of |
Nloss |
matrix with a number of rows corresponding to simulation time steps (number of rows of |
Nmin |
matrix with a number of rows corresponding to simulation time steps (number of rows of |
Nmin.sink.1 , ... , Nmin.sink.n |
matrices with a number of rows corresponding to simulation time steps (number of rows of |
Nbalance |
matrix with a number of rows corresponding to simulation time steps (number of rows of |
The SORCERING
code was written in C++ using the R packages Rcpp
\insertCiteEddelbuettel2021sorcering
and RcppArmadillo
\insertCiteEddelbuettel2021asorcering.
This documentation was built with the help of the R packages mathjaxr
\insertCiteViechtbauer2021sorcering
and Rdpack
\insertCiteBoshnakov2021sorcering.
Marc Scherstjanoi marc.scherstjanoi@thuenen.de, Rene Dechow
#1 RothC application with fictional input for a single site
#1.1 Input
data(RothC_Cin_ex, RothC_Nin_ex, RothC_N0_ex, RothC_C0_ex, RothC_xi_ex,
RothC_site_ex, RothC_env_in_ex) #fictional data
#1.2 Simulations
#In the following two methods are presented, one with a RothC as a predefined
#model (1.2.1), one where the RothC rate modifying factors must be calculated
#beforehand (1.2.2). Both methods lead to the same results.
#1.2.1 Simulation with predefined model
out_rothC <- sorcering( model="RothC", site=RothC_site_ex, env_in=RothC_env_in_ex,
Cin=RothC_Cin_ex, Nin=RothC_Nin_ex, N0=RothC_N0_ex, C0=RothC_C0_ex,
calcN=TRUE, tsteps="monthly")
#1.2.2 Simulation with own model definition and rate modifying factor definition
A_RothC <- fget_A_RothC(clay=30) #create transfer matrix for RothC
out_rothC_own <- sorcering(A=A_RothC , xi=RothC_xi_ex, Cin=RothC_Cin_ex,
Nin=RothC_Nin_ex, N0=RothC_N0_ex, C0=RothC_C0_ex, calcN=TRUE, tsteps="monthly")
#Note that RothC_xi_ex contains site and model specific rate modifying factors that
#are only valid in this specific example. Generally, xi must be calculated by the
#user for different environmental conditions and SOC models used.
#1.3 Results
#output structure summary
summary(out_rothC)
#show that results of 1.2.1 and 1.2.2 differ negligibly
all( abs(out_rothC$C-out_rothC_own$C) < 1e-14)
all( abs(out_rothC$N-out_rothC_own$N) < 1e-14)
#example plot
oldpar <- par(no.readonly=TRUE) #save old par
par(mfrow=c(1,1),mar=c(4,4,1,4))
plot(rowSums(out_rothC$N),axes=FALSE, col=1, cex.lab=2,xlab="",ylab="",ylim=c(0,9),
pch=20)
par(new=TRUE)
plot(rowSums(RothC_Cin_ex)/rowSums(RothC_Nin_ex),
axes=FALSE,col=2, cex.lab=2,xlab="",ylab="",ylim=c(0,60),pch=20)
axis(side=2, pos=0,
labels=(0:6)*1.5, at=(0:6)*10, hadj=0.7, padj=0.5, cex.axis=2,las=1,col.axis=1)
axis(side=4, pos=60,
labels=(0:6)*10, at=(0:6)*10, hadj=0, padj=0.5, cex.axis=2, las=1,col.axis=2)
axis(side=1, pos=0,
labels= (0:6)*10 , at=(0:6)*10, hadj=0.5, padj=0, cex.axis=2)
title(ylab=expression("total N [t ha"^-1*"]"), line=2, cex.lab=2)
mtext("C input / N input", side=4, line=2, cex=2,col=2)
title(xlab="time", line=2, cex.lab=2)
par(oldpar) #back to old par
#2 RothC application with fictional input for a multiple site application
#2.1 Input
data(RothC_Cin_ex_sl, RothC_Nin_ex_sl, RothC_N0_ex, RothC_C0_ex, RothC_site_ex,
RothC_env_in_ex) #fictional data
#2.2. Simulation
out_multi_rothC <- sorcering( model="RothC", site_sl=rep(list(RothC_site_ex),3),
env_in_sl=rep(list(RothC_env_in_ex),3), Cin_sl=RothC_Cin_ex_sl,
Nin_sl=RothC_Nin_ex_sl, N0_sl=rep(list(RothC_N0_ex),3),C0_sl=rep(list(RothC_C0_ex),3),
calcN=TRUE, tsteps="monthly", multisite=TRUE,
sitelist=list("normal","half_input","double_Cin"))
#2.3 Results
#output structure summary
summary(out_multi_rothC$normal)
summary(out_multi_rothC$half_input)
summary(out_multi_rothC$double_Cin)
#example plot
oldpar <- par(no.readonly=TRUE) #save old par
par(mfrow=c(1,1),mar=c(4,4,1,4))
for (listelement in c(1:3))
{
lwidth<-1
if (listelement==2)lwidth<-3
plot(rowSums(out_multi_rothC[[listelement]]$N),axes=FALSE, col=1,type="l", lwd=lwidth,
lty=listelement+2,cex.lab=2,xlab="",ylab="",ylim=c(0,18))
par(new=TRUE)
plot(rowSums(RothC_Cin_ex_sl[[listelement]])/rowSums(RothC_Nin_ex_sl[[listelement]]),
type="l", lwd=lwidth, lty=listelement+2,axes=FALSE,col=2, cex.lab=2,xlab="",
ylab="",ylim=c(0,120))
par(new=TRUE)
}
axis(side=2, pos=0,
labels=(0:6)*3, at=(0:6)*20, hadj=0.7, padj=0.5, cex.axis=2,las=1,col.axis=1)
axis(side=4, pos=60,
labels=(0:6)*20, at=(0:6)*20, hadj=0, padj=0.5, cex.axis=2, las=1,col.axis=2)
axis(side=1, pos=0,
labels= (0:6)*10 , at=(0:6)*10, hadj=0.5, padj=0, cex.axis=2)
title(ylab=expression("total N [t ha"^-1*"]"), line=2, cex.lab=2)
mtext("C input / N input", side=4, line=2, cex=2,col=2)
title(xlab="time", line=2, cex.lab=2)
legend(x=40,y=100,legend=c("normal","half_input","double_Cin"),lty=c(3,4,5),
lwd=c(1,3,1))
par(oldpar) #back to old par
#3 RothC application with fictional input
#and fictional measurement data to calculate C0 and N0
#3.1 Input
#fictional data
data(RothC_Cin_ex_sl, RothC_Nin_ex_sl, RothC_site_ex, RothC_env_in_ex, meas_data_ex)
#3.2. Simulation
out_rothC_C0<-sorcering( model="RothC", site=RothC_site_ex, env_in=RothC_env_in_ex,
Cin=RothC_Cin_ex, Nin=RothC_Nin_ex, calcC0=TRUE, calcN=TRUE, calcN0=TRUE,
tsteps="monthly", meas_data=meas_data_ex)
#3.3 Results
#output structure summary
summary(out_rothC_C0)
#example plot
oldpar <- par(no.readonly=TRUE) #save old par
par(mfrow=c(1,1),mar=c(4,4,1,4))
plot(rowSums(out_rothC_C0$N),axes=FALSE, col=1, cex.lab=2,xlab="",ylab="",ylim=c(0,9),
type="l",lwd=1)
par(new=TRUE)
plot(rowSums(out_rothC_C0$C),axes=FALSE, col=2, cex.lab=2,xlab="",ylab="",ylim=c(0,90),
type="l",lwd=1)
par(new=TRUE)
plot(x=meas_data_ex[,1],y=meas_data_ex[,3],axes=FALSE, col=1, cex.lab=2,xlab="",ylab="",
xlim=c(0,length(rowSums(out_rothC_C0$N))),ylim=c(0,9),pch=4,cex=3)
par(new=TRUE)
plot(x=meas_data_ex[,1],y=meas_data_ex[,2],axes=FALSE, col=2, cex.lab=2,xlab="",ylab="",
xlim=c(0,length(rowSums(out_rothC_C0$N))),ylim=c(0,90),pch=4,cex=3)
par(new=TRUE)
axis(side=2, pos=0,
labels=(0:8)*1, at=(0:8)*10, hadj=1, padj=0.5, cex.axis=2,las=1,col.axis=1)
axis(side=4, pos=60,
labels=(0:8)*10, at=(0:8)*10, hadj=0, padj=0.5, cex.axis=2, las=1,col.axis=2)
axis(side=1, pos=0,
labels= (0:8)*10 , at=(0:8)*10, hadj=0.5, padj=0, cex.axis=2)
title(ylab=expression("SON [t ha"^-1*"]"), line=2, cex.lab=2)
mtext(expression("SOC [t ha"^-1*"]"), side=4, line=3, cex=2,col=2)
title(xlab="time", line=2, cex.lab=2)
legend(x=30,y=30,legend=c("model result","measurement"),lwd=c(1,0))
legend(x=31,y=30,legend=c("",""),pch=4,pt.cex=c(0,3),bty="n")
par(oldpar) #back to old par
#4 Yasso15 application using multiple sites and
#input values of different wood diameters which take uncertainties into account
#4.1 Input
data(Yasso_Cin_ex_wood_u_sl, Yasso_Nin_ex_wood_u_sl, Yasso_C0_ex_sl, Yasso_N0_ex_sl,
RothC_env_in_ex) #fictional data
#show last entries of C input for 3rd site, 2nd wood layer, 4th uncertainty layer
tail(Yasso_Cin_ex_wood_u_sl[[3]][[2]][[4]])
#diameter of wood input: 2 classes of 0 cm and 10 cm for each of the 3 sites
wood_diam_ex_sl<-list(c(0,10),c(0,10),c(0,10))
#environmental variables
Yasso_env_in_ex<-RothC_env_in_ex[,1:2]
#4.2 Simulation
out_multi_yasso_wood_unc <- sorcering( model="Yasso15", C0_sl=Yasso_C0_ex_sl,
env_in_sl=rep(list(Yasso_env_in_ex),3), wood_diam_sl=wood_diam_ex_sl,
Cin_wood_sl=Yasso_Cin_ex_wood_u_sl,Nin_wood_sl=Yasso_Nin_ex_wood_u_sl,
N0_sl=Yasso_N0_ex_sl, calcN=TRUE, tsteps="monthly", multisite=TRUE,
sitelist=list("a","b","c"))
#4.3 Results
#show the last C results for 3rd site, 4th uncertainty layer
tail(out_multi_yasso_wood_unc[[3]][[4]]$C)
#5 RothC application using stochastically varying parameters
#and multiple sites
#5.1 fictional data
data(RothC_Cin_ex_sl, RothC_Nin_ex_sl, RothC_C0_ex, RothC_N0_ex,
RothC_site_ex, RothC_env_in_ex)
#standard deviations [%] used for each of the 7 RothC theta parameters
RothC_theta_unc <- c(0,0,1,1,1,1,2)
#5.2 Simulation
out_sl <- sorcering( model="RothC", site_sl=rep(list(RothC_site_ex),3),
env_in_sl=rep(list(RothC_env_in_ex),3), Cin_sl=RothC_Cin_ex_sl,
Nin_sl=RothC_Nin_ex_sl, C0_sl=rep(list(RothC_C0_ex),3),
N0_sl=rep(list(RothC_N0_ex),3),calcN=TRUE,theta_n_unc=10,
theta_unc=RothC_theta_unc, multisite=TRUE,
sitelist=list("normal","half_input","double_Cin"))
#5.3 Means and standard deviation
#60 time steps, 5 pools, 9 output types, 10 theta_n_unc, 3 sites
out_sl_arr <- array(unlist(out_sl),c(60,5,9,10,3))
out_sl_arr_N <- out_sl_arr[,,2,,] #only output type 2: N
#mean over all uncerts
out_sl_arr_N_mean <- apply( out_sl_arr_N , c(1,2,4), na.rm=TRUE, FUN=mean )
#standard deviation
out_sl_arr_N_sd<-
array(0, dim=c(dim(out_sl_arr_N)[1],dim(out_sl_arr_N)[2],dim(out_sl_arr_N)[4]))
for (dim3 in c(1:dim(out_sl_arr_N)[4]))
out_sl_arr_N_sd[,,dim3]<-apply(out_sl_arr_N[,,,dim3],c(1:2),sd)
#5.4 Results
#show the last N means for stand 1
tail(out_sl_arr_N_mean[,,1])
#show the last N standard deviations for stand 1
tail(out_sl_arr_N_sd[,,1])
#6 How to create input lists for a RothC application using stochastically
#varying inputs and input scenarios
#6.1 Input
#fictional data
data(RothC_Cin_ex_sl, RothC_C0_ex, RothC_site_ex, RothC_env_in_ex)
#create input list of 3 scenarios, 100 uncertainties each
set.seed(17) #to make 'random' results reproducible
f1<-1
for (no in c(1:3)) #loop over 3 input scenarios
{
#normal, half and double input
Cin <- switch (no, RothC_Cin_ex, RothC_Cin_ex/2, RothC_Cin_ex*2)
f2 <- 1
#create fictional uncertainties
for (unc in c(1:100)) #loop over 100 uncertainties
{
randnum<-max(0,rnorm(1,1,0.5)) #out of normal dist. with 50% sd.
if (f2==1) Cin_u <- list(Cin*randnum) else
Cin_u[[length(Cin_u)+1]] <- Cin*randnum
f2 <- 0
}
if (f1==1) Cin_u_sl <- list(Cin_u) else
Cin_u_sl[[length(Cin_u_sl)+1]] <- Cin_u
f1 <- 0
}
#show input of scenario 3, uncertainty 51
head(Cin_u_sl[[3]][[51]])
#6.2 Simulation
out_sl <- sorcering( model="RothC", site_sl=rep(list(RothC_site_ex),3),
env_in_sl=rep(list(RothC_env_in_ex),3),Cin_sl=Cin_u_sl,
C0_sl=list(RothC_C0_ex,RothC_C0_ex,RothC_C0_ex), tsteps="monthly",
multisite=TRUE, sitelist=list("normal","half_input","double_Cin"))
#6.3 Means and standard deviation
#60 time steps, 5 pools, 1000 uncertainties, 3 sites
out_sl_arr <- array(unlist(out_sl),c(60,5,100,3))
#means
out_sl_arr_mean <- apply( out_sl_arr , c(1,2,4), na.rm=TRUE, FUN=mean )
#standard deviation
out_sl_arr_sd<-
array(0, dim=c(dim(out_sl_arr)[1],dim(out_sl_arr)[2],dim(out_sl_arr)[4]))
for (dim3 in c(1:dim(out_sl_arr)[4]))
out_sl_arr_sd[,,dim3]<-apply(out_sl_arr[,,,dim3],c(1:2),sd)
#6.4 Results
#C-pool sums of means for the 3 scenarios
totalC_m1<-rowSums(out_sl_arr_mean[,,1])
totalC_m2<-rowSums(out_sl_arr_mean[,,2])
totalC_m3<-rowSums(out_sl_arr_mean[,,3])
#C-pool sums of standard deviations for the 3 scenarios
totalC_s1<-rowSums(out_sl_arr_sd[,,1])
totalC_s2<-rowSums(out_sl_arr_sd[,,2])
totalC_s3<-rowSums(out_sl_arr_sd[,,3])
#example plot
oldpar <- par(no.readonly=TRUE) #save old par
par(mfrow=c(1,1),mar=c(4,4,1,4))
plot(totalC_m1,axes=FALSE, col=2, cex.lab=2,xlab="",ylab="",ylim=c(0,100),
type="l",lwd=1)
par(new=TRUE)
plot(totalC_m2,axes=FALSE, col=3, cex.lab=2,xlab="",ylab="",ylim=c(0,100),
type="l",lwd=1)
par(new=TRUE)
plot(totalC_m3,axes=FALSE, col=4, cex.lab=2,xlab="",ylab="",ylim=c(0,100),
type="l",lwd=1)
par(new=TRUE)
polygon(c(1:60,60:1),c(totalC_m1+totalC_s1, rev(totalC_m1-totalC_s1)),
border=NA,col=rgb(1,0,0,0.27),density=40,angle=180,xlab="",ylab="")
par(new=TRUE)
polygon(c(1:60,60:1),c(totalC_m2+totalC_s2, rev(totalC_m2-totalC_s2)),
border=NA,col=rgb(0,1,0,0.27),density=30,xlab="",ylab="")
par(new=TRUE)
polygon(c(1:60,60:1),c(totalC_m3+totalC_s3, rev(totalC_m3-totalC_s3)),
border=NA,col=rgb(0,0,1,0.27),density=25,angle=90,xlab="",ylab="")
par(new=TRUE)
axis(side=2, pos=0,
labels=(0:10)*1, at=(0:10)*10, hadj=1, padj=0.5, cex.axis=2,las=1,col.axis=1)
axis(side=1, pos=0,
labels= (0:6)*10 , at=(0:6)*10, hadj=0.5, padj=0, cex.axis=2)
title(ylab=expression("SOC [t ha"^-1*"]"), line=2, cex.lab=2)
title(xlab="time", line=2, cex.lab=2)
legend(x=20,y=30,fill=c(0,0,0,4,2,3),density=c(0,0,0,25,40,30),angle=c(0,0,0,90,0,45),
border=c(0,0,0,1,1,1),legend=c("mean double input scenario",
"mean regular input scenario", "mean half input scneario",
"uncertainty range double input scenario", "uncertainty range regular input scenario",
"uncertainty range half input scenario"))
legend(x=20,y=30,lty=c(1,1,1,0,0,0),seg.len=c(1,1,1,0,0,0), col=c(4,2,3,0,0,0),
legend=c("","","","","",""),bty="n")
par(oldpar) #back to old par
#7 RothC application with fictional input for a spin-up application
#7.1 Input
#fictional data
data(RothC_Cin_ex_sl_spin, RothC_Nin_ex_sl_spin, RothC_site_ex, RothC_env_in_ex)
#7.2. Simulation
out_multi_rothC <- sorcering( model="RothC", site_sl=rep(list(RothC_site_ex),3),
env_in_sl=rep(list(RothC_env_in_ex[1:12,]),3), Cin_sl=RothC_Cin_ex_sl_spin,
Nin_sl=RothC_Nin_ex_sl_spin, calcN=TRUE, tsteps="monthly", multisite=TRUE,
sitelist=list("normal","half_input","double_Cin"), spinup=TRUE, t_spin_sl=36000,
C0=c(0,0,0,0,20), N0=c(0,0,0,0,2), CN_spin=c(100,100,50,50,10))
#7.3 Results
#example plot
oldpar <- par(no.readonly=TRUE) #save old par
par(mfrow=c(1,1),mar=c(4,4,1,4))
for (listelement in c(1:3))
{
lwidth<-1
if (listelement==2)lwidth<-3
printN<-rowSums(out_multi_rothC[[listelement]]$N)
printseq<-seq.int(1L,length(printN),100L)
printC<-rowSums(out_multi_rothC[[listelement]]$C)
plot(printN[printseq],axes=FALSE, col=1,type="l", lwd=lwidth,
lty=listelement+2,cex.lab=2,xlab="",ylab="",ylim=c(0,30))
par(new=TRUE)
plot(printC[printseq],axes=FALSE, col=2,type="l", lwd=lwidth,
lty=listelement+2,cex.lab=2,xlab="",ylab="",ylim=c(0,180))
par(new=TRUE)
}
axis(side=2, pos=0,
labels=(0:6)*5, at=(0:6)*30, hadj=0.7, padj=0.5, cex.axis=2,las=1,col.axis=1)
axis(side=4, pos=360,
labels=(0:6)*30, at=(0:6)*30, hadj=0, padj=0.5, cex.axis=2, las=1,col.axis=2)
axis(side=1, pos=0,
labels= (0:6)*6000 , at=(0:6)*60, hadj=0.5, padj=0, cex.axis=2)
title(ylab=expression("total N [t ha"^-1*"]"), line=2, cex.lab=2)
mtext(expression("total C [t ha"^-1*"]"), side=4, line=2, cex=2,col=2)
title(xlab="time", line=2, cex.lab=2)
legend(x=120,y=140,legend=c("normal","half_input","double_Cin"),lty=c(3,4,5),
lwd=c(1,3,1))
par(oldpar) #back to old par
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.