An inverse method to relate organic carbon reactivity to isotope composition from serial oxidation

Serial oxidation coupled with stable carbon and radiocarbon analysis of sequentially evolved CO2 is a promising method to characterize the relationship between organic carbon (OC) chemical composition, source, and residence time in the environment. However, observed decay profiles depend on experimental conditions and oxidation pathway. It is therefore necessary to properly assess serial oxidation kinetics before utilizing decay profiles as a measure of OC reactivity. We present a regularized 5 inverse method to estimate the distribution of OC activation energy (E), a proxy for bond strength, from ramped temperature pyrolysis/oxidation (RPO). This method directly compares reactivity to isotope composition by determining the E range for OC decaying within each temperature interval over which CO2 is collected. By analyzing a decarbonated test sample at multiple masses and oven ramp rates, we show that OC decay during RPO analysis follows a superposition of parallel first-order kinetics and that resulting E distributions are independent of experimental conditions. We therefore propose the E distribution 10 as a novel proxy to describe OC reactivity and suggest that E vs. isotope relationships can provide new insight into the compositional controls on OC source and residence time. This manuscript is accompanied by an open-source Python package for performing all analyses. Copyright statement. ©Author(s) 2017. CC BY 3.0 License

remineralization exerts a major control on the global carbon cycle and on atmospheric CO 2 levels (Lasaga et al., 1985). However, OC remineralization rates are spatiotemporally heterogeneous, leading to decay timescales that range from minutes to millions of years (Boudreau and Ruddick, 1991;Forney and Rothman, 2012a;Middelburg, 1989). To explain this variability, it has been hypothesized that remineralization depends on multiple chemical and environmental factors such as OC molecular structure (Burdige, 2007;Tegelaar et al., 1989), microbial community composition (Pedler et al., 2014;Schmidt et al., 5 2011), secondary chemical interactions (Schmidt et al., 2011), and physical protection by particles (Mikutta et al., 2006;Keil and Mayer, 2014). The relative importance of these governing mechanisms remains actively debated and is thought to vary depending on environmental setting (Hedges et al., 2001;Rothman and Forney, 2007;Schmidt et al., 2011), thus limiting our mechanistic understanding of OC decay.
While bulk measurements include all OC contained within a sample, they offer no information on the distribution of chemical structure or reactivity within a complex mixture. In contrast, biomarker analysis is highly specific but individual compounds nonetheless still represent the average of multiple sources. Furthermore, biomarkers typically constitute ≤1 % of total OC and can be subject to production, transport, and preservation biases . 15 To bridge the information gained by these methods, a novel class of analytical techniques, termed "serial oxidation," has emerged. Such analyses separate carbon within a bulk sample based on its susceptibility to decomposition by chemical hydrolysis (Helfrich et al., 2007), uv light (Beaupré et al., 2007;Follett et al., 2014), heat (Rosenheim et al., 2008), or microbial respiration (Beaupré et al., 2016) and measure the stable carbon ( 13 C/ 12 C, expressed as δ 13 C) and radiocarbon ( 14 C/ 12 C, here expressed as fraction modern or Fm) content of evolved CO 2 . By separating CO 2 into multiple lability intervals, isotope ratios 20 are obtained for carbon atoms exhibiting similar physical and/or chemical properties. Because δ 13 C provides information on the source of organic matter while Fm reflects the amount of time that has passed since organic compounds were initially syntehsized, serial oxidation is a promising method to directly probe the compositional controls on OC source and residence time.
Still, a theoretical treatment of serial oxidation kinetics is lacking, hindering our ability to correlate measured isotope distri- 25 butions with intrinsic chemical properties and reactivity. In this study, we relate OC thermal recalcitrance to its corresponding δ 13 C and Fm values using ramped-temperature pyrolysis/oxidation (RPO). This method involves heating OC at a controlled rate while continuously quantifying and collecting evolved CO 2 , which is binned over user-defined time intervals (termed "fractions") and analyzed for δ 13 C and Fm (Rosenheim et al., 2008;Hemingway et al., 2017). We describe non-isothermal OC decay rates as a function of E, the Arrhenius activation energy, using a novel inverse solution to the distributed activation 30 energy model (Braun and Burnham, 1987;Burnham and Braun, 1999;Cramer, 2004;White et al., 2011). By conducting a set of kinetic experiments, we show that the E distribution within a given OC mixture does not depend on experimental conditions and is thus a reliable proxy for bond strength and OC chemical composition.
We begin in Section 3 by deriving the governing equations to describe a parallel superposition of first-order, non-isothermal decay. Then, in Section 4, we describe a method to solve for the distribution of E using a regularized inverse approach.     Table 1 for symbol definitions).
Finally, in Section 5, we determine the subset of E that is contained within each RPO fraction and directly relate OC reaction energetics to corresponding isotope values. All calculations were performed using the accompanying 'rampedpyrox' Python package (Hemingway, 2017).

5
As a representative sample, we analyzed particulate organic carbon (POC) contained in suspended sediments from the surface of the Narayani River. This sample (PB-60) was collected at the base of the Himalayas (27.70°N, 84.43°E) and has been analyzed for bulk OC and plant-wax carbon isotopes (Galy et al., 2008;Galy and Eglinton, 2011;Galy et al., 2011).
Aliquots were taken for RPO analysis from freeze-dried, archived material and acidified under HCl fumes at 60 • C for 72 h to remove carbonates (Whiteside et al., 2011). Because residual chloride has been shown to interact with the RPO catalyst wire 10 , acidified aliquots were rinsed 3× in 18.2 MΩ MilliQ water and freeze-dried overnight at −40 • C prior to analysis. For consistency and to properly calculate RPO isotope mass balance, bulk %OC, δ 13 C, and Fm values were re-measured using fumigated and rinsed material (McNichol et al., 1994a, b Tables 2-3 for values).
than that for un-rinsed aliquots (Galy et al., 2008), reflecting a minor loss of acid soluble OC for this sample during the rinsing step.
To test if the presence of inorganic carbon (IC; e.g. CaCO 3 ) affects decay kinetics, we analyzed a pure CaCO 3 laboratory working standard (Icelandic spar; Hemingway et al., 2017) as well as carbonate-rich sediment from the Southern Ocean (60.24°S , 170.19°W) collected for the Joint Global Ocean Flux Study (JGOFS; Sayles et al., 2001). JGOFS aliquots were taken from 5 archived core-top material (0 cm to 0.5 cm, stored at −80 • C), freeze-dried overnight at −40 • C, and homogenized prior to RPO analysis. IC content, OC content, and bulk 13 C composition were re-quantified at NOSAMS (McNichol et al., 1994a).

Instrumental setup
RPO analysis has been described in detail previously (Rosenheim et al., 2008;Hemingway et al., 2017). In summary, a solid sample containing ≈150 µg C to 250 µg C is loaded into a pre-combusted (850 • C, 5 h) quartz reactor and placed into a two-10 stage oven, as shown in Fig. 1. The reactor is then sealed and the sample is exposed to an atmosphere of 92:8 He:O 2 with a flow rate of 35 mL min −1 (oxidation mode). During analysis, the oven surrounding the sample is programmed to heat at a userdefined ramp rate, termed β (see Table 1 for symbol descriptions). Instantaneous temperature within the oven is measured using two thermocouples separated by ≈1 cm to monitor temperature heterogeneity, which is typically <5 • C. Following standard practice (Rosenheim et al., 2008), a ramp rate of 5 • C min −1 was used for all experiments in which CO 2 gas was collected for isotope analysis. In the second (downstream) oven, eluent gas is passed over a Cu, Pt, and Ni catalyst wire held at 800 • C to facilitate oxidation of reduced carbon-containing gases to CO 2 .
After exiting the second oven, eluent gas is distilled through a water trap and passed into a flow-through infrared gas analyzer (IRGA) to measure CO 2 concentration (in parts per million by volume; ppm CO 2 ) with 1-s temporal resolution. Resulting 5 ppm CO 2 vs. temperature plots are referred to as "thermograms" (Fig. 2). IRGA measurements were calibrated using a twopoint calibration curve before each analysis to account for instrument drift and are precise to ±5 ppm CO 2 . Downstream of the IRGA, eluent gas is passed into one of two switchable traps and CO 2 is cryogenically frozen while He and O 2 are vented to the atmosphere. Traps are switched at user-defined time points and CO 2 is further distilled, quantified, transferred into glass tubes packed with ≈100 mg CuO and ≈10 mg Ag, and flame sealed. Finally, CO 2 was recombusted at 10 525 • C for 1 h to remove trace contaminant gases.

Isotope measurement, blank correction, and data analysis
Radiocarbon compositions of all bulk samples and RPO fractions were measured at NOSAMS following standard graphitization methods (McNichol et al., 1994b). All radiocarbon results are expressed in fraction modern notation (Fm). We note that Fm used here is corrected for 13 C fractionation and is thus identical to the F 14 C notation of Reimer et al. (2004). Bulk and RPO 15 fraction stable carbon isotope compositions were measured on CO 2 gas using a dual-inlet IRMS located at NOSAMS (McNichol et al., 1994a), with resulting 13 C content expressed in δ 13 C per mille (‰) notation relative to Vienna Pee Dee Belemnite (VPDB). RPO fraction masses, δ 13 C values, and Fm values were corrected for blank carbon contribution, and δ 13 C was additionally corrected to ensure 13 C mass balance as incomplete oxidation to CO 2 has been shown to impart a small fractionation effect . Analytical uncertainty was propagated throughout all corrections. Thermograms and isotope 20 results for both Narayani POC and JGOFS sediment are plotted in Fig. 2, while temperature ranges, carbon masses, and isotope values are additionally reported in Tables 2-3.

Deriving a model of decay kinetics
We derive the distributed activation energy model by first considering the case where OC is separated into a set of discrete components with unique E values. We then generalize this description to allow for a continuous E distribution (Braun and 25 Burnham, 1987;Burnham and Braun, 1999;Cramer, 2004).

Discrete model
During OC remineralization, the decay rate of carbon contained in a particular component i is often described as as a first-order process with respect to g i (t), the mass of carbon remaining in component i at any time t (Westrich and Berner, 1984;Braun and Burnham, 1987), according to where k i is the first-order rate coefficient associated with component i. Total OC decay is then treated as the sum over all components. Although it is possible that OC decay in the natural environment additionally depends on oxidant concentration, we omit this dependency here since O 2 is provided in excess in our experimental setup (Fig. 1). In contrast to the "multi-G" and 5 "reactive continuum" models that are often used to describe environmental OC degradation rates (Westrich and Berner, 1984;Boudreau and Ruddick, 1991;Forney and Rothman, 2012a, b), here we allow k i to vary with time. Because rate coefficients are related to temperature and activation energy, k i can be written as a time-dependent function of E following the Arrhenius equation:

10
where ω is the empirically derived Arrhenius pre-exponential ("frequency") factor, R is the ideal gas constant, E i is the activation energy of carbon contained in component i, and T (t) is the measured temperature of the system at time t. For nonisothermal systems, time-dependent decay coefficients can therefore be described by the static property E i and the observed variable T (t). Although T (t) is related to t by a constant ramp rate β during RPO analysis, we leave this written as is to emphasize that our model is valid for any measured time-temperature history. Substituting Eq.
(1), we write the 15 first-order decay at time t during a non-isothermal process as The mass of carbon remaining in component i at time t can be determined by integrating Eq. (3) from an initial time t = 0: is the time integrated decay coefficient at time t and g i (0) is the initial mass of carbon contained in component i. Equation (5) states that g i (t) depends on the entire time-temperature history of the experiment. That is, the evolution of dg i (t)/dt reflects both a decrease in g i (t) as OC is remineralized and an increase in k i (t) with increasing T (t) as the experiment progresses.
This results in a predictable shift in RPO thermograms toward higher elution temperatures with increasing β (Miura and Maki,25 1998), as shown in Fig. 3.
Following Boudreau (1997) and Westrich and Berner (1984), an environmental sample containing a complex OC mixture can be described as a superposition of a finite set of n components, each decaying according to a unique k i (t) and thus corresponding to a unique E i value. G(t), the total carbon mass remaining at t, is then the sum of the mass remaining in each component: Substituting Eq. (4) into Eq. (6), this can be written as 5 We then define G 0 , the initial OC mass present in the entire sample, as the sum of initial mass contained in each component: Finally, we define p i (0), the fraction of total carbon initially contained in component i, as and note that Substituting Eq. (9) into Eq. (7) yields which describes the evolution of the fraction of initial carbon remaining at any time. The fraction of OC initially present within each component, p i (0), can be determined by fitting Eq. (11) to the observed G(t) profile measured by the RPO instrument.
While informative, this discrete description of the model suffers from two major limitations: (i) n must be set a priori or 5 determined empirically (Boudreau and Ruddick, 1991) and (ii) any noise recorded in the data will result in large uncertainty in best-fit p i (0) and E i values (Forney and Rothman, 2012b). To circumvent the first of these issues, we derive a more general description of non-isothermal first-order decay that does not assume a finite set of components with unique E i , but rather allows E to vary continuously (Boudreau, 1997;Braun and Burnham, 1987;Burnham and Braun, 1999;Cramer, 2004). The second problem is then solved using Tikhonov regularization (Section 4.2; Forney and Rothman, 2012b;Hansen, 1994).

Continuous model
In the continuous model, discrete components g i (t), κ i (t) and E i are replaced by continuous variables g(t, E), κ(t, E) and E, respectively (Table 1). Analogous to Eq.
(3), we calculate the decay of OC associated with an infinitesimal range dE about any non-negative value of E following first-order Arrhenius kinetics as 15 The mass of carbon associated with any value of E that remains unreacted at time t is then calculated by integrating Eq. (12) to obtain where g(0, E) is the initial mass of carbon associated with activation energy value E and 20 The total carbon remaining at time t can now be written as the integral of g(t, E) over all possible values of E as Substituting Eq. (13) into Eq. (15), we obtain Analogous to Eq. (9), we then define the fraction of total carbon initially associated with any value of E as Substituting Eq. (17) into Eq. (16) yields: The distribution of p(0, E) over all values of E describes the initial probability density function (pdf) of E that will lead to the observed OC decay profile when a sample is analyzed in the RPO instrument. As RPO analysis proceeds, this pdf must evolve with time to reflect the fact that some carbon has been remineralized to CO 2 . Like g(t, E), p(t, E) follows first-order Arrhenius kinetics according to where p(t, E) is the fraction of initial carbon mass that remains associated with E at time t. This can be obtained by integrating Eq. (20) from an initial time t = 0: This implies that the carbon initially remineralized to CO 2 must be associated with the lowest E values, as low E will lead to a 15 κ(t, E) term in Eq. (21) that approaches zero most rapidly. Put differently, OC that is described by higher E values will resist remineralization until more time has passed and, therefore, higher temperatures have been reached -i.e. it is more thermally recalcitrant.

First-order verification
Because our model is a specific case of n th -order non-isothermal kinetic models (Braun and Burnham, 1987;White et al., 20 2011), we must verify that carbon degradation within the RPO instrument behaves according to a superposition of parallel first-order reactions rather than higher-order processes. By replacing g(t, E) with G 0 p(t, E) on the right-hand side of Eq. (12), it can be seen that By integrating over all non-negative values of E and utilizing the definition of G(t) from Eq. (15), this can be written as The first-order model describes dG(t)/dt as a linear function of G 0 multiplied by an integral term that depends on p(t, E) but is independent of G 0 . In contrast, if carbon decomposition within the RPO instrument were to follow a higher-order process, then the relationship between dG(t)/dt and G 0 would be nonlinear and would evolve as a function of time (Follett et al., 2014).
If we define 5 then the carbon decay at time t as predicted by parallel first-order kinetics simplifies to Therefore, similar to the isothermal case (Follett et al., 2014), a superposition of parallel first-order decay reactions will result in a linear relationship between dG(t)/dt and G 0 with a zero intercept and a time-dependent slope. Thus, m(t) can be interpreted as the G 0 -normalized OC decay rate at time t. 10 We verify that OC remineralization within the RPO instrument follows parallel first-order kinetics by assessing the linearity between Narayani POC dG(t)/dt and G 0 at any time t across a range of G 0 values. For 4 arbitrarily chosen time points, this relationship is linear with an ordinary least squares R 2 ≥ 0.999, resulting in identical G 0 -normalized thermograms within analytical uncertainty (Fig. 4a-b). Thus, the decay of complex OC mixtures contained in carbonate-free samples during RPO analysis can indeed be accurately described as a superposition of parallel first-order reactions.

A note of caution on carbonates
While most RPO studies to date have isolated OC by acidifying to remove carbonates (e.g. Rosenheim et al., 2008;Rosenheim and Galy, 2012;Rosenheim et al., 2013;Schreiner et al., 2014;Bianchi et al., 2015), it has been argued that acid hydrolysis and/or dissolution of short-range-order minerals during acid treatment can alter the OC chemical bonding environment and therefore affect thermal stability (Plante et al., 2013). While analyzing samples without acid treatment can circumvent these 20 issues, the effect of carbonates on decay kinetics has not yet been considered. To test if carbonate-rich samples follow parallel first-order kinetics, we analyzed JGOFS sediment for a range of G 0 values (Fig. 4c-d). Prior to t ≈ 4500 s, when δ 13 C values of eluted CO 2 indicate a predominantly OC source (Table 3; Fig. 2b), dG(t)/dt can be accurately described as a linear function of G 0 (R 2 ≥ 0.999). However, as carbonate begins to decompose above t ≈ 4500 s, the relationship between dG(t)/dt and G 0 becomes nonlinear and the carbonate peak shifts toward higher t with increasing G 0 (Fig. 4d). 25 To investigate if non-first-order decomposition is an intrinsic property of CaCO 3 or if this is due to interactions with other materials within the sample (so-called "matrix effects"), we analyzed a pure Icelandic spar CaCO 3 labroatory standard at multiple masses (G 0 = 258 µg C, 492 µg C and 1014 µg C; β = 5 • C min −1 ; not shown). Results indicate that pure carbonate, unlike JGOFS sediment, does follow first-order kinetics with a maximum decomposition rate occurring at (700 ± 6) • C independent of G 0 . Interaction with reduced organic carbon, corresponding hetero-atoms (e.g. N, P, S), or trace metals contained within  samples. Thus, while avoiding the issues of acid treatment, the presence of carbonate will result in thermograms that cannot be accurately described by the model presented here, and we therefore argue in favor of acid treatment when using the RPO instrument to determine reaction energetics of carbonate-containing samples. In contrast to previous solutions (Braun and Burnham, 1987;Burnham and Braun, 1999;Cramer, 2004), this approach does not require an a priori assumption about the form of p(0, E) (e.g. Gaussian). Because this problem is sensitive to noise at the level of our analytical uncertainty (Forney and Rothman, 2012b), we seek a smooth solution using Tikhonov regularization 5 (Section 4.2; Forney and Rothman, 2012b;Hansen, 1994).
To numerically calculate p(0, E), we discretize the continuous variable t over the time course of the experiment into a vector t containing n t nodes such that For j = 1 and j = n t , t j−1 and t j+1 in Eq. (26) are, respectively, replaced by t j since t is undefined outside of this range.

10
For generality, and because the distributed activation energy model is frequently applied over geologic timescales with nonuniformly distributed time measurements (Braun and Burnham, 1987;Burnham and Braun, 1999;Cramer, 2004), Eq. (26) does not require a uniform time step (i.e. it is possible that ∆t j = ∆t i =j ). Similarly, we generate a vector E containing n E nodes over the range values considered for the model solution such that 15 noting that E is uniformly spaced since this vector is not constrained by observations. We constrain E to be within 50 kJ mol −1 to 350 kJ mol −1 based on published biomass and petroleum E ranges (Braun and Burnham, 1987;Burnham and Braun, 1999;Cramer, 2004;White et al., 2011 ∆t u ∆E, j = 1, . . . , n t , l = 1, . . . , n E . 25 The A matrix is often termed the model "design matrix." We then calculate the fraction of initial carbon remaining at each time point as where α(t) is the G 0 -normalized, integrated RPO thermogram at time t. We generate a discretized vector g by interpolating G(t)/G 0 onto each t j in t (j = 1, . . . , n t ). Our model can then be written in matrix form as where p is an unknown, discretized vector of p(0, E) with length n E such that 5 While Eq. (30) can be solved by multiplying g by the computed inverse of A, if g contains noisy data this may result in negative values of p l that are mathematically possible but physically unreasonable (Forney and Rothman, 2012b). Here, we find the solution that satisfies subject to the constraints 10 p l ≥ 0, l = 1, . . . , n E , and n E l=1 p l = 1, l = 1, . . . , n E .
Eqs. (32)-(34) describe the model solution that minimizes the norm of the residual error (i.e. the root mean square error, or RMSE) while fulfilling the constraints that p is non-negative and sums to unity.

Choice of frequency factor
In order to construct the A matrix and solve for p, our method requires that the Arrhenius frequency factor ω is prescribed a priori. There exists significant discussion in the literature on the best choice of ω, as multiple values can describe laboratory results equally well but will result in drastically different predictions of OC degradation rates over geologic timescales (Braun and Burnham, 1987;Dieckmann, 2005). Furthermore, it has been argued that ω represents a variable change in entropy asso-20 ciated with the decay of specific organic compounds and should therefore be parameterized as a function of E (the so-called "kinetic compensation effect" or KCE; Dieckmann, 2005;Lakshmanan et al., 1991;Tang et al., 2000). For example, a linear ω increase with E from ≈10 8 s −1 (E = 175 kJ mol −1 ) to ≈10 26 s −1 (E = 400 kJ mol −1 ) has been utilized to better predict petroleum formation rates (Dieckmann, 2005). To circumvent the issue of multiplicity, and to account for the KCE, Miura and Maki (1998) developed a method to estimate the best-fit ω for each E value by comparing the shift in elution temperatures 25 when a sample is analyzed at multiple ramp rates. However, because this approach is based on large extrapolations in 1/T vs.
β/T 2 space, it is highly sensitive to noise in temperature and β measurements (Burnham and Braun, 1989). To select a best-fit ω, we construct the A matrix for a range of ω values and determine the residual error norm between measured G(t)/G 0 and that predicted by the resulting p vector determined by Eqs. (32)-(34). We consider the KCE by calculating ω as a function of E according to Resulting residual errors for Narayani POC using a range of KCE slopes and intercepts are shown in Fig. 5 (β = 5 • C min −1 , 5 E ranging from 50 kJ mol −1 to 350 kJ mol −1 ). By setting an "acceptable" residual error norm cutoff of ≤10 −4 , it can be seen that there exist multiple KCE slope and intercept combinations that can equally fit the observed data. Additionally, we estimate the best-fit ω using a range of ramp rates (β = 2 • C min −1 , 5 • C min −1 and 10 • C min −1 ) following the method of Miura and Maki (1998) (Fig. 5, white circle). While this estimate falls outside of the cutoff range, likely due to noise in temperature and β measurements, it results in a KCE slope near zero and suggests that ω is constant during RPO oxidation of this sample. To 10 accurately compare RPO results between samples, we therefore select a constant ω value of 10 10 s −1 , well within the cutoff range, for all samples analyzed herein (Fig. 5, red star). While a different choice of ω will shift p(0, E) to higher or lower absolute values of E, we emphasize that it will not affect the distribution of E and that only relative changes in E between RPO fractions should be interpreted.
For example, a shift in ω from a constant value of 10 7 s −1 to 10 12 s −1 results in an increase in the mean of the pdf of E, 15 termed E and calculated as from 150 kJ mol −1 to 224 kJ mol −1 for Narayani POC. However, the relative standard deviation of the pdf of E, calculated as 20 remains constant at 20%. A higher ω value therefore results in a broader p(0, E) distribution that is centered at a higher mean E value but has no effect on the shape of the distribution.

Tikhonov regularization
In principle, after choosing ω and constructing the A matrix, the pdf of E that best describes an RPO thermogram can be determined by solving Eqs. (32)-(34). However, the inverse solution is sensitive to noise at the level of RPO instrument 25 precision (±5 ppm CO 2 , ±5 • C; Hemingway et al., 2017), and is therefore ill-posed (Hansen, 1994;Lakshmanan et al., 1991).
We address this sensitivity to data uncertainty using Tikhonov regularization (Hansen, 1994;Forney and Rothman, 2012a, b).
This approach finds an optimal solution that minimizes p(0, E) complexity (as determined by the intensity of fluctuations, or "roughness") while maximizing solution accuracy. Following Forney and Rothman (2012b), we calculate roughness as the residual error norm, ||g -A·p|| first derivative of the solution vector: where R is the bi-diagonal first derivative operator with an additional first row equal to [1 0] and last row equal to [0 − 1] to account for p being equal to zero outside of the range E min < E < E max . The regularized inverse solution can then be determined by including this roughness term when solving the constrained least squares: where λ is a scalar that determines how much to weight the roughness R · p relative to the residual error g − A · p . The best choice of λ is considered to be the value that optimizes this balance. As described in Hansen (1994), a common approach is to choose the value corresponding to the point of maximum curvature in a log − log plot of residual error and roughness while allowing λ to range over many orders of magnitude (i.e. the so-called "L-curve"). From this point, increasing λ strongly 10 increases residual error but has little effect on solution roughness, while decreasing λ greatly increases roughness but has little 10 -5 10 -4 10 -3 10 -2 10 -1 residual error norm, ||g -A·p|| effect on residual error. Thus, here we choose the λ value corresponding to the corner of the L-curve for each sample, as exemplified in Fig. 6.

p(0, E): A novel proxy for chemical bond strength
In order to interpret p(0, E) as an intrinsic property of OC contained within a sample, we must show that results do not depend on experimental conditions such as ramp rate β and initial carbon mass G 0 . To test this, we analyzed Narayani POC using a 5 range of masses (G 0 = 268 µg C, 533 µg C and 828 µg C) and ramp rates (β = 2 • C min −1 , 5 • C min −1 and 10 • C min −1 ). Fig.   7 shows that the regularized pdfs of E are nearly identical across all experimental conditions. This supports the hypothesis that the p(0, E) distribution is an intrinsic property of a given sample when exposed to a particular oxidation pathway. Although there exist small differences between individual analyses due to measurement uncertainty and variability in best-fit λ values (ranging from 0.044 to 0.448, n = 5), the main features of p(0, E) are robust across all conditions. 10 We propose p(0, E) as a novel proxy to describe the distribution of carbon bond strength (Braun and Burnham, 1987;Burnham and Braun, 1999;Cramer, 2004). For example, Narayani POC is known to integrate recently fixed biomass, pre-aged soils, and eroded rock-derived material (Galy et al., 2008(Galy et al., , 2011Galy and Eglinton, 2011;Rosenheim and Galy, 2012). Such integration should lead to large chemical diversity and a broad, complex E distribution, as is observed (Fig. 7). Furthermore, slow environmental turnover has been shown to diversify the distribution of chemical bonds due to a combination of microbial 15 alterations (Schmidt et al., 2011), OC aggregation (Keil and Mayer, 2014), and stabilization by mineral surfaces (Keil and  Mayer, 2014; Mikutta et al., 2006). Thus, OC reactivity within the RPO instrument and the resulting E distribution likely reflects both the strength of covalent bonds between carbon atoms as well as interactions with mineral surfaces (Keil and Mayer, 2014;Mikutta et al., 2006). We therefore propose that combining p(0, E) with serial oxidation isotope measurements is an ideal method to test the effects of mineral interactions and selective preservation on OC turnover time.
5 Relating E and isotope composition 5

Determining the distribution of E within each RPO fraction
To relate p(0, E) distributions to RPO isotope measurements, we calculate the subset of the pdf of E that is contained within each RPO fraction. Because we can predict the evolution of p(t, E) at any time t following Eq. (21), this can be calculated as where n f is the number of RPO fractions collected for a given sample, Π f (E) is the subset of p(0, E) contained in RPO along the y axis for visual clarity. Π f (E) distributions do not follow any predictable functional form and are highly overlapping due to the fact that OC associated with a given E value decays over a wide time interval (Cramer, 2004).
we calculate the mean E value contained in each RPO fraction as and the standard deviation of E contained in each RPO fraction as σ f , where Resulting E f and σ f values are reported in Table 2. It can be seen in Fig. 8 that Π f (E) distributions do not follow any particular 5 form (e.g. Gaussian) and are highly overlapping, reflecting the fact that the CO 2 isotope composition for each RPO fraction is itself a weighted average of multiple sources.

Kinetic isotope fractionation
While not necessary for Fm because it is fractionation-corrected by definition (Reimer et al., 2004), it is important to correct for any kinetic isotope effects occurring within the RPO instrument before interpreting δ 13 C as a carbon source tracer (Hemingway 10 et al., 2017). If kinetic fractionation is large, as has been observed both during thermogenic methane formation (Tang et al., Cramer, 2004) and dissolved OC oxidation by uv light (Oba and Naraoka, 2008), then this effect could overprint carbon source δ 13 C signals. However, when directly measured using single-compound standards, Hemingway et al. (2017) concluded that 13 C fractionation within the RPO instrument must be ≤2 ‰. Still, we correct the measured δ 13 C values of each RPO fraction using the ratio of carbon-normalized 13 C and 12 C decomposition rates at each time point according to 13/12 where we have added a preceding 12 or 13 superscript to specify isotope-specific variables. Following the Arrhenius equation, 13/12 r(t) can be written as a function of the difference in E between 13 C-and 12 C-containing molecules: Although 13 -12 ∆E is likely not identical for all compounds due to differences in the entropy and enthalpy of isotope substitution (Tang et al., 2000), the estimated range of values for RPO analysis is small (0.3 × 10 −3 kJ mol −1 to 1.8 × 10 −3 kJ mol −1 ; 10 Hemingway et al., 2017). We therefore assume a 13 -12 ∆E value of 1.8 × 10 −3 kJ mol −1 for all RPO fractions, noting that a choice of 0.3 × 10 −3 kJ mol −1 would result in δ 13 C values that are identical to those calculated here within analytical uncertainty.
13 C-containing molecules decay at rates governed by a pdf of E that is identical to p(0, E) but has been shifted by 1.8 × 10 −3 20 kJ mol −1 . We then correct the measured δ 13 C values of each RPO fraction f according to where 13/12 r(t) av f is the average 13/12 r(t) f value over the time of collection for RPO fraction f . For the samples analyzed here, 13/12 r(t) is initially ≈0.999, indicating slightly faster 12 C decay at low temperatures, and gradually increases to ≈1.002 when 25 G(t) 0.01G 0 , as has been described previously (Cramer, 2004;Hemingway et al., 2017). Resulting kinetic fractionation corrections are near or within analytical uncertainty, with absolute δ 13 C values for all RPO fractions shifted by <0.2 ‰.  Galy et al., 2011;Galy and Eglinton, 2011). Note that plant-wax fatty acids are known to contain less 13 C (lower δ 13 C values) than corresponding biospheric OC. Each point is plotted at E = E f . Error bars in E are equal to σ f , while δ 13 C and Fm analytical uncertainty is always smaller than point marker and is therefore not shown.

Comparing E to δ 13 C and Fm
Finally, we describe a framework to directly relate OC reactivity and isotope distributions by plotting E f for each RPO fraction vs. the corresponding measured δ 13 C and Fm values ( acid isotope values (Galy et al., 2011;Galy and Eglinton, 2011), are shown in Fig. 9. Within this framework, it can be seen that Narayani POC must contain at least two end members with drastically different isotope compositions and unique yet overlapping E distributions. Previous studies have shown that ≈(20 ± 5) % of OC contained in this sample is derived from the erosion of carbon-rich bedrock (Galy et al., 2008;Rosenheim and Galy, 2012). Rock-derived OC is the likeliest source of high-E, low-Fm material, as this end member is 14 C-free by definition. Plant-wax FA δ 13 C and Fm values are similar to those 5 for low-E RPO fractions (Fig. 9), suggesting that vascular plant OC is the source of low-E material. Narayani POC isotope trends are thus consistent with predominantly biospheric carbon below ≈150 kJ mol −1 , a mixed region from ≈150 kJ mol −1 to ≈200 kJ mol −1 , and exclusively rock-derived OC above ≈200 kJ mol −1 . This result exemplifies the utility of RPO E vs.
isotope relationships to directly relate the distribution of OC sources, environmental turnover times, and chemical bonding environments.

Conclusions
In this study, we present a regularized, inverse method to determine the distribution of E, a measure of OC reactivity, when natural organic matter is exposed to serial oxidation. We show that OC decay follows parallel, first-order kinetics. In contrast, the kinetics of carbonate oxidation cannot be constrained due to matrix effects. We propose that p(0, E), the distribution of E contained within a sample, is a useful proxy to describe the range of OC bonding environments. Importantly, our method 15 does not require a priori assumptions about the distributional form of p(0, E). Finally, we determine the subset of E contained within each RPO fraction in order to directly relate reaction energetics with the distribution of carbon isotope compositions within a complex OC mixture. We suggest that E vs. isotope relationships can provide new insight into understanding the compositional controls on OC source and residence time. This manuscript is accompanied by an open-source Python package for performing all analyses. *δ 13 C f is additionally corrected following Hemingway et al. (2017) to ensure that the mass-weighted mean matches the measured bulk value.
**Assuming L-curve best-fit λ value and ω = 10 × 10 10 s −1 . *δ 13 C f is additionally corrected following Hemingway et al. (2017) to ensure that the mass-weighted mean matches the measured bulk value.