Technical note: A Bayesian mixing model to unravel isotopic data and quantify trace gas production and consumption pathways for time series data – Time-resolved FRactionation And Mixing Evaluation (TimeFRAME)

Isotopic measurements of trace gases such as N₂O, CO₂, and CH₄ contain valuable information about production and consumption pathways. Quantification of the underlying pathways contributing to variability in isotopic time series can provide answers to key scientific questions, such as the contribution of nitrification and denitrification to N₂O emissions under different environmental conditions or the drivers of multiyear variability in atmospheric CH₄ growth rate. However, there is currently no data analysis package available to solve isotopic production, mixing, and consumption problems for time series data in a unified manner while accounting for uncertainty in measurements and model parameters as well as temporal autocorrelation between data points and underlying mechanisms. Bayesian hierarchical models combine the use of expert information with measured data and a mathematical mixing model while considering and updating the uncertainties involved, and they are an ideal basis to approach this problem.

Here we present the Time-resolved FRactionation And Mixing Evaluation (TimeFRAME) data analysis package. We use four different classes of Bayesian hierarchical models to solve production, mixing, and consumption contributions using multi-isotope time series measurements: (i) independent time step models, (ii) Gaussian process priors on measurements, (iii) Dirichlet–Gaussian process priors, and (iv) generalized linear models with spline bases. We show extensive testing of the four models for the case of N₂O production and consumption in different variations. Incorporation of temporal information in approaches (i)–(iv) reduced uncertainty and noise compared to the independent model (i). Dirichlet–Gaussian process prior models have been found to be most reliable, allowing for simultaneous estimation of hyperparameters via Bayesian hierarchical modeling. Generalized linear models with spline bases seem promising as well, especially for fractionation estimation, although the robustness to real datasets is difficult to assess given their high flexibility. Experiments with simulated data for δ¹⁵N^bulk and δ¹⁵N^SP of N₂O showed that model performance across all classes could be greatly improved by reducing uncertainty in model input data – particularly isotopic end-members and fractionation factors. The addition of the δ¹⁸O additional isotopic dimension yielded a comparatively small benefit for N₂O production pathways but improved quantification of the fraction of N₂O consumed; however, the addition of isotopic dimensions orthogonal to existing information could strongly improve results, for example, clumped isotopes.

The TimeFRAME package can be used to evaluate both static and time series datasets, with flexible choice of the number and type of isotopic end-members and the model setup allowing simple implementation for different trace gases. The package is available in R and is implemented using Stan for parameter estimation, in addition to supplementary functions re-implementing some of the surveyed isotope analysis techniques.

1 Introduction

Analysis of isotopic signatures is frequently used in environmental sciences to infer production and consumption pathways for trace gases. For example, N₂O isotopic composition reflects the production via different pathways (including microbial denitrification, nitrification, fungal denitrification), mixing within the soil airspace, and consumption via complete denitrification. The different production pathways have distinct isotopic end-members, which describe the isotopic composition of emitted N₂O. Following emission, N₂O from different pathways mixes in the soil airspace, described using the approximated mixing equation (Ostrom et al., 2007; Fischer, 2023):

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}_{\text{mix}}\approx \sum _{k=\mathrm{1}}^K{\mathit{f}}_k{\mathit{\delta }}_k,\end{array}$$

where δ_mix is the isotopic composition of a mixture of two or more sources enumerated by $k=\mathrm{1},\mathrm{\dots },K$ with isotopic compositions designated by δ_k and fractional contributions to the mixture designated by f_k. This approximation assumes that the light isotope has a much greater concentration than the heavy isotope, which is valid for common trace gases such as CO₂, CH₄, and N₂O.

N₂O is consumed during complete denitrification to N₂, which favors the light isotope and thus leads to progressive enrichment of the remaining N₂O pool. The isotopic effect of consumption can be approximated using the Rayleigh equation (Mariotti et al., 1981; Ostrom et al., 2007; Fischer, 2023):

$$\begin{array}{}\text{(2)}& {\mathit{\delta }}_{\text{substr},r}\approx {\mathit{\delta }}_{\text{substr},r=\mathrm{1}}+\mathit{ϵ}\ln \left(r\right),\end{array}$$

where ${\mathit{\delta }}_{\text{substr},r=\mathrm{1}}$ and δ_substr,r are the isotopic composition of the initial substrate prior to consumption (r=1) and when a certain fraction (1−r) has been consumed, ϵ is the fractionation factor for the reaction in per mill (‰), and r is the fraction of substrate remaining where r=0 represents a complete reaction.

We can combine Eqs. (1) and (2) for a full model of mixing and fractionation of the subsequent mixture – for example, mixing of N₂O from different sources within the soil airspace, followed by complete reduction of a certain fraction of N₂O, before measurement of N₂O isotopic composition (Fischer, 2023):

$$\begin{array}{}\text{(3)}& \mathit{\delta }=\sum _{k=\mathrm{1}}^K{\mathit{f}}_k{\mathit{\delta }}_k+\mathit{ϵ}\ln \left(r\right),\end{array}$$

where δ is the measured isotopic composition. In this equation, we assume that mixing occurs before fractionation, when in reality mixing and fractionation are likely occurring simultaneously depending on the soil pore size distribution and connectivity, the availability of different substrates, and the microbial community present (Denk et al., 2017; Yu et al., 2020; Lewicka-Szczebak et al., 2020). A detailed discussion of the implications of this assumption is given in Sect. 1 in the Supplement. Further uncertainty in the model equation relates to open- vs. closed-system fractionation, describing renewal of the N₂O pool relative to the rate of N₂O consumption (Yu et al., 2020; Lewicka-Szczebak et al., 2020). However, the largest uncertainties in evaluation of this equation to interpret the measured isotopic composition δ relate to the end-members for different sources δ_k and the fractionation factor ϵ.

A commonly used approach to interpret trace gas isotopic measurements is the application of dual-isotope mapping, which utilizes the relationship between two isotopic parameters to infer pathways, for example, δ¹⁵N^bulk and δ¹⁵N^SP in the case of N₂O. The mapping approach can be used to roughly estimate the dominance of different pathways and the importance of fractionation during consumption (Wolf et al., 2015; Lewicka-Szczebak et al., 2017; Wu et al., 2019; Yu et al., 2020; Rohe et al., 2021). However, these approaches fail to provide a quantitative determination of different pathways or to estimate uncertainty for individual samples. Moreover, mapping approaches are limited to mixing scenarios involving only two sources, which – for example – does not allow for the differentiation of contributions from the nitrification and fungal denitrification pathways which have similar δ¹⁵N^SP signatures. In addition, there are no statistical packages available to implement these mapping approaches, calling into question the reproducibility among studies using this approach.

Bayesian approaches to solve isotopic mixing models have been successfully implemented in several well-known frameworks (R packages MixSIAR and simmr; Parnell et al., 2013; Stock et al., 2018). These advanced models are used to resolve the contribution of multiple sources to a mixture using a range of Bayesian statistical techniques, and they are widely used for applications such as animal diet partitioning (Stock et al., 2018). However, these packages do not offer the capability to deal with pool consumption and Rayleigh fractionation, and thus they are not suitable for the interpretation of trace gas isotopic measurements where consumption/destruction plays a key role (Fischer, 2023; Lewicka-Szczebak et al., 2020).

The FRAME (Fractionation And Mixing Evaluation) model provided the first Bayesian tool to include both mixing and fractionation for the interpretation of isotopic data (Lewicka-Szczebak et al., 2020; Lewicki et al., 2022). FRAME applies a Markov chain Monte Carlo (MCMC) algorithm to estimate the contribution of individual sources and processes, as well as the probability distributions of the calculated results (Lewicki et al., 2022). However, the FRAME model can only be applied to time series data by solving separately for single or aggregated points. Although the contribution of different pathways may vary strongly on short timescales, model parameters – such as isotopic end-members and fractionation factors – are expected to show minimal variability between subsequent points in a time series. Time series information can be added to isotopic models (i) through statistical approaches using smoothing and other techniques to account for temporal autocorrelation and measurement noise or (ii) through the application of dynamic approaches incorporating differential equations (Bonnaffé et al., 2021). In the Time-resolved FRactionation And Mixing Evaluation (TimeFRAME), we use the statistical approach as a natural extension to the implementation of FRAME; investigation of dynamical approaches may be challenging due to high uncertainties in all inputs and should be a focus of further research.

Here we present the TimeFRAME extension to the FRAME model to allow for efficient analysis of time series data. TimeFRAME uses one independent time step model in which points in a time series are treated independently and three classes of models to fully incorporate time series information: (i) independent time step models, (ii) Gaussian process priors on measurements, (iii) Dirichlet–Gaussian process priors, and (iv) generalized linear models with spline bases. The models are solved for the contribution of different pathways, end-members, and fractionation factors within a MCMC framework, and the full posterior distributions of parameters are reported. The isotopes, end-members, fractionation factors, and model setup are defined by the user, allowing flexible application to many isotopic problems.

2 Methodology

2.1 Inference of source contributions

One objective of studying isotopic signatures is to determine the source contributions f₁⋯f_K from measurements of the mixture. However, measuring one single isotopic species will only be efficient in distinguishing between a maximum of two sources. For additional sources, or if consumption of the mixture needs to be accounted for, multiple isotopic species are necessary. Analysis of N₂O sources and pathways, for instance, can include analysis of δ¹⁵N^bulk, δ¹⁵N^SP, and δ¹⁸O. The vector of d different isotopic species shall be denoted by X∈ℝ_d. Measurements of the isotopic end-member for each individual source enumerated by $k=\mathrm{1},\mathrm{\dots },K$ are assumed to be known and denoted by $S_{\mathrm{1}},\mathrm{\dots },S_K\in {\mathbb{R}}_d$ together with the fractionation factor ϵ∈ℝ_d. Using vector and matrix notation they can later be used to state the mixing equation in vector form:

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }\mathbf{f}:=[{\mathit{f}}_{\mathrm{1}}\mathrm{\cdots }{\mathit{f}}_K{]}^T\in {\mathbb{R}}^K,\text{(5)}& {\displaystyle }& {\displaystyle }\mathbf{S}:=[S_{\mathrm{1}}\mathrm{\cdots }{S}_K{]}^T\in {\mathbb{R}}^{d\times K}.\end{array}$$

The case of Rayleigh fractionation as expressed in Eq. (3) can be similarly expressed in vectorized form:

$$\begin{array}{}\text{(6)}& X=\mathit{\mu }(\mathbf{f},r):=\mathbf{Sf}+\mathit{ϵ}\ln \left(r\right).\end{array}$$

In a simple example with two sources and measurements of two isotopic species $K=d=\mathrm{2}$, the mixing equation can be solved (assuming convergence is possible) for the parameters of interest using linear algebra (see Fischer, 2023, for details). More generally, mixing and fractionation according to Eq. (6) can be solved for an arbitrary number of sources as long as an equal number of isotopic species is available, i.e., $K=d\ge \mathrm{2}$. In this case the linear system of equations can be written in matrix terms and augmented with the sum constraint on f:

$$\begin{array}{}\text{(7)}& \stackrel{\mathrm{̃}}{X}:=\left[\begin{array}{c}X\\ \mathrm{1}\end{array}\right]=\left[\begin{array}{cc}\mathbf{S}& -\mathit{ϵ}\\ {\mathrm{1}}^{\mathrm{T}}& \mathrm{0}\end{array}\right]\left[\begin{array}{c}\mathbf{f}\\ \ln \left(r\right)\end{array}\right]=:\stackrel{\mathrm{̃}}{\mathbf{S}}\stackrel{\mathrm{̃}}{\mathbf{f}}.\end{array}$$

This d+1-dimensional linear system of equations can be addressed with decomposition techniques, and its solution can be expressed as $\stackrel{\mathrm{̃}}{\mathbf{f}}={\stackrel{\mathrm{̃}}{\mathbf{S}}}^{-\mathrm{1}}\stackrel{\mathrm{̃}}{X}$. A unique solution exists if $\stackrel{\mathrm{̃}}{\mathbf{S}}$ is invertible or, equivalently, if none of the mixing lines and the consumption line are co-linear. Only non-negative solutions $\stackrel{\mathrm{̃}}{\mathbf{f}}\ge \mathrm{0}$ are feasible to ensure that the source contributions f correspond to mixing weights and that $\mathrm{0}<r\le \mathrm{1}$.

A flaw of the isotope mapping approach as presented above is that it does not take measurement uncertainty into account. However, this can easily be added by formulating the measurements X as random variables with expected value given by the mixing equation $\mathbb{E}\left[X\right]=\mathit{\mu }(\mathbf{f},r)$. Most commonly, measurements are modeled using the Gaussian distribution with independent components and variance η²∈ℝ, thus allowing the mixing-fractionation equation to be stated as

$$\begin{array}{}\text{(8)}& X\sim {\mathcal{N}}_d\left(\mathit{\mu }\right(\mathbf{f},r),{\mathit{\eta }}^{\mathrm{2}}\mathbb{1})={\mathcal{N}}_d(\mathbf{Sf}+\mathit{ϵ}\ln (r),{\mathit{\eta }}^{\mathrm{2}}\mathbb{1}).\end{array}$$

A maximum-likelihood solution to this mixing-fractionation equation can be pursued (see Fischer, 2023, for details); however, this limits the framework to parameters that can be approximated with a Gaussian distribution. Often, the epistemic uncertainty of source isotopic end-members is modeled as a uniform distribution to best account for the combination of measurement uncertainty and true variability in end-member values (Lewicki et al., 2022). Bayesian statistics is useful to incorporate all assumptions and constraints into the model, as well as to employ numerical inference methods for source contribution estimation and uncertainty estimation.

2.2 Stationary inference

Stationary inference involves inference for the source contributions f∈𝒮_K, where 𝒮_K is the K-simplex, and the fraction of the pool remaining $r\in [\mathrm{0},\mathrm{1}]$ for one single measurement independent of time. This can be accomplished with the original FRAME model (Lewicki et al., 2022). FRAME constructs a prior and likelihood structure where the isotopic species measurements X∈ℝ^d are independently normally distributed with variance vector ${\mathit{\eta }}^{\mathrm{2}}\in {\mathbb{R}}_{+}^d$ around a mean given by an arbitrary mixing equation μ(f,r). The source contributions f are then equipped with a flat Dirichlet prior and the fraction remaining r with a uniform prior:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& \mathbf{f}\sim \text{Dir}\left(\mathrm{1}\right),r\sim \text{Uni}(\mathrm{0},\mathrm{1})\\ & X\mathrm{|}\mathbf{f},r\sim {\mathcal{N}}_d\left(\mathit{\mu }\right(\mathbf{f},r),{\mathit{\eta }}^{\mathrm{2}}).\end{array}\end{array}$$

The auxiliary parameters for the mixing equation S and ϵ are understood to be random variables as well, with predetermined fixed priors that are omitted from the model description above. Choosing those priors is dependent on the origin of the data and thus not subject to the further engineering of inference models in the following sections. The likelihood of X is understood to be implicitly conditional on these auxiliary parameters. This means that a joint posterior $\mathit{\pi }(\mathbf{f},r,\mathbf{S},\mathit{ϵ}\mathrm{|}X)$ is fit by the model, and the reported posterior π(f,r|X) is simply its marginalization.

The FRAME model can be extended by taking different choices of prior distributions for the parameters of interest (the source contributions f and the fraction remaining r). The Jeffreys prior for source contributions is constructed by computing the Fisher information matrix and choosing the probability distribution proportional to the square root of its determinant. For the source contributions $\mathbf{f}=\left(\begin{array}{c}\mathrm{1}-f\\ f\end{array}\right)\in {\mathcal{S}}_{\mathrm{2}}$ of two sources $S_{\mathrm{1}},S_{\mathrm{2}}\in \mathbb{R}$ the computation can be done by omitting the influence of fractionation, leading to a uniform prior over the domain of f. For multiple source contributions this is equivalent to the flat Dirichlet distribution used in the original FRAME model.

Taking now the fraction remaining with regards to Rayleigh fractionation $r\in [\mathrm{0},\mathrm{1}]$ independently of the mixing weights f into account, the Jeffreys prior can be computed relative to the pure mixing solution $M=X-S_{\mathrm{1}}(\mathrm{1}-f)-S_{\mathrm{2}}f\in \mathbb{R}$ with fractionation factor ϵ∈ℝ (Fischer, 2023):

$$\begin{array}{}\text{(10)}& {\mathcal{I}}_r\left(r\right)\propto -\mathbb{E}\left[{\displaystyle \frac{d^{\mathrm{2}}r}{d{r}^{\mathrm{2}}}}{\displaystyle \frac{(M-\mathit{ϵ}\ln (r){)}^{\mathrm{2}}}{{\mathit{\eta }}^{\mathrm{2}}}}|r\right]={\displaystyle \frac{\mathrm{2}{\mathit{ϵ}}^{\mathrm{2}}}{{\mathit{\eta }}^{\mathrm{2}}{r}^{\mathrm{2}}}}\propto {\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}.\end{array}$$

Therefore, the objective prior for the fraction remaining r is given by $\mathit{\pi }\left(r\right)\propto \frac{\mathrm{1}}{r}$ for $r\in [\mathrm{0},\mathrm{1}]$, which is also known as the logarithmic prior; since it cannot be normalized it is an improper prior. Additionally, even though this prior can be considered uninformative for r individually according to the Jeffreys criterion, a joint prior could lead to different results, although priors are typically chosen as independent distributions.

In the case of Rayleigh fractionation, however, it might be more reasonable to use a different distribution that is somewhere in between the uniform and logarithmic prior and incorporates the bounds to the interval [0,1] as well, which is a similar idea to prior averaging (Berger et al., 2015). The beta distribution offers a functional form that is similar to the Jeffreys prior but can be normalized. Parameterizing the distribution with a restricted concentration parameter $\mathit{\alpha }\in [\mathrm{0},\mathrm{1}]$, the form beta(α, 1) is the uniform distribution for α=1 and converges to the Jeffreys prior for α→0, thus expressing a generalized approach.

2.3 Time series inference

In order to incorporate time series information in the inference procedure, the model can be extended to work with multiple measurements at different points in time. The source contribution f and fraction remaining r are assumed to be functions with respect to time f_τ and r_τ, and the measurements correspond to samples in time X_t=X(τ_t) at n discrete time points ${\mathit{\tau }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\tau }}_n$.

Now the measurements can be grouped into a measurement matrix $\mathbf{X}:=\left[X_{\mathrm{1}}\mathrm{\dots }{X}_N\right]\in {\mathbb{R}}^{d\times N}$, where the time dimension is along the matrix rows. Inference of the parameters can be done at the identical time points f_t=f(τ_t) and r_t=r(τ_t), so that they can be grouped into similar matrices as well: $\mathbf{F}:=\left[{\mathbf{f}}_{\mathrm{1}}\mathrm{\dots }{\mathbf{f}}_N\right]\in {\mathbb{R}}^{K\times N}$ and $\mathbf{r}:=\left[r_{\mathrm{1}}\mathrm{\dots }{r}_N\right]\in {\mathbb{R}}^{\mathrm{1}\times N}$. This grouping has the advantage that the mixing equation can be expressed in vectorized form over all time points without changing its general layout (Fischer, 2023):

$$\begin{array}{}\text{(11)}& \mathbb{E}[\mathbf{X}\mathrm{|}\mathbf{F},\mathbf{r}]=\mathit{\mu }(\mathbf{F},\mathbf{r})=\mathbf{SF}+\mathit{ϵ}\ln \left(\mathbf{r}\right).\end{array}$$

2.3.1 Independent time steps

The simplest method to extend the stationary model is to assume complete independence between all points in time. This reduces the time series problem to a set of N independent stationary problems with one single measurement point each; thus, the same stationary FRAME model can be used for each point. The vector of measurement errors η∈ℝ^d is now also allowed to vary in time as ${\mathit{\eta }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\eta }}_N$:

$$\begin{array}{}\text{(12)}& \begin{array}{rlr}& {\mathbf{f}}_t\sim \text{Dir}\left(\mathrm{1}\right),r_t\sim \text{Uni}(\mathrm{0},\mathrm{1})& \forall t\\ & X_t\mathrm{|}{\mathbf{f}}_t,r_t\sim {\mathcal{N}}_d\left(\mathit{\mu }({\mathbf{f}}_t,r_t),{\mathit{\eta }}_t^{\mathrm{2}}\right)& \forall t.\end{array}\end{array}$$

The prior on the series of source contributions f_t and pool fraction remaining r_t is now fully independent in time, and the information that could be contained in the fact that some measurements are closer in time than others is ignored.

Prior information can be encoded into the prior distribution for f_t by introducing a concentration parameter $\mathit{\sigma }\in {\mathbb{R}}_{+}^K$ as well as a parameter $\mathit{\alpha }\in (\mathrm{0},\mathrm{1})$ that interpolates between the uniform prior for r_t and the Jeffreys prior using the beta distribution as described in Sect. 2.2. This allows for the inclusion of information that is universally true for all time points simultaneously. If no information is available, the model can be extended by adding an additional hierarchical layer for these parameters, with weakly informative hyperpriors being the gamma distribution Γ(2,2) on the positive real axis for concentrations σ and the uniform Uni(0,1) for α:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& \mathit{\sigma }\sim \mathrm{\Gamma }(\mathrm{2},\mathrm{2}),\mathit{\alpha }\sim \text{Uni}(\mathrm{0},\mathrm{1})\\ & {\mathbf{f}}_t\sim \text{Dir}\left(\mathit{\sigma }\right),r_t\sim \mathrm{B}(\mathit{\alpha },\mathrm{1})& \forall t\\ & X_t\mathrm{|}{\mathbf{f}}_t,r_t\sim {\mathcal{N}}_d\left(\mathit{\mu }({\mathbf{f}}_t,r_t),{\mathit{\eta }}_t^{\mathrm{2}}\right)& \forall t.\end{array}\end{array}$$

2.3.2 Gaussian process priors

Time series information can be incorporated by various methods in the case of unconstrained random variables. An obvious method is to apply the FRAME model with independent time steps on a preprocessed measurement series. The time series preprocessing can be done before model application, without consideration of the Bayesian mixing model. Candidate preprocessing algorithms are kernel smoothing, spline smoothing, and local polynomial regression. While the latter can offer uncertainty estimates of the smoothed time series, a holistic treatment of estimation with uncertainty can be offered by Gaussian process regression.

Despite the possibility of running the simple time-independent model on preprocessed measurement time series, it is beneficial to combine both steps into an advanced model; for example, the problem-specific geometry could influence the feasibility of a region in measurement space. A combined model will include a Gaussian process prior on the measurements X_t such that posterior means W_t can be estimated and used to drive the FRAME model with independent time steps (Sect. 12). The Gaussian process is shifted and scaled to align with the empirical mean ${\widehat{\mathit{\mu }}}_X$ and standard deviation ${\widehat{\mathit{\sigma }}}_X$ of the measurements X_t and controlled by a kernel function G (Fischer, 2023):

$$\begin{array}{}\text{(14)}& \begin{array}{rlr}& {\mathbf{f}}_t\sim \text{Dir}\left(\mathrm{1}\right),r_t\sim \text{Uni}(\mathrm{0},\mathrm{1})& \forall t\\ & W_t\mathrm{|}{\mathbf{f}}_t,r_t\sim {\mathcal{N}}_d\left(\mathit{\mu }({\mathbf{f}}_t,r_t),{\displaystyle \frac{{\mathit{\eta }}_t^{\mathrm{2}}}{\mathrm{2}}}\right)& \forall t\\ & W\sim \mathcal{G}{\mathcal{P}}_d\left({\widehat{\mathit{\mu }}}_X,{\widehat{\mathit{\sigma }}}_X,\mathbf{G}\right)\\ & X_t\mathrm{|}{W}_t,{\mathbf{f}}_t,r_t\sim {\mathcal{N}}_d\left(W_t,{\displaystyle \frac{{\mathit{\eta }}_t^{\mathrm{2}}}{\mathrm{2}}}\right)& \forall t.\end{array}\end{array}$$

The distribution on the latent estimates W_t is the product of the Gaussian process prior as well as the independent normal distribution around the mixing estimate. Ideally, this model does not need to include sampling of f_t and r_t because if the mixing equation can be expressed as a linear system of equations (Eq. 7), then the smooth measurement series W_t is sufficient to solve for the source contribution and fractionation parameters directly. In practice, this approach reduces to applying isotopic mapping techniques to the time series that is preprocessed using Gaussian process smoothing.

If the mixing equation is not explicitly inverted but evaluated by sampling the parameters f_t and r_t, then the latent variables W can be marginalized over and eliminated from the model. The product density of W can also be expressed using known identities (Pedersen and Petersen, 2012) for each separate isotopic measurement dimension $j=\mathrm{1},\mathrm{\dots },d$ in terms of its empirical mean ${\widehat{\mathit{\mu }}}_{X,j}$, empirical standard deviation ${\widehat{\mathit{\sigma }}}_{X,j}$, and noise variance ${\mathit{\eta }}_j^{\mathrm{2}}$. Using the Cholesky decomposition, these distribution parameters can efficiently be computed and used for sampling; thus, the latent parameters W can be eliminated from the model, and the likelihood of each row X_j: can be directly computed (Fischer, 2023):

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& \mathbf{F}\sim \text{Dir}\left(\mathrm{1}\right),\mathbf{r}\sim \text{Uni}(\mathrm{0},\mathrm{1})\\ & \mathit{\mu }(\mathbf{F},\mathbf{r})\sim \mathcal{G}{\mathcal{P}}_d\left({\widehat{\mathit{\mu }}}_X,{\widehat{\mathit{\sigma }}}_X,\mathbf{G}\right)\\ & {\mathbf{X}}_{j:}^T\mathrm{|}\mathbf{F},\mathbf{r}\sim {\mathcal{N}}_N\left({\stackrel{\mathrm{̃}}{\mathit{\mu }}}_j,{\stackrel{\mathrm{̃}}{\mathrm{\Sigma }}}_j+{\displaystyle \frac{{\mathit{\eta }}_t^{\mathrm{2}}}{\mathrm{2}}}\mathbb{1}\right)& \forall j=\mathrm{1},\mathrm{\dots },d.\end{array}\end{array}$$

Gaussian process priors on measurements use only one single hyperparameter, which is the correlation length ρ used to compute the kernel matrix ${\mathbf{G}}_{ij}={\mathit{\kappa }}_{\mathit{\rho }}({\mathit{\tau }}_i,{\mathit{\tau }}_j)$. The scale of the Gaussian process is always set to the empirical standard deviation of the data and is thus considered fixed. In order to compile a fully hierarchical Bayesian model, an inverse gamma distribution $\frac{\mathrm{1}}{\mathit{\rho }}\sim \mathrm{\Gamma }(\mathrm{2},\mathrm{2})$ can be used as hyperprior for the correlation length, assuming that the timescales are properly normalized.

Figure 1 Examples of a generalized Gaussian process prior with the radial basis function kernel using different values for correlation length (ρ) and concentration (σ). (a) Prior observation for three sources (f₁, f₂, f₃) mapped to the simplex using the centered log ratio transform, shown over an arbitrary time axis using different values of ρ with σ=1. (b) Estimated marginal densities (left-hand side) of transformed Gaussian process priors for different values of the concentration parameters σ with ρ=0.1, with a prior observation for three sources (f₁, f₂, f₃) shown over an arbitrary time axis on the right-hand side.

2.3.3 Dirichlet–Gaussian process priors

To make use of time series information in the source contributions f and the fraction reacted r, direct priors are desired. These priors can be constructed by sampling auxiliary variables from multiple independent Gaussian processes Z∼𝒢𝒫_K(G) and at each point in time inverting the log ratio transformations on the simplex in order to create a time series of simplex-valued variables f_t. The fraction reacted r_t is constrained to the interval [0,1] and can thus be linked for instance by applying the logit transform logit(r) = $\ln \frac{r}{\mathrm{1}-r}$ at each point in time, which maps it to the entire real axis. Hyperparameters for correlation length $\mathit{\rho }\in {\mathbb{R}}_{+}$ and concentration $\mathit{\sigma }\in {\mathbb{R}}_{+}$ are used to compute the kernel matrix ${\mathbf{G}}_{ij}={\mathit{\sigma }}^{\mathrm{2}}{\mathit{\kappa }}_{\mathit{\rho }}({\mathit{\tau }}_i,{\mathit{\tau }}_j)$ for the Gaussian process. The general shape of these priors is visualized in Fig. 1. Working with the matrix of source contributions $\mathbf{F}=[{\mathbf{f}}_{\mathrm{1}},\mathrm{\dots },{\mathbf{f}}_n]\in {\mathbb{R}}^{K\times N}$ and of the fraction reacted $\mathbf{r}=[r_{\mathrm{1}},\mathrm{\dots },r_n]\in {\mathbb{R}}^{\mathrm{1}\times N}$, the model can be stated in vectorized form, where the link functions are understood to be columnwise (Fischer, 2023):

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& \text{clr}\left(\mathbf{F}\right)\sim \mathcal{G}{\mathcal{P}}_K\left(\mathbf{G}\right)\\ & \text{logit}\left(\mathbf{r}\right)\sim \mathcal{G}\mathcal{P}\left(\mathbf{G}\right)\\ & X_t\mathrm{|}{\mathbf{f}}_t,r_t\sim {\mathcal{N}}_d\left(\mathit{\mu }({\mathbf{f}}_t,r_t),{\mathit{\eta }}_t^{\mathrm{2}}\right)& \forall t.\end{array}\end{array}$$

Both link functions used can easily be inverted once random variables Z∼𝒢𝒫_K(G) and Y∼𝒢𝒫(G) are sampled from Gaussian processes over time $t=\mathrm{1},\mathrm{\dots },N$. The inverse of the clr transform is given by the softmax function, and the inverse of the logit link is given by the sigmoid function (see Fischer, 2023). The prior on the source contribution parameters G is known as a generalized Gaussian process prior, and techniques such as Taylor expansion can be used to derive analytic approximations (Chan and Dong, 2011). Its marginal is a softmax transformed multivariate Gaussian, which is also known as a logistic normal distribution and serves as an approximation to the Dirichlet distribution (Aitchison and Shen, 1980; Devroye, 1986; Fischer, 2023). Mapping using centered log ratio (clr) transforms thus creates a time series of random variables with approximate Dirichlet marginals, which is referred to as a Dirichlet–Gaussian process (DGP; Chan, 2013).

The marginals are controlled by the parameter σ of the Gaussian process, which now acts as the concentration parameter of the Dirichlet distribution. Since the covariance kernel G is scaled to generate Gaussian random variables with unit variance if σ=1, the marginal distribution in that case is approximately the uniform Dir(1). This can be seen by sampling from the generalized Gaussian process priors and estimating the marginals, as shown in Fig. 1b.

Using the isometric log ratio transform ilr instead of clr reduces the number of Gaussian processes that need to be sampled to K−1 for the source contributions. Inverting this link function is accomplished by applying an orthonormal base transform U to the random variables and then applying the softmax function. Since interpretability of the sampled Gaussian process variables is not required, any orthonormal basis is suitable, and a simple construction using Gram–Schmidt orthogonalization is chosen (Nesrstová et al., 2022; Fischer, 2023):

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& \text{ILR}\left(\mathbf{F}\right)\sim \mathcal{G}{\mathcal{P}}_{K-\mathrm{1}}\left(\mathbf{G}\right)\\ & \text{logit}\left(\mathbf{r}\right)\sim \mathcal{G}\mathcal{P}\left(\mathbf{G}\right)\\ & X_t\mathrm{|}{\mathbf{f}}_t,r_t\sim {\mathcal{N}}_d\left(\mathit{\mu }({\mathbf{f}}_t,r_t),{\mathit{\eta }}_t^{\mathrm{2}}\right)& \forall t.\end{array}\end{array}$$

Shorthand notation $\mathcal{G}\mathcal{P}\left(\mathbf{G}\right)=\mathcal{G}\mathcal{P}(\mathrm{0},\mathrm{1},\mathbf{G})$ is used with correlation length ρ and scale σ included in the kernel computation ${\mathbf{G}}_{ij}={\mathit{\sigma }}^{\mathrm{2}}{\mathit{\kappa }}_{\mathit{\rho }}({\mathit{\tau }}_i,{\mathit{\tau }}_j)$. Both kernel parameters ρ and σ can be set in advance or given weak hyperpriors. The inverse gamma distribution $\frac{\mathrm{1}}{\mathit{\rho }}\sim \mathrm{\Gamma }(\mathrm{2},\mathrm{2})$ and the regular gamma distribution $\mathit{\sigma }\sim \mathrm{\Gamma }(\mathrm{2},\mathrm{2})$ are chosen under the assumption that the time variables ${\mathit{\tau }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\tau }}_N$ are scaled appropriately. This hierarchical model benefits especially from the reduced number of Gaussian processes sampled when using the ilr transform, since the kernel covariance matrix must be reconstructed in every sampling step. The number of hyperparameters could be increased by using separate concentrations and correlation lengths for the source contributions F and the fraction remaining r.

Figure 2 Examples of the spline prior for three source contributions (f₁, f₂, f₃, where ${\mathit{f}}_{\mathrm{1}}+{\mathit{f}}_{\mathrm{2}}+{\mathit{f}}_{\mathrm{3}}=\mathrm{1}$) transformed to the simplex with the clr transformation using different degrees of freedom M that can be used to control the covariance of source contributions at separate points in time.

2.3.4 Spline-based priors

An alternative to Gaussian process priors is spline basis functions, which can be used to construct a linear fitting operation that is then mapped to simplex space. This allows for the addition of exogenous variables as predictors of source contributions or fractionation. A cubic spline basis of M basis functions (Fig. 2) is evaluated at the measurement points ${\mathit{\tau }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\tau }}_N$ to form the evaluation matrix $\mathbf{H}\in {\mathbb{R}}^{N\times M}$ with H_ij=s_j(τ_i) for polynomial basis functions $s_{\mathrm{1}}(\cdot ),\mathrm{\dots },s_M(\cdot )$. The time series of source contributions in simplex space is reconstructed with the basis coefficients b_k∈ℝ^M for each source $k=\mathrm{1},\mathrm{\dots },K$ arranged to the matrix $[{\mathbf{b}}_{\mathrm{1}}\mathrm{\cdots }{\mathbf{b}}_K{]}^T=\mathbf{B}\in {\mathbb{R}}^{L\times M}$ and coefficients for fractionation $\mathbf{c}\in {\mathbb{R}}^{\mathrm{1}\times M}$. This type of model is therefore part of the generalized linear model class (Nelder and Wedderburn, 1972) and allows for easy extension with fixed effects relating to measurement dimensions as well as random effects for experiment replication. It will thus further be referred to as the generalized linear model with spline basis (spline GLM; Fischer, 2023):

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& \mathbf{B},\mathbf{c}\sim \mathcal{N}(\mathrm{0},\mathrm{1})\\ & \text{clr}\left(\mathbf{F}\right)={\mathbf{BH}}^T\\ & \text{logit}\left(\mathbf{r}\right)={\mathbf{cH}}^T\\ & X_t\mathrm{|}{\mathbf{f}}_t,r_t\sim {\mathcal{N}}_d\left(\mathit{\mu }({\mathbf{f}}_t,r_t),{\mathit{\eta }}_t^{\mathrm{2}}\right)& \forall t.\end{array}\end{array}$$

In consequence, the distribution of the basis coefficient vector ${\mathbf{b}}_k={\mathbf{B}}_{k:}^T\in {\mathbb{R}}^M$ before transformation has distribution ${\mathbf{b}}_k\sim {\mathcal{N}}_M(\mathrm{0},\mathbb{1})$ for source $k=\mathrm{1},\mathrm{\dots },K$. After application of the spline basis transform it is thus still Gaussian ${\mathbf{Hb}}_k\sim {\mathcal{N}}_N(\mathrm{0},{\mathbf{HH}}^{\mathbf{T}})$, although with a modified covariance matrix ${\mathbf{HH}}^{\mathbf{T}}\in {\mathbb{R}}^{N\times N}$. Since the inverse centered log ratio transform maps Gaussian random variables with unit variance approximately to a uniform Dirichlet distribution, it makes sense to scale the basis transform such that $\frac{\mathrm{1}}{N}\text{Tr}\left({\mathbf{HH}}^T\right)=\mathrm{1}$, as the spline basis vectors are not semi-orthogonal in general.

2.4 Prior distribution for the fraction remaining r

Sampling the prior distribution of the fractionation weight r for closed-system Rayleigh fractionation is challenging, because it is connected through the non-linear logarithm to the effect on measurements. Although a uniform prior usually does not inform the posterior about anything except the boundaries, the effect of the logarithm on the posterior is much more unclear. Different choices for prior distributions of r were thus tested for their effect on the generated posterior sample. The simulated dataset used 17 different values for the fractionation index d ranging from 0.05 to 0.95. Source contributions were fixed to ${\mathit{f}}_{\mathrm{1}}^{\ast }=\mathrm{0.7}$ and ${\mathit{f}}_{\mathrm{2}}^{\ast }=\mathrm{0.3}$. Each value of r is used to generate Q=64 data points $X^{\left(\mathrm{1}\right)},\mathrm{\dots },X^{\left(Q\right)}$ with measurement error η=4 for a total of 1088 data points (Fig. 4a). The stationary inference model given in Eq. (9) was fitted to each point individually, which makes this setting analogous to the inference procedure used in the original FRAME model (Lewicki et al., 2022).

2.5 Model comparison

2.5.1 Data simulation

No time series datasets with known N₂O production and consumption pathway contributions are available; therefore, simulated data must be used for thorough comparison of models. The four models were compared by simulating the data-generating process multiple times and then comparing the resulting posterior sample with fixed truth input values. The time series of source contributions f and pool fractions reacted r used to simulate the data are denoted by ${\mathbf{F}}^{\ast }=\left[{\mathbf{f}}_{\mathrm{1}}^{\ast }\mathrm{\cdots }{\mathbf{f}}_N^{\ast }\right]$ and ${\mathbf{r}}^{\ast }=\left[{\mathbf{r}}_{\mathrm{1}}^{\ast }\mathrm{\cdots }{\mathbf{r}}_N^{\ast }\right]$, and the mixing equation with Rayleigh fractionation is used (Eq. 3). Measurement generation is then repeated Q times by sampling the source isotopic signature $S^{\left(q\right)}\in {\mathbb{R}}^{d\times K}$ and fractionation factor ϵ^(q)∈ℝ^d from their respective priors and then adding independent Gaussian measurement errors $E_t^{\left(q\right)}\sim {\mathcal{N}}_d(\mathrm{0},{\mathit{\eta }}^{\mathrm{2}})$ with noise variance η²∈ℝ^d for $q=\mathrm{1},\mathrm{\dots },Q$ (Fischer, 2023):

$$\begin{array}{}\text{(19)}& X_t^{\left(q\right)}={\mathbf{S}}^{\left(q\right)}{\mathbf{f}}_t^{\ast }+{\mathit{ϵ}}^{\left(q\right)}\ln \left(r_t^{\ast }\right)+E_t^{\left(q\right)}.\end{array}$$

Table 1 Prior distribution parameters for N₂O source isotopic signatures and the fractionation factor for consumption used to simulate datasets for model testing (Yu et al., 2020). Ranges for sources indicate the full range of previous data used to construct the uniform distribution, whereas for reduction the range indicates the standard deviation of the Gaussian distribution.

Auxiliary data for source isotopic signatures and fractionation factors are taken from Yu et al. (2020) (Table 1). These values correspond to the major N₂O sources, nitrification (S₁) and bacterial denitrification (S₂). Priors are uniform for the sources $S_j\sim \text{Uni}(b_j,{\mathrm{\Delta }}_j),j=\mathrm{1},\mathrm{2}$ and Gaussian for the fractionation factor with variance matched to the reported bounds $\mathit{ϵ}\sim \mathcal{N}(c,\frac{{\mathrm{\Delta }}_{\mathit{ϵ}}^{\mathrm{2}}}{\mathrm{2}})$.

Figure 3 True parameter series used to simulate datasets with different properties. The production–consumption scenario corresponds to production by nitrification (S₁, red) and denitrification (S₂, green) followed by mixing and subsequently reduction (r is the fraction remaining following reduction, blue) in complete denitrification. (a) Six parameter series to test time series modeling, which illustrate fast-changing source contributions (FastS), slow-changing source contributions (SlowS), extremal source contributions (ExtrS), high fractionation (HighF), average and variable fractionation (AvgF), and low fractionation (LowF). (b) The general example (GenE) and the resultant isotopic time series show measurement values simulated accordingly together with LOESS estimates on the right.

The simulations are done on fixed parameter sets that intend to be illustrative for six given cases that might occur in reality. The focus is on investigation of temporal patterns; other fractionation scenarios have been explored previously using a stationary model setup (Lewicka-Szczebak et al., 2020; Lewicki et al., 2022). The true parameter time series ${\mathbf{f}}_t^{\ast }$ and $r_t^{\ast }$ are shown in Figure 3 and are sampled at N=32 equally spaced time steps. One additional general example (GenE) is used for simulation with properties being less extreme than for the other six, which may be more representative of average datasets that would be encountered in practice. For GenE, a Gaussian error with magnitude η=5 is used to sample N=64 measurements $X_{\mathrm{1}},\mathrm{\dots },X_N$. The fixed parameter values and the simulated data are shown in Fig. 3.

Bayesian parameter estimation is then tested on each generated dataset ${\mathbf{X}}^{\left(q\right)}=\left[X_{\mathrm{1}}^{\left(q\right)}\mathrm{\cdots }{X}_N^{\left(q\right)}\right]\in {\mathbb{R}}^{d\times N}$ for $q=\mathrm{1},\mathrm{\dots },Q$ individually, and a total of S posterior samples of all parameters is produced each time. The posterior samples shall be denoted by ${\mathbf{F}}^{(q,s)}=\left[{\mathbf{f}}_{\mathrm{1}}^{(q,s)}\mathrm{\cdots }{\mathbf{f}}_N^{(q,s)}\right]$ and ${\mathbf{r}}^{(q,s)}=\left[r_{\mathrm{1}}^{(q,s)}\mathrm{\cdots }{r}_N^{(q,s)}\right]$, respectively, for $s=\mathrm{1},\mathrm{\dots },S$.

2.5.2 Measuring quality of inference

Sampling from the posterior distribution does not give unique point estimates for the parameters involved, and multiple ways of computing final parameter estimates exist. Most commonly the posterior mean is used as point estimate, although using the median could, for example, be a useful strategy for posterior distributions that are highly dissimilar to a Gaussian distribution.

The accuracy of the estimation can be assessed by computing the distance between these pointwise estimates and the true value using root mean squared error (RMSE) and mean average error (MAE). Although it would be possible to compute the metrics at each point in time, they are averaged for simpler model comparison:

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }{\text{RMSE}}_k^{\left(q\right)}:=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{t=\mathrm{1}}^N{\left({\widehat{\mathbf{F}}}_{kt}^{\left(q\right)}-{\mathbf{F}}_{kt}^{\ast }\right)}^{\mathrm{2}}},\text{(21)}& {\displaystyle }& {\displaystyle }{\text{MAE}}_k^{\left(q\right)}:={\displaystyle \frac{\mathrm{1}}{N}}\sum _{t=\mathrm{1}}^N\mathrm{|}{\widehat{\mathbf{F}}}_{kt}^{\left(q\right)}-{\mathbf{F}}_{kt}^{\ast }\mathrm{|}.\end{array}$$

Computations for r^(q,s) are analogous.

Since the parameters to be estimated are interpreted as a time series, it makes sense to also compare specific time series information across the model estimates. The rate of change can be significantly confounded in the measurement time series, since the measurement errors follow a white noise distribution that introduces high-frequency changes. The ability of models to filter this noise can be measured by comparing the rate of change, which is approximated using first differences: $\mathrm{\Delta }{\mathbf{F}}_{t,j}^{\ast }={\mathbf{F}}_{t,j+\mathrm{1}}^{\ast }-{\mathbf{F}}_{t,j}^{\ast }$ and $\mathrm{\Delta }{r}_t^{\ast }=r_{t+\mathrm{1}}^{\ast }-r_t^{\ast }$ for $t=\mathrm{1},\mathrm{\dots },N-\mathrm{1}$. The magnitude of changes is not necessarily relevant, since poor estimates of the magnitude of changes would also lead to poor pointwise comparison metrics. Therefore, the ratio of variances of first differences shall serve as comparison metric for time series information, which can be understood as comparing a notion of curvature or acceleration in the time series (Fischer, 2023):

$$\begin{array}{}\text{(22)}& {\text{VarFD}}_k^{\left(q\right)}={\displaystyle \frac{{\sum }_{t=\mathrm{1}}^{N-\mathrm{1}}{\left(\mathrm{\Delta }{\mathbf{F}}_{kt}^{\left(q\right)}-{\sum }_{j=\mathrm{1}}^{N-\mathrm{1}}{\mathbf{F}}_{kj}^{\left(q\right)}\right)}^{\mathrm{2}}}{{\sum }_{t=\mathrm{1}}^{N-\mathrm{1}}{\left(\mathrm{\Delta }{\mathbf{F}}_{kt}^{\ast }-{\sum }_{j=\mathrm{1}}^{N-\mathrm{1}}{\mathbf{F}}_{kj}^{\ast }\right)}^{\mathrm{2}}}}.\end{array}$$

2.5.3 Metrics for comparison of Bayesian posterior distributions

Bayesian models are mainly used to derive pointwise estimates, but their advantage is the creation of a sample from the posterior distribution. It is thus also important to take distribution properties into account. Posterior interval coverage is a useful metric to evaluate simulated data from the full Bayesian model, meaning in particular that F^∗ and r^∗ are sampled from their priors as well. Since the data simulation for model comparison used here has fixed parameter values, interpreting interval coverage becomes less meaningful. It is still practical to use the size of the credible interval as measure of uncertainty, and thus the interval span can be compared, which is desired to be as small as possible given otherwise accurate estimates. The credible interval with level $\mathit{\gamma }\in [\mathrm{0},\mathrm{1}]$ is then the set excluding the tails with a proportion of $\frac{\mathrm{1}-\mathit{\gamma }}{\mathrm{2}}$ of the most extreme observations on either side (Fischer, 2023):

$$\begin{array}{}\text{(23)}& I\left(\mathit{\gamma }\right)=\left[q\left({\displaystyle \frac{\mathrm{1}-\mathit{\gamma }}{\mathrm{2}}}\mathrm{|}\mathbf{X}\right),q\left({\displaystyle \frac{\mathrm{1}+\mathit{\gamma }}{\mathrm{2}}}\mathrm{|}\mathbf{X}\right)\right].\end{array}$$

Further metrics for the quality of the entire posterior distributions can be taken into consideration. Posterior predictive checks are typically used in cases where no true values for the parameter estimates are available, in order to assess the models' capability of representing the input data well (Rubin, 1984). The posterior predictive distribution is the likelihood of hypothetical future measurements calculated using the posterior distribution over parameter values in the presence of the actually available data:

$$\begin{array}{}\text{(24)}& p\left(\stackrel{\mathrm{̃}}{\mathbf{X}}\mathrm{|}\mathbf{X}\right)=\int p(\stackrel{\mathrm{̃}}{\mathbf{X}}\mathrm{|}\mathbf{F},\mathbf{r})\mathit{\pi }(\mathbf{F},\mathbf{r}\mathrm{|}\mathbf{X})\mathrm{d}\mathbf{F}\mathrm{d}\mathbf{r}.\end{array}$$

This posterior predictive distribution is not unique due to the auxiliary parameters S and ϵ. It is unclear whether the posterior predictive distribution should be proportional to the marginalized likelihood $p(\stackrel{\mathrm{̃}}{\mathbf{X}}\mathrm{|}\mathbf{F},\mathbf{r})$ or rather the likelihood conditioned on the auxiliary parameters $p(\stackrel{\mathrm{̃}}{\mathbf{X}}\mathrm{|}\mathbf{F},\mathbf{r},\mathbf{S},\mathit{ϵ})$. This discrepancy renders comparison of predictive density values across dataset simulations X^(q) ineffective, since the models fit a joint posterior and thus assume that future data must be sampled using identical auxiliary parameter values, whereas the simulation resamples their S^(q) and ϵ^(q) values every time.

The log pointwise predictive density (LPPD) is a metric for the quality of the posterior predictive density and thus by proxy of the Bayesian model. It is computed by evaluating the predictive density at the original data points and can be approximated using the posterior samples (Gelman et al., 2014):

$$\begin{array}{}\text{(25)}& \begin{array}{rl}& {\text{LPPD}}^{\left(q\right)}=\ln \int p(\mathbf{X}\mathrm{|}\mathbf{F},\mathbf{r})\mathit{\pi }(\mathbf{F},\mathbf{r}\mathrm{|}\mathbf{X})\\ & \mathrm{d}\mathbf{F}\mathrm{d}\mathbf{r}\approx {\displaystyle \frac{\mathrm{1}}{S}}\sum _{s=\mathrm{1}}^{S}p\left(\mathbf{X}\mathrm{|}{\mathbf{F}}^{\left(s\right)}\right).\end{array}\end{array}$$

2.6 Model implementation

TimeFRAME is implemented in R R Core Team (2017). Bayesian modeling in TimeFRAME uses Stan for Hamiltonian Monte Carlo sampling (see Sect. 2 in the Supplement). TimeFRAME can be installed using the links provided in the “Code and data availability” section. All experiments were run on an Intel Core i9-10900K CPU. The reported run times in Sect. 2 in the Supplement are for a single sampling chain, and in other sections the reported times are the maximum of four simultaneously run chains.

Table 2 Model configurations used to compare performance on time series data across different edge cases for the four model classes described in Sect. 2.3. The “abbreviation” column refers to the abbreviations for each configuration used in subsequent discussion.

3 Results

In this section, we will first discuss the influence of the prior distribution for r and then present a detailed comparison of different aspects of model performance. The four classes of models were compared for performance on time series data across a wide range of scenarios representative of edge cases that might occur in reality. For these experiments, fixed parameter values for source contributions F^∗ and fraction remaining r^∗ were sampled as described in Sect. 2.5.1. The model configurations to be compared are summarized in Table 2.

3.1 Influence of the chosen prior distribution for the fraction remaining r

Using a uniform prior for the fraction reacted π(r)∝1 is the natural choice used as standard by all models. This was compared to the Jeffreys prior, which was derived to be $\mathit{\pi }\left(r\right)\propto \frac{\mathrm{1}}{r}$ and thus is an improper prior (see Sects. 2.2 and Fischer, 2023, for details). A middle ground between these two choices was given by the beta prior $r\sim B(\frac{\mathrm{1}}{\mathrm{2}},\mathrm{1})$, which has $\mathit{\pi }\left(r\right)\sim \frac{\mathrm{1}}{\sqrt{r}}$. The first argument of the beta distribution can also be used as a free parameter $\mathit{\alpha }\in (\mathrm{0},\mathrm{1}]$ alongside a uniform hyperprior to construct the hierarchical model $r\sim B(\mathit{\alpha },\mathrm{1}),\mathit{\alpha }\sim \text{Uni}(\mathrm{0},\mathrm{1})$.

Figure 4 Comparison of different priors for the fraction remaining r. (a) Dual-isotope plot for the simulated data, constructed with two sources and different values of r. Data points are colored by the true value of r. (b) Posterior densities of the fractionation weight r averaged across simulations for true r values of 0.05, 0.5, and 0.95 using different prior distributions. (c) Mean absolute error of Bayesian models using different prior distributions of r, shown for different true values of r and f. Vertical lines indicate the standard deviation over Q=64 repetitions. The performance on source contributions is identical for both sources (Ni and bD), since they are perfectly correlated, so only one panel is shown for f.

Each data point was supplied to the model individually, and S=5000 posterior samples were generated. The samples were combined for each value of r to marginalize over the distributions of the auxiliary parameters (Fig. 4b). The source S₂ (bD) has a wide prior distribution (Table 1) which confounds with the effect of fractionation and introduces high uncertainty. Since this uncertainty is proportional to the logarithm of r, the effect on the posterior distribution is much more pronounced when r is high. The Jeffreys prior introduces a shift towards lower values compared to the uniform prior, with the beta prior being in between the two. Inference performance was compared by taking the posterior means ${\widehat{\mathbf{f}}}^{\left(q\right)}=\frac{\mathrm{1}}{S}{\sum }_{s=\mathrm{1}}^S{\mathbf{f}}^{(q,s)}$ and ${\widehat{\mathbf{r}}}^{\left(q\right)}=\frac{\mathrm{1}}{S}{\sum }_{s=\mathrm{1}}^S{\mathbf{r}}^{(q,s)}$ and comparing them against the true values (Fig. 4c). Clearly, the Jeffreys prior performs worst for source contributions f, with the uniform prior performing best, closely followed by the beta and hierarchical prior. Regarding the fraction remaining r, the Jeffreys prior performs best for low values of r, with the uniform prior performing worst; however, this relationship switches at about r=0.4 to the opposite. Therefore, choosing any prior can be justified if one expects the fractionation index to be in a certain range. The standard deviation between different repetitions, however, clearly shows that the effect of prior choice is overwhelmed by the variation introduced through the distribution of the sources and consumption, due to the large uncertainty in the source and fractionation factor priors.

Figure 5 Model performance for different configurations (see Table 2 for details: fast-changing source contributions (FastS), slow-changing source contributions (SlowS), extremal source contributions (ExtrS), high fractionation (HighF), average and variable fractionation (AvgF), and low fractionation (LowF)). (a) Mean absolute errors for all examples averaged over Q=64 dataset simulations with standard deviations shown as vertical lines. The results for the two source contributions were identical since they are perfectly correlated, so the plots only show one result for f in addition to the fraction remaining r. (b) Log variance ratio of first differences of the estimated time series against the true values for source contributions f. The time series for fraction remaining r are linear or constant in most examples and are thus not suitable to be used for model comparison with this metric.

3.2 Comparison of overall model performance

The examples described in Sect. 2.5.1 were sampled for Q=64 repetitions. Measurements were simulated with a Gaussian measurement error of magnitude η=5. The posteriors were sampled for a total of S=10 000 steps using four parallel chains. Goodness of estimation was quantified with estimates computed as the posterior means ${\widehat{\mathbf{f}}}^{\left(q\right)}=\frac{\mathrm{1}}{S}{\sum }_{s=\mathrm{1}}^S{\mathbf{f}}^{(q,s)}$ and ${\widehat{\mathbf{r}}}^{\left(q\right)}=\frac{\mathrm{1}}{S}{\sum }_{s=\mathrm{1}}^S{\mathbf{r}}^{(q,s)}$ and taking the mean absolute error (Eq. 21) to the ground truth (Fig. 5a). Conclusions using the root mean squared error or with posterior medians were similar. All models using fixed hyperparameters use default values that are not specifically tuned for the examples at hand. Therefore, the reported performance is not indicative of best-case performance and only shows the quality of the chosen values. Hierarchical models do not have this problem since they can estimate the hyperparameters for each example specifically.

Overall performance was best for hierarchical DGP models and spline GLMs (Fig. 5a). Estimation of fraction reacted r is relatively poor in most cases (MAE > 0.4); however, spline models perform well for all cases (MAE usually < 0.3), particularly for extremal fractionation amounts (all examples except AvgF). The default number of degrees of freedom that spline models use seems surprisingly robust in all examples, whereas the default correlation length of Gaussian processes does not. Gaussian process priors on measurements appear to be slightly worse than DGP priors and spline-based priors, especially SlowS and ExtrS. Independent time step models have worse performance than the rest for all examples, and the hierarchical extension to it only has good performance in examples SlowS and ExtrS, which represent cases where the contribution of the sources is either very low or very high. In other cases, using flat priors as default values seems to work best.

Spline GLMs have very small errors on fraction consumed r whenever the value is slowly changing and close to 0 or 1. This could suggest that the chosen hyperparameters are suitable for all examples. Another plausible explanation is the fact that the spline bases used have an intercept term, which allows the center of estimation to freely move, whereas DGP models do not. For the source contributions this likely does not matter, but allowing the Gaussian process to have non-zero mean could be beneficial for estimating r. An additional spline model was added with Laplace priors on the coefficients, which seems to be beneficial in cases where parameters are close to their boundaries since the prior allows for values farther from zero. The model using Gaussian process priors on measurements was implemented using both formulations with latent variables and with analytically computed likelihood. Performance is identical in all cases, strongly indicating that the formulations are equivalent for parameter estimation. Additionally, using clr or ilr transformations for DGP models and Spline GLMs does not make a difference in estimation accuracy, as is to be expected from their derivations.

Time series models are not only expected to give accurate estimates of mean source contributions and fractionation, but the resulting time series should also have similar properties to the ground truth. The variance ratio of first differences (Eq. 25) measures accuracy of the estimated curvature (Fig. 5b). Clearly, the hierarchical Gaussian process and DGP models estimate the correlation length well, resulting in a time series with similar rates of change to the true values. Spline GLMs perform well, especially for the examples FastS, SlowS, and ExtrS. Independent time step models result in high variation, which is due to the fact that the measurement noise is not adequately filtered. The fixed correlation length Gaussian processes seem to have misspecified hyperparameters, since they also overestimate the rates of change in the time series.

Fitting times ranged from 20 to 72 s on average for the S=10 000 posterior samples generated by each model, split into 2500 over four chains (see Fischer, 2023, for details). Hierarchical models tend to be slowest due to the additional parameters and repeated matrix decompositions that need to be computed, whereas fixed parameter models, especially independent time steps and spline GLMs, are sampled fastest. Spline models using the Laplace prior have long fitting times, which could indicate that the high parameter values – which are allowed due to weaker regularization of parameter ranges far from zero – are not sufficiently identifiable, resulting in slow sample generation.

3.3 Influence of fractionation extent on model performance

Models can have varying performance at different reaction extents and thus different levels of fractionation. For this reason, a time series of source contribution values ${\mathbf{f}}_t^{\ast }$ was taken and paired with different constant fractionation values $r_t^{\ast }=r^{\ast }$ to generate measurements and monitor performance. In total, 17 equally spaced fractionation values ranging from $r_{\ast }=\mathrm{0.02}$ to 0.98 were used, and a total of Q=32 datasets were generated per value. True values for source contributions were taken from the general example GenE (Sect. 2.5.1), and a measurement error magnitude of η=5 was used. This experiment was done only with representative models of the four main model classes in order to reduce the number of comparisons. The model configurations used were the independent time step model, the Gaussian process prior on measurements with hierarchical estimation of correlation length, DGP prior with hierarchical estimation of correlation length, and spline-based GLM with fixed hyperparameters for degrees of freedom.

Figure 6 Impact of the fraction remaining r on model performance. Each reported value is the average over Q=32 dataset simulations, with vertical lines indicating standard deviations. (a) Mean absolute error of the four main model classes over different values of the fraction remaining r. (b) The 95 % credible interval spans of the four main model classes over different values of the fraction remaining r.

Overall performance is best when r^∗ is close to 0.5 (Fig. 6a). Estimation is also more accurate with very low values of r^∗, because small r^∗ values have large impacts on isotopic measurements, and thus estimation can become more accurate. Spline models are expected to perform well here because the time series of fraction remaining is constant, which can be reflected in the low degrees of freedom used. The different model classes seem to be equally affected by changes in r^∗ otherwise, showing that the choice of hyperparameters to reflect the situation of interest is more important than selecting a particular model.

Spans of 95 % credible intervals can give additional insight into the pattern observed for the estimation accuracy over different values of r^∗ (Fig. 6b). If parameter estimation is good, then a smaller credible interval span shows a narrow posterior around the correct mean. DGP prior models have the smallest credible interval span for source contributions and fraction remaining. All other models have an interval width of over 0.75 for large values of fraction remaining, thus spanning over half of the possible domain. Clearly, due to the Rayleigh fractionation equation being non-linear in r, it is difficult to estimate large values of r (resulting in low fractionation) with high accuracy. If the amount of remaining substrate is larger than 0.1, the data do not give sufficient information regarding r, so estimates of r group around the prior mean of 0.5 and show large credible interval spans.

3.4 Influence of measurement noise on model performance

The main advantage that smooth models such as Gaussian processes and splines have over the independent time step assumption is that they promise to filter measurement noise and thus produce estimates that are more accurate and have a narrower posterior distribution. For this reason an experiment was conducted using values of source contributions ${\mathbf{f}}_t^{\ast }$ and fraction remaining $r_t^{\ast }$ from the general example GenE (Sect. 2.5.1) to simulate datasets with different levels of measurement noise. Noise values range from η=0.5 to η=20, and for each separate value of η, a total of Q=32 datasets were generated.

Figure 7 Impact of measurement noise on model performance. Each reported value is the average over Q=32 dataset simulations, with vertical lines indicating standard deviations. (a) Mean absolute error of estimation of the four main model classes over different measurement noise levels η. (b) Log variance ratio of first differences for the estimated time series against the true values for source contributions over different measurement noise levels η.

Performance of all models generally decreases with increased measurement noise as expected (Fig. 7a). However, below η=5, reduction in noise does not lead to further improvement in performance: most of the estimation error at this noise level comes from the source end-member uncertainty rather than the measurement noise. The variance ratio of first differences (Fig. 7b) can be used to assess the quality of the estimated time series in the presence of high-frequency changes due to measurement noise. Variance ratios of the independent time series model gradually increase with increasing measurement noise magnitude. All other models seem to filter the noise well, having much lower overestimations of the first difference variance. The hierarchical DGP model seems to be less equipped to deal with very high noise, which could simply be due to the fact that the weakly informative hyperprior on the correlation length is not suitable here. The spline GLM appears to have constant low values for the ratio of first difference variance, due to the fact that the fixed degrees of freedom predetermine the smoothness of the estimates independent of measurement noise. Moreover, the model run time for the spline model was close to the run time for the independent model and less than half the run time for the GP and DGP models across all noise levels. These results show that use of the spline model can be particularly advantageous for data featuring high statistical uncertainty.

Figure 8 Model performance considering different improvements to the input data. The original model performance is shown in gray. (a) Mean absolute error for all models on Q=64 generated datasets using one additional isotopic measurement (δ¹⁸O). (b) Mean absolute error for all models on Q=64 generated datasets using the idealized sources and fractionation factor with 10 % uncertainty in isotopic end-members.

3.5 Potential impact of improvements in data quality

Estimation accuracy can be improved not only by choosing the right model but also by improving the quality of the data available to the model. Several ways of adding more or higher-quality data exist, and the effects on model performance were studied in order to find what would be most beneficial. The above examples use two sources with two isotopic measurements, making the system well determined. If additional isotopic measurements are available, they can be added to make the system overdetermined and thus eliminate some noise. To investigate this approach, the additional isotopic measurement δ¹⁸O was added with source end-members and uncertainty from Yu et al. (2020) (Table 1). The same dataset generation procedure as in Sect. 3.2 was used with a total of Q=64 datasets generated. Resulting improvement in estimation accuracy for the same two sources is shown in Fig. 8a. It is worth noting that the additional measurement is not ideal in quality, with large uncertainty for the fractionation factor ϵ. An improvement can clearly be seen for estimation of the fraction reacted r, especially in the examples HighF and AvgF. Very little improvement is seen for estimation of f. Spline GLMs improved the most, especially in their already good ability to estimate the fraction reacted. The addition of δ¹⁸O to this model did not strongly improve results, due to large uncertainty in source end-members and fractionation factors. However, the addition of isotopic dimensions with low uncertainty or strong differences compared to existing information could improve results, for example, clumped isotopes or δ¹⁸O for determination of fungal denitrification.

Instead of adding additional measurements, efforts could be put into determining the end-members of the sources and the fractionation factors more accurately, thus reducing uncertainty in the input data. To study this case an idealized set of sources and fractionation factors were selected to have the mixing and reduction line exactly perpendicular, with an uncertainty of 10 % in each dimension with respect to the mean. This renders the mixing and fractionation components independent, since they cannot confound each other. We therefore set the end-members for S₁ to $-\mathrm{1}\pm \mathrm{0.2}$ and 1±0.2 ‰ for δ¹⁵N^bulk and δ¹⁵N^SP, respectively, and for S₂ to 1±0.2 and $-\mathrm{1}\pm \mathrm{0.2}$ ‰ for δ¹⁵N^bulk and δ¹⁵N^SP, respectively. The fractionation factor was set to 1±0.1 ‰ for both isotopes. Measurements sampled in this setting follow exactly the same procedure as in Sect. 3.2 but only use a Gaussian measurement error with magnitude η=0.1. For each example, Q=64 datasets were generated, and the mean absolute error of estimation is shown in Fig. 8b.

Improving all uncertainties involved at a minimum has a great impact on model performance. Almost all mean absolute errors of estimation are below an error margin of 0.05 for source contributions and below 0.1 for fraction remaining. Furthermore, model choice becomes less relevant, and even the independent time step models perform similarly to the other more sophisticated ones. Interestingly, the DGP model with fixed hyperparameters and the spline GLM with fixed spline basis underperform in source contribution estimation, for example, HighF. This could be evidence that the default parameters become less robust when noise is removed, and they should be selected more carefully. This experiment is an extremely idealized case, and natural variability likely precludes this level of precision in end-members for microbial N₂O production in soil; however, it shows the high potential for improvements in input data to enhance results and moreover to make results more robust towards model configuration. Currently, the level of uncertainty in direct anthropogenic N₂O and CH₄ source end-members (e.g., industrial production, energy and transport emissions) is very high due to the scarcity of measurements (Eyer et al., 2016; Röckmann et al., 2016; Harris et al., 2017) – further investigation of the isotopic range of these sources, as well as consideration of end-members for novel isotopes such as clumped species, may lead to the level of uncertainty reduction required to achieve accurate source partitioning.

4 Application of TimeFRAME to real and simulated datasets

In this section, we will present use of the TimeFRAME package for the analysis of simulated and real datasets to illustrate aspects of model configuration choice and output data under different scenarios.

4.1 Model selection and application

TimeFRAME allows different models to be applied with minimal effort, meaning that data can be analyzed with several different model setups to investigate the robustness of results. The independent time step model does not incorporate time series information; thus, it is recommended only for datasets with independent measurements. The DGP and spline models both perform well, reproducing the input data values and time series properties – the spline model was better able to estimate r. All models estimate f of different sources across the full range with similar accuracy; however, when the fraction remaining r is very low or high, the results show much larger error (Fig. 8). This is compounded by the difference between mixing followed by reduction (MR) and reduction followed by mixing (RM) models at low values of r (Sect. 1 in the Supplement). We therefore recommend users test both DGP and spline models for time series data and treat results with caution when these models differ strongly. Estimates of very low fraction remaining should also be treated with caution. Despite these points, we find that TimeFRAME offers a strong improvement on previously available methods: accounting for information contained within time series significantly reduces the uncertainty in estimates of f and r, and the package application is simple and fast, as well as easy to document and reproduce.

The testing here focuses on interpretation of N₂O isotope data to unravel production and consumption pathways. TimeFRAME can also be applied to other scenarios, for example, trace gases such as CO₂ or CH₄, or datasets with many more isotopic dimensions through clumped isotope measurements. The number of sources is indefinite as the model can be extended by the user; however, when the number of sources is larger than the number of isotopic dimensions, the model will be poorly constrained. The model can currently only include one consumption pathway applied after mixing – future versions will include more complex setups; however, the uncertainty in input data currently precludes this level of complexity. The examples shown here use time as the dimension of autocorrelation, as time series are the most common kind of data. However, other dimensions could be used, such as temperature in the case of measurements across a gradient of temperature-dependent processes or soil moisture for a set of incubations across a moisture gradient. TimeFRAME's setup allows simple adaptation to different user-defined mixing and fractionation models and fast and reproducible interpretation of these models.

4.2 Application to the general simulated time series

The main purpose of the time series models is to provide estimates of source contribution and fractionation weights with uncertainty. In the sections above, only the performance metrics aggregated over many simulations have been shown. To illustrate the modeling capabilities, the representative general example (GenE) was simulated from fixed parameter values, and the inference results are shown in comparison to the true values.

Figure 9 Posterior means of the three model types compared to the true parameter values. Shaded areas indicate 95 % credible intervals, and the true parameter values used to simulate the measurements are shown as black lines.

In order to run the Bayesian models and estimate source contributions and fractionation over time, the auxiliary distributions of the source isotopic signatures S and the fractionation factor ϵ as well as the noise magnitude η must be supplied in addition to the dataset. Three different model classes were run to illustrate the computed output: (i) the independent time step model described in Eq. (12), (ii) the spline GLM described in Eq. (18), and (iii) the hierarchical DGP prior model described in Eq. (17). From the output that the models produce, either summary statistics of the posterior (such as its mean and quantiles) can be gathered or the mean and credible intervals from all posterior time series sampled can be extracted as shown in Fig. 9.

The independent time step model clearly shows poorer performance due to the large effect of measurement error on the estimated parameters. Nevertheless, the credible interval covers the true parameter values well and is reasonably narrow. Estimation of the fraction remaining r seems to be biased toward higher values, which could be due to overlap with the variation in source isotopic signatures. The B-spline basis for the GLM seems to have default values for degrees of freedom that are fairly optimal in this case. The time series of parameter estimates is now smooth similarly to the actual parameter series. Estimation using the hierarchical DGP prior model gives the best results: the time series are adequately smooth, and estimates are close to the true values with narrow credible intervals.

4.3 Comparison of TimeFRAME and a dual-isotope mapping approach on a stationary dataset

The samples used to compare the TimeFRAME model with the traditional dual-isotope mapping approach were taken from Kenyan livestock enclosures (bomas in Kiswahili, also known as Swahili) at the Kapiti Research Station and Wildlife Conservancy of the International Livestock Research Institute (ILRI) located in the semi-arid savanna region (1°35.8^′–1°40.9^′ N, 36°6.4^′–37°10.3^′ E). The different samples represent the isotopic composition of N₂O taken from boma clusters of varying age classes (0–5 years after abandonment). At Kapiti, bomas are set up in clusters of three to four, which are used for the duration of approximately 1 month before setting up a new cluster. The sampling campaign was conducted in October 2021 in order to understand the underlying mechanisms of the huge N₂O emissions observed in these systems (Butterbach-Bahl et al., 2020), and the findings will be published in a separate paper (Fang et al., submitted). The dataset contains measurements of δ¹⁵N^bulk, δ¹⁵N^SP, and δ¹⁸O of N₂O as well as δ¹⁵N of soil. A stable isotope analysis of these samples was done using isotope ratio mass spectrometry (IRMS) as described in Verhoeven et al. (2019) and Gallarotti et al. (2021). A dual-isotope mapping approach was used to extract the production pathways and reduction extent for this dataset based on δ¹⁵N^SP and δ¹⁸O using the scenarios of mixing followed by reduction (MR) and reduction followed by mixing (RM), described in detail in Fang et al. (submitted) and Ho et al. (2023).

As the boma dataset is a set of independent measurements, TimeFRAME was applied to the data using the independent time step model. The time axis was replaced by numbering of points in the dataset. The standard deviation of all isotope measurements was set to 1 ‰ as the measurement uncertainty was not quantified; however, this may be a low estimate given the many sources of uncertainty from measurement error to international scale calibration. TimeFRAME was applied in two configurations:

1. using only δ¹⁵N^SP and δ¹⁸O as well as Ni and bD pathways to closely mimic the configuration of the dual-isotope MR approach and

2. using δ¹⁵N^bulk, δ¹⁵N^SP, and δ¹⁸O with Ni, bD, and fD (fungal denitrification) pathways, which is the most detailed configuration available with these data.

Figure 10 Comparison of TimeFRAME results with pathway estimates from a dual-isotope plot approach. TimeFRAME was applied in two configurations: (i) using only δ¹⁵N^SP and δ¹⁸O as well as Ni and bD pathways and (ii) using δ¹⁵N^bulk, δ¹⁵N^SP, and δ¹⁸O with Ni, bD, and fD pathways. The dual-isotope approach used δ¹⁵N^SP and δ¹⁸O as well as Ni and bD pathways in configurations MR (mixing followed by reduction) and RM (reduction followed by mixing). (a) A one-to-one comparison of estimates from the two TimeFRAME configurations and the dual-isotope plot MR and RM methods. (b) A plot of the estimates for each pathway from the two TimeFRAME configurations and from the dual-isotope MR method.

The agreement between the dual-isotope method and the analogous TimeFRAME two-isotope implementation is very good (mean absolute deviation of 8 % and 15 % for bD/Ni and reduction, respectively, for the MR method, Fig. 10). The dual-isotope method results are not significantly different for the MR and RM implementations, supporting the assumptions made in Eq. (3) as the basis for the TimeFRAME package: these models only deliver significantly different results in cases where N₂O reduction is very high (see Sect. 1 in the Supplement). TimeFRAME offers the major advantage that uncertainty bounds for the prior are incorporated and thus calculated for the posterior. Moreover, using the TimeFRAME package functions, the calculations are reproducible, can be performed in just two lines of code, and can be easily adapted to consider different end-members and model setups. The TimeFRAME three-isotope implementation shows very different results to the other estimates due to the inclusion of δ¹⁵N^bulk and the fungal denitrification pathway. This pathway has high δ¹⁵N^SP (Table 1) and thus strongly impacts the model estimates of nitrification and reduction, which are also evidenced by high δ¹⁵N^SP.

The pathway estimates were used to reconstruct the isotope measurements, and the RMSE between true measurements and reconstructed measurements was found as an estimate of model performance in the absence of true knowledge of pathways (Table 3). The TimeFRAME two-isotope implementation is able to reproduce the isotopic data more accurately than the dual-isotope plot due to the Bayesian optimization of the fit. The TimeFRAME three-isotope implementation shows poorer RMSE due to the additional challenge of fitting δ¹⁵N^bulk as well as the fD pathway. The difference between MR and RM implementations of the dual-isotope approach is minimal, showing that model configuration and uncertainty in end-members is far more important for results than the specific formulation of the fractionation and mixing equation.

Table 3 RMSE (‰) for isotopic data estimated using pathway contributions from the two TimeFRAME and the two dual-isotope model configurations, compared to the measured isotopic data used as input for the models.

Figure 11 Using TimeFRAME to understand the impact of water-filled pore space (WFPS) on N₂O production and consumption. (a) Comparison of TimeFRAME pathway estimates (independent, spline (M = 6 and M_r = 3), and DGP (ρ = 0.2 and ρ_r = 0.5) implementations) with estimates from the FRAME model, from an SP-δ¹⁸O mapping approach, and from a ¹⁵N-labeled gas flux approach. Comparisons are only shown for bD (bacterial denitrification) and reduction as other pathways are not estimated by the mapping and gas flux approaches. The gas flux approach does not strictly estimate bD but the proportion of N₂O arising from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ substrate. The lines show the linear regression, and the legends show the mean absolute deviation for each comparison. (b) Impact of WFPS on bD and reduction estimated using TimeFRAME spline and DGP results. The error bars show the estimated standard deviation at each point from the TimeFRAME fit. The legend shows details of the weighted linear regression (Rlm() function weighted by $\frac{\mathrm{1}}{\mathit{\sigma }}$) for each dataset.

These results show the importance of considering different pathways and model configurations. TimeFRAME users should aim to include as much isotopic data as possible and to use other complementary approaches such as microbial ecology to constrain potential production and consumption pathways, for example, to decide whether it is appropriate to include fungal denitrification. Users should consider both the estimated uncertainty for a particular model setup, provided by the TimeFRAME package, and the variation between estimates for different plausible scenarios.

4.4 Comparison of time series analysis with existing approaches

TimeFRAME was applied to two irregularly time-spaced datasets from soil incubations at different soil moisture levels (L1 = drier = 55 %–66 % WFPS; L2 = wetter = 69 %–82 % WFPS). The L1 and L2 incubations were sampled on 8 and 11 dates, respectively, with between 1 and 7 duplicate measurements taken on each sampling date and a total of 41 and 24 measurements made, respectively. The incubations are described in detail in Lewicka-Szczebak et al. (2020). TimeFRAME was compared to results from the 3DIM/FRAME model (Lewicka-Szczebak et al., 2020; Lewicki et al., 2022), with both models considering four pathways (bacterial denitrification, nitrifier denitrification, fungal denitrification, nitrification) as well as N₂O reduction using the end-members and fractionation factors reported in Lewicka-Szczebak et al. (2020). FRAME solves the isotopic equations independently for each sampling date, whereas TimeFRAME spline and DGP implementations are able to consider temporal correlations between sampling times. Additionally, dual-isotope mapping and ¹⁵N labeling approaches were compared, as described in Lewicka-Szczebak et al. (2020).

The agreement was good between pathway estimates from TimeFRAME and FRAME, although the spline implementation estimated lower reduction than other methods (Fig. 11). Agreement with the mapping approach was very poor for bD and good for reduction, reflecting the low ability of the mapping approach to unravel pathway contributions with similar end-members. Agreement with the ¹⁵N gas flux method was good for reduction and acceptable for bD, considering the denitrification contribution being quantified is not identical. The results clearly showed the influence of WFPS on bD and reduction, with both pathways increasing by 2 % for every 1 % increase in WFPS (Fig. 11).

5 Conclusions

The TimeFRAME data analysis package uses Bayesian hierarchical modeling to estimate production, mixing, and consumption pathways based on isotopic measurements. The package was particularly developed for the analysis of N₂O isotopic data and contains default isotopic end-members and fractionation factors for N₂O, but the flexible implementation means it can also be applied to other trace gases such as CH₄ and CO₂. TimeFRAME provides a simple and standardized method for analysis of pathways based on isotopic datasets, which has previously been lacking. The package will contribute strongly to reproducibility and uncertainty quantification in the analysis of these datasets.

TimeFRAME has four main classes of models, which have been extensively tested in a range of scenarios. For time series data, the Dirichlet–Gaussian process and spline GLM priors show very good results. These models are able to smooth time series to reduce the impact of noisy data and deliver good pathway quantification compared to the ground truth for simulated datasets. The independent time step model is strongly impacted by measurement noise but delivers good performance compared to the dual-isotope mapping approach, with simpler and more reproducible implementation as well as uncertainty quantification.

Model application and testing showed that uncertainty in end-members and fractionation factors was the major source of uncertainty in pathway quantification. Reduction of uncertainty in these parameters will strongly improve the insights that can be gained from isotopic data. Model setup is also critical: the sources/pathways chosen in the model strongly affect results and should be informed based on any other relevant sources of information, for example, profiling of the microbial community present at a measurement site. TimeFRAME provides a robust and powerful analysis tool, but the accuracy of results gained from TimeFRAME depends on careful definition of the model setup and configuration by the user.