A generic hierarchical model of organic matter degradation and preservation in aquatic systems

Organic matter degradation and preservation are crucial components of Earth’s carbon cycle. Empirical and phenomenological models usually contain parameters determined by site-specific data and focus on different aspects of the decay and accretion characteristics. To investigate more fundamental mechanisms, this study suggests a hierarchical model that links microscopic physical quantities to macroscopic degradation and preservation patterns. This mechanistic model predicts several commonly observed phenomena, including the lognormal distribution of degradation rate constants, the recalcitrance-dependent sensitivity to temperature, the dependence of a heterogeneous organic-matter system’s persistence on its complexity, logarithmic-time decay, and power-law degradation behavior. The theoretical predictions of this model are consistent with the observational data from marine and lake environments. This hierarchical model may provide a step towards a fundamental theory of organic matter degradation and preservation in aquatic and other ecosystems. Organic matter in aquatic systems is broken down in a hierarchical ensemble, involving a decrease in effective kinetic barriers, according to a mechanistic model using decay kinetics that predicts degradation and preservation characteristics consistent with observations.

E arth's carbon cycle provides us a window through which to investigate the interactions between life and the environment in different spatial and temporal scales 1,2 . As two main components of the carbon cycle, oxygenic photosynthesis and organic degradation play important roles in controlling these interactions 1,2 . In oxygenic photosynthesis, photoautotrophs (e.g., plants, algae, and cyanobacteria) use water (H 2 O) and carbon dioxide (CO 2 ) as raw materials and energy in the form of photons to produce organic matter (CH 2 O) and molecular oxygen (O 2 ) via the forward reaction in the following chemical equation: energy + H 2 O + CO 2 ⇌ CH 2 O + O 2 , where CH 2 O is a schematic chemical formula representing a large variety of organic compounds produced by photosynthetic activities. The backward reaction (i.e., aerobic respiration) burns organic matter with O 2 and provides energy to heterotrophic organisms.
The loop between photosynthesis and respiration drives the material and energy exchanges between the biological and geological components of Earth's carbon cycle. However, there exists a tiny leak in this loop. A small amount of organic matter originating from photosynthetic activities escapes degradation and is buried in sedimentary rocks, leading to the accumulation of O 2 in Earth's atmosphere and oceans 1,2 . As a result, organic matter degradation and preservation determine the net CO 2 flux out of Earth's surface and the net O 2 flux into it.
Organic matter is highly heterogeneous 3 ; so, too, is the environment in which it is deposited 4 . The complex interplay between the intrinsic molecular composition and extrinsic environmental variables influences the degradability of organic matter. In aquatic (e.g., marine and lake) systems, the chemical properties and degradation/preservation mechanisms of organic matter are altered as it sinks through the water column and accumulates in sediments 5,6 . Correspondingly, different interpretations have been proposed for the recalcitrant fraction in dissolved and sedimentary organic matter. For example, the dilution hypothesis [7][8][9] suggests that low concentrations of labile components in dissolved organic matter may significantly reduce the effective energy gains from microbial respiration, resulting in the recalcitrance of these organic compounds in the water column. Instead, the physical protection hypothesis 10,11 proposes that the association of sedimentary organic matter with minerals can protect the former from enzymatic attacks of microorganisms when it passes through the thermodynamic ladder 12 of redox zones in sediments, leading to the persistence of organic matter. These hypotheses offer qualitative understandings of organic matter degradation and preservation in marine and lake environments; investigating the fundamental microscopic mechanisms of degradation/preservation and predicting the long-term fates of organic matter, however, require a more quantitative theoretical framework.
Different empirical and phenomenological models [13][14][15][16][17] have been proposed to interpret the spatial and temporal variations of organic matter. Nevertheless, the mathematical formulations of existing models usually are not well constrained and contain parameters determined by site-specific data 18 . Moreover, studies on organic matter decay and accretion often focus on thermodynamics 12,19 ; biogeochemical models considering kinetics are generally dictated by simple rate constants [13][14][15][16][17] . Understanding the underlying microscopic mechanisms of organic matter degradation and preservation is necessary for developing a more fundamental theory 4 .
This work proposes a hierarchical model that predicts several phenomena commonly observed in previous studies, including the lognormal distribution of degradation rate constants 3,20 , the recalcitrance-dependent sensitivity to temperature [21][22][23] , the dependence of a heterogeneous organic-matter system's persistence on its complexity 24-26 , logarithmic-time decay [27][28][29] , and power-law degradation behavior 14,[30][31][32] . The plausible consistency of this model's predictions with the field observations in marine and lake environments suggests that the theory presented here potentially provides a step toward a fundamental theoretical framework of organic degradation and preservation in aquatic and possibly other ecosystems.

Results
A discrete hierarchical ensemble of organic matter. Previous studies have suggested a variety of physical, chemical, and biological factors determining the persistence of organic matter 4-6 . Many factors influence the decay and burial of organic compounds by affecting their reactivities and activation barriers 10,11 . For instance, the physical protection hypothesis suggests that the long-term preservation of some types of organic compounds in marine and lake sediments is a consequence of their strong association with mineral surfaces, which elevates the activation energy for degradation and prevents these organic compounds from the attacks of microbial enzymes 10,11 . Here, I suggest a hierarchical model for organic matter degradation from the perspective of kinetics (i.e., reactivity and activation energy).
I consider a hierarchical ensemble of organic compounds {C i } with degradation rate constant k i , where i ∈ [0, N], N is the highest level in the hierarchy, and C i represents the organic matter at level i. Figure 1(a) shows a series of decay reactions occurring along a degradation coordinate: C N ! k N C NÀ1 ; ::: ; C i ! k i C iÀ1 ; ::: In this hierarchical ensemble, C N represents the most recalcitrant organic matter, which has the longest persistence, while C 0 is the most labile organic matter, which is preferentially utilized by organisms in natural environments. The activation energy for the final step on the degradation coordinate is denoted by E 0 , and the activation energy for the degradation reaction occurring i steps prior to the final step is denoted by E i . One should note that "hierarchical" here does not mean "sequential"; the degradation of organic compounds at different hierarchical levels occurs in parallel simultaneously. The objective of this hierarchical model is to investigate the connection of organic degradation patterns to kinetics (e.g., k i and E i ) rather than to thermodynamics (e.g., the Gibbs free energy). Thus, I assume that each reaction at level i along the degradation coordinate has a negative Gibbs free energy change (i.e., ΔG i < 0) and therefore is thermodynamically feasible. Moreover, in this hierarchical model, organic compounds having the same activation energy E i are treated as identical component C i located at the same hierarchical level (i.e., the i-th level) regardless of their compositions and structures.
According to the transition state theory 33,34 , the rate constant for the i-th reaction is expressed by where k B is the Boltzmann constant, ℏ is the Planck constant, T is thermodynamic temperature, and R is the ideal gas constant. Organic compounds at different hierarchical levels start to decay from distinct points on the degradation coordinate ( Fig. 1(a)). Their eventual fates, however, are the same: being completely degraded to CO 2 and removed from the ensemble. I define the total activation energy for the complete degradation of organic matter at the i-th level as its effective activation energy: organic matter at the i-th level is defined as I substitute Eqs. (1) and (2) into Eq. (3), apply the product rule for exponents, and obtain Two direct predictions of the hierarchical model. From the perspective of stochastic reaction kinetics 35,36 , k j represents the probability per unit time that the i-th reaction on the degradation coordinate (i.e., the transformation of organic matter from the j-th level to the (j − 1)-th level) will occur. Therefore, κ i represents the probability per unit time that organic matter at the i-th level will completely decay to CO 2 . The above relation between the effective rate constant (κ i ) and conventional rate constant (k j ) can be rewritten as ln κ i $ ∑ i j¼0 ln k j . Because (i) organic matter degradation usually consists of plenty of steps to complete (i.e., 1 ≪ i) and (ii) k m and k n are independent for m ≠ n, the Central Limit Theorem implies that κ i satisfies a lognormal distribution 37,38 : where N is a normal distribution with mean μ i $ ∑ i j¼0 Eðln k j Þ and variance σ 2 i $ ∑ i j¼0 Varðln k j Þ. This result (Eq. (4)), a particular case of which appears in Fig. 1(b), is consistent with previous theories and observations 3,20 and suggests an alternative interpretation for the origin of the lognormal distribution of organic matter degradation rate constants.
The other direct prediction of the hierarchical model is the recalcitrance-dependent sensitivity to temperature. The relative sensitivity of the organic matter degradation rate at the i-th hierarchical level to temperature can be expressed as This equation, also illustrated in Fig. 1(c), shows that, when organic matter has higher effective activation energy, the dependency of its degradation rate on temperature is stronger. This relation is consistent with a phenomenon observed in previous studies [21][22][23] : the decay of more recalcitrant organic matter responds more sensitively to changes in temperature.
Organic matter degradation and preservation on a continuum. As discussed above, in the hierarchical model, organic compounds having the same activation energy E i are treated as identical component C i ; they are initially located at the i-th level and start to decay from this position. I assume that the activation barrier for the i-th decay reaction consists of n i barrier units and denote the activation energy of each barrier unit by Δu. The activation energy at the i-th reaction then can be written as I denote the initial total amount of organic matter in the whole hierarchical ensemble {C i } by G 0 , and the fraction of organic matter at the i-th level (i.e., with n i activation barriers) by f i ; the initial amount of component C i then is expressed as G 0 f i . I assume that barrier units are uniformly distributed at each hierarchical level along the degradation coordinate, implying that n i is proportional to C i amount and therefore f i ' n i ∑ N j¼0 n j . Following the conventional first-order kinetics [13][14][15] , the amount of leftover C i after t is expressed by g i ðtÞ ¼ G 0 f i expðÀκ i tÞ. The ensemble consists of N hierarchical levels; thus, there are N different degradation processes occurring in parallel simultaneously. Therefore, the total amount of remaining organic compounds in the whole ensemble at time t is GðtÞ ¼ ∑ N i¼0 g i ðtÞ, and the preservation efficiency is defined as B(t) = G(t)/G 0 . I substitute Eqs. Fig. 1 A profile of activation energy along the discrete degradation coordinate and two direct predictions of the hierarchical model. a Activation energy profile for organic matter degradation, in which C i represents the organic matter at level i, and E i represents the activation energy for the i-th reaction. The energy profile illustrated here appears along a discrete degradation coordinate; its continuous form is introduced later. b Lognormal distribution of the effective degradation rate constant (κ i ) predicted by Eq. (4). The mean and variance of this lognormal distribution depends on the type of organic matter and the environment in which it decays. The specific distribution presented in panel b is for the degradation rate constants of plant biomass 3 . c Recalcitrance-dependent sensitivity to temperature predicted by Eq. (5). The orange contours represent the temperature sensitivity of rate constants (i.e., ∂κ i /∂T). When organic matter is more recalcitrant (i.e., the effective activation energy, ϵ i , is higher), its degradation responds more sensitively to the changes in temperature.

Energy Degradation Coordinate
(1)-(3) and (6) into this definition of B(t) and express the preservation efficiency of the whole ensemble as a linear superposition of the fractions of leftover organic matter at each hierarchical level: where ϵ max ¼ ∑ N i¼0 E i is the effective activation energy at the highest hierarchical level. Eq. (7) links the macroscopic quantity B to the microscopic quantities n i and ϵ max on a discrete degradation coordinate.
Theoretical studies 4,14,15 and laboratory/field observations 39,40 have suggested that organic matter in natural environments can be represented by a continuum of species with a wide range of reactivities. Compared to the discrete coordinate ( Fig. 1(a)) from which Eq. (7) is derived, a continuous coordinate, which is henceforth called a degradation continuum, would provide a more realistic depiction of the degradation path 4 . In such a scenario, the summation in Eq. (7) can be rewritten as an integral. To do so, I first assume that the number of activation barrier units (n i ) changes smoothly with i along the degradation coordinate ( Fig. 1(a)). This assumption allows us to introduce an (almost) continuous variable x i = n i δ, where δ ≪ 1. I define a functionf ðx i Þ to represent the fraction of the continuous variable x i over the summation of x i on all N levels: x j , which is analogous to the function f i for the discrete variable n i . The fraction of organic matter that has effective activation energy up to level i on the degradation continuum then can be expressed by R x i x 0 e f ðx 0 Þdx 0 . To obtain the mathematical expression of B(t) on the degradation continuum from Eq. (7), x 0 e f ðx 0 Þdx 0 (the Methods section), where x max ¼ Nδ. This variable s is the normalized, continuous form of ð∑ i j¼0 n j Þ ð∑ N j¼0 n j Þ; the denominator x 0 e f ðx 0 Þdx 0 in s is a normalization factor that guarantees 0≤s≤1 (the Methods section). Eq. (7) then can be rewritten as the following integral form (the Methods section): BðtÞ ¼ R 1 0 exp Àκ 0 t exp Àðϵ max sÞ=ðRTÞ Â Ã À Á ds. Further calculation (the Methods section) transforms this integral to an expression for B(t) in terms of integral exponential functions: where τ N = 1/κ N and τ 0 = 1/κ 0 are the length of time for the most recalcitrant and the most labile organic compounds, respectively, in the ensemble to completely decay to CO 2 , and (8) is characterized by four parameters: τ 0 and τ N determine the time scales, and ϵ max and T determine the functional shape of B. It is important to note that τ 0 and τ N are functions of T and ϵ max (Eq. (3)); therefore Eq. (8) does not simply claim that B is positively correlated with T while negatively correlated with ϵ max . In contrast, as demonstrated in the Methods section, B depends negatively on T and positively on ϵ max . The negative correlation between B and T is consistent with a broadly observed phenomenon: an increase in temperature would promote the remineralization of organic matter, resulting in less preservation 41,42 . In the following sections, I show that several commonly observed phenomena of organic matter degradation and preservation can be derived from Eq. (8) and that the theoretical predictions are consistent with observational data from aquatic systems.
Logarithmic oxygen-exposure-time decay. In oxic environments, t in Eq. (8) is conventionally interpreted as oxygen exposure time (t O 2 ) 27,28 , that is, the length of time during which organic matter is exposed to O 2 . Oxygen exposure time has been suggested as a proxy for the aerobic degradation of organic matter 27,28 . Figure 2 shows the data on "preservation efficiency versus oxygen-exposure time" in sediments (circles) compiled from field observations [27][28][29]43,44 , which are used to test the theoretical relation between B and t O 2 in Eq. (8). These data are measured from the sediments deposited on the surface of continental shelves or lake bottoms, where most of the organic matter burial on Earth takes place 1,10 . The temperature of such environments usually varies from 273 K to 303 K with an average of 288 K 41,45 . The timescale of the compiled data ranges from 10 −2 to 10 3 years; I assume that τ 0 and τ N in Eq. (8) equal 10 −2 and 10 3 years, respectively. Supplementary Table 1 summarizes the activation energy for organic matter degradation in marine and lake sediments. According to these data, I take 25 kJ/mol and 140 kJ/mol as the lower and upper limits of ϵ max , respectively. Figure 2 displays the theoretical relation between B and t O 2 predicted by Eq. (8). The red and purple curves represent the predictions with the lower and upper limits of ϵ max , respectively; more than 85% of the compiled data fall in the region between them (Fig. 2). These curves are obtained without fitting free parameters, which is different from the traditional approaches in which models with free parameters (e.g., linear regression) 27,28,43,44 are used to fit the empirical relation between B and t O 2 .
Based on the data in Fig. 2, previous studies [27][28][29]43,44 have suggested that the preservation efficiency of organic matter in marine and lake sediments is a monotonically decreasing function of the logarithm of oxygen-exposure time. Here, I show that logarithmic-t O 2 decay of organic matter can be derived from Eq. (8). In natural environments, highly labile organic compounds decay in several minutes or even shorter, while extremely recalcitrant organic compounds can be stored for millions of years or even longer 4,46 . Asymptotic analysis of Eq. (8) (the Methods section) shows that, when degradation time t O 2 falls between its lower and upper limits (i.e., τ 0 ( t O 2 ( τ N ), preservation efficiency in sediments can be written as the following asymptotic form: which suggests that the slope of preservation efficiency in sediments versus oxygen-exposure time depends on ϵ max and T.
The logarithmic-t O 2 decay shown in Eq. (9) effectively ceases at time τ N because the Bðt O 2 Þ ' 0 when t O 2 ) τ N (the Methods section). This theoretical logarithmic-t O 2 decay pattern (Eq. (9)) is consistent with the empirical relations suggested by previous studies [27][28][29]43,44 , indicating that the hierarchical model provides insights into the underlying mechanisms of organic matter decay and accretion in marine and lake sediments.
Dependence of an organic-matter system's persistence on its complexity. By definition, the preservation efficiency function B(t) represents the fraction of the organic matter at age t (i.e., the waiting time before degradation occurs), where t = 0 is set at the time point when organic matter initially arrives in ocean or lake.
Field studies have shown that the age of organic matter ranges from a few minutes (or even shorter) to millions of years (or even longer) 4,46 . Here, I assume that the length of the degradation time for organic compounds in the hierarchical ensemble varies from 0 to ∞. As shown in the Methods section, B(t) → 1 as t → 0 while B(t) → 0 as t → ∞. These two limits are consistent with the physical definition of B(t) as a ratio between the remaining and initial amounts of organic matter: this ratio can be no higher than 1 because secondary production always comes at the expense of degradation of other organic compounds and will reach 0 when all organic matter is consumed. Since B(t) ranges from 0 to 1, it also represents the cumulative distribution function (starting at t = ∞) of the ages of organic compounds in the ensemble, and the corresponding probability density function p(t) is pðtÞ ¼ dBðtÞ=dt . The average age of all organic compounds in the hierarchical ensemble, hti ¼ R 1 0 tpðtÞdt, is then expressed as in which 〈 ⋅ 〉 represents the average of a quantity. The derivation of Eq. (10) is provided in the Methods section. The relation in Eq. (10) suggests that the average age of an organic-matter system increases with its maximum effective activation energy, which depends on the complexity of this system. When the complexity of composition (e.g., molecular diversity) and the environment (e.g., spatial heterogeneity and temporal variability) rises, the number of parallel degradation processes and thus the total number of hierarchical levels (N) will also increase. This growth of N would lead to a higher ϵ max because ϵ max ¼ ∑ N i¼0 n i Δu. Consequently, the ensemble of organic matter rewards less energy to microbial metabolisms 47 and therefore has a longer average age 〈t〉. In other words, Eq. (10) suggests that an organic-matter ensemble with higher complexity tends to be more persistent, which is consistent with previous observations in aquatic systems [24][25][26] .
Power-law degradation behavior. To investigate how the degradation kinetics change with time, I express the degradation of the whole hierarchical ensemble in terms of the conventional first-order kinetics [13][14][15] : dG(t)/dt = − K(t)G(t), where K(t) is the average degradation rate constant of all organic compounds in the ensemble. When the length of time (i.e., age) t is within the two extreme timescales, that is, τ 0 ≪ t ≪ τ N , I have The derivation of Eq. (11) is provided in the Methods section. On a log-log graph, the relation between K and t in Eq. (11) is represented by a straight line with a slope of -1 and an intercept of b, which is consistent with the power-law behavior shown by the previous observations 14,30-32 : log 10 K ¼ Àa log 10 t þ b. Table 1 summarizes the values of a and b obtained by empirical studies from marine 14,30,31 and lake 32 environments and predicted by the hierarchical model in this study. The theoretical value of a equals 1; this result is consistent with the empirical values of a, which approximate 1, from field studies 14,30-32 . The value of b, however, varies in a relatively large range. The systemdependent variation of b values demonstrated in Table 1 may be explained by the dependence of b on τ N shown in Eq. (11): the factors influencing τ N and therefore b, such as organic matter composition and environmental conditions, change across the observation sites in these studies 14,30,32 .
Data (blue circles and orange squares in Fig. 3) measured from the marine 30 and lake 32 systems are compiled to test the predictions of Eq. (11). To obtain the theoretical estimate of intercept b over these data, I first determine the value of the characteristic degradation time τ N . The minimum degradation rate constant in the compiled dataset (Fig. 3) is K min ' 2 10 À7 year −1 , implying that Table 1 Parameter values in the power-law relation, log 10 K ¼ Àa log 10 t þ b, in published studies.  Fig. 3 Theoretical prediction and field data for the power-law relation between organic carbon reactivity (K) and age (t) in marine and lake systems. Blue circles and orange squares represent data from marine 30 and lake systems 32 , respectively. Red and grey straight lines are the theoretical approximation of Eq. (11) and linear regression, respectively, for the data from both marine and lake environments. R 2 represents the coefficient of determination of the theoretical approximation or linear regression. τ N ¼ 1=K min ' 5 10 6 year. I substitute this τ N value into Eq. (11) and calculate the average of log 10 exp Àt=τ N À Á =Ei t=τ N À Á Â Ã over t of the data in Fig. 3; this average, − 0.94, is the theoretical estimate of b for the compiled dataset. The red line (log 10 K ¼ Àlog 10 t À 0:94; R 2 = 0.98) in Fig. 3 shows this theoretical relation, while the grey line (log 10 K ¼ Àlog 10 t À 0:64; R 2 = 0.96) in Fig. 3 is obtained via linear regression, which is the conventional method used in the previous studies 30,32 . The theoretical line has a coefficient of determination (R 2 ) that is as good as the linear-regression line. Given that the value of b is system-dependent, I apply the same analyses to marine data or lake data individually to test if the predictions by Eq. (11) hold for separate, smaller datasets. The results, which are presented in Supplementary Figure 1 and Supplementary Table 2, show that the theoretical formula still performs as well as or even better than the traditional linear-regression method.
Relating organic preservation hypotheses to the hierarchical model. The results presented in previous sections (Eqs. (9), (10), and (11)) show that the patterns of organic matter decay and accretion are primarily constrained by the lower and upper extremes of effective activation energy in the hierarchical ensemble: ϵ 0 and ϵ max (note that κ 0 , κ N , τ 0 and τ N are all functions of ϵ 0 or ϵ max ). This suggests that the degradation and preservation of organic matter are influenced by three microscopic quantities in the hierarchical model, N, Δu, and n (Eqs. (2) and (6)), which characterize an organic-matter ensemble's complexity, the energy of each barrier unit, and the number of unit barriers at each hierarchical level, respectively. These physical quantities potentially shed light on the underlying mechanisms of organic matter decay and persistence. Figure 4 illustrates the possible connections of the three microscopic quantities in the hierarchical model with several hypotheses of organic matter preservation. Each of these hypotheses can be characterized by either one quantity or a combination of multiple quantities. For example, the complexity hypothesis [24][25][26] is represented in hierarchical model by N. As the complexity (e.g., the molecular diversity, spatial heterogeneity, and temporal variability) of an organic-matter ensemble rises, its retention time (i.e., age) becomes longer. Other explanations for the long-term preservation of organic matter, including intrinsic recalcitrance (e.g., inherent qualities of organic matter) 4-6 , external environments (e.g., physical and chemical factors shaping microbial communities) [4][5][6] , and physical protection (e.g., mineral-organic matter association) 10,11 , are characterized by nΔu in the hierarchical model. As these factors become stronger, nΔu increases, leading to a more persistent organicmatter ensemble.
The conceptual diagram in Fig. 4 provides a primitive but plausible approach to link microscopic quantities to the existing hypotheses of organic preservation. However, this framework requires further theoretical validation and experimental tests. Moreover, the hierarchical model is based on the classical theory of chemical kinetics; incorporating quantities/principles from thermodynamics (e.g., Gibbs free energy 12,19 ) and quantum statistics (e.g., Fermi-Dirac statistics 48 ) into the current model may offer deeper insights into more fundamental mechanisms of organic matter degradation and preservation.

Discussion
The degradation and preservation of organic matter appears to be a complex process because its underlying mechanisms have not been well understood 4 . Occam's razor 49 and the simplicity postulate 50 both point in the direction of developing simple models and evaluating their predictions. From the perspective of kinetics, this study develops a hierarchical model to depict organic matter decay and accretion and predicts several commonly observed phenomena; in particular, this model provides excellent predictions that correspond to logarithmic-t O 2 decay 27-29 and power-law degradation behavior 14,[30][31][32] . The consistency between theoretical predictions and observational data in marine and lake environments suggests that this hierarchical model potentially offers a framework for a fundamental theory of organic matter degradation and preservation in aquatic and possibly other ecosystems.
There are a variety of valid approaches for modeling organic matter decay [13][14][15][16][17] . Particularly, the mathematical form of Eq. (11) in this study is identical to that of Eq. (5) in the theory of Rothman and Forney 16 , which presents a reaction-diffusion model for the degradation of mineral-associated organic matter in the porous marine sediments. Nevertheless, these two formulas are obtained with different physical assumptions and via distinct mathematical derivations. The model of Rothman and Forney 16 made an assumption that reactivity k of organic matter satisfies a special probability distribution p(k)~1/k and obtained its Eq. (5) from the connection of reactivitties to reaction-diffusion dynamics. These differ from the assumptions of hierarchical structures (the Results section) and the mathematical approach of deriving Eq. (11) (the Methods section) in this study.
I do not claim that the hierarchical model is superior to the previous models. This is not only because models for any open natural system are nonunique, but also because models, even partially confirmed, have their own heuristic values 51 . However, I have found that, compared to previous empirical and phenomenological models 13-17 , the hierarchical model provides an alternative approach to depict organic degradation and Quantity N represents an organic-matter ensemble's complexity, Δu is the energy of each barrier unit, and n characterizes the number of unit barriers at each hierarchical level. Quantity N is related to the complexity hypothesis; nΔu is linked to several hypotheses, including intrinsic recalcitrance, external environments, and physical protection. When N, n, and Δu of an organic matter system are all relatively small (i.e., inside the orange tetrahedral region), it tends to be labile. In contrast, an ensemble of organic matter is recalcitrant when one or more of N, n, and Δu are higher than some thresholds (i.e., outside the orange tetrahedral region).
preservation from a mechanistic perspective. Although some quantities, such as activation energy, have been widely applied to investigate organic matter's fate in natural environments, this work mathematically connect these quantities to logarithmic-time decay (Eq. (9)) and power-law degradation behavior (Eq. (11)) and shows plausible consistency of theory with observations. This study mainly focuses on the geochemical aspects of organic matter degradation/preservation and develops the hierarchical model on the basis of the traditional first-order kinetics, in which the biological factors are implicitly integrated into the first-order rate constants (i.e., κ i 's). Although the hierarchical model offers a potential theoretical framework to link the chemically-oriented hypotheses for the persistence of recalcitrant organic matter (Fig. 4) to microscopic physical quantities, providing insights into the biological-oriented explanations for recalcitrance, such as the dilution hypothesis [7][8][9] , is beyond the ability of the hierarchical model. As introduced above, the dilution effect, which derives from the imbalance between microbial respiration and its energy rewards from dilute substrates, is basically independent of the chemical properties of organic matter and therefore is not covered by the geochemically-based theoretical framework in this work.
Similar to the hierarchical model presented here, studies on organic matter decay in the water column and sediments [13][14][15][16][17] often formulate models using the first-order kinetics with simplistic assumptions that microbial consumers are present in sufficient numbers and thus degradation processes are independent of the abundance of microorganisms. Microbial populations and metabolisms, however, significantly influence the fates of organic matter in the water column 24,52,53 and sediments 5,6,18 in aquatic systems and therefore deserve separate treatment in future work. Explicitly incorporating biological factors into the theoretical framework proposed in this work may provide a better understanding and justification of the biologically-oriented interpretations for the persistence of recalcitrant organic matter.
Due to the limited availability of data, I tested only the predictions of logarithmic-t O 2 decay and power-law degradation behavior using the measurements from marine and lake environments. The fundamental mechanisms in Earth's carbon cycle, however, should be universal 4,54 ; the results here (Eqs. (9), (10), and (11)) therefore are expected to apply to organic matter degradation and preservation in other ecosystems. Examining these theoretical predictions in other environments, such as soils 55,56 and rivers 57,58 , may provide further support for the theory presented in this study. Moreover, several measurable indices, such as amino acid and amino sugar compositions, have been suggested as proxies for the degradation states of organic matter [59][60][61] ; these degradation states are analogous to the hierarchical levels in the degradation coordinate ( Fig. 1(a)). Researching the relations between activation energy/reactivity and decay patterns of organic matter in different degradation states (measured using those indices) would offer valuable validation for the hierarchical model and its theoretical predictions.
Hierarchical structures are ubiquitous in nature 62 , implying potentially wide applications of the theory in this work. For example, many synthetic polymers (e.g., plastics) possess similar molecular structures and compositions, backbone linkages, and degradation mechanisms to those of natural organic matter 63,64 . Investigating synthetic polymers using the hierarchical model would yield insights into their decay and accretion patterns and offer clues to develop more effective, environmentally friendly management of polymer wastes. Furthermore, this study focuses on organic matter degradation and preservation in modern environments; their underlying mechanisms, however, are expected to have been the same or similar on the ancient Earth 1,2 .
Incorporating the theoretical results presented here, especially the preservation efficiency function, with geochemical identification of the specific types of buried organic compounds in deep time and laboratory measurements of their effective activation energy, into global-scale biogeochemical models 65,66 may provide new perspectives on significant events in the evolution of Earth's carbon and oxygen cycles, such as carbon isotope excursions 67,68 and atmosphere-ocean oxygenation 69,70 , over geologic timescales.

Methods
Preservation efficiency on a degradation continuum. In the Results section, I defined x i = n i δ and e f ðx i Þ ¼ On the degradation continuum, the ratio between the two summations in Eq. (12) can be written as a ratio between two integrals: in which the integral over ½x 0 ; x max , i.e., R x max x 0 e f ðx 0 Þdx 0 , is a normalization factor and the " ≡ " sign means that the variable s ∈ [0, 1] is defined as the ratio between the two integrals.
According to the first fundamental theorem of calculus 71 , I express the first derivative of s with respect to x i as With Eqs. (13) and (14), I rewrite Eq. (7) as ð15Þ Finally, I derive Eq. (8) from Eq. (15): where Ei( ⋅ ) is defined as in which R Ã represents non-zero real numbers. Eq. (17) is an integral exponential function. In the above derivation, two variables ν ¼ exp Àðϵ max sÞ=ðRTÞ Â Ã and ω = κ 0 tv are introduced for "integration by substitution". With the definitions of effective activation energy and effective degradation rate constant (Eqs. (2) and (3)), the last line of Eq. (16) can be rewritten as Eq. (8).
The dependence of B on T and ϵ max . Rewriting the discrete form of the preservation efficiency (Eq. (7)) to its continuous form (Eq. (8)) does not change the intrinsic characteristics of the preservation efficiency function itself. Here, I use the discrete form (Eq. (8)) to investigate the dependence of the preservation efficiency B on temperature T and the maximum effective activation energy ϵ max . With Eqs.
(1) and (3), the expression of B in Eq. (7) can be rewritten as Then I take the partial derivatives of B with respect to T and ϵ max : Since the exponential function expðÁÞ is always positive, Eqs. (19) and (20) imply These two inequalities show that B is negatively correlated with T but positively correlated with ϵ max .
Asymptotic analyses of B(t). There are two timescales in Eq. (8): the long timescale determined by τ N and the short timescale determined by τ 0 , which are corresponding to the organic compounds that have the highest and lowest effective activation energy, respectively. Here, I perform asymptotic analyses for B(t) in three regimes of degradation time t: short (i.e., t ≪ τ 0 ), intermediate (i.e., τ 0 ≪ t ≪ τ N ), and long (i.e., τ N ≪ t). In the short-time regime, t/τ N ≪ t/τ 0 ≪ 1, under which an asymptotic expansion 72 of Eq. (8) leads to This indicates that the preservation efficiency slightly decreases as degradation time increases in this short-time regime. When degradation time t is very short, Eq. (22) implies BðtÞ ' 1 as t ! 0: In the intermediate-time regime, t/τ N ≪ 1 ≪ t/τ 0 , under which an asymptotic expansion 72