In regions of sediment accumulation, the rates of deposition can vary in space and time even when factors like climatic conditions are kept constant. This type of "random" variability complicates the interpretation of the geologic record because many important chemical signals are affected by changes in sediment deposition rates. To explore the specific effects of random variability, we constructed a model with realistic, time-varying sedimentation rates that tracks the burial and oxidation of organic matter. We chose organic matter because it is important to the global carbon cycle and is also commonly used to infer past environmental changes. We find that random changes in sedimentation rates can meaningfully impact organic matter burial and create patterns that look conspicuously like field observations interpreted as reflecting changes in environmental conditions. By predicting the patterns that can form randomly, our approach can help improve the interpretations of the geologic record.
In the existing diagenetic models, sedimentation is often assumed to be constant, or, if allowed to vary, is parameterized using step functions where the sedimentation rate is held constant over millions of years or tens to hundreds of meters of sediment cores (Arndt et al., 2009;Wehrmann et al., 2013). However, this is likely insufficient as natural sedimentation rate variability is much more complex. On monthly to decadal timescales, sedimentation could vary as a function of productivity (Collins et al., 2011), river discharge (Syvitski & Morehead, 1999), and strength of tidal and coastal currents (Nittrouer et al., 1986). Within a depositional environment, sedimentation rates also vary due to self-organization wherein deposition is locally focused at any given moment, but the position of this locus varies in time. These factors, which arise due to the internal dynamics of sedimentary systems, are termed "autogenic" behavior (Hajek & Straub, 2017). They are prevalent in all depositional environments (Sadler & Jerolmack, 2015) and give rise to timescale dependence in measurements of apparent sedimentation rates (Sadler, 1981).
Although long-term, basin-wide deposition must be deterministic due to mass balance constraints (Hajek & Straub, 2017;Sadler & Jerolmack, 2015), sedimentation in 1D, the dimension in which chemostratigraphic records are often generated and interpreted, has a stochastic component due to these autogenic dynamics (Hajek & Straub, 2017). Stochastic sedimentation variability should influence diagenetic reactions by introducing variations in the amount of time that sediments remain in diffusive contact with oxidant-laden overlying (pore) water. Consequently, we hypothesize that autogenic variations in preservation efficiency are an inherent component of geochemical signals found in sedimentary records, and to test this hypothesis, incorporate autogenic sedimentation dynamics into a reactive-transport model.
To simulate OC diagenesis influenced by autogenic sedimentary dynamics, we constructed a reactive-transport model with the following diagenetic equation (Berner, 1980):
where x is the depth below the sediment water interface, t is time, G is the OC concentration in sediments, ω is the sedimentation rate, and k is the reaction rate constant of OC oxidation, describing the reactivity of OC.
We parameterize the reaction term as a first-order decay with respect to OC concentration (-k ⋅ G). While more complicated kinetic models exist (e.g., Monod kinetics (Soetaert et al., 1996), multi-G approach (Archer et al., 2002;Boudreau, 1996)), we opt to use this simple linear model initially in order to be able to understand the fundamental interactions between OC degradation and sedimentary dynamics. Likewise, we lump all oxidants together and consider a single oxidant that controls the reactions. For conceptual and interpretative simplicity, we refer to this oxidant as oxygen, but acknowledge that other terminal electron acceptors are important (Bowles et al., 2014).
Both k and 𝐴𝐴 𝐴𝐴 are highly variable in natural systems (Freitas et al., 2021;Middelburg, 1989;Sadler, 1981) and can be uncertain. We therefore use the Dämkohler number (Da), which is a dimensionless parameter that describes the relative timescales of reaction (OC degradation) to transport (sedimentation), to characterize individual systems. Da can be calculated using Equation 2, where 𝐴𝐴 𝐴 𝐴𝐴 is the mean sedimentation rate of the system and z oxygen depletion is the depth at which virtually all oxygen is consumed by OC (i.e., depth of oxygen depletion), which can be derived from steady-state, analytical solutions to Equation 1 (Berner, 1980) (see Text S1 in Supporting Information S1).
Da encapsulates the combined effects of sedimentation, oxidant availability and intrinsic reactivity on OC diagenesis at the steady state. Thus, by performing simulations over the relevant range of Da (see Text S1 and Table S3 in Supporting Information S1), we can circumvent the need to experiment with all possible combinations of k and 𝐴𝐴 𝐴 𝐴𝐴 while still capturing the range of behaviors expected for most natural systems and facilitating the comparison between our multivariate model simulations. Although the depth of oxygen depletion is expected to vary temporally as a result of varying sedimentation rates within each of our model simulations, Da defined by the mean sedimentation rate is still a reliable descriptor of the system. At any given time, the system does not deviate significantly from the steady-state scenario because the depth of oxygen depletion is set by sedimentation rates averaging over the timescale of OET rather than responding to changes in sedimentation rate instantaneously (see Text S1 in Supporting Information S1).
Variability in sediment accumulation can arise from transient periods of nondeposition (i.e., stasis) and/or periods of erosion (Ganti et al., 2011;Tipper, 2015), but the relative importance of these two drivers of rate variability remains incompletely understood (Paola et al., 2018) and likely varies between environments (Sadler & Jerolmack, 2015).
With the ultimate goal of comparing the effects of nondeposition and erosion on OC preservation, we start with model simulations with only sediment accumulation and no erosion. Adopting the approach in Trampush and Hajek (2017), time series of sedimentation rates with annual resolution are generated through random sampling from a set probability distribution. For deposition-only simulations, we use two distributions: exponential and generalized Pareto (Ganti et al., 2011;Straub et al., 2012), both of which are right-skewed such that most random samples are at near-zero sedimentation rates, which, in effect, produces periods of nondeposition (Figure 2a, Table S4 in Supporting Information S1). The sedimentation rate time series are then used for time-dependent numerical solutions to the reactive transport model (Equation 1, see Text S1 in Supporting Information S1), which keeps track of temporal evolution of OC content within the sediment column.
In Figure 2c, we show the stratigraphic variability in OC preservation corresponding to the sedimentation rate time series in Figure 2a with the same mean sedimentation rate 𝐴𝐴 ( ω𝜔) . Stochasticity in deposition alone introduces temporal variations in the amount of OC preserved in the sedimentary records. And the magnitude of this variability scales with the variability in sedimentation rates (Figure 3b). These autogenic signals are robust as they persist when smoothing is applied (Figure S1 in Supporting Information S1). A smoothing window of 50 cm should be sufficient to simulate the expected effects of bioturbation and compaction, which are not included in our model. The mean value of OC burial efficiency calculated from the models with rate variability is the same as the OC burial efficiency for a system with the same mean sedimentation rate, but no rate variability (i.e., steady state, Figure 3a). So, in this case, accounting for sedimentary dynamics does not impact the magnitude of OC burial, but instead generates temporal patterns that might affect proxy record interpretation. The noise in OC (Algeo et al., 2008;Sageman et al., 2003;Seki et al., 2019;Tada et al., 1999). Error bars represent one standard deviation. content would also translate into spurious signals with amplitudes of 1-2 permil in δ 13 C isotope records (Figure S2 in Supporting Information S1) if we consider isotopic fractionation during OC respiration (Wynn, 2007).
To simulate sediment accumulation in systems with transient erosion, we use a truncated double Pareto distribution (Ganti et al., 2011;Straub et al., 2012;Trampush & Hajek, 2017). This heavy-tailed distribution allows for negative sedimentation rates (i.e., erosion), but is centered on a positive mean sedimentation rate such that there is net sediment accumulation in the long term. Two distributions of sedimentation rates with different shape and truncation parameters (Ganti et al., 2011;Trampush & Hajek, 2017) but the same mean are shown in Figure 2b and the corresponding simulated stratigraphic variability in OC preservation in Figure 2d. In contrast to simulations with deposition only, we find that the mean OC burial efficiency is greater in systems with erosion than predicted based on a constant mean sedimentation rate (Figure 3a). This occurs as erosion substantially skews the distribution of OET: the mean OET is much shorter compared with the deposition-only scenarios (Figure 3c). The distributions are also heavy-tailed (Figure 3c), giving rise to the greater stratigraphic variability in OC preservation.
It may be counterintuitive that erosion increases OC burial efficiency as postdepositional mobilization of sediments is considered to promote OC degradation (Aller, 1998;Aller & Blair, 2006). However, this discrepancy may result from a difference in the scale of the system under consideration. Specifically, OC burial can be considered at the scale of a single core or at the scale of an entire depositional basin. These measures can be different as the material eroded from a single point must be redeposited elsewhere to conserve mass. The fate of this eroded material is not considered in the type of 1-D models we use here, which represent a single sediment core. However, eroded sediment must be taken into account when considering OC burial at the scale of an entire depositional basin.
Locally enhanced burial (i.e., 1-D core-scale) and regionally enhanced oxidation (basin scale) are not mutually exclusive. The model results we show in Figure 3 display higher burial efficiencies compared to the constant sedimentation rate case because erosion shortens the transit time of OC through the diagenetic zone. If the eroded sediment is redeposited elsewhere within the same basin, then the effect at the basin scale is to increase transit times as material eroded from depth is redeposited at the surface and further exposed to oxidants. However, if the eroded material is exported from the system (e.g., from the shelf to slope), then OET at the basin scale will be similar to OET simulated at the core scale given that there is little to no deposition of re-eroded material. Likely, different systems will display different degrees of these two end-member scenarios. Moreover, the exact behavior will depend upon how the naturally blurry boundaries between depositional systems are defined. Accordingly, to apply any inferences from our modeling approach to natural systems, it is necessary to consider how the specific depositional processes control the distribution of OET, which may be highly variable between systems.
As our model is generic and not aimed at reproducing any particular sedimentary OC record, we use spectral analysis to compare our model predictions with field observations instead of trying to find matching patterns in OC profiles. Model simulations and real records of sedimentary OC variability are transformed to signals in the frequency domain (i.e., spectral analysis, see Text S1 in Supporting Information S1), which describe the distribution of energy among different frequencies. Sedimentary OC records from IODP cores (Seki et al., 2019;Tada et al., 1999) and sedimentary rock sections (Algeo et al., 2008;Sageman et al., 2003) are chosen to compare with model-simulated stratigraphic OC variability as they preserve relatively long or high-resolution records. Additional model simulations that consider an OC reactive continuum (Boudreau & Ruddick, 1991) that in principle captures the complex and potentially more realistic kinetics of OC oxidation were performed to avoid model artifacts (see Text S1 in Supporting Information S1). The power spectra of both model simulations and observations exhibit power law scaling with negative exponents commonly found in nature, and the estimated spectral slopes of simulations are in the same range with those of field data (Figure 3d and4b). Moreover, there appears to be a better match between observations and model simulations with transient erosion (Figure 3d), although more definitive interpretation could only be made with additional, higher-resolution data sets. Nevertheless, the fact that our model simulations can reproduce similar spectral slopes of observations highlights the relevance of stochasticity in early diagenesis and the manifestation of sedimentary dynamics in stratigraphic OC records.
Similar to purely physical models of autogenic behavior (Jerolmack & Paola, 2010), a spectral slope break is present at a high frequency corresponding to the inverse of the mean OET in our model simulations (Figure 4a). This is not observed in the power spectra of observed OC records, likely because field data with comparable resolution is effectively absent when considering analytical noise (Figure 4b). The spectral slope break does not affect the spectral slope estimation as it appears at a very high frequency where power spectral density is vanishingly small, but might be a potential proxy on past diagenetic conditions given field data with high enough resolution.
Real depositional systems are modulated by a combination of internal dynamics and external forcings. As a result, periodic signals from climatic/environmental forcings may be modified by the internal dynamics we simulate, leading to the preservation of convoluted signals in the rock record. Oscillating boundary conditions with different periodicities imitating external forcings (see Text S1 in Supporting Information S1) produce sharp spectral peaks at corresponding frequencies in simulations with constant sedimentation rates, whereas the spectral peaks are dampened or concealed in models with stochastic sedimentation variability (Figures 4c andS3 in Supporting Information S1).
Mixed evidence exists for the origin of sedimentation rate variability and hence the incompleteness of stratigraphic records. Sedimentary structures in fine-grained strata (Trabucho-Alexandre, 2014) and highly variable focusing factors in pelagic environments (Anderson, 2003) suggest erosion to be a potential driver for the incompleteness (Dutkiewicz & Mülle, 2022), while other lines of evidence points to stasis being the primary cause (Tipper, 2015). Sequence stratigraphy has limited utility in resolving this issue due to the preservation of mostly ordinary conditions in strata (i.e., "strange ordinariness," (Paola et al., 2018)). Mechanistically, the "strange ordinariness" has been suggested to be caused by hierarchies in the forms of surface topography (Ganti et al., 2020;Paola et al., 2018), and this mechanism is not inconsistent with either erosion-or hiatus-driven gaps. Our approach may provide a new constraint: we simulated the OC variability given the same age-depth model but different sedimentation rate time series where gaps result solely from nondeposition and erosion, respectively (Figure 5a). The distinct patterns of OC burial efficiencies (and potentially other chemical tracers) near unconformities (Figure 5b) may help distinguish between the physical mechanism(s) responsible for the missing time in the rock record.
Our findings suggest that stochastic variability in sedimentation as a result of the internal dynamics of sedimentary systems coupled with diagenetic reactions can generate autogenic signals. Such signals may be prevalent and persistent in sedimentary OC records, obscuring signals produced by external environmental forcings. Therefore, care needs to be taken when interpreting sedimentary OC records, and likely, any other geochemical proxy record that is sensitive to diagenetic transformations, as patterns can emerge without environmental perturbation. Accordingly, our approach could be used to construct mechanistic null hypotheses for improved interpretations of sedimentary proxy records and to infer past depositional conditions. (Seki et al., 2019). (c) Power spectra of model simulations with oscillating boundary conditions as periodic external forcings. Spectral peak prominent in the steady-state simulations with a constant sedimentation rate is concealed when stochasticity in sedimentation rates is incorporated in the model.