Inverse method for estimating respiration rates from decay time series

. Long-term organic matter decomposition experiments typically measure the mass lost from decaying organic matter as a function of time. These experiments can provide information about the dynamics of carbon dioxide input to the atmosphere and controls on natural respiration processes. Decay slows down with time, suggesting that organic matter is composed of components (pools) with varied lability. Yet it is unclear how the appropriate rates, sizes, and number of pools vary with organic matter type, climate, and ecosystem. To better understand these relations, it is neces-sary to properly extract the decay rates from decomposition data. Here we present a regularized inverse method to identify an optimally-ﬁtting distribution of decay rates associated with a decay time series. We motivate our study by ﬁrst evaluating a standard, direct inversion of the data. The direct inversion identiﬁes a discrete distribution of decay rates, where mass is concentrated in just a small number of discrete pools. It is consistent with identifying the best ﬁtting “multi-pool” model, without prior assumption of the number of pools. However we ﬁnd these multi-pool solutions are not robust to noise and are over-parametrized. We therefore in-troduce a method of regularized inversion, which identiﬁes the solution which best ﬁts the data but not the noise. This method shows that the data are described by a continuous distribution of rates, which we ﬁnd is well approximated by a lognormal distribution, and consistent with the idea that decomposition results from a continuum of processes at different rates. The ubiquity of the lognormal distribution suggest that decay may be simply described by just two parameters: a mean and a variance of log rates. We conclude by describ-ing a procedure that estimates these two lognormal parameters from decay data. Matlab codes for all numerical methods and procedures are provided.


Introduction
Over 10 7 different types of organic substances (Wackett, 2006) comprise the roughly 1800 Gt of carbon in soils (Denman et al., 2007) and biomass, and 750 Gt of dissolved and particulate organic carbon in the oceans (Benner and Herndl, 2011) and marine sediments (Denman et al., 2007). Diverse decomposer communities (up to 6000 bacterial species per gram soil and up to 11 000 bacterial species per gram marine sediments, Horner-Devine et al., 2004) respire these compounds, converting 75 Gt of terrestrial carbon (Schlesinger and Andrews, 2000) and 50 Gt of oceanic carbon to CO 2 every year (Denman et al., 2007). During degradation, organic tissues are broken down to particulate or dissolved organic matter (Eijsackers and Zehnder, 1990), which are then processed microbially to carbon dioxide or converted to microbial compounds and by-products (Madigan et al., 2005) that are subsequently transformed and broken down again (Agren and Bosatta, 1998), eventually resulting in a complete conversion of organic carbon to carbon dioxide. Quantitatively estimating the rates of these processes is difficult for many reasons: the variety of components initially found in tissues vary in lability (Tenney and Waksman, 1929;Burdige, 2006;Lutzow et al., 2006;Minderman, 1968;Berg and McClaugherty, 2007); as decomposition proceeds, organic 3602 D. C. Forney and D. H. Rothman: Inverse method for estimating respiration rates molecules interact chemically (Berg and Laskowski, 2006;Lee et al., 2004;Lutzow et al., 2006), forming humus and other hard-to-degrade compounds (Berg and McClaugherty, 2007;Paul, 2007;Eijsackers and Zehnder, 1990); microbial processes produce compounds which vary in lability (Andrén and Paustian, 1987;Agren and Bosatta, 1998); and particulate and dissolved carbon bond and sorb to clays and minerals (Oades, 1988;Hedges and Oades, 1997;Mayer, 1994;Vetter et al., 1998;Nieder and Benbi, 2008), forming organomineral complexes that also affect decomposability.
This wide range of decomposition processes and mechanisms lead to heterogeneous kinetics, with rates of decay ranging from weeks to thousands of years (Janssen, 1984;Yang and Janssen, 2000;Trumbore, 2000) or greater (Middelburg, 1989). Knowledge of rate heterogeneity is important because it tells us how long carbon resides in organic matter. Rates of decay can be related to other dynamic properties such as the turnover times and ages of soil carbon (Jenkinson et al., 1990;Forney and Rothman, 2012;Bolin and Rodhe, 1973;Rodhe, 1992;Feng, 2009b). However, it is difficult to quantitatively model decomposition because we lack fundamental constitutive relations between rates, mechanisms, composition, and environment. Our primitive understanding of decomposition dynamics is evident in state-of-the art ecosystem models (Cox, 2001;Moorcroft et al., 2001;Medvigy et al., 2009;Sitch et al., 2003;Krinner, 2005). The treatment of primary production in those models is more sophisticated and mechanistic than the treatment of organic matter decomposition. Empirically identifying constitutive relations requires (1) mathematical relations underlying the dynamics of decay and (2) estimating dynamic parameters of the model from decomposition data. This paper focuses on the latter problem. Because decomposition dynamics involves multiple time-scales and is highly complex, heuristic models are used. Because degradation can take place over such a wide range of timescales, decay experiments sample only a portion of the decomposition history. As a compromise, we investigate plant litter decay phases and early transformations to young soil organic matter by analyzing a long-term litter decay study that spans up to 10 yr.
Current degradation models differ in the way they account for kinetic heterogeneity. They vary in the number of rate pools, mass flow partitioning, and complexity ). Models of organic matter decomposition can be classified in terms of three types. The first and most commonly applied model is the multi-pool model (Adair et al., 2008;Harmon et al., 2009;Currie et al., 2010), otherwise known as multi-compartment (Jenkinson, 1977;Nieder andBenbi, 2008), multi-component (van Keulen, 2001;Andrén and Paustian, 1987), or multi-G (Berner, 1980) models. In these models organic matter is partitioned into one (Olson, 1963) or more (Minderman, 1968;Berner, 1980;Jenkinson, 1977;Harmon et al., 2009;Currie et al., 2010) pools, decaying exponentially at different rates. Pools are suggested to be associated with different compounds present in plant matter (Minderman, 1968). However, the number of rates present, amount of material at each rate, and the value of each rate may vary across different types of organic matter, ecosystems and climates, and the relations between these parameters are not well understood (Adair et al., 2008;Feng, 2009b). The second class of models are continuous parallel models, also called reactive continuum (Boudreau and Ruddick, 1991) or disordered kinetic models (Forney and Rothman, 2012). In these models, organic matter is described as a continuum of qualities, and decomposition proceeds by a continuous distribution of exponential decay rates (Boudreau and Ruddick, 1991;Feng, 2009b). These have only more recently been applied Rothman and Forney, 2007;Feng, 2009a;Forney and Rothman, 2012). The third and most detailed type of models are transformational models, which incorporate transformations to other types of soil organic matter, decomposer biomass, humus, etc. These models can be discrete, consisting of a network of pools exchanging carbon with one another (van Veen and Paul, 1981;Parton et al., 1987Parton et al., , 1993Eijsackers and Zehnder, 1990;Beare et al., 1992), or continuous (Carpenter, 1981;Bosatta, 1985;Agren and Bosatta, 1998). Some transformational models also couple the rates of organic carbon decomposition with nutrient dynamics (Manzoni and Porporato, 2007;Agren and Bosatta, 1998) and may be nonlinear Porporato, 2007, 2009). Transformations in organic matter can also be described by a first-order decay constant which decreases according to a function of time (Janssen, 1984;Yang and Janssen, 2000;Bosatta, 1995).
The three types of models are similar. The multi-pool model can be written in terms of the disordered kinetic model when the continuous rate distribution is comprised of delta functions (Forney and Rothman, 2012). The disordered kinetic model can be represented by a transformational model when quality is preserved during transformation (Bosatta, 1995). When linear, transformational network models can sometimes be mathematically represented by a multi-pool model (Bolker et al., 1998). This simplification also occurs when the timescale of transformation is short with respect to the decay timescale of slowly degrading products (Forney and Rothman, 2012).
Using the disordered kinetic approach, we have recently found that plant matter decay is well described by a lognormal distribution of decay rates (Forney and Rothman, 2012). We also determined how the rate distribution and the two parameters of the lognormal are qualitatively related to climatic and compositional variables (Forney and Rothman, 2012).
In this paper, we elaborate on the formation and implementation of the regularized inversion method that we employed to identify the best fitting rate distribution (Forney and Rothman, 2012). To do so, we first consider a simpler direct inversion of the data, which forms the basis for the regularization method. The direct inversion tends to find discrete rate distributions with mass contained at only a few rates. In this sense, the direct inversion is a new procedure for finding the best fitting multi-pool solution without needing to assume how many pools are present. We use the resulting multi-pool solutions to motivate the use of a regularized inverse. We then describe the regularization method in detail and provide an example of its implementation. Because the regularization method identifies a common lognormal pattern (Forney and Rothman, 2012), we also present a simple procedure for estimating the lognormal parameters from the decay data. This procedure provides a more precise estimate of the best fitting lognormal, resulting in slightly better fits to the data than the regularized solution. The remainder of the paper proceeds as follows: In Sect. 2 we discretize the disordered kinetic model for numerical manipulation. In Sect. 3 we calculate the model inverse. This basic inversion easily and rapidly identifies the appropriate multi-pool solution without prior assumption regarding the number and type of pools. Later in that section we show that the best fitting multi-pool solution is sensitive to noise in the data as both the number and type of best fitting pools fluctuate with the noise levels found in the data. Thus, in Sect. 4 we address the noise sensitivity by employing a method of regularization to invert the data and provide details of its implementation. Finally, in Sect. 5, we analyze litter decay data by assuming rates are distributed lognormally and fit just the two lognormal parameters to the data. Matlab codes for all numerical procedures and approaches are provided online in the Supplement.

The parallel model of decay
As organic matter is composed of many different compounds, some are more resistant to degradation and break down more slowly than others (Tenney and Waksman, 1929;Burdige, 2006;Lutzow et al., 2006;Minderman, 1968;Berg and McClaugherty, 2007). Substrate heterogeneity suggests that degradation proceeds at different rates in parallel. This heterogeneity can result in either a discrete or a continuous distribution of decay rates.
We identify the best fitting rate distribution by inverting a continuous parallel model of decay. In this model, the fraction g(t) of original mass remaining is described by a continuous superposition of exponential decays (Boudreau and Ruddick, 1991), where p(k) is the probability distribution of components with rate k and ∞ 0 p(k)dk = 1. We have previously found that litter decay rate distributions are well characterized by a lognormal p(k), although other possible forms have been hypothesized (Boudreau and Ruddick, 1991;Bolker et al., 1998;Rothman and Forney, 2007;Feng, 2009b).
Mathematically, Eq. (1) is the Laplace transform of p(k). We obtain the distribution p(k) by taking the inverse Laplace transform of the data g(t) (Forney and Rothman, 2012). To do so, we transform the model to facilitate calculations. We first make a change of variables from k to ln k as we expect to find a wide range of decay rates. Because probability is conserved, p(k)dk = ρ(ln k)d ln k, and Eq. (1) becomes where ρ(ln k) is the probability distribution in ln k space.
(2) can then be written in matrix form as g = Aρ. (3) g is the vector of predicted time series data points g i = g(t i ) having length m. ρ is a vector of length n, representing the average value of ρ(ln k) over a discrete interval ln k. A is an m × n matrix representing the discrete Laplace transform operator. See Appendix A1 for more details. We proceed to identify the underlying rate distribution ρ(ln k) from an observed decay time seriesĝ(t) by inverting the model (Eq. 3) and solving for ρ.

Direct inversion with constraints
Our approach in this section directly calculates the inverse Laplace transform of data with non-negative constraints. This approach identifies a handful of separate and distinct pools which are present during decay, having values ρ j > 0; the remaining pools are not present and have values ρ j = 0. This technique therefore provides a direct estimate of the best fitting multi-pool model, expressed as a discrete distribution ρ(ln k). We find however that the multi-pool solution is very sensitive to noise in the data and is overparametrized. We therefore proceed in Sect. 4 to refine this approach using a regularization technique to invert the data (Forney and Rothman, 2012), providing a continuous rate distribution ρ(ln k) which is both simpler and less sensitive to noise.

Calculating the constrained direct inverse
The distribution ρ can be directly calculated by computing the inverse of (Eq. 3) from the measured decay dataĝ which includes the first data pointĝ 1 = 1 at t 1 = 0, Because solutions to Eq. (4) fit the noisy data exactly, negative values of ρ j are possible. In order to find solutions with non-negative ρ i , we instead solve the constrained least squares problem such that where x ≡ x 2 i is the norm of vector x. The elements A 1j corresponds to t 1 = 0 and constraint Eq. (7) ensures that ρ is a probability distribution, which sums to one. We use the Matlab function lsqnonneg.m to calculate the nonnegative solution to the least squares problem (5). Constraint (7) is met by weighing the first data point at g(0) = 1 more heavily than the others. Details of the solver can be found in Appendix A2.
We apply this method to data from the Long-Term Intersite Decomposition Experiment Team (LIDET) study (Gholz et al., 2000;Harmon et al., 2009;Adair et al., 2008;Currie et al., 2010;Harmon, 2007). The LIDET study monitored the decomposition of 27 different types of litter, including needles, leaves, roots, wood, grass, and wheat distributed amongst 28 different locations across North America ranging from Alaskan tundra to Panamanian rainforests. Litter was collected and then re-distributed in litter bags at different sites in order to investigate the effect of composition, ecosystem, and climatic parameters on decomposition. Litter bags were collected and analyzed each year for up to ten years, with four replicates for each site, litter type, and removal time. We call a data point the average mass fraction remaining of all replicates of a given plant matter type, site, and duration. A data set is the time series of all data points associated with a particular combination of litter type and site. Bags at tropical and sub-tropical sites were more frequently collected at three to six month intervals.
An example of a decay dataset is shown in Fig. 1a. The rate distribution ρ corresponding to the solution of Eqs. (5-7) is shown in Fig. 1b. Three pools are associated with the decay shown in Fig. 1b: a very rapid pool, a moderately labile pool, and an extremely slow pool. However, by varying the search domain of decay rates k min and k max we find that the slow pool always takes the smallest value, k min , which in this case is 10 −6 yr. This suggests that the pool at k min represents an inert or constant mass fraction. The fast pool on the other hand is located at k = 38.5 yr −1 (1.4 wk −1 ). The amount remaining of this fast pool at the first measurement t = 1 yr is e −38.5 = 1.9 × 10 −17 , which is just past the double precision limit of 16 significant figures. Therefore this pool is numerically zero for every measured data point, and is therefore indistinguishable from an instantaneous decay. The inversion therefore suggests that this dataset contains three types of pools: a rapid, instantaneous pool; an exponentially decaying, active pool; and a constant, inert pool. Because the inert and instantaneous pools do not have rates associated with them, this dataset is described by three parameters.
Models with various combinations of active, inert and instantaneous pools have been previously suggested to describe litter decay (Harmon et al., 2009;Nieder and Benbi, 2008). In 1945, Henin and Dupuis (Nieder and Benbi, 2008) suggested using an inert pool to represent the transformation of incoming carbon to stabilized soil carbon. However, because the turnover time is the mean inverse rate k −1 of each individual pool (Feng, 2009b;Forney and Rothman, 2012), decomposition models containing an inert pool have an infinite turnover time and result in the unphysical situation where soil organic matter grows indefinitely. These models require additional parameters and heuristics to calculate an effective rate of turnover (Currie et al., 2010;Harmon et al., 2009).
The direct inversion technique provides an approach for identifying the simplest multi-pool solution which best fits the data. It generally provides solutions consistent with fitting different multi-pool models which contain various types of pools and choosing the model having the least error. This approach therefore gives a simple, direct estimate of the appropriate multi-pool solution and answers the question "how many pools?". The multi-pool approach however exhibits some serious shortcomings, which we discuss in the following sections.

Analysis of the multi-pool solutions
We use the same inversion technique to find the best multipool solution associated with 191 datasets chosen (Forney and Rothman, 2012) from the LIDET study. We then determine whether pools are active, inert, or rapid as described above. If any pool contains a mass fraction of less than 1 %, we neglect that pool as it requires additional model parametrization but has a negligible effect on the residual error. We distribute the mass of the small pool proportionately amongst the remaining pools, as these pools may be Table 1. This table shows the wide and seemingly unpredictable variation in the types of pool models that best fit 191 LIDET datasets. Each model was obtained by solving the non-negative least squares problem (5) and analyzing the resulting ρ to determine whether the pools were active (decaying exponentially), inert, or rapid, as discussed in the text. The first column lists the model type, the second column indicates kinds of pools present in that model, the third column provides the formula for the model, the fourth column gives the total number of pools in the model, the fifth column indicates the total number of parameters in the model, and the sixth column indicates how many LIDET datasets are best fit by that model. The sizes of the constant, rapid, and active pools in the models are c, r, and a, respectively. The size of one of the pools is determined by the others since mass fractions must sum to one.

Model
Para-# Datatype Type of pools Formula Pools meters sets artifacts due to noise. 1 Table 1 shows the results of the inversion and indicates that eight different types of multi-pool models are needed to describe all 191 of these LIDET data sets. 17 % of the datasets are described by one pool, 56 % are described by two pools, and 27 % are described by three or more pools. The maximum number of pools found was four and all datasets are characterized by 5 or less parameters. The data from Fig. 1a are of type 4 in Table 1; 24 other LIDET datasets also exhibit the same behavior. This information can be used to determine some patterns underlying litter decomposition. For example, we find that 29 % of all datasets involving either needles or wood appear to decay exponentially, while roots on the other hand decay exponentially only for 2 % of their datasets. Roots tend to have 3 or 4 pools for 48 % of their datasets, while leaves, needles, and wood each have 3 or 4 pools less than 25 % of the time. However, extracting further information from the inversion results in order to determine which model and parameters are appropriate for different combinations of tissue types, substrates, or environments is difficult if not impossible since there is too much variation in ρ and model type. Indeed this is one of the shortcomings of the multi-pool model. It is too overparametrized to be useful in identifying fundamental constitutive relations between composition and dynamics.

Sensitivity to noise
In this section we test the direct inversion method and evaluate the appropriateness of multi-pool models for describing noisy decay data. The standard shortcoming of a multipool solution is its inherent sensitivity to noise in the datâ g (Kroeker and Henkelman, 1986;Kleinberg, 1996), as models with differing numbers of pools can well approximate one another to a high degree of accuracy (Yeramian and Claverie, 1987). We proceed to investigate the solution's sensitivity to noise levels present in litter decay data. If we assume that the multi-pool model is correct, then the line in Fig. 1a is the true decay g(t). We then estimate the variance of the data noise σ 2 by calculating its maximum likelihood σ 2 = iĝ 2 i /m (Konishi and Kitagawa, 2010), assuming the noise is Gaussian with zero mean. We then generate numerical trials of decay data by adding Gaussian noise of the same strength, in this case σ = 0.024, to the predicted decay g(t). Then, for each trial we recalculate the associated inverse solution ρ. The results of the inversion of 60 000 numerical trials is shown in Fig. 2a-g. There is wide variation in the active pools just due to the noise in the data. Figure 2b shows how the model type varies due to noise. The "true" model, type 4, is only identified 39 % of the time. The other LIDET datasets show similar results. These results highlight why it is difficult to gain fundamental understanding of decomposition processes from fitting a multi-pool model to the data.

Regularized inversion method
Because the inverse Laplace transform Eq. (4) is sensitive to noise in the data, it is ill-posed (Hansen, 1987(Hansen, , 1994. In the previous section, we have demonstrated that the inverse Laplace transform with non-negative constraints Eqs. (5-7) is also ill-posed (Yeramian and Claverie, 1987) at the level of noise present in our decay functions g(t). There are a number of approaches to solve the inverse Laplace transform robustly without being sensitive to small changes in the data. These include Laplace-Padé approximation (Yeramian and Claverie, 1987), the phase function method (Zhou and Zhuang, 2006), regularization (Hansen, 1994(Hansen, , 1987Istratov and Vyvenko, 1999;Press et al., 1992)

Fig. 2.
Distributions ρ(ln k) for 60 000 trials of adding random noise to the predicted decay in Fig. 1a. For each trial, ρ(ln k) is estimated using the method discussed in Sect. 3.1. The pools in each solution are then analyzed to determine the model type as described in the text and Table 1, except for model type 9 which is a three active-pool model and type 10 which is any combination of pools not described by models 1-9. (A-G) shows the mean ρ(k) for each model type, calculated by summing all of the observed distributions for each type and dividing by the total number of times each type was observed. (H) Fraction of trials that each model type was identified.
others (Istratov and Vyvenko, 1999). We choose regularization because of its straightforward implementation. Regularization is commonly applied elsewhere, for example in the analysis of NMR spin relaxation (Lamanna, 2005) in biological tissue (Kroeker and Henkelman, 1986) and in porous media (Gallegos and Smith, 1988;Kleinberg, 1996). Like most other robust inverse Laplace transform methods, regularization works by determining the solution that, in principle, fits the data without fitting the noise. This solution represents a tradeoff between simplicity and accuracy. Here we choose a specific type of regularization called Tikhonov regularization (Hansen, 1994(Hansen, , 1987Press et al., 1992) to calculate a solution to the constrained inverse problem (Eqs. 5-7) because this method handles the constraints on ρ(ln k) well. The goal of Tikhonov regularization is to minimize both the residual error and the complexity of the solution ρ(ln k). Solution complexity is assumed to be associated with the roughness, or the intensity of fluctuations, in ρ(ln k). Here, we measure roughness by the norm of the first derivative of the solution vector, where R is the bi-diagonal first derivative operator, with an additional first row [1 0] and additional final row [0 −1] to account for ρ being zero outside the domain k min < k < k max . The regularization method proceeds by finding the solution ρ(ln k) that minimizes a sum of the residual error and the roughness: where ω is the regularization parameter which controls the relative weight of the solution roughness to the residual error. Because it is unclear how much to weigh the roughness, ω is typically varied many orders of magnitude from ω 1, which emphasizes the importance of residual error, to ω 1, which emphasizes the smoothest solution. The regularized ρ(ln k) is determined by identifying the value of ω that sets an optimal balance between the residual error and the roughness.

The L-curve
An approach for finding the optimal ω is to use the "L-curve" technique (Hansen, 1987(Hansen, , 1994). An L-curve is generated by parametrically varying ω and solving Eq. (9) for ρ, obtaining values for the residual error norm ĝ − Aρ and roughness norm Rρ for each value of ω. The L-curve is made by plotting Rρ vs. ĝ − Aρ on a log-scale, resulting in a characteristic "L" shape. Figure 3a shows the L-curve for the dataset shown in Fig. 1a. While each point on the curve is associated with an optimal solution for a certain value of ω, the corner of the Lcurve is associated with the regularized solution. The corner represents the point where increasing the information contained in the solution no longer improves the residual error. Solutions at the upper left of the L-curve contain high information as n elements ρ j are free and independent (Pierce,  Fig. 3. Determining the optimal value of regularization parameter ω via the "L-curve" method (Hansen, 1994(Hansen, , 1987.  Fig. 1a. The decay predicted from the solutions in B,C,D is also plotted. The corner with ω = 0.66 is chosen as the optimal solution since it predicts the data with significantly less residual error and has no trend in the residual, unlike the corner having ω = 20.1 The optimal solution is associated with a corner of the curve. 1980), whereas fewer modes are active in the smooth solution at the corner. The corner is also the point where relatively little decrease in information results in large increases in the residual error. Often, the corner is associated with a change in the number of maxima in ρ(ln k); solutions above the corner may have more modes than solutions below.

Biogeosciences
Generally, the L-curve method identifies a smooth solution which predicts a decay having a residual error approximately equal to the noise in the dataset. Solutions for three different values of ω are shown in Fig. 3b-d with the unregularized solution shown in Fig. 3b. Decays predicted from these three solutions are shown in Fig. 3e.
When there are multiple corners, as shown in Fig. 3a, the optimal corner needs to be chosen. We apply the approach discussed in Appendix A3 in order to choose the appropriate corner in Fig. 3a. The distribution corresponding to the lower corner (Fig. 3d) predicts a decay that has a large error compared to the noise in the data and has a significant trend in the residual error as seen in Fig. 3e. However, the distribution corresponding to the higher corner (Fig. 3c) predicts a decay which is hardly distinguishable from the decay predicted by the unregularized solution, shown in green, after t = 1. Therefore, the distribution shown in Fig. 3c is considered the optimal, regularized solution. Figure 4 shows the overall variation in the inversions due to noise by summing all 60 000 multi-pool solutions from Fig. 2a-g. We also plot the solution ρ(ln k) from the regularized inversion. The regularization provides a simpler solution that captures the variation of the rate heterogeneity associated with uncertainty of the multi-pool solutions.

Results of regularization
More specifically, we find that this heterogeneity takes a certain form: the regularized solution presented in Fig. 3c appears lognormal, i.e. Gaussian in ln k space. We have shown in a separate publication (Forney and Rothman, 2012) that the entire LIDET dataset is lognormal on average. This lognormal description is considerably simpler than the multipool approach as the rate heterogeneity and all dynamic information are now contained within just two parameters, a mean µ and a variance σ 2 of the order of magnitude of decay rates. This approach also emphasizes the view that the quality of litter and young soil organic matter quality is continuously distributed.  Fig. 2 plotted alongside the regularized solution (red). ρ(ln k) is calculated by summing ρ from all 60 000 noisegenerated trials and dividing the sum by 60 000.

Fitting lognormal parameters to decay data
The ubiquity of the the lognormal distribution suggests that litter decay can be well described by the two lognormal parameters, µ and σ . In this section, we assume a model of lognormally distributed rates and present a procedure to estimate the values of the mean of log rates, µ, and standard deviation of log rates, σ , associated with a decay dataset. We then compare its results to the results of the regularized inversion.
We proceed by substituting the Gaussian distribution, for ρ(ln k) in Eq.
(2), providing a prediction of the mass fraction g(t) remaining when decay rates are lognormally distributed. g(t) as a function of the parameters µ and σ is We then identify the values of µ and σ which best fit the data by solving the non-linear least squares minimization problem whereĝ i are the measured data points and g(t) is the decay predicted from Eq. (11). We solve Eq. (12) using Matlab's non-linear least squares solver nlinfit.m. This is repeated for all 191 LIDET datasets that are appropriately described by a superposition of exponential decays (Forney and Rothman, 2012). sion (Eqs. 5-7), regularized inversion Eq. (9), and the twoparameter fit Eq. (12). The root mean square error (RMSE) of the fit to the data is calculated for each dataset and the cumulative density function of the RMSE rescaled by the total number of datasets (191) is plotted; the vertical axis shows the number of datasets having an RMSE smaller than the value on the horizontal-axis. The RMSE of the unregularized inversion technique (green) is clearly smaller than that of the other two techniques, but these solutions suffer from the problems discussed in Sect. 3.3. Mean RMSE values for all 191 datasets are 0.048, 0.053, and 0.052 for the direct inversion, regularized inversion, and lognormal model, respectively. Surprisingly, solutions from the two-parameter fit (magenta) appear to have slightly smaller residual error than solutions identified from the regularization technique (blue). Thirty-two of the datasets are predicted by both the two-parameter fit and regularization method to have a single rate; the two-parameter lognormal model however predicts the decay from 126 of the remaining 159 datasets better than the distribution obtained by regularization. Therefore, the lognormal model fits the data equally or better than the regularized solution for 158/191 or 83 % of the datasets. This is unexpected because the lognormal model has only two degrees of freedom. The marginally poorer fit of the regularization method is due to the emphasis of the method on smoother (wider) solutions. We proceed to investigate this effect by comparing the fitted parameters σ to the standard deviation of the regularized inversions ρ(ln k). Figure 6 compares the fitted values µ and σ to the mean, ln k , and standard deviation (ln k − ln k ) 2 associated with the regularized solutions ρ. Figure 6a shows that values of µ and the mean of the inversion ρ(ln k) are consistent with one another. There is a departure at high ln k because some of these inversions are bi-modal: one mode is active whereas the other mode is distributed over extremely fast rates associated with instantaneous decays. These bi-modal datasets are effectively described by an instantaneous mass loss followed by decay that proceeds with a distribution of rates. Figure 6b, on the other hand, tells another story. Values of σ tend to be less than the standard deviation of ρ(ln k), indicating that the regularized solution is wider than the best fitting lognormal distribution. This is a consequence of weighing the solution roughness during regularization. A narrower lognormal solution with smaller residual error exists, but the weight of the roughness is large enough that the regularization method chooses a slightly wider solution that fits the data almost as well. Collectively, Figs. 5 and 6b indicate that when rates are heterogeneous, the transition from unimodal to multimodal solutions near the corner of the L-curve occurs before finding the unimodal solution with the smallest residual error.
These results indicate that the regularization method is useful for identifying general trends and shapes of solutions. If regularization suggests that rate distributions are lognormal, then fitting µ and σ to the data identifies more precisely the specific lognormal distribution that best fits the data.

Conclusions
Direct calculation of the inverse Laplace transform with a non-negativity constraint provides the best fitting multipool solution without specification of the number of pools a priori. However, this multi-pool solution is very sensitive to small changes in the decay function g(t) (Yeramian and Claverie, 1987), as the number of pools and the rates of each pool vary widely due to the level of noise associated with litter bag data. The regularization method described here is robust to noise and commonly indicates that a lognormal distribution provides a concise representation of the rates associated with decay data. If the lognormal distribution may be assumed, a marginally better fit to the data is found by directly estimating the lognormal parameters µ and σ from the decay g(t). However, we cannot recommend applying the lognormal model without first using the regularization procedure.
Matlab codes for all numerical procedures (direct inversion, multi-pool estimation, regularized inversion, and the two-parameter fitting procedure) are provided online in the Supplement.

A1 Discretizing the Laplace transform
Various methods can be used to discretize Eq. (2), such as quadrature (Lamanna, 2005;Hansen, 1994), linear or logspaced discretization. Here, we choose to discretize (2) at n nodes λ j spaced uniformly along the ln k axis between the limits ln k min < ln k < ln k max . The nodes are therefore spaced at intervals of width λ = ln k max − ln k min n .
This discretization is chosen in order to provide resolution over the appropriate wide range of k. 2 Eq.
The m×n matrix A is the discrete Laplace transform operator with elements A ij = e e λ j t i λ.
g is the vector of predicted time series data points g i = g(t i ) having length m. The Matlab function lsqnonneg.m employs Lagrange multipliers (Strang, 1986) to calculate the non-negative solution to the least squares problem (5). Constraint (7) is met by weighing the first data point at g(0) = 1 more heavily than the others, although an additional Lagrange multiplier could be used. Because m < n, A is underdetermined and has rank m. In Matlab 2008a and older, the function lsqnonneg.m uses the algorithm mldivide.m (Mathworks, 2009) to calculate A −1ĝ . This algorithm utilizes a rank-revealing QR factorization with column pivoting which calculates ρ only from the m most linearly independent orthogonal components of A. As a result, it returns at most m non-zero components in the vector ρ. In newer versions of Matlab, lsqnonneg.m uses the pseudo-inverse, pinv.m (Mathworks, 2011) to calculate A −1ĝ . The pseudo-inverse identifies the solution ρ with minimum norm. For our problem, we find that results from both the newer and older versions of lsqnonneg.m are the same.

A3 Procedure for choosing the corner of the L-curve
When there are multiple corners, as shown in Fig. 3a, one corner is typically more appropriate than the others. We use the following method to identify the corner associated with the simplest solution: 1. For each value of ω, calculate the rank correlation (Kendall and Gibbons, 1990) of the residual error and the corresponding P value.
2. Identify the region of the L-curve where the rank correlation of solutions have P > 0.1, suggesting no significant trend in the residual.
3. Choose the corner within the region of the L-curve having P > 0.1.
4. If a corner seems too rounded or undefined, the L-curve can also be plotted on linear Aρ −ĝ and Rρ axes. The corner collapses onto a much tighter volume in linear space and can sometimes be more readily chosen if the corner in log space has large curvature.
5. If there are multiple corners with P > 0.1, compare more closely the residual trend in all corners, as values P > 0.1 do not necessarily indicate a lack of trend in the residual. If there is a noticeable difference in the quality of the fit associated with each corner, choose the corner that does not appear to have a trend in the residual. An additional test is to look at the error bias. Here, the error bias is simply the sum of residual error divided by the error norm, | (Aρ−ĝ)| Aρ−ĝ . Experience with this data suggests choosing solutions with error bias < .25. If multiple corners still have a bias < .25 and a similar strength in trend in their residual error, choose the corner with lowest roughness. Usually the corner having lower residual error and higher roughness tends to be not biased and to have no trend in the residual.