Full documentation and tutorials can be found on the isobxr website.
The isobxr package is a set of R tools designed to perform and explore stable isotope box modelling of open or closed systems. It provides a ready-to-use tool allowing users to develop and test isotopic box models of their system of interest. It allows the user to explore the behavior of these systems in both static (e.g., at steady state) or dynamic modes (e.g., in reaction to a perturbation), build complex scenarios, and sweep the space of parameters in both static and dynamic modes.
The global definition of the systems considered for the stable isotope box modelling by isobxr is presented in the cheatsheet below.
run_isobxr is the main function of the isobxr package.
run_isobxr is used to perform single runs of stable isotope box models of all types of systems.
The general workflow of the run_isobxr function is summarized in the cheatsheet below.
In this section we present the mathematical formalism used in isobxr stable isotope box modelling.
The theory and formalism shown in this section is derived from the chapter 7 (Dynamic systems) of
: Albarède, F., 1995. Introduction to Geochemical Modelling. Cambridge University Press, Cambridge.
https://doi.org/10.1017/CBO9780511622960
We consider a system consisting of \(n\) boxes \(j\), exchanging an element X that has at least two stable isotopes \(M\) and \(N\). Any box \(i\) is defined by:
We call \(R_i\) the \(M/N\) isotope abundance ratio: \[R_i = \left(\dfrac{n^M}{n^N}\right)_i = \left(\dfrac{m^M}{m^N}\right)_i \times \left(\dfrac{M^N}{M^M}\right) = \left(\dfrac{m^M}{m^N}\right)_i A \] where
Finally, we define \(D_{i \mapsto j}^{M/N}\) the relative fractionation coefficient between M and N from \(i\) to \(j\):
\[D_{i \mapsto j}^{M/N} = \dfrac{{J_{i \mapsto j}^M}/{J_{i \mapsto j}^N}}{m^M_i/m^N_i}\]
The principle of modelling the system consists now in predicting the evolution of box sizes (\(m^X_i\)) during time and the evolution of the isotope ratios (\(R_i\)), knowing the fluxes, initial box sizes and fractionation coefficients.
Considering X has no radioactive nor radiogenic component, mass conservation for element X reads :
\[\dfrac{dm^X_i}{dt} = - \sum \limits_{j\neq i} J_{i \mapsto j}^X + \sum \limits_{j\neq i} J_{j \mapsto i}^X\]
When considering constant fluxes, this equation can be integrated in its linear form as follows:
\[m^X_i = \left(- \sum \limits_{j\neq i} J_{i \mapsto j}^X + \sum \limits_{j\neq i} J_{j \mapsto i}^X\right) \times t + m^X_{i, t_0}\]
where \(m^X_{i, t_0}\) is the mass of element X in box \(i\) at t = 0.
Considering \(N\) and \(M\) are stable and non radiogenic, mass conservation reads :
\[\dfrac{dm^N_i}{dt} = - \sum \limits_{j\neq i} J_{i \mapsto j}^N + \sum \limits_{j\neq i} J_{j \mapsto i}^N = - \sum \limits_{j\neq i} \dfrac{J_{i \mapsto j}^N}{m^N_i} m^N_i + \sum \limits_{j\neq i} \dfrac{J_{j \mapsto i}^N}{m^N_j} m^N_j\]
Likewise,
\[\dfrac{dm^M_i}{dt} = - \sum \limits_{j\neq i} \dfrac{J_{i \mapsto j}^M}{m^M_i} m^M_i + \sum \limits_{j\neq i} \dfrac{J_{j \mapsto i}^M}{m^M_j} m^M_j\]
We can meanwhile derive the ratio \(R_i\) as :
\[\dfrac{dR_i}{dt} = \dfrac{d}{dt}\left(\dfrac{m^M}{m^N}\right)_i A = \dfrac{A}{{m^N_i}^2} \left(m^N_i\dfrac{dm^M_i}{dt} - m^M_i \dfrac{dm^N_i}{dt}\right) = \dfrac{A}{m^N_i} \left(\dfrac{dm^M_i}{dt} - \dfrac{R_i}{A} \dfrac{dm^N_i}{dt}\right)\]
We can insert mass conservation derivative into previous relationship:
\[\dfrac{dR_i}{dt} = \dfrac{A}{m^N_i} \left(- \sum \limits_{j\neq i} \dfrac{J_{i \mapsto j}^M}{m^M_i} m^M_i + \sum \limits_{j\neq i} \dfrac{J_{j \mapsto i}^M}{m^M_j} m^M_j + \dfrac{R_i}{A} \sum \limits_{j\neq i} \dfrac{J_{i \mapsto j}^N}{m^N_i} m^N_i - \dfrac{R_i}{A} \sum \limits_{j\neq i} \dfrac{J_{j \mapsto i}^N}{m^N_j} m^N_j \right)\]
Modified as such :
\[\dfrac{dR_i}{dt} = A \left(- \sum \limits_{j\neq i} \dfrac{J_{i \mapsto j}^M}{m^M_i} \dfrac{m^M_i}{m^N_i} + \sum \limits_{j\neq i} \dfrac{J_{j \mapsto i}^M}{m^M_j} \dfrac{m^M_j}{m^N_j} \dfrac{m^N_j}{m^N_i} + \dfrac{R_i}{A} \sum \limits_{j\neq i} \dfrac{J_{i \mapsto j}^N}{m^N_i} - \dfrac{R_i}{A} \sum \limits_{j\neq i} \dfrac{J_{j \mapsto i}^N}{m^N_j} \dfrac{m^N_j}{m^N_i} \right)\]
\[\dfrac{dR_i}{dt} = A \left(- \sum \limits_{j\neq i} \dfrac{J_{i \mapsto j}^M}{m^M_i} \dfrac{R_i}{A} + \sum \limits_{j\neq i} \dfrac{J_{j \mapsto i}^M}{m^M_j} \dfrac{R_j}{A} \dfrac{m^N_j}{m^N_i} + \dfrac{R_i}{A} \sum \limits_{j\neq i} \dfrac{J_{i \mapsto j}^N}{m^N_i} - \dfrac{R_i}{A} \sum \limits_{j\neq i} \dfrac{J_{j \mapsto i}^N}{m^N_j} \dfrac{m^N_j}{m^N_i} \right)\]
\[\dfrac{dR_i}{dt} = R_i \sum \limits_{j \neq i} \left(\dfrac{J^N_{i \mapsto j}}{m^N_i} - \dfrac{J^M_{i \mapsto j}}{m^M_i} \right) + \sum \limits_{j \neq i} \dfrac{J^M_{j \mapsto i}}{m^M_j} \dfrac{m^N_j}{m^N_i} R_j - \sum \limits_{j \neq i} \dfrac{J^N_{j \mapsto i}}{m^N_j} \dfrac{m^N_j}{m^N_i} R_i\]
\[\dfrac{dR_i}{dt} = R_i \sum \limits_{j \neq i} \left(\dfrac{J^N_{i \mapsto j}}{m^N_i} - \dfrac{J^M_{i \mapsto j}}{m^M_i} \right) + \sum \limits_{j \neq i} \left( \dfrac{J^M_{j \mapsto i}}{m^M_j} R_j - \dfrac{J^N_{j \mapsto i}}{m^N_j} R_i \right) \dfrac{m^N_j}{m^N_i}\]
\[\dfrac{dR_i}{dt} = R_i \sum \limits_{j \neq i} \dfrac{J_{i \mapsto j}^N}{m^N_i} \left(1- \dfrac{{J_{i \mapsto j}^M}/{J_{i \mapsto j}^N}}{{m^M_i}/{m^N_i}}\right) + \sum \limits_{j \neq i} \dfrac{J_{j \mapsto i}^N}{m^N_j} \left(\dfrac{{J_{j \mapsto i}^M}/{J_{j \mapsto i}^N}}{m^M_j / m^N_j} R_j - R_i \right) \dfrac{m^N_j}{m^N_i}\]
Resulting in :
\[\dfrac{dR_i}{dt} = R_i \sum \limits_{j \neq i} \dfrac{J_{i \mapsto j}^N}{m^N_i} \left(1- D_{i \mapsto j}^{M/N} \right) + \sum \limits_{j \neq i} \dfrac{J_{j \mapsto i}^N}{m^N_j} \left(D_{j \mapsto i}^{M/N} R_j - R_i \right) \dfrac{m^N_j}{m^N_i}\]
In order to solve this equation, we assume that:
\[\begin{array}{ccccc} \dfrac{J_{i \mapsto j}^N}{m^N_i} \simeq \dfrac{J_{i \mapsto j}^X}{m^X_i} & , & \dfrac{J_{j \mapsto i}^N}{m^N_j} \simeq \dfrac{J_{j \mapsto i}^X}{m^X_j} & , & \dfrac{m^N_j}{m^N_i} \simeq \dfrac{m^X_j}{m^X_i} \end{array}\]
These assumptions are valid for amplitude of stable isotope fractionations of the order of tens of ‰.
These assumptions find their limitation in the case of ample stable isotope fractionations such as sometimes observed for isotopic systems of light elements (such as D/H). In these cases, the differences in reactivity can be too great and amplitudes of stable isotope fractionations too marked for the assumptions to be applicable.
We thus have the following system of first order differential equations used thereafter:
\[ \dfrac{dR_i}{dt} = R_i \sum \limits_{j \neq i} \dfrac{J_{i \mapsto j}^X}{m^X_i} \left(1- D_{i \mapsto j}^{M/N} \right) + \sum \limits_{j \neq i} \dfrac{J_{j \mapsto i}^X}{m^X_j} \left(D_{j \mapsto i}^{M/N} R_j - R_i \right) \dfrac{m^X_j}{m^X_i} \]
This system can be solved in order to reconstruct the evolution of stable isotope ratios \(R_i(t)\) in each box over time.
Two main situations can be distinguished here.
The inward and outward fluxes of element X for any given box are balanced
In this case, all box sizes (\(m^X_i\)) are constant, which means that \({dm^X_i}/{dt} = 0\) for all boxes.
The reconstruction of \(R_i(t)\) only requires the solving of the system of ordinary first order differential equations.
This can be done analytically via the algebraic inversion of the equation system matrix. This is what the ana_slvr core functions does.
The inward and outward fluxes of element X for at least one box is unbalanced
In this case, at least one box displays a varying mass of element X with time (\(m^X_i(t)\)) which means that \({dm^X_i}/{dt} \neq 0\) for at least one box.
The reconstruction of \(R_i(t)\) requires the simultaneous solving of the system of ordinary first order differential equations describing \(R_i(t)\) and of the derivative of \(m^X_i(t)\) (with a linear solution).
This is performed numerically via an incremental integration as done by the num_slvr core function.
All \(R_i(t)\) stable isotope ratios can be also expressed as \(\delta_i(t)\) values expressed in the ‰ unit, defined as follows and where \(R_{std}\) is the stable isotope ratio of the standard or reference material:
\[\delta^{M/N}X_i(t) = \left(\dfrac{R_i(t)}{R_{std}} - 1\right) \times 1000 \]
The isobxr package uses two core solver functions to solve the system of differential equations: num_slvr and ana_slvr.
These functions are not intended to be directly called by users themselves, although this remains possible.
These functions are called within broader functions (e.g., run_isobxr).
Because these functions are the actual solvers of the stable isotope box model equation systems, you will find here a description of their functioning and applicability.
You can also refer to the documentation for these two functions by typing the following code into your R console:
?ana_slvr
?num_slvrThe INPUT file contains the set of commands required to run either ana_slvr or num_slvr.
Its format and name are standardized and should not be altered.
Its direct modification by user is not required because INPUT files are edited by the broader and user-friendly run_isobxr function, described in another section of this vignette.
The INPUT file contains at least the 4 following objects strictly named as follows:
The INPUT file is an output from the run_isobxr function and does not require the user intervention.
The numerical solver num_slvr iteratively integrates the derivatives of the box sizes (\(\frac{dm^X_i}{dt}\)) and isotope ratios (\(\frac{dR_i}{dt}\)) over the whole run duration.
The num_slvr uses the ode ordinary differential equation integrating function from the deSolve package.
This function solves the evolution of both the isotope ratios and the sizes of each box defined by user.
The num_slvr function therefore can be used in two distinct cases:
In the case of balanced inward and outward fluxes of element X between all involved boxes.
The derivative of box sizes will be null (\(\frac{dm^X_i}{dt} = 0\)), in other words all box sizes will be constant.
In the case of unbalanced inward and outward fluxes of element X for at least one box of the system.
The derivative of box sizes will be non null in at least one case (\(\frac{dm^X_i}{dt} \neq 0\)), in other words at least on box loses and accumulates the element X.
The drawbacks of this function is that
The use of this function is thus preferred in the case of unbalanced fluxes, where at least one box loses or accumulates the element X over time (case #2).
It then allows to study the dynamic responses of a system to a perturbation of the flux balance for at least one box.
Remark: The direct use of num_slvr function by user is not intended because this function is integrated in the global run_isobxr function described thereafter.
The num_slvr function reads an INPUT file (file name structure: [ RUN name + INPUT.Rda ]) containing all commands for the required model run, namely:
The num_slvr function returns the numerically determined evolution of the stable isotope ratios and masses of element X in all boxes over the run duration.
The num_slvr OUTPUT contains the following:
By default (save_run_outputs = FALSE), all OUTPUT data is stored in a temporary directory. If save_run_outputs = TRUE, the outputs are saved in the same directory as the input file.
By default, all OUTPUT data is saved in a Rda output file (RUN name + OUT.Rda).
The num_slvr function optionally allows user to edit the OUTPUT data as csv files stored in a DIGEST folder ([ RUN name + DIGEST ]) with the following name structures:
The analytical solver ana_slvr calculates the analytical solutions of the system of ordinary differential equations (ODES) of stable isotope ratios for an element X in all boxes.
It determines the eigenvalues, eigenvectors and constants defined by initial conditions associated to the ODES, and allows to determine the functions that describe the evolution of stable isotope ratios as a function of time.
In this case, the box sizes are considered as constant because the solving of the system of differential equations only considers the isotope ratios as varying with time.
In other words, in the case of balanced inward and outward fluxes for all finite boxes (\(\frac{dm^X_i}{dt} = 0\)), this function allows to accurately determine the expression of the stable isotope ratio in a given box as a function of time.
The function also edits a dataset of the evolution of isotope compositions in each box over the time range, with the time resolution determined by the user.
The use of this function is thus not applicable in the case of unbalanced fluxes for one finite box or more.
However, the ana_slvr function can be used in balanced open systems, with boxes defined as infinite and unbalanced sources or sinks, as described in the model design section thereafter.
Remark: The direct use of ana_slvr function by user is not intended because this function is integrated in the global run_isobxr function described thereafter.
The ana_slvr function reads an INPUT file (file name structure: [ RUN name + INPUT.Rda ]) containing all commands for the required model run, namely:
The ana_slvr function returns the analytical solutions of the ODES which describes the stable isotope ratios in each box. The ana_slvr function also returns the evolution of the delta values in each box. (in ‰, relative to the reference material).
The ana_slvr OUTPUT contains the following:
ODE_SOLNs data file summarizing outputs of the analytical solutions of the ODES
(eigenvalues, eigenvectors, relaxation times, constants according to initial conditions)
evD data file storing the evolution of the delta values in each box
By default (save_run_outputs = FALSE), all OUTPUT data is stored in a temporary directory. If save_run_outputs = TRUE, the outputs are saved in the same directory as the input file.
By default, all OUTPUT data is saved in a Rda output file (RUN name + OUT.Rda).
The ana_slvr function optionally allows user to edit the OUTPUT data as csv files stored in a DIGEST folder ([ RUN name + DIGEST ]) with the following name structures: