Modeling radiocarbon dynamics in soils: SoilR version 1.1

Radiocarbon is an important tracer of the global carbon cycle that helps to understand carbon dynamics in soils. It is useful to estimate rates of organic matter cycling as well as the mean residence or transit time of carbon in soils. We included a set of functions to model the fate of radiocarbon in soil organic matter within the SoilR pack-5 age for the R environment for computing. Here we present the main system equations and functions to calculate the transfer and release of radiocarbon from di ﬀ erent soil organic matter pools. Similarly, we present functions to calculate the mean transit time for di ﬀ erent pools and the entire soil system. This new version of SoilR also includes a group of datasets describing the amount of radiocarbon in the atmosphere over time, 10 data necessary to estimate the incorporation of radiocarbon in soils. Also, we present examples on how to obtain parameters of pool-based models from radiocarbon data using inverse parameter estimation. This implementation is general enough so it can also be used to trace the incorporation of radiocarbon in other natural systems that can be represented as linear dynamical systems. 15


Introduction
To study the global carbon cycle and its interaction with climate, it is necessary to develop models that can accurately represent the size and the amount of transfers among different C reservoirs within the Earth system.Soils are one of the most important C reservoirs, storing between 800 to 1700 Pg C in the first 1 m, and exchanging between 53-57 Pg C yr −1 with the atmosphere in the form of heterotrophic respiration (Schlesinger and Andrews, 2000;Lal, 2004;Bond-Lamberty and Thomson, 2010;Todd-Brown et al., 2013).However, there are large uncertainties in these estimations, which are related to uncertainties in C stocks of arctic peatlands, coarse woody debris, and C stocks below topsoil (Jobbágy and Jackson, 2000;Harmon et al., 2011;Todd-Brown et al., 2013).It is also highly debated whether climate change may destabilize Figures current soil C stocks (Trumbore, 1997;Schlesinger and Andrews, 2000;Kirschbaum, 2006;Davidson and Janssens, 2006;von Lützow and Kögel-Knabner, 2009;Conant et al., 2011;Sierra, 2012).
Radiocarbon can be used as a tracer of the interactions between terrestrial ecosystems and the atmosphere, and provides information about the rates of carbon inputs and losses from soils (Trumbore, 2009).Radiocarbon is a cosmogenic radionuclide that is constantly produced in the upper layers of the stratosphere.In the lower atmosphere, the amount of radiocarbon at any given time is given by the balance between cosmogenic production, radioactive decay, and sources and sinks from oceans, and the terrestrial biosphere.Atmospheric concentrations of radiocarbon are well known for the past 7000-1000 years, and the continuous record even extends to 50 000 years into the past (Reimer et al., 2009(Reimer et al., , 2013)).Therefore, it is possible to know with good precision when a C atom entered the terrestrial biosphere and for how long it has been stored in a terrestrial reservoir.
Radiocarbon is also used in tracer studies in which known amounts of radiocarbon label are introduced in vegetation or soils and its fate is followed as it moves among different compartments and subsequently leaves the system.During the late 1950s and early 1960s nuclear weapon tests considerably increased the amount of radiocarbon in the atmosphere, creating a global-scale labeling experiment that allows researchers to follow the fate of this spike in atmospheric radiocarbon concentrations across many different reservoirs of the biosphere.In soils, radiocarbon studies have proved useful for estimating the residence times of carbon in organic matter that cycles on time-scales ranging from years to millennia (Trumbore, 2009).Organic matter is subject to different transformation processes in soils, it can be quickly consumed by microorganisms once it enters the soil, it can be transformed into different compounds as a result of microbial-mediated reactions, or it can also react with soil mineral surfaces (Sollins et al., 1996;Schmidt et al., 2011;Gleixner, 2013).These different processes create a heterogeneity of rates of organic matter decomposition that are of fundamental importance in determining long-term Introduction

Conclusions References
Tables Figures
In this manuscript, we present the implementation of the radiocarbon component within the SoilR package, a software tool developed for modeling soil organic matter dynamics (Sierra et al., 2012a).First, we present the mathematics behind the new implementation.Then, we present some details about the numerical implementation in R and the particular functions implemented in SoilR.At the end of the manuscript, we present some particular examples about its use.

General radiocarbon model
Previously, we have defined a general model of soil organic matter decomposition as a linear dynamical system of the form (Sierra et al., 2012a) where the amount of carbon in different pools is represented as a vector C(t), with total inputs of carbon represented by the vector I(t).Similarly, the dynamical system for radiocarbon in soil organic matter can be represented as where the amount of radiocarbon in each pool i is represented by the vector 14 C(t), with radiocarbon inputs represented by I14 C (t), and λ as the radioactive decay constant.
Both I14 C (t) and 14 C(t) represent the total amount of radiocarbon in a sample in relation to an international standard (Stuiver and Polach, 1977).
The fate of radiocarbon in soils can also be described in fractional form as where F (t) is a vector of length m and • represents the entry-wise product between the two vectors.The fraction F (t) represents the activity ratio of a sample with respect to a reference material (see Sect. 2.2 for details, and Stuiver and Polach, 1977;Mook and Van Der Plicht, 1999).The system of equations can therefore be expressed as where F a (t) is a scalar value that represents the fraction of radiocarbon in the atmosphere, which is not constant and has changed considerably over time due to the action of cosmic rays, the storage and release of carbon from oceans and the biosphere, and human activities (Reimer et al., 2009;Levin et al., 2010;Reimer, 2012).
In SoilR, we compute the time-dependent solution of Eq. ( 4), solving for F (t) using standard numerical methods (see Sect. 2.4.1).F (t) contains the radiocarbon fraction for each pool i for a given time (t).
We are also interested in calculating the total radiocarbon in soil organic matter weighted by its mass F C (t), and the total amount of released radiocarbon weighted Figures

Back Close
Full by the total amount of released carbon F R (t).These weighted averages, or expectations, can be related to the average radiocarbon content of a soil sample and the average radiocarbon content of the released (respired) carbon from a sample, respectively.Mathematically, both concepts can be expressed as and respectively.In both equations the sum is over all pools at each time t.

Reporting radiocarbon
In reporting radiocarbon, there are different ways to refer to the proportion of radiocarbon in a sample.Atmospheric radiocarbon data for the pre-bomb period is commonly reported as ∆ 14 C (Reimer et al., 2013), which is defined according to Stuiver and Polach (1977) as with where A SN represents the activity of a sample normalized for 13 C fractionation, and A ABS the activity of the oxalic acid standard normalized for 13 C fractionation and corrected for decay since 1950.Introduction

Conclusions References
Tables Figures

Back Close
Full For post-bomb applications, radiocarbon is better expressed as F 14 C, which according to Reimer et al. (2004) is expressed as where A ON is the activity of the oxalic acid standard with 13 C normalization, but without decay correction; i.e. y−1950) .( 10) Hua et al. (2013) report atmospheric radiocarbon values for the post-bomb period as F 14 C and as ∆ 14 C, the later expressed as i.e., the activity of the standard does not change with time during the post-bomb period.
As both representations of ∆ 14 C (Eqs. 7 and 11) are algebraically similar, we take both types of ∆ 14 C values and treat them equally in our calculations.
We define an absolute fraction modern F value as where ∆ 14 C is expressed as Eq. ( 7) for radiocarbon data previous to 1950, and as Eq. ( 11) after 1950.The system of differential equations of Eq. ( 4) is solved using the values of F as previously described.

Definitions and assumptions
A commonly used metric to compare different compartment models is the concept of mean transit time, also known as mean residence time (Eriksson, 1971; Bolin and Introduction

Conclusions References
Tables Figures

Back Close
Full  , 1973;Nir and Lewis, 1975;Thompson and Randerson, 1999;Manzoni et al., 2009).In previous studies, the mean transit time of a system has been defined as the average time a particle of carbon spends in the system from entry to exit.This definition however, has been proposed for linear time invariant (LTI) systems in which the solution does not change over time and the system is in steady-state.This contrast with the more general models that SoilR can solve (Eqs. 1 and 2) that allow time dependent input fluxes and decomposition rates.In addition, this definition of transit times does not specify the set of particles whose transit times contribute to the average, suggesting an average over all particles in the system.

Rodhe
Here we provide a more general definition of mean transit time that takes into account the more general models that SoilR can solve and specifies the set of particles used for calculating the average.Our formal definition states: Given a system described by the complete history of inputs I(t) for t ∈ (t start , t 0 ) to all pools until time t 0 and the cumulative output O(t 0 ) of all pools at time t 0 the mean transit time Tt 0 of the system at time t 0 is the average of the transit times of all particles leaving the system at time t 0 .
Accordingly, we define the related density distribution: Given a system described by the complete history of inputs I(t) for t ∈ (t start , t 0 ) to all pools until time t 0 and the cumulative output O(t 0 ) of all pools at time t 0 the transit time density ψ t 0 (T ) of the system at time t 0 is the probability density with respect to T implicitly defined by Methods for calculating the mean transit time and transit time density for the general case and the models of the form of Eqs. ( 1) or (4) will be described in a forthcoming more detailed publication.Here we will limit to describe the most common calculation of mean transit time for the LTI case, i.e. for models in steady-state (total inputs are equal to total outputs), constant coefficients, and constant inputs.The general form of Introduction

Conclusions References
Tables Figures

Back Close
Full these LTI models, a special case of Eq. ( 1), is given by

Implementation
For the LTI case, it has been shown previously that the transit time density distribution ψ(T ) for a transit time T is identical to the output O(T ) observed at time T of a different system which started with a normalized impulsive input I I at time T = 0 (Nir and Lewis, 1975;Manzoni et al., 2009); where I represents the sum of all elements of the vector I. Translated to the language of an ODE solver, an impulsive input becomes a vector of initial conditions I I at time T = 0, and S r the release flux of the solution of the initial value problem observed at time T Note that from the perspective of the ode solver, S r depends only on the decomposition operator A and the distribution of the input among the pools (Eq.14).It is therefore possible to implement the transit time distribution as a function only of the decomposition operator and the fixed input flux distribution.To insure steady-state conditions the decomposition operator is not allowed to be a true function of time.We therefore implement the method only for the subclass ConstantDecompositionOperator, a new native class of SoilR objects for the time invariant decomposition operator A.
To compute the mean transit time for the distribution we need to compute the integral However, to avoid issues with numerical integration, we do not use ∞ as upper limit of integration, but cut the integration interval prematurely.For this purpose we calculate Introduction

Conclusions References
Tables Figures

Back Close
Full a maximum response time of the system as (Lasaga, 1980) where λ i are non-zero eigenvalues of the matrix A. The upper limit of integration in Eq. ( 16) is replaced by τ cycle in or calculations.
In future versions of SoilR, it will be possible to compute a dynamic, time-dependent transit-time distribution for objects of class Model with a time argument specifying for which time the distribution is sought.

Implementation of the general radiocarbon model
The implementation of the general model of radiocarbon is similar to the implementation of the general decomposition model presented in version 1.0 of SoilR (Sierra et al., 2012a).The system of ordinary differential equations is solved using the deSolve package of Soetaert et al. (2010).
In this new version, we introduced a new set of R classes to distinguish between the time-dependent (Eq. 1) and time-invariant (Eq.14) versions of our general models.In particular, we use the virtual super class DecompOp for different types of decomposition operators, and the virtual super class InFlux for different types of input fluxes.For radiocarbon related objects, we use the classes ConstFc and BoundFc to represent the radiocarbon fractions of time-invariant and time-bounded vectors, respectively.These classes must include an argument about the format of the radiocarbon values, either Delta14C or AbsoluteFractionModern.

Model initialization
All models that include radiocarbon dynamics are initialized in SoilR by the function GeneralModel_14().The arguments for this function are t: a vector containing the points in time where the solution is sought.

Conclusions References
Tables Figures

Back Close
Full -ivList: a vector containing the initial amount of carbon for the m pools.
-initialValF an object of class ConstFc containing a vector with the initial values of the radiocarbon fraction for each pool and a format string describing in which format the values are given (Delta14C or AbsoluteFractionModern).
-inputFluxes: an object of class InFlux consisting of a vector valued function describing the inputs to the pools.
-inputFc: an object of class BoundFc consisting of a function describing the fraction of 14 C in per mille of the input fluxes.
lambda: a scalar with the radiocarbon decay constant.By default, we use −0.0001209681 yr −1 .
solverfunc: the function used to solve the ODE system.This can be SoilR.euleror deSolve.lsoda.wrapperor any other user provided function with the same interface.
pass: if set to TRUE it forces the constructor to create the model even if it violates mass balance principles.By default, it is set ot FALSE.
Once a model of class Model14 has been initialized, it can be queried with one of the functions described in Table 1.The model can also be queried by the functions getC, getReleaseFlux, and getAccumulatedReleaseFlux.
For models with constant coefficients, the mean transit time can be calculated with the function getMeanTransitTime() applied to an object of class Introduction

Conclusions References
Tables Figures

Back Close
Full

Radiocarbon datasets
We introduced five new datasets in SoilR to facilitate the representation and analysis of soil radiocarbon dynamics.These datasets contain information on the atmospheric radiocarbon concentration over time for different spatial and temporal domains.For the pre-bomb period, IntCal09 (Reimer et al., 2009) and IntCal13 (Reimer et al., 2013)  For the post-bomb period (after 1950 AD) two additional datasets were included.The dataset C14Atm_NH was assembled for the Northern Hemisphere using data provided by Levin et al. (2010) and other measurements from North America.This dataset contains the atmospheric radiocarbon concentration in ∆ 14 C for 111 years, form 1900 to 2010 AD.
We also included the dataset compiled by Hua et al. (2013) for four different zones in the northern and Southern Hemispheres (Table S3 therein).This dataset, Hua2013 in SoilR, was implemented as an R list containing 5 data.frame,each representing an atmospheric zone with 5 variables.The variables are: the year AD, mean ∆ 14 C value, its standard deviation, mean F14 value, and its standard deviation.
We also included a dataset of observations of the ∆ 14 C value of respired CO 2 from soils of the Harvard Forest, MA, USA (Sierra et al., 2012b).This dataset, HarvardForest14CO2, was implemented as a data.framewith the variables: year of observation, ∆ 14 C value of respired CO 2 , and the site of measurement within the Harvard Forest.

GMDD Figures Back Close
Full

Auxiliary functions
A few functions were also introduced in this version of SoilR to help with processing of radiocarbon data.These are: bind.C14curves: binds pre-and a post-bomb ∆ 14 C curves together.The result can be expressed in years BP or AD.
-turnoverFit: finds the turnover times of a soil sample using the ∆ 14 C value measured at a particular year, the amount of litter inputs to soil, and an initial amount of C.
-PlotC14Pool: plots the output from a call to getF14 along with a radiocarbon curve.
For more details see the documentation of each function.

Model structure and transit times
To interpret radiocarbon observations in soil organic matter, it is common to use models with two or three pools that capture different cycling rates of carbon (O' Brien and Stout, 1978;Jenkinson and Rayner, 1977;Bruun et al., 2004;Gaudinski et al., 2000;Trumbore, 2000).However, a multi-pool model may have different connections among pools representing processes related to the stabilization and destabilization of organic 3173 Introduction

Conclusions References
Tables Figures

Back Close
Full matter (Sierra et al., 2011).In this example, we show how the connections among the pools may yield very different outcomes for interpreting soil radiocarbon data.
We will look at three different model structures of a three-pool model (Fig. 1), which are special cases of the general model of Eqs. ( 1) and (4).In this example we will ignore external environmental effects on decomposition rates, therefore we assume In the first case, carbon enters the soil and it is split among the three pools in different proportions (γ i ).Decomposition occurs in each pool independently without any transfer of carbon to other compartments.We call this model three-pool parallel, and can be written as In the second case, carbon enters only one of the reservoirs and it is transferred to other reservoirs in a cascade or series structure in which the residues of decomposition from one compartment may transfer to other compartments with lower decomposition rates (Swift et al., 1979;Manzoni and Porporato, 2009;Manzoni et al., 2009).This three-pool series model can be expressed mathematically as The third model structure considers a return of carbon residues to pools that decompose faster, mimicking processes of carbon destabilization from slowly cycling pools (Manzoni et al., 2009).Mathematically, the model can be expressed as To model radiocarbon dynamics under these three different assumptions of model structure, we transform C(t) in Eqs. ( 18), ( 19), and (20) to F (t) • C(t) and add a radiodecay term similarly as in the general models of Eqs. ( 1) and (4).In SoilR, these models are implemented by the functions ThreepParallel-Model14, ThreepSeriesModel14, and ThreepFeedbackModel14.We can run simulations for the period between the years 1901 and 2009 incorporating the atmospheric radiocarbon record of the Northern Hemisphere in the provided dataset C14Atm_NH.Using some arbitrary initial conditions and similar decomposition rates for all model structures (Table 2), we can observe differences between the radiocarbon content of the different pools as well as the radiocarbon content in the bulk soil and the respired CO 2 (Fig. 2).
Code to run these simulation is provided in the example of the function ThreepFeedbackModel14 of SoilR.To see the example simply type ?ThreepFeedbackModel14 in the R command shell.To run the example type example("ThreepFeedbackModel14").
The simulations show that even with the same amount of inputs and decomposition rates for the three pools, the temporal behavior of radiocarbon may change significantly (Fig. 2) posing challenges for the interpretation of measured data.
Furthermore, the mean transit times of carbon obtained from these three different model structures differ significantly among them.For the parallel model structure the mean residence time is 21 years, for the series model structure 29 years, and for the feedback model structure 79 years.The higher the complexity of the model (number of connections among pools), the longer carbon stays in the system (Bruun et al., 2004;Manzoni et al., 2009), which has a direct effect on the radiocarbon signature of the different pools, the bulk soil, and the respired CO 2 (Fig. 2).Introduction

Conclusions References
Tables Figures

Back Close
Full

Inverse parameter estimation: fitting a one pool model to a radiocarbon sample
Soil radiocarbon data is commonly used to estimate the turnover time (τ = 1/k) of a one-pool model.However, this is generally an ill-defined parameter estimation problem because the objective is to estimate the value of one parameter from one radiocarbon value.The problem gets exacerbated by the fact that there are always two possible solutions given the nature of the bomb-radiocarbon curve.
We introduced a function to estimate the two possible values of turnover time that can be obtained from one radiocarbon sample.This function, turnoverFit, takes as arguments the ∆ 14 C value of the soil sample and the year of measurement, the annual amount of litter inputs to soil either as a constant value or as a data.frame of inputs by year.It also requires an initial amount of carbon for the first year of the simulation, and a radiocarbon hemispheric zone according to Hua et al. (2013).
The function runs an optimization algorithm that minimizes the squared difference between the observation and the output of OnepModel14.It returns the two possible values of turnover time (τ = 1/k) that minimizes this difference between predictions and observations and a plot that illustrates the problem (Fig. 3).An example on how to run this function for a radiocarbon sample taken at a temperate forest soil is presented below.turnoverFit(obsC14=115.22, obsyr=2004.5, C0=2800, yr0=1900, In=473, Zone="NHZone2") The function runs much faster if not plot is produced, i.e. with the argument plot=FALSE.
One important limitation of this algorithm is the lack of uncertainty estimation for the predicted turnover times.We do not recommend this function for formal scientific analyses and reporting, but rather for preliminary exploration of laboratory results.A formal

GMDD Figures Back Close
Full estimation of turnover times can be achieved by performing inverse parameter estimation, which is described in the following example.

Inverse parameter estimation: Harvard Forest example
The assumption that soil organic carbon can be represented as a single, homogeneous pool is generally not supported by theory and observations of soil organic matter cycling (Swift et al., 1979;Bosatta and Agren, 1991;Trumbore, 2009;Manzoni and Porporato, 2009;Sierra et al., 2011), therefore the use of turnoverFit is not recommended for heterogenous organic matter.To account for this heterogeneity, it is necessary to use multi-pool models such as those in Fig. 2 or even more complex models with more pools and connections among them (e.g.O'Brien and Stout, 1978;Jenkinson and Rayner, 1977;Bruun et al., 2004;Gaudinski et al., 2000;Trumbore, 2000;Braakhekke et al., 2014).Parameters for these models can be objectively obtained using inverse parameter estimation (Schädel et al., 2013;Ahrens et al., 2014;Braakhekke et al., 2014).SoilR can be coupled with R package FME (Soetaert and Petzoldt, 2010) to obtain parameter values for a specific model.We will present an example on how to integrate both packages and use Markov chain Monte Carlo to obtain parameter values for a simple model of soil organic matter dynamics derived from measured radiocarbon data from the Harvard Forest, USA.Radiocarbon measurements of respired CO 2 have been collected at this site for the past decade as well as data on soil carbon stocks and proportions of organic matter in different fractions (Gaudinski et al., 2000;Sierra et al., 2012b).These radiocarbon data are provided in SoilR as HarvardForest14CO2.In a previous study, we found that a six-pool model can reproduce very well the observed patterns of soil radiocarbon over time (Sierra et al., 2012b).However, we are interested here in finding whether a simpler three-pool model containing roots, organic, and mineral carbon can reproduce the temporal behavior observed over time.This three pool model is expressed Introduction

Conclusions References
Tables Figures

Back Close
Full To implement this model in SoilR it is necessary to provide the arguments described in Sect.2.4.1 to the function GeneralModel_14.The code for this implementation is presented in the supplementary material as well as the code for creating a cost function using package FME with the function modCost, and fitting a preliminary model to data using the function mofFit.The mean squared residuals and the covariance matrix of the estimated parameters from this optimization are used to run a Markov chain Monte Carlo estimation procedure using the function modMCMC.
The results from this inverse parameter estimation procedure show that the model agrees well with the observed data (Fig. 4).Similarly, the distribution of the parameters seem to indicate unimodal posterior distributions of the parameters and some degree of correlation among them (Fig. 5).

Extrapolation of the atmospheric radiocarbon time series
Atmospheric radiocarbon data are only released at irregular intervals to the scientific community (e.g.Levin et al., 2010;Hua et al., 2013).For forward modeling of soil radiocarbon it is sometimes necessary to extrapolate existing data for some time into the future.There are a large number of tools in R for time series analyses and forecasting.For our specific problem, the forecast package (Hyndman and Khandakar, 2008) offers a simple and powerful extrapolation routine.
The function ets in package forecast automatically finds the best possible model for the given time series using exponential smoothing state-space modeling.Based on the fitted model, the function forecast produces predictions forward for a given number of periods for forecasting.

Conclusions References
Tables Figures

Back Close
Full Applying this procedure to the northern-hemisphere zone 1 series in Hua et al. (2013), we can forecast for example the concentration of radiocarbon in the atmosphere from 2010 to 2020 for this region (Fig. 6).The results from this forecast can be subsequently merged with the original dataset and run simulations using SoilR as described before.However, care must be taken with the interpretation of results using forecasted atmospheric radiocarbon data.

Conclusions
We introduced a number of functions and datasets within SoilR to model radiocarbon dynamics in soil organic matter.With this tool it is possible to model the temporal dynamics of radiocarbon in soils and respired CO 2 using models with any number of pools and connections among them.These models are generalizable to other systems where the incorporation of bomb radiocarbon is used to infer turnover or transit timesincluding human tissues, plants, sediments, etc. Radiocarbon data and other auxiliary information can also be used for model identification; i.e. to obtain parameter values of decomposition and transfer rates in models of soil organic matter decomposition.This is accomplished in SoilR with an interface to R package FME, but other inverse parameter estimation methods could also be used.
Depth profiles of radiocarbon cannot be simulated with this current implementation, but this dimension will be added in a future version of SoilR.Introduction

Conclusions References
Tables Figures

Back Close
Full  Full  Full   .It also requires an initial amount of carbon for the first do not recommend this function for formal scientific analyses       et al., 2013), including prediction intervals, for the period 2010-2020 using the forecast package (Hyndman and Khandakar, 2008).
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | The decomposition operator A(t), a square matrix of dimension m × m, contains in its main diagonal the decomposition rates k i for each pool i , and coefficients representing the proportion of carbon transferred from one pool to another in the off-diagonals.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | -A: a DecompOp object consisting of a matrix valued function describing the whole model decay rates for the m pools, connection and feedback coefficients as functions of time, and a time range for which this function is valid.The dimensions of this matrix must be equal to the number of pools.The time range must cover the times given in the t argument.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Fig. 1.Possible structures for a three-pool model.Each box represent a pool with a specific decomposition rate, and arrows represent inputs to or outputs from the pools.In the first case, carbon enters the system and it is split among the three pools in different propor-

Fig. 1 .
Fig.1.Possible structures for a three-pool model.Each box represent a pool with a specific decomposition rate, and arrows represent inputs to or outputs from the pools.In the first case, carbon enters the system and it is split among the three pools in different proportions without any transfer between pools.In the second case, carbon enters the system through one reservoirs and it is transferred serially between compartments.In the third case, carbon is returned back to donor pools.

Fig. 3 .Fig. 4 .
Fig. 3. Output of the function turnoverFit for a radiocarbon sample taken at a temperate forest soil subject to annual inputs of 473 Mg C ha −1 yr −1 .The upper panel shows the two possible curves that can match the observed radiocarbon value.The bottom curve shows the squared residuals between predictions and observations for different values of k in a one pool model.See documentation of function turnoverFit for additional details.

Fig. 4 .Fig. 5 .Fig. 6 .
Fig. 4. Predictions of respired radiocarbon values from the model of Eq. (21) vs. observations.Model predictions include uncertainty range for the mean ± standard deviation, and the minimum-maximum range.Radiocarbon concentration in the atmosphere is depicted in blue.

Fig. 6 .
Fig.6.Forecast of the atmospheric radiocarbon data of the northerhemisphere zone 1(Hua et al., 2013), including prediction intervals, for the period 2010-2020 using the forecast package (Hyndman
provide global-scale atmospheric radiocarbon data on an annual time-scale for the period 0-50 000 years BP.In SoilR, these datasets are called IntCal09 and IntCal13, respectively.They are implemented as data.framewith 5 variables: calibrated age in years BP, 14 C age in years BP, ∆ 14 C value in per mil, and corresponding uncertainty values for each.For additional details, see ?IntCal09 and ?IntCal13 in SoilR.

Table 1 .
Main functions implemented in SoilR version 1.1 to calculate the radiocarbon fraction in soil organic matter.Calculates the radiocarbon fraction for each pool at each time step.It returns a matrix of dimension n × m; i.e. n time steps as rows and m pools as columns.Calculates the average radiocarbon fraction weighted by the amount of carbon release at each time step.It returns a vector of length n with the value of F R for each time step.Introduction

Table 2 .
Parameter values and initial conditions used to simulate a three-pool model with different structures as represented in Fig.1and Eqs.(18)-(20).

Table 2 .
Parameter values and in three-pool model with different s 1 and equations18, 19, and 20.
the example simply type ?T in the R command shell.example(''ThreepFeedThesimulations show that 440 inputs and decomposition rate