Modeling Organic Carbon Accumulation Rates and Residence Times in Coastal Vegetated Ecosystems

Modeling Organic Carbon Accumulation Rates and Residence Times in Coastal Vegetated Ecosystems E. Fay Belshe1, Jose Sanjuan2, Carmen Leiva-Dueñas3, Nerea Piñeiro-Juncal3,4, Oscar Serrano5, Paul Lavery5, and Miguel Angel Mateo3,5 1Leibniz Center of Tropical Marine Research, Fahrenheitstraße, Bremen, Germany, 2Deutsches Zentrum Für Luft- und Raumfahrt, Bremen, Germany, 3Centro de Estudios Avanzados de Blanes, Consejo Superior de Investigaciones Cientificas, Blanes, Spain, 4Departamento de Edafoloxía e Química Agrícola, Universidade de Santiago de Compostela, Santiago de Compostela, Spain, 5School of Science and Centre for Marine Ecosystems Research, Edith Cowan University, Joondalup, Western Australia, Australia Abstract Coastal vegetated “blue carbon” ecosystems can store large quantities of organic carbon (OC) within their soils; however, the importance of these sinks for climate change mitigation depends on the OC accumulation rate (CAR) and residence time. Here we evaluate how two modeling approaches, a Bayesian age-depth model alone or in combination with a two-pool OC model, aid in our understanding of the time lines of OC within seagrass soils. Fitting these models to data from Posidonia oceanica soil cores, we show that age-depth models provided reasonable CAR estimates but resulted in a 22% higher estimation of OC burial rates when ephemeral rhizosphere OC was not subtracted. This illustrates the need to standardize CAR estimation to match the research target and time frames under consideration. Using a two-pool model in tandem with an age-depth model also yielded reasonable, albeit lower, CAR estimates with lower estimate uncertainty, which increased our ability to detect among-site differences and seascape-level trends. Moreover, the two-pool model provided several other useful soil OC diagnostics, including OC inputs, decay rates, and transit times. At our sites, soil OC decayed quite slowly both within fast cycling (0.028 ±0.014 yr−1) and slow cycling (0.0007 ±0.0003 yr−1) soil pools, resulting in OC taking between 146 and 825 yr to transit the soil system. Further, an estimated 85% to 93% of OC inputs enter slow-cycling soil pools, with transit times ranging from 891 to 3,115 yr, substantiating the importance of P. oceanica soils as natural, long-term OC sinks. 1. Introduction Coastal vegetated ecosystems (seagrasses, tidal marshes, and mangroves) occupy a very small area of Earth's surface (<1%), yet play an outsized role in the global carbon cycle by burying an estimated 0.08–0.22 Pg C yr−1, which is equivalent to approximately 10% of the residual land sink (Duarte et al., 2013; National Academies of Sciences, 2019; Spivak et al., 2019). Moreover, these “blue carbon” ecosystems are widely reported to store substantial quantities of organic carbon (OC) within their soils (Brindgham et al., 2006; Donato et al., 2011; Fourqurean et al., 2012; Lavery et al., 2013; Macreadie et al., 2014; Rovai et al., 2018). However, to better understand the potential of these ecosystems as sinks for atmospheric CO2, robust estimations of global OC accumulation rates are needed. Further, the ability of a reservoir to act as a relevant sink for climate change mitigation depends on (a) the time scale OC resides within it, with the number of years that carbon is stored away from the atmosphere equaling the number of years green house gas emissions are delayed, and (b) the time horizon chosen for mitigation accounting (Levasseur et al., 2012). The longer carbon is stored relative to the period of time its impact is considered, the greater the benefits of delaying its release (Levasseur et al., 2012). Under the most common time horizon (100 yr) used for carbon reporting under the United Nations Framework Convention on Climate Change (UNFCCC), OC stored for decades or longer becomes increasingly relevant (Levasseur et al., 2012; Marland et al., 2001; Trumbore, 2000). This necessitates the robust estimation of accumulation rates and residence times of soil OC within pools that persists over these longer time frames. Currently, the accumulation rate of OC within soils of blue carbon ecosystems is widely estimated using age-depth models (Chmura et al., 2003; Breithaupt et al., 2012; Jankowska et al., 2016; Mazarrasa et al., 2017; Rozaimi et al., 2016; Ouyang & Lee, 2014; Samper-Villarreal et al., 2018; Serrano, Ricart, et al., RESEARCH ARTICLE 10.1029/2019JG005233 Key Points: • Conservative long-term carbon burial rates can be estimated from soils below the rhizosphere • Estimate uncertainty with pool models increases the ability to detect seascape-level trends • Pool model estimates of soil carbon diagnostics help to unveil dynamics governing carbon storage Supporting Information: • Supporting Information S1 •DataSetS1 Correspondence to: E. F. Belshe, [email protected] Citation: Belshe, E. F., Sanjuan, J., Leiva-Dueñas, C., Piñeiro-Juncal, N., Serrano, O., Lavery, P., & Mateo, M. A. (2019). Modeling organic carbon accumulation rates and residence times in coastal vegetated ecosystems. Journal of Geophysical Research: Biogeosciences,124, https://doi.org/10.1029/2019JG005233 Received 1 MAY 2019 Accepted 30 SEP 2019 Accepted article online 1 NOV 2019 Author Contributions Conceptualization: E. Fay Belshe, Oscar Serrano, Paul Lavery, Miguel Angel Mateo Data curation: E. Fay Belshe, Jose Sanjuan Funding Acquisition: Miguel Angel Mateo Methodology: E. Fay Belshe, Jose Sanjuan Validation: Oscar Serrano, Paul Lavery Writing - Original Draft: E. Fay Belshe Formal Analysis: E. Fay Belshe, Jose Sanjuan ©2019. The Authors. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made. BELSHE ET AL. 3652 3652 3671. – Published online 28 NOV 2019 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 2016; Wilkinson et al., 2018; Wright et al., 2017). In general, these approaches model sedimentation rates (or mass accumulation rates) either implicitly or explicitly from dated layers within the soil profile (Appleby & Oldfield, 1978; Blaauw & Christen, 2011; Blaauw, 2010; Haslett & Parnell, 2008; Ramsey, 2008; Trachsel & Telford, 2016). Once soil (or mass) accumulation rates are estimated from a given age-depth model, the OC accumulation rate (CAR) is calculated by simply multiplying the OC density (or %OC) by the rate from the corresponding depth of the soil profile. The mean CAR for a particular site is strongly dependent on the thickness (time frame) of the soil profile used for its estimation. Depending on the research target, CARs can be estimated from, or including, surface soil layers to provide relevant information about recent (or new) rates of OC accumulation within an ecosystem (Marba et al., 2015; Mazarrasa et al., 2017; Samper-Villarreal et al., 2018). Alternatively, CARs adjusted to account for ephemeral OC mineralization or estimated from deeper soil layers where enough time has elapsed since deposition for faster-cycling OC to decay, can serve as a proxy for OC burial rates, which represent the rate OC is removed from the short-term, biosphere-atmosphere carbon cycle (Gälman et al., 2008; Sobek et al., 2009). Without subtracting the contribution of the faster-cycling OC pool, burial rates within blue carbon ecosystems will be systematically overestimated when using age-depth models. Again, the research project must clarify the target and time frame used to subjectively differentiate CARs from longer-term OC burial rates. Additionally, age-depth models estimated from dating the mineral components of the soil (such as 210Pb) or using closed-system materials for 14C dating, do not directly provide information about the age of soil OC (Trumbore et al., 2016). They simply provide the time frame when the mineral portion of the soil was deposited. Although this provides useful information for reconstructing past events (López-Merino et al., 2015; Mazarrasa et al., 2017; Serrano et al., 2011), the age of OC contained within a particular envelope of the soil profile can contain a large uncertainty when roots are present (Trumbore, 2009). Depending on the site-specific sedimentation rates, rooting depths, and OC residence times, OC ages within a soil layer can range from very recent (the current year's root exudates and litter) to the time of original deposition of the soil layer, which can be centuries ago. In efforts to achieve a better understanding of the temporal OC dynamics within soils of blue carbon ecosystems, an alternative modeling approach commonly utilized in terrestrial systems can be coupled with age-depth models, albeit with some assumptions and simplifications. These models are based on a simplified conceptualization of OC as a suite of pools or compartments and can be used to estimate two diagnostic characteristics that define the soil's effectiveness to retain OC, the transit time and age (Manzoni et al., 2009). The transit time (T) is defined as the time a soil OC molecule remains in the soil since deposition until it is remineralized and can be interpreted as the time it takes for OC to transit the soil system (Sierra et al., 2017). The age (𝜏) describes the age distribution of OC molecules within the soil system at a given time (Eriksson, 1971; Manzoni et al., 2009; McGuire & McDonnell, 2006). The simplest pool model describes the change in OC over time of a single pool as a result of constant inputs (I) and first-order linear outputs: dx/dt=I−kx(t), where kis the first-order kinetic rate constant, in units of inverse time (Manzoni et al., 2009). To adapt pool models to data from soil cores, they can be fit to the cumulative OC density (x) down the soil profile versus time (t) estimated from an age-depth model, with the assumptions that the system is in steady state, OC is accreted vertically, and recent OC dynamics (top of core) are representative of initial OC dynamics of the system (Clymo, 1984; Trumbore & Harden, 1997). The goal of this work is to evaluate how these two modeling approaches (age-depth models alone or in combination with a two-pool model) aid in our understanding of the time lines of OC within soils of coastal vegetated ecosystems. Specifically, we explore the use of the Bayesian age-depth model Bacon to estimate the mean CAR for different sections of the soil profile (whole core, rhizosphere, and subsoil), and a pool model to estimate the age, transit time, and CAR of two discrete OC pools decaying in parallel. By critically analyzing and comparing model estimates, along with their associated uncertainty, we aim to provide a strategy on how to better estimate the CAR, burial rate, age, and transit time of OC within soils of blue carbon ecosystems. In addition, we provide recommendations for auxiliary data collection to aid in decisions for delineating the soil profile for OC burial rate estimates. In this work, we present analysis and data from seagrass soils and compare our results to findings from tidal marshes and mangroves. In particular, we analyze soils of Posidonia oceanica, an endemic, long-lived seagrass from the Mediterranean Sea, which forms stable and continuous meadows (Duarte et al., 2006). P. oceanica ranks among the most productive of all marine ecosystems and accumulates organic-rich soils underneath the meadows over thousands of years (Mateo et al., 1997). Investigation: Oscar Serrano, Paul Lavery, Miguel Angel Mateo Resources: Jose Sanjuan, Miguel Angel Mateo Supervision: Miguel Angel Mateo Visualization: E. Fay Belshe Writing - review & editing: Jose Sanjuan, Oscar Serrano, Paul Lavery, Miguel Angel Mateo BELSHE ET AL. 3653 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 Figure 1. Sites where soil cores were extracted within Santa Maria Bay and Es Port Bay, located within the Cabrera Archipelago of the Balearic Islands of Spain. The benthic habitat characteristics within the Bays is indicated, including the distribution of seagrass (P. oceanica and Cymodocea nodosa)meadows. 2. Materials and Methods 2.1. Soil Core Collection and Processing For this study, we analyzed soil cores collected within seagrass meadows from two Bays of Cabrera Island: Santa Maria Bay and Es Port Bay (Figure 1). Cabrera is the largest island of the Cabrera Archipelago, located a few kilometers south of Mallorca (Balearic Islands, western Mediterranean). Soil cores were collected at sites along a water depth gradient (5, 10, 15, and 25 m) in Santa Maria and at one site (5 m water depth) in Es Port (Figure 1). A single, hand-driven pvc core (7 cm diameter) was collected at each site within P. oceanica meadows during June 2015 on SCUBA. The length of extracted cores was dictated by the ability to penetrate the soil and ranged from 90 to 161 cm. All cores were fitted with a “core catcher” to minimize soil loss during extraction, and corers were slowly hammered into the soil with gentle rotation to minimize soil compression. After full insertion, soil depth was measured inside and outside of the core to calculate soil compression during coring. Once extracted, cores were sealed, transported to the lab, and stored frozen in a vertical position until processed. In the laboratory, cores were opened by cutting the core tubes on each side with an electric saw and then split longitudinally into two halves using a ceramic knife. Opened cores were photographed with a high-resolution color-line scan camera and then one hemicore was subdivided into 1cm slices. Soil subsections (slices) were weighted before and after drying at 60 ◦C until constant weight, and soil dry bulk density (DBD; g cm−3) was calculated as the dry weight of soil per volume of each sample. The volume of each 1 cm slice was calculated for each hemicore's specific dimensions, to account for any uneven longitudinal cutting of the cores. After drying, soils from the 1 cm slices were homogenized and further split into three subsamples for 14C and 210Pb chronology and OC determination. Every centimeter for the first 20 cm of each core, and every other centimeter along the core thereafter, was selected for OC analysis. These subsamples were fully homogenized by grinding in an agate mortar (Mortar Grinder RM 200). 2.2. OC Content Determination, and OC Density and Stock Calculations Homogenized soil samples were sent to the University of Hawai'i at Hilo's Analytical Laboratory where the %OC and 𝛿13C was determined, after acidification with 1 M HCL to remove carbonates, on an elemental BELSHE ET AL. 3654 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 analyzer (Costech) coupled to an isotope ratio mass spectrometer. OC density (g OCcm−3) of each 1 cm slice was calculated from the measured %OC and the dry bulk density (g cm−3). An exponential decompression function was applied to correct for the gradual shortening (compression) of the cores during coring (Morton & White, 1997). All core depths reported throughout this study are decompressed depths (cm). Because OC was measured in every other centimeter after the top 20 cm of the core and compression corrections created small gaps in data along the core, missing OC density values were linearly interpolated from adjacent values. In addition, one core only extended 90 cm into the soil. To gap fill these final 10 cm, a negative exponential model was used, which allowed us to estimate OC stocks for all cores to at least 1 m. The amount of OC stored within the soil was calculated by summing the OC density in each depth increment (cm) down to a standardized depth of 1m. 2.3. Radio-Isotope Dating To construct a chronology of the soil profile for each site, both radiocarbon (14C) and lead-210 (210Pb) were measured. 14C-based chronology was built from samples of P. oceanica sheath debris, which were selected from four to nine depths spaced along each core. The materials used for 14C dating were first rinsed in Milli-Q water to remove fine soil particles, then placed into an ultrasonic bath of Mili-Q water for 5 min to further remove attached allochthonous particles, and finally inspected under a binocular microscope for any remaining particles. If particles were present, the ultrasonic bath was repeated until no particles could be seen. Samples were then sent to DirectAMS (https://www.directams.com) and measured with an accelerator mass spectrometer following standard procedures (Karlen et al., 1964; Stuiver et al., 1998). The top 30 cm of soil from each site were also dated using the 210Pb method. Total 210Pb activity was determined from its granddaughter 210Po, assuming radioactive equilibrium between both radionuclides, using alpha spectrometry (Appleby & Oldfield, 1978) measured on a PIPS detector (CANBERRA, Mod. PD-450.18 AM) at the Laboratori de Radioactivitat Ambiental de la Universitat Autònoma de Barcelona. The age of soil layers were estimated using the Constant Rate of Supply (CRS) model (Appleby & Oldfield, 1978; 2004). This model assumes that the flux of excess 210Pb onto the soil surface is constant over time and is exactly fit (“mapped”) to the measured profile of excess 210Pb to estimate the chronology along the soil profile (Mabit et al., 2014). 2.4. Age-Depth Models: Soil Accumulation Times, Rates, and Ages To model the age-depth relationship of each core, we utilized the Bayesian modeling approach Bacon, which explicitly models the sediment accumulation rate to establish a coherent evolution of deposition along the soil profile (Blaauw & Christen, 2011). Specifically, the accumulation rates of linear segments along the profile are controlled with a gamma autoregressive semiparametric model, and a t-walk (Christen & Fox, 2010), self-adjusting, Markov Chain Monte Carlo (MCMC) sampling algorithm is used to estimate the accumulation rate (yr cm−1) for each segment (Blaauw & Christen, 2011). From each Bacon model run, three data sets are generated, which we used in subsequent analysis. First, calibrated 14C values along with their error distribution are given for each 14C date entered. Second, the data forming the age-depth model, consisting of the chronology (ages in years) and their estimated probability density function (PDF) for each depth (cm) along the soil core is output. Third, the accumulation rate (yr cm−1) and its associated PDF for each segment of the core is output. We fed data consisting of both 210Pb calendar dates (±1𝜎) estimated with the CRS model and measured 14C dates (±1𝜎) from each core into Bacon (version 2.3.4) in R (Blaauw & Christen, 2011). Input dates were expressed as years before present (BP), where present (year 0 BP) corresponds to 1950. To calibrate 14C dates we used the marine 13 calibration curve (Reimer et al., 2013). Information about the reservoir residence time (𝛥R) for the water bodies surrounding our study sites were also incorporated into the 14C calibration, to correct for the marine reservoir effect (Stuiver & Reimer, 1993). As no 𝛥Rvalues are reported for the Balearic Islands, we assumed that water masses from the east coast of Spain and from the Algerian coast (brought by mesoscale eddies) are representative of our study sites (Jordi et al., 2006; 2009; Robinson et al., 2001). Therefore, we used a reservoir residence time correction of 𝛥R=76 ±59 yr BP, which is the mean from waters off the coast of Northern Spain and Algeria (Reimer & McCormac, 2016; Siani et al., 2016). For each core, we provided information about (a) the minimum year of each core, (b) the accumulation rate prior mean and shape, (c) the memory prior mean and shape, (d) the thickness of the segments, and (e) core length (Table S1 in the supporting information). Knowledge about the hydrology, distance from the coast and past climatic events in the region aided in our choosing of the priors. BELSHE ET AL. 3655 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 2.4.1. Age-Depth Model CARs The Bacon-estimated sediment accumulation rates were used to calculate the CAR for each core. Since the Bacon accumulation rate is given in yrcm−1(which is actually the accumulation times), we converted from yr cm−1to cm yr−1. First, the posterior PDF of the accumulation times were converted from a Gamma distribution to an inverse Gamma distribution, that is, 𝑓i(r)= 𝛽𝛼i i Γ(𝛼i)r−𝛼i−1e−𝛽i∕r(1) where 𝛤(·) denotes the Gamma function, istands for the segments along the core, and ris the sediment accumulation rate in units of cmyr−1. The parameters 𝛼and 𝛽are the shape and scale parameters, respectively, which are provided by the Bacon MCMC outputs. However, we are interested in CAR in units of gOCm −2yr−1. This is calculated by multiplying the accumulation rate (r) by the OC density (𝜌iin units of gOCcm −3) of the segment and by applying a factor of 104to convert from cm3to m2(for 1 cm height segments), so all rates are reported for the areal extent of a square meter. This is equivalent to scaling 𝛽i,r by 104𝜌iin equation (1). In order to estimate the confidence intervals (CI) of the CAR for each segment, the cumulative density function (CDF) is computed, which for an inverse-Gamma distribution is Fi(r)=Γ(𝛼i,𝛽 i∕r) Γ(𝛼i)(2) where 𝛤(· ,·) is the upper incomplete Gamma function and now rcorresponds to the CAR instead of the sediment accumulation rate. The mean CAR (𝜇i,), the median, and 95% CIs for every segment are calculated from the CDF. The CAR uncertainties only include those associated with the rates and not the ones associated with the OC density (𝜌) since replicate cores for each site were not taken. 2.4.2. Age-Depth Model Mean CARs We estimated the mean CAR (g OC m−2yr−1) for the entire length of the core (whole core) and for different sections of the core, the rhizosphere, and the subsoil. This delineation was done to compare rates with and without separating, and subtracting, the contribution of the fast-cycling OC from the slower cycling pools in the rhizosphere. Two approaches were used to delineate the rhizosphere from the subsoil. First, we used the top 25 cm of the soil profile to estimate the rhizosphere mean CAR and used the remaining portion of the soil (>25 cm) for the subsoil mean CAR. This standard depth delineation was chosen because the maximum root biomass of seagrasses is often found between 10 and 20 cm depth (Duarte et al., 2005), which we adjusted slightly downward to 25cm to encapsulate the portion of the rhizosphere where the majority of roots occur, since we did not have explicit data on the vertical distribution of root density. Second, we divided each core based on the maximum rooting depth of the seagrass at the site. To determine the maximum rooting depth, we utilized the high-resolution color images of each core, along with data on the %OC, 𝛿13C, and DBD of the core. Although there were not always clear trends in %OC or DBD, in general we found that the rooting depth determined from the images corresponded to the depth of the core where DBD stabilized (Figure 2). The delineation of the maximum rooting depths were further supported by core-specific trends in 𝛿13C, where enriched values indicated the presence of living roots, which for our sites fell between −12.4 to −14% 0 (Figure S1 in the supporting information). The area from the soil surface to the maximum rooting depth was considered the rhizosphere, and all soil below was considered the subsoil. The Bacon estimated mean CAR ( 𝜇) for each core section (whole core, rhizosphere, and subsoil) was calculated as 𝜇=1 N N ∑ i 𝜇i,(3) where Nis the total number of centimeters within a section. The PDF of 𝜇is the convolution of the PDFs of each segment i, which are an inverse-Gamma distribution (equation (1)) with 𝛽iscaled by 1∕N. The CDF of 𝜇is calculated from its PDF in order to estimate the median and the 95% CIs. 2.5. Pool Models: CARs, Transit Times, and Ages To estimate time lines of OC along with CARs, we utilized a pool model consisting of two parallel compartments, where each pool is fed by a fraction of OC inputs and decay independently according to first-order kinetics (Manzoni et al., 2009), which can be described by the following system of ordinary differential equations, dx1(t) dt=𝛼I(t)−k1x1(t)(4) BELSHE ET AL. 3656 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 Figure 2. Depth-specific soil properties (% organic carbon, OC; dry bulk density; OC density) and timespan (years) estimate with the Bacon age-depth model, with year zero corresponding to the year the cores where extracted (2015). Arrows indicate the centimeter where dated materials were sampled and are colored according to dating method, with gray arrows corresponding to 210Pb dates, the black arrows to 14C dates, the hollow gray arrow to the core extraction date for site SM25. The shaded green region corresponds to the maximum rooting depth of seagrass at each site. Note the difference in the Timespan yaxis for core SM05. BELSHE ET AL. 3657 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 dx2(t) dt=(1−𝛼)I(t)−k2x2(t)(5) which have the following solution: xglobal(t)= x1(t)+x2(t) =I0[𝛼 k1 (1−e−k1t)+1−𝛼 k2 (1−e−k2t)](6) where xis OC density (g OC cm−3), 𝛼is the portion of inputs I(t)entering the pools, and kiis the linear decay rate of each pool. Time (t) is the median chronology (in years) estimated by the Bacon age-depth model. For each core, the parameters 𝚯(𝛼,k1,k2, and I0) were estimated by minimizing the negative log-likelihood function: (Θ) = 𝛜(Θ)Σ−1𝛜(Θ)T(7) with 𝛜=xmeas −xglobal (8) where xmeas and xglobal are the data and the model, respectively. The bold fonts indicate that the variables are vectors or matrices. xmeas is the cumulative sum of 𝜌(OC density): xmeas(tk)= k ∑ i=0 𝜌(ti).(9) with k=0,1,…,N−1, where Nis the total number of data points. The covariance matrix, ˝, is not a diagonal one because the PDFs are correlated since xmeas is generated by the cumulative sum of OC density (see equation (9)). ˝is composed of two terms: (i) the OC density uncertainty (𝜎2 𝜌) and (ii) the time uncertainty (𝜎2 t). Because we did not have replicate cores for each site, we defined the OC density standard deviation as 10% of its value at each centimeter along the core and assumed it to be Gaussian. We choose 10% based on the variability in OC density or %OC for replicate cores taken from P. oceanica soils reported in other studies (Mazarrasa et al., 2017; Serrano et al., 2016) 𝜎𝜌(tk)=0.1𝜌(tk).(10) The second term, the time uncertainty, is converted to OC uncertainty as 𝜎2 OC,t(tk)=(dxmeas(t) dt||||t=tk)2 𝜎2 t(tk),(11) where the derivative is calculated numerically and smoothed by applying a low-pass filter. 𝜎2 tis the variance of the chronology provided by the Bacon age-depth model, which follows a Gaussian PDF with different variances for each core segment. Thus, the covariance matrix is Σ=Σ 𝜌+Σ OC,t(12) where ˝𝜌is the symmetric matrix: .Σ𝜌= ⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ 0 ∑ i=0 𝜎2 𝜌(ti) 0 ∑ i=0 𝜎2 𝜌(ti)… 0 ∑ i=0 𝜎2 𝜌(ti) 1 ∑ i=0 𝜎2 𝜌(ti)… 1 ∑ i=0 𝜎2 𝜌(ti) … 2 ∑ i=0 𝜎2 𝜌(ti) ⋱⋮ N−1 ∑ i=0 𝜎2 𝜌(ti) ⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (13) BELSHE ET AL. 3658 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 and ΣOC,t=⎛⎜⎜⎜⎜⎝ 𝜎2 OC,t(t0)0…0 𝜎2 OC,t(t1)… 0 ⋱⋮ 𝜎2 OC,t(tN−1) ⎞⎟⎟⎟⎟⎠ .(14) The uncertainties in the parameter estimates were determined using Monte Carlo simulations, where 10,000 realizations were performed. Each realization sampled the random variables just described. Our goal was to model OC dynamics for two parallel pools to conceptually reflect a pool that decays quickly (fast pool) within the biologically active rhizosphere, and a pool with a slower decay rate (slow pool) that escapes remineralization in the short-term and continues to reside within the subsoil. However, for comparison, we also fit data using the same methods described above to the simpler one-pool model: x(t)=I0 k(1−e−kt).(15) For each core, one- and two-pool models were compared using Akaike Information Criterion (AIC), and a model was deemed to be superior if its AIC score was at least two AIC units lower than the other candidate model (Anderson et al., 1998). 2.5.1. Pool Model Transit Times and Model Ages Using the parameters (I0,𝛼,k1, and k2), we can calculate the mean transit time of OC through soil system from the transient time distribution, that is, the distribution of times for OC molecules to reach the system exit, which is equivalent to the output flow resulting from an impulse input. The mean transit time of the whole soil system (global) is (Manzoni et al., 2009)  T=1 k1k2 [(1−𝛼)k1+𝛼k2].(16) Likewise, the parameters (I0,𝛼,k1, and k2) can be used to estimate the age distribution, which is defined as the amount of OC entered at time t−𝜏that have transit times greater or equal to its age (𝜏) over the total mass at time t. For the two compartments in parallel, the mean age is (Manzoni et al., 2009) 𝜏 =1 k1 +1 k2 +1 (1−𝛼)k1+𝛼k2 .(17) For the case of individual systems, that is, when considering only one of the pools (either the fast or slow), the mean transit time and mean age are  T=𝜏 =1 k.(18) ThePDFsofthe mean transienttime and themeanage werecalculated with the 10,000parameters estimated during the Monte-Carlo realizations. 2.5.2. Pool Model CARs Using the parameter estimates (I0,𝛼,k1, and k2), the CAR is calculated as the time derivative of equation (6), that is, 𝜇(t)=I0[𝛼e−k1t+(1−𝛼)e−k2t],(19) which, in turn, can be split into the two respective pools: 𝜇fast (t)=𝛼I0e−k1t,(20) 𝜇slow (t)=(1−𝛼)I0e−k2t(21) BELSHE ET AL. 3659 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 The mean CARs for the global system and for the individual pools are calculated as 𝜇=1  T∫ T 0 𝜇(t)dt =I0 (1−𝛼)k1+𝛼k2[𝛼k2(1−e(𝛼−1)k1∕k2−𝛼)+(1−𝛼)k1(1−e𝛼(1−k2∕k1)−1)](22) 𝜇fast =1  Tfast ∫ Tfast 0 𝜇fast (t)dt=I0𝛼(1−e−1)(23) 𝜇slow =1  Tslow ∫ Tslow 0 𝜇slow (t)dt=I0(1−𝛼)(1−e−1)(24) where we only take into account the predicted data until the mean transient time (  T) is reached for each respective pool. The PDFs of the mean CARs are derived from the 10,000 Monte-Carlo simulations. For the one-pool model case the CAR reduces to equation (20) and the mean CAR to equation (23), both with 𝛼=1. 2.5.3. Modeling Procedures and Estimate Comparisons All Bacon age-depth models and subsequent analysis to estimate sedimentation times, sedimentation rates, CARs, and mean CARs were done in R (version: R-3.5.0-2018-04-23; R Core Team, 2018) using the packages rbacon (Blaauw & Christen, 2018), zipfR (Evert & Baroni, 2007), and zoo (Zeileis & Grothendieck, 2005), and graphics were produced using the ggplot2 package (Wickham, 2016). Pool models were fit and Monte Carlo simulations were run using Matlab R2015a (8.5.0.197613; MATLAB, 2015). To compare estimates from age-depth models (differences in mean CARs among sites or core sections within a site) and pool models (differences among sites in parameter estimates, mean CARs, ages, and transit times) we utilized the estimates' confidence intervals (Hector, 2015). Differences were determined to be significant when there were no overlap in 83% confidence intervals, which is an approximation method for assessing the significant difference among means at the p =0.05 level (Austin & Hux, 2002). 3. Results 3.1. Soil Properties and OC At our sites within Santa Maria Bay and Es Port Bay, the percentage of soil OC ranged from a maximum of 10.37% to a minimum of 0.38% and in general followed the trend of decreasing %OC with soil depth (Figure 2). The exception to this trend was SM05, in which %OC was variable with depth and exhibited no clear trend. This may be due to the short length of the SM05 core (85 cm), which was not long enough to extend below the maximum rooting depth of the seagrass at this site. At SM25, %OC within the rhizosphere was comparatively lower than those in other sites, with a maximum of only 2.3%. The seagrass maximum rooting depth decreased with increasing water depth within Santa Maria Bay, ranging from 85 cm at SM05 (5 m water depth) to 26cm at SM25 (25m water depth). Rooting depth was 64 cm at EP05 (Figure 2). Bulk density generally increased with soil depth and ranged from 0.98 to 1.93 g cm−3(Figure 2). OC density ranged from 0.018 to 0.037 g OC cm−3, and was on average highest in soils at EP05 (Figure 2). The largest OC stocks in the top 1 m of soil were located at EP05, which contained an estimated 329 Mg OC ha−1. The second highest stock was found at SM25 with 227 Mg OC ha−1, followed by SM05 with 179 Mg OC ha−1. Soils at SM10 and SM15 stored similar amounts of OC, storing 99 and 98 Mg OC ha−1, respectively, which were comparatively lower than the other sites. 3.2. Age-Depth Models: Core Chronologies, Timespans, and CARs For our Bacon age-depth models, we input both 210Pb and 14C dates, except at site SM25 where sediment mixing produced an excess 210Pb profile that did not allow a CRS age-depth model to be calculated. Bacon-estimated core chronologies were nonlinear with respect to depth and site specific (Figure 2). The median timespan captured within each core (whole core), ranged from a maximum of 1,917 yr (95% CI: 1,515–2,156 yr) at SM15 to a minimum of 164 yr (95% CI: 136–198 yr) at SM05, with cores from EP05, SM10, and SM15 encompassing intermediate timespans (Table 1). These differences reflect both site-level differences in sedimentation rates and core lengths. The timespan captured in the top meter of soil did not exhibit any clear trend in relation to OC stocks in the same volume. Sites with similar OC stocks (SM10 =99 Mg BELSHE ET AL. 3660 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 Table 4 Pool Model Estimated OC Decay Rates for Fast and Slow-Cycling Soil Pools, Along With Literature Values of Soil OC and Aboveground (AG) and Belowground (BG) Litter Decay Within Blue Carbon Ecosystems Ecosystem/Type Decay rate k(yr−1) This study: mean; range Fast-pool (k1) 0.028; 0.0012–0.079 Slow-pool (k2) 0.0007; 0.0003–0.0013 Other studies: mean; range (citation) Slow-cycling soil organic matter Seagrass 0.0003; 0.00008–0.0005 (a) Tidal marsh 0.023; 0.005–0.042 (b) Mangrove 0.001 (c) Litter Seagrass AG 3.77; 1.06–7.3 (d) Seagrass BG 0.82; 0.11–2.37 (e) Tidal marsh AG 1.45; 0.40–5.91 (f) Tidal marsh BG 0.65; 0.29–2.55 (g) Mangrove AG 2.33; 0.98–5.48 (h) Mangrove BG 0.77; 0.36–1.27 (i) Note.Corresponding sources: (a) Serrano et al. (2012) and Mateo et al. (1997); (b) Rybczyk et al. (1998); (c) Chen and Twilley (1999); (d) Fourqurean and Schrlau (2003), Harrison (1989), Nicastro et al. (2012), Chiu et al. (2013), and Mateo and Romero (1996); (e) Fourqurean and Schrlau (2003), Kenworthy and Thayer (1984), and Romero et al. (1992); (f) Kirwan and Blum (2011), Hemminga and Buth (1991), and Christian (1984); (g) Hemminga and Buth (1991) and Kirwan et al. (2014); (h) Middleton and McKee (2001); (i) Middleton and McKee (2001), Morrisey et al. (2010), and Poret et al. (2007). the mixed-oxic conditions (0.042yr−1; Lovelock et al., 2017) typically found within the seagrass rhizosphere (Brodersen et al., 2017), and the deviation of our k1estimates from decay rates of above- and below-ground litter from seagrass and other blue carbon ecosystems (Table 4). The slow decay of soil OC reported here helps to explain the high OC stocks found in P. oceanica soils (Serrano, Lavery, et al., 2016; Fourqurean et al., 2012) and calls for further research quantifying OC decay rates under in situ soil conditions. Slow rates of decay translated into long residence times of OC within the soil. We estimated that OC within the global soil pool takes on average between 146 and 825 yr to transit the soil system at our sites (Table 2 and Figure 6). Further, an estimated 85% to 93% of OC inputs remain within slow-cycling pools, with transit times ranging from on average 791 to 3,115 yr (Table 2). These century to millennial estimates bolster previous estimates of residence times of OC within P. oceanica soils (Mateo et al., 1997; Serrano et al., 2012), and substantiate the importance of seagrass ecosystems for long-term OC storage. Our estimates of fast-pool OC transit times averaged between 13 and 44 yr, which is 10 times longer than the typical transit time (<1 yr) estimated for litter from seagrasses, mangroves, and tidal marshes (Enriquez et al., 1993). Although this further supports a decoupling between litter decay processes and soil OC dynamics, it points out a bias in our modeling approach toward longer time scales. Here, pool models rely on data of the remaining soil OC and age-depth model estimates of time, with a minimum 1 yr time step, which cannot capture the contribution of faster processes (<1yr). Therefore, our estimates of mean age and transit time are skewed toward longer time scales and should be interpreted with this limitation in mind. In addition, we found that data from soil cores did not support the use of the more complex two-pool model in two out of our five sites. This hindered us from differentiating global from slow pool CARs, which could be approximated using Bacon whole-core and subsoil estimates for all sites. BELSHE ET AL. 3667 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 In spite of these limitations, the ability to estimate transit times in addition to CARs enabled us to unveil some of the underlying dynamics potentially governing OC storage at our sites. First, the site with the highest OC storage (EP05) had a combination of fast CARs and long transit times. But a site can achieve high OC storage even with a very slow CAR if the time it takes OC to transit the system is long (SM25); and the fastest CAR does not necessarily equate to the highest stocks if OC transit times are short (SM05; Figure 6). These counteracting OC dynamics were further illustrated at SM10 and SM15, which contained very similar OC stocks (99 and 98 Mg OC ha−1, respectively) with contrasting CARs and transit times. These results fall in line with ecological equilibrium theory, which posits that soil OC persistence depends on the dynamic equilibrium between continuous losses (decay) and gains (inputs), and the extent to which protection mechanisms affect the rate of OC turnover (Caruso et al., 2018). Although this may be intuitive, it highlights the need to quantify OC transit times in blue carbon ecosystems because CARs alone only provide part of the picture. Finally, we found contrasting trends in CARs and OC transit times along the water depth gradient in Santa Maria Bay (Figure 6). The decrease in CARs along the depth-induced light gradient in tandem with the reduction in rhizosphere extent with depth (Figure 2) suggests that declining seagrass productivity (Serrano et al., 2014) and/or lower below-ground OC inputs reduces CARs within seagrass soils. On the other hand, the shortening of the biologically active rhizosphere may increase the likelihood for OC to escape remineralization (Ekschmitt et al., 2008) and perhaps leads to the long OC transit times at deep water sites. Although further investigation is required to pinpoint the factors driving the trends revealed here, these results could explain the lack of clear patterns in OC stocks along the depth-gradient in this study and others (Mazarrasa et al., 2017; York et al., 2018), and points to the potential importance of deep-water seagrasses for long-lasting OC storage. 5. Conclusion We put forth here a simple solution of using the average subsoil CAR as a conservative estimate of long-term OC burial rates and point out the need to standardize the estimation and reporting of CARs to fit the research target and time frame under consideration. We demonstrate the utility of pool models used in tandem with age-depth models, which can yield reasonable estimates of CARs with lower estimate uncertainty than age-depth models alone, thereby increasing the ability to detect among-site differences and seascape-level trends. In addition, pool models provide several other diagnostics of the soil OC system, including OC inputs, decay rates, and transit times. These additional estimates provided insights into the dynamics of OC within the soils of our sites, including (i) the importance of in situ soil stabilization mechanisms for maintaining the slow rates of soil OC decay needed to retain OC for long timespans, and (ii) how the interplay of CARs and OC transit times can result in differences in soil OC stocks. Our estimates of centuries long OC transit times further substantiate the role of seagrass ecosystems in climate change mitigation and emphasize the importance of conserving blue carbon ecosystems to preserve the efficacy of their soil OC protection mechanisms. Finally, although we show here that robust modeling approaches can help to diminish estimation uncertainty, there are other important contributors to the uncertainty in global CAR estimates. Further efforts to standardize sampling methods, increase our understanding of variability within sites and among geomorphological settings, and improve estimates of habitat extent would substantially improve our ability to quantify the rate OC accumulates globally within soils of blue carbon ecosystems. References Anderson, D. R. K., Burnham, K. P., & White, C. C. (1998). Comparison of Akaike information criterion and consistent Akaike information criterion for model selection and statistical inference from capture-recapture studies. Journal of Applied Statistics,25, 263–282. Appleby, P. G., & Oldfield, F. (1978). The calculation of lead-210 dates assuming a constant rate of supply of unsupported 210Pb to the sediment. Catena,5(1), 1–8. Appleby, P. G., & Oldfield, F. (2004). The assessment of 210Pb data from sites with varying sediment accumulation rates. Hydrobiologia, 103, 29–35. Austin, P. C., & Hux, J. E. (2002). A brief note on overlapping confidence intervals. Journal of Vascular Surgery,36(1), 194–195. Blaauw, M. (2010). Methods and code for ‘classical’ age-modelling of radiocarbon sequences. Quaternary Geochronology,5(5), 512–518. Blaauw, M., & Christen, J. A. (2011). Flexible paleoclimate age-depth models using an autoregressive gamma process. Bayesian Analysis, 6(3), 457–474. Blaauw, M., & Christen, J. A. (2018). rbacon: Age-depth modelling using Bayesian statistics. R package version 2.3.4 [Computer software manual]. Breithaupt, J. L., Smoak, J. M., Smith III, T. J., Sanders, C. J., & Hoare, A. (2012). Organic carbon burial rates in mangrove sediments: Strengthening the global budget. Global Biogeochemical Cycles,26, GB3011. https://doi.org/10.1029/2012GB004375 Acknowledgments For this work E. Fay Belshe was supported by the German Academic Exchange Service (DAAD) with funds from the German Federal Ministry of Education and Research (BMBF) and the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme (FP7/2007-2013) under REA Grant Agreement 605728 (P.R.I.M.E. - Postdoctoral Researchers International Mobility Experience). O. Serrano was supported by an ARC DECRA (DE170101524). The work leading to this publication was supported by the projects SUMILEN (CTM2013-47728-R) and PALEOPARK (1104/2014), funded by the Spanish National Parks and State Research Schemes. We also would like to acknowledge the Coastal Carbon Research Coordination Network CCRCN (NSFDEB 1655622) for helping us submit our data to the CCRCN database at the Smithsonian Institute and providing insights into coastal soil carbon dynamics during the December 2018 Soil Carbon Working Group workshop. All data used in this publication can be found on the Smithsonian Institution's Figshare data repository (https://doi. org/10.25573/data.9856769.v1) and can also be accessed through the Coastal Carbon Atlas (https://ccrcn. shinyapps.io/CoastalCarbonAtlas/). In addition, in the supporting information we have uploaded a folder containing the data sets, R and Matlab scripts, and user information to allow for the reproduction of our results. Please acknowledge these data products, using the following citation: Belshe, E. F., Sanjuan, J., Leiva-Dueñas, C., Piñeiro-Juncal, N., Serrano, O., Lavery, P. S., Mateo, M. A. (2019) Dataset: Modeling organic carbon accumulation rates and residence times in coastal vegetated ecosystems. figshare. Dataset. (https:// doi.org/10.25573/data.9856769). BELSHE ET AL. 3668 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 Brindgham, S. D., Megonigal, J. P., Keller, J. K., Bliss, N. B., & Trettom, C. (2006). The carbon balance of North American wetlands. Wetlands, 26(4), 889–916. Brodersen, K. E., Koren, K., Moßhammer, M., Ralph, P. J., Kühl, M., & Santner, J. (2017). Seagrass-mediated phosphorus and iron solubilization in tropical sediments. Environmental Science & Technology,51(24), 14,155–14,163. Callaway, J. C., Borgnis, E. L., Turner, R. E., & Milan, C. S. (2012). Carbon sequestration and sediment accretion in San Francisco Baytidal wetlands. Estuaries and Coasts,35(5), 1163–1181. Caruso, T., de Vries, F. T., Bardgett, R. D., & Lehmann, J. (2018). Soil organic carbon dynamics matching ecological equilibrium theory. Ecology and Evolution,2, 1–10. Chapman, S. K., Hayes, M. A., Kelly, B., & Langley, JA. (2019). Exploring the oxygen sensitivity of wetland soil carbon mineralization. Biology Letters,15, 20180407. Chen, R., & Twilley, R. R. (1999). A simulation model of organic matter and nutrient accumulation in mangrove wetland soils. Biogeochemistry,44(1), 93–118. Chiu, S.-H., Huang, Y.-H., & Lin, H.-J. (2013). Carbon budget of leaves of the tropical intertidal seagrass Thalassia hemprichii.Estuarine Coastal and Shelf Science,125, 27–35. Chmura, G. L., Anisfeld, S. C., Cahoon, D. R., & Lynch, J. C. (2003). Global carbon sequestration in tidal, saline wetland soils. Global Biogeochemical Cycles,17(4), 1111. https://doi.org/10.1029/2002GB001917 Christen, J. A., & Fox, C. (2010). A general purpose sampling algorithm for continuous distributions (the t-walk). Bayesian Analysis,5(2), 263–281. Christian, R. R. (1984). A life-table approach to decomposition studies. Ecology,65(5), 1693–1697. Clymo, R. S. (1984). The limits to Peat Bog growth. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 303(1117), 605–654. Donato, D. C., Kauffman, J. B., Murdiyarso, D., Kurnianto, S., Stidham, M., & Kanninen, M. (2011). Mangroves among the most carbon-rich forests in the tropics. Nature Geoscience,4, 293–297. Donnelly, A. P., & Herbert, R. A. (1999). Bacterial interactions in the rhizosphere of seagrass communities in shallow coastal lagoons. Journal of Applied Microbiology Symposium Supplement,85, 151S–160S. Duarte, C. M., Fourqurean, J. W., Krause-Jensen, D., & Olesen, B. (2006). Dynamics of seagrass stability and change, Seagrasses: Biology, ecology and conservation (pp. 1–24). The Netherlands: Springer Dordrecht Duarte, C. M., Holmer, M., & Marba, N. (2005). Plant-microbe interactions in seagrass meadows, Interactions between macro- and microorganisms in marine sediments (pp. 31–60). Washington, D. C.: American Geophysical Union. Duarte, C. M., Losada, I. J., Hendriks, I. E., Mazarrasa, I., & Marba, N. (2013). The role of coastal plant communities for climate change mitigation and adaptation. Nature Publishing Group,3(11), 961–968. Ekschmitt, K., Kandeler, E., Poll, C., Brune, A., Buscot, F., Friedrich, M., et al. (2008). Soil-carbon preservation through habitat constraints and biological limitations on decomposer activity. Journal of Plant Nutrition and Soil Science,171, 27–35. Enriquez, S., Duarte, C. M., & Sand-Jensen, K. (1993). Patterns in decomposition rates among photosynthetic organisms: The importance of detritus C:N:P content. Oecologia,94, 457–471. Eriksson, E. (1971). Compartment models and reservoir theory. Annual Review of Ecology and Systematics,2, 67–84. Evert, S., & Baroni, M. (2007). zipfR: Word frequency distributions in R, Proceedings of the 45th Annual Meeting of the Association for Computational Linguistics, Posters and Demonstrations Sessions (pp. 29–32). Prague, Czech Republic. Fourqurean, J. W., Duarte, C. M., Kennedy, H., Marba, N., Holmer, M., Mateo, M. A., et al. (2012). Seagrass ecosystems as a globally significant carbon stock. Nature Geoscience,5(7), 505–509. Fourqurean, J. W., & Schrlau, J. E. (2003). Changes in nutrient content and stable isotope ratios of C and N during decomposition of seagrasses and mangrove leaves along a nutrient availability gradient in Florida Bay, USA. Chemistry and Ecology,19(5), 373–390. Gälman, V., Rydberg, J., de Luna, S. S., Bindler, R., & Renberg, I. (2008). Carbon and nitrogen loss rates during aging of lake sediment: Changes over 27 years studied in varved lake sediment. Limnology and Oceanography,53(3), 1076–1082. Gentile, R., Vanlauwe, B., & Six, J. (2011). Litter quality impacts shortbut not long-term soil carbon dynamics in soil aggregate fractions. Ecological Applications,21(3), 695–703. Harrison, P. G. (1989). Detrital processing in seagrass systems: A review of factors affecting decay rates, remineralization and detritivory. Aquatic Botany,23, 263–288. Haslett, J., & Parnell, A. (2008). A simple monotone process with application to radiocarbon-dated depth chronologies. Journal of the Royal Statistical Society: Series C (Applied Statistics),57(4), 399–418. Hector, A. (2015). The new statistics with R, An introduction for biologists. Oxford: Oxford University Press. Hemminga, M. A., & Buth, G. J. C. (1991). Decomposition in salt marsh ecosystems of the S.W. Netherlands: The effects of biotic and abiotic factors. Vegetatio,92, 73–81. Jankowska, E., Michel, L. N., Zaborska, A., & Włodarska-Kowalczuk, M. (2016). Sediment carbon sink in low-density temperate eelgrass meadows (Baltic Sea). Journal Of Geophysical Research: Biogeosciences,121, 1–17. https://doi.org/10.1002/2016jg003424 Johannessen, S. C., & Macdonald, R. W. (2016). Geoengineering with seagrasses: Is credit due where credit is given? Environmental Research Letters,11, 113001. Jordi, A., Basterretxea, G., Orfila, A., & Tintore, J. (2006). Analysis of the circulation and shelf-slope exchanges in the continental margin of the northwestern Mediterranean. Ocean Science,2, 173–181. Jordi, A., Basterretxea, G., & Wang, D.-P. (2009). Evidence of sediment resuspension by island trapped waves. Geophysical Research Letters, 36, L18610. https://doi.org/10.1029/2009GL040055 Kaldy, J. E., Eldridge, P. M., Cifuentes, L. A., & Jones, W. B. (2006). Utilization of DOC from seagrass rhizomes by sediment bacteria: 13C-tracer experiments and modeling. Marine Ecology Progress Series,317, 41–55. Karlen, I., Olsson, I. U., Killburg, P., & Kilici, S. (1964). Absolute determination of the activity of two 14C dating standards. Arkiv Geofysik, 4, 465–471. Keil, R. G., & Mayer, L. M. (2014). Mineral matrices and organic matter. In H. D. Holland, & K. K. Turekian (Eds.), Treatise on geochemistry (second edition) (pp. 337–359). Oxford: Elsevier Ltd. Kenworthy, J. W., & Thayer, G. W. (1984). Production and decomposition of the roots and rhizomes of segrasses, Zostera marina and Thalassia testudinum in temperate and subtropical marine ecosystems. Bulletin of Marine Science,35(3), 364–379. Kirwan, M. L., & Blum, L. K. (2011). Enhanced decomposition offsets enhanced productivity and soil carbon accumulation in coastal wetlands responding to climate change. Biogeosciences,8, 987–993. Kirwan, M. L., Guntenspergen, G. R., & Langley, J. A. (2014). Temperature sensitivity of organic-matter decay in tidal marshes. Biogeosciences,11, 4801–4808. BELSHE ET AL. 3669 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 Kuzyakov, Y. (2010). Priming effects: Interactions between living and dead organic matter. Soil Biology and Biochemistry,42, 1363–1371. López, N. I., Duarte, C. M., Vallespinos, F., Romero, J., & Alcoverro, T. (1995). Bacterial activity in NW Mediterranean seagrass (Posidonia oceanica) sediments. Journal of Experimental Marine Biology and Ecology,187, 39–49. López-Merino, L., Serrano, O., Adame, M. F., Mateo, M. A., & Cortizas, A. M. (2015). Glomalin accumulated in seagrass sediments reveals past alterations in soil quality due to land-use change. Global and Planetary Change,133, 87–95. Lavery, P. S., Mateo, M. A., Serrano, O., & Rozaimi, M. (2013). Variability in the carbon storage of seagrass habitats and its implications for global estimates of blue carbon ecosystem service. PLOS ONE,8(9), e73748. Levasseur, A., Brandão, M., Lesage, P., Margni, M., Pennington, D., Clift, R., & Samson, R. (2012). Valuing temporary carbon storage. Nature Climate Change,2,6–8. Lovelock, C. E., Fourqurean, J. W., & Morris, J. T. (2017). Modeled CO2emissions from coastal wetland transitions to otherland uses: Tidal marshes, mangrove forests, and seagrass beds. Frontiers in Marine Science,4, 143. MATLAB (2015). version 8.5.0.197613 (R2015a). Natick, Massachusetts: The MathWorks Inc. Mabit, L., Benmansour, M., Abril, J. M., Walling, D. E., Meusburger, K., Iurian, A. R., et al. (2014). Fallout 210-Pb as a soil and sediment tracer in catchment sediment budget investigations: A review. Earth Science Reviews,138, 335–351. Macreadie, P. I., Baird, M. E., Trevathan-Tackett, S. M., Larkum, AWD, & Ralph, P. J. (2014). Quantifying and modelling the carbon sequestration capacity of seagrass meadows—A critical assessment. Marine Pollution Bulletin,83, 430–439. Macreadie, P. I., Ewers-Lewis, C. J., Whitt, A. A., Ollivier, Q., Trevathan-Tackett, S. M., Carnell, P., & Serrano, O. (2018). Comment on ‘Geoengineering with seagrasses: Is credit due where credit is given?’ Environmental Research Letters,13, 028002. Manzoni, S., Katul, G. G., & Porporato, A. (2009). Analysis of soil carbon transit times and age distributions using network theories. Journal of Geophysical Research,114, G04025. https://doi.org/10.1029/2009JG001070 Marba, N., Arias Ortiz, A., Masqué, P., Kendrick, G. A., Mazarrasa, I., Bastyan, G. R., et al. (2015). Impact of seagrass loss and subsequent revegetation on carbon sequestration and stocks. JournalofEcology,103, 296–302. Marland, G., Fruit, K., & Sedjo, R. (2001). Accounting for sequestered carbon: The question of permanence. Environmental Science & Policy, 4, 259–268. Mateo, M. A., & Romero, J. (1996). Evaluating seagrass leaf litter decomposition an experimental comparison between litter-bag and oxygen-uptake methods. Journal of Experimental Marine Biology and Ecology,202, 97–106. Mateo, M. A., Romero, J., Pereza, M., Littler, M. M., & Littler, D. S. (1997). Dynamics of millenary organic deposits resulting from the growth of the Mediterranean seagrass Posidonia oceanica.Estuarine Coastal and Shelf Science,44, 103–110. Mazarrasa, I., Marba, N., Garcia Orellana, J., Masqué, P., Arias Ortiz, A., & Duarte, C. M. (2017). Effect of environmental factors (wave exposure and depth) and anthropogenic pressure in the C sink capacity of Posidonia oceanica meadows. Limnology and Oceanography, 62(4), 1436–1450. McGuire, K. J., & McDonnell, J. J. (2006). A review and evaluation of catchment transit time modeling. Journal of Hydrology,330, 543–563. Middleton, B. A., & McKee, K. L. (2001). Degradation of mangrove tissues and implications for peat formation in Belizean Island forests. Journal of Ecology,89, 818–828. Miyajima, T., Hori, M., Hamaguchi, M., Shimabukuro, H., Adachi, H., Yamano, H., & Nakaoka, M. (2015). Geographic variability in organic carbon stock and accumulation rate in sediments of East and Southeast Asian seagrass meadows. Global Biogeochemical Cycles,29, 1–19. https://doi.org/10.1002/2014GB004979 Morris, J. T., Barber, D. C., Callaway, J. C., Chambers, R., Hagen, S. C., Hopkinson, C. S., et al. (2016). Contributions of organic and inorganic matter to sediment volume and accretion in tidal wetlands at steady state. Earth's Future,4, 110–121. https://doi.org/10.1002/ 2015EF000334 Morris, J. T., & Bowden, W. B. (1986). A mechanistic, numerical model of sedimentation, mineralization, and decomposition for marsh sediments. Soil Science Society of America Journal,50, 96–105. Morrisey, D., Swales, A., Dittmann, S., Morrison, M., Lovelock, C., & Beard, C. (2010). The ecology and management of temperate mangroves. Oceanography and Marine Biology,48, 43–160. Morton, R. A., & White, W. (1997). Characteristics of and corrections for core shortening in unconsolidated sediments. Journal of Coastal Research,13(3), 761–769. National Academies of Sciences, Engineering, and Medicine (2019). Negative emissions technologies and reliable sequestration: A research agenda. Washington, DC: The National Academies Press. Nicastro, A., Onoda, Y., & Bishop, M. J. (2012). Direct and indirect effects of tidal elevation on eelgrass decomposition. Marine Ecology Progress Series,456, 53–62. Ouyang, X., & Lee, S. Y. (2014). Updated estimates of carbon accumulation rates in coastal marsh sediments. Biogeosciences,11, 5057–5071. Piñeiro-Juncal, N., Mateo, MA, Holmer, M., & Martínez-Cortizas, A. (2018). Potential microbial functional activity along a Posidonia oceanica soil profile. Aquatic Microbial Ecology,81, 189–200. Poret, N., Twilley, R. R., Rivera-Monroy, V. H., & Coronado-Molina, C. (2007). Be lowground decomposition of mangrove roots in Florida coastal everglades. Estuaries and Coasts,30(3), 491–496. R Core Team (2018). R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria. Ramsey, C. B. (2008). Deposition models for chronological records. Quarternary Science Reviews,27(1-2), 42–60. Reimer, P. J., Bard, E., Bayliss, A., Beck, J. W., Blackwell, P. G., Ramsey, C. B., et al. (2013). IntCal13 and Marine13 radiocarbon age calibration curves 0–50,000 years cal BP. Radiocarbon,55(4), 1869–1887. Reimer, P. J., & McCormac, F. G. (2016). Marine radiocarbon reservoir corrections for the Mediterranean and Aegean Seas. Radiocarbon, 44(1), 159–166. Robinson, A. R., Leslie, W. G., Theocharis, A., & Lascaratos, A. (2001). Mediterranean Sea circulation. In S. A. Thorpe (Ed.), Ocean currents (pp. 283–298). San Diego. Romero, J., Pergent, G., Pergent-Martini,C., Mateo, M. A., & Regnier, C. (1992). The detritic compartment in a Posidonia oceanica meadow: Litter features, decomposition rates, and mineral stocks. Marine Ecology,13(1), 69–83. Rovai, A. S., Twilley, R. R., Castañeda-Moya, E., Riul, P., Cifuentes-Jara, M., Manrow-Villalobos, M., et al. (2018). Global controls on carbon storage in mangrove soils. Nature Climate Change,8, 534–538. Rozaimi, M., Fairoz, M., Hakimi, T. M., Hamdan, N. H., Omar, R., Ali, M. M., & Tahirin, S. A. (2017). Carbon stores from a tropical seagrass meadow in the midst of anthropogenic disturbance. Marine Pollution Bulletin,119, 253–260. Rozaimi, M., Lavery, P. S., Serrano, O., & Kyrwood, D. (2016). Long-term carbon storage and its recent loss in an estuarine Posidonia australis meadow (Albany, Western Australia). Estuarine, Coastal and Shelf Science,171, 58–65. Rybczyk, J. M., Callaway, J. C., & Day, J. W. (1998). A relative elevation model for a subsiding coastal forested wetland receiving wastewater effluent. Ecological Modelling,112, 23–44. BELSHE ET AL. 3670 Journal of Geophysical Research: Biogeosciences 10.1029/2019JG005233 Rybczyk, JM, Garson, G., & Day, JW (1996). Nutrient enrichment and decomposition in wetland ecosystems: Models, analyses and effects. Current Topics in Wetland Biogoechemistry,2, 52–72. Samper-Villarreal, J., Mumby, P. J., Saunders, M. I., Barry, L. A., Zawadzki, A., Heijnis, H., et al. (2018). Vertical accretion and carbon burial rates in subtropical seagrass meadows increased following anthropogenic pressure from European colonisation. Estuarine Coastal and Shelf Science,202, 40–53. Serrano, O., Lavery, P. S., López-Merino, L., Ballesteros, E., & Mateo, M. Á. (2016). Location and associated carbon storage of erosional escarpments of seagrass Posidonia Mats. Frontiers in Marine Science,3(42), 1–19. Serrano, O., Lavery, P. S., Rozaimi, M., & Mateo, M. A. (2014). Influence of water depth on the carbon sequestration capacity of seagrasses. Global Biogeochemical Cycles,28, 950–961. Serrano, O., Mateo, M. A., Dueñas-Bohórquez, A., Renom, P., López-Sáez, J. A., & Cortizas, A. M. (2011). The Posidonia oceanica marine sedimentary record: A Holocene archive of heavy metal pollution. Science of The Total Environment,409(22), 1–10. Serrano, O., Mateo, MA, Renom, P., & Julià, R. (2012). Characterization of soils beneath a Posidonia oceanica meadow. Geoderma,185-186, 26–36. Serrano, O., Ricart, A. M., Lavery, P. S., Mateo, M. A., Arias-Ortiz, A., Masqué, P., et al. (2016). Key biogeochemical factors affecting soil carbon storage in Posidonia meadows. Biogeosciences,13, 4581–4594. Siani, G., Paterne, M., Arnold, M., Bard, E., Métivier, B., Tisnerat, N., & Bassinot, F. (2016). Radiocarbon reservoir ages in the Mediterranean Sea and Black Sea. Radiocarbon,42(2), 271–280. Sierra, C. A., Müller, M., Metzler, H., Manzoni, S., & Trumbore, S. E. (2017). The muddle of ages, turnover, transit, and residence times in the carbon cycle. Global Change Biology,23, 1763–1773. Sobek, S., Durisch-Kaiser, E., Zurbrügg, R., Wongfun, N., Wessels, M., Pasche, N., & Wehrli, B. (2009). Organic carbon burial efficiency in lake sediments controlled by oxygen exposure time and sediment source. Limnology and Oceanography,54(6), 2243–2254. Spivak, A. C., Sanderman, J., Bowen, J. L., Canuel, E. A., & Hopkinson, C. S. (2019). Global-change controls on soil-carbon accumulation and loss in coastal vegetated ecosystems. Nature Geoscience,12, 685–962. Stuiver, M., & Reimer, P. J. (1993). Extended 14C data base and revised CALIB 3.0 14C age calibration program. Radiocarbon,35(1),215–230. Stuiver, M., Reimer, P. J., Bard, E., Beck, J. W., Burr, G. S., Hughen, K. A., et al. (1998). INTCAL98 radiocarbon age calibration, 24,000-0 cal BP. Radiocarbon,40, 1041–1083. Trachsel, M., & Telford, R. J. (2016). All age–depth models are wrong, but are getting better. The Holocene,27(6), 860–869. Trevathan-Tackett, S. M., Thomson, AlexandraC. G., Ralph, P. J., & Macreadie, P. I. (2018). Fresh carbon inputs to seagrass sediments induce variable microbial priming responses. Science of the Total Environment,621, 663–669. Trumbore, S. E. (2000). Age of soil organic matter and soil respiration: Radiocarbon constraints on belowground C dynamics. Ecological Applications,10(2), 399–411. Trumbore, S. E. (2009). Radiocarbon and soil carbon dynamics. Annual Review of Earth and Planetary Sciences,37, 47–66. Trumbore, S. E., & Harden, J. W. (1997). Accumulation and turnover of carbon in organic and mineral soils of the BOREAS northern study area. Journal of Geophysical Research,102, 28,817–28,830. Trumbore, S. E., Sierra, C. A., & Hicks Pries, C. E. (2016). Radiocarbon nomenclature, theory, models, and interpretation: Measuring age, determining cycling rates, and tracing source pools. In E. A. G. Schuur, E. R. M. Druffel, & S. E. Trumbore (Eds.), Radiocarbon and climate change (pp. 45–82). Switzerland. von Lützow, M., Kögel-Knabner, I., Ludwig, B., Matzner, E., Flessa, H., Ekschmitt, K., et al. (2008). Stabilization mechanisms of organic matter in four temperate soils: Development and application of a conceptual model. Journal of Plant Nutrition and Soil Science,171(1), 111–124. Wickham, H. (2016). ggplot2: Elegant graphics for data analysis. Verlag, New York: Springer. Wilkinson, G. M., Besterman, A., Buelo, C., Gephart, J., & Pace, M. L. (2018). A synthesis of modern organic carbon accumulation rates in coastal and aquatic inland ecosystems. Scientific Reports,8, 15736. Wright, A. J., Edwards, R. J., van de Plassche, O., Blaauw, M., Parnell, A. C., van der Borg, K., et al. (2017). Reconstructing the accumulation history of a saltmarsh sediment core: Which age-depth model is best? Quaternary Geochronology,39, 35–67. York, P. H., Macreadie, P. I., & Rasheed, M. A. (2018). Blue carbon stocks of Great Barrier Reef deep-water seagrass. Biology Letters,14, 20180529. Zeileis, A., & Grothendieck, G. (2005). zoo: S3 infrastructure for regular and irregular time series. Journal of Statistical Software,14(6), 1–27. Zhang, D., Hui, D., Luo, Y., & Zhou, G. (2008). Rates of litter decomposition in terrestrial ecosystems: Global patterns and controlling factors. Journal of Plant Ecology,1(2), 85–93. BELSHE ET AL. 3671