Coastal marshes and their carbon pools adapt to sea level rise largely through ecogeomorphic feedbacks in which increased flooding stimulates plant growth, mineral sediment deposition, and vertical soil development (D'Alpaos et al., 2007;Kirwan and Megonigal, 2013;Marani et al., 2007;Morris et al., 2002). Because plant productivity and organic matter accumulation are inherently linked in anaerobic soils, these ecogeomorphic feedbacks also determine the amount of carbon that marsh soils accumulate through time (Gonneea et al., 2019;Kirwan & Mudd, 2012). Recent work indicates that coastal wetlands are an important global carbon sink (Chmura, 2013;Hopkinson et al., 2012;Mcleod et al., 2011), and that carbon accumulation rates increase with accelerated sea level rise (Rogers et al., 2019;Wang et al., 2019). Coastal wetlands therefore potentially represent a unique negative carbonclimate feedback where carbon emissions lead to faster rates of sea level rise and enhanced carbon sequestration, making them an important tool towards mitigating changes in the Earth's climate (Crooks et al., 2011;Holmquist et al., 2018).
Nevertheless, the areal extent of coastal marshes has declined worldwide (Duarte, 2008), and there are concerns over the stability of coastal carbon pools in the face of interacting components of global change. Previous numerical modeling and stratigraphic observations suggest that the elevation of marshes and the size of their carbon pools increase with sea level rise until some threshold rate, beyond which they drown (Kirwan et al., 2010;Morris et al., 2002). Elevated CO2 (eCO2) increases marsh elevation gain in both short-term field experiments (Langley et al., 2009;Reef et al., 2017) and longterm modeling efforts (Ratliff et al., 2015), but how the interacting effects of SLR and eCO2 influence marsh resilience and carbon accumulation over decades to centuries is poorly understood. Empirical studies suggest that eCO2 increases soil carbon in terrestrial systems up to a saturation point (Heimann & Reichstein, 2008;van Groenigen et al., 2014) and may make coastal marshes more resilient to SLR by increasing soil elevation via C3 plant productivity (Langley et al., 2009;Reef et al., 2017). However, increases in plant productivity may in turn lead to increased decay of older carbon via root-derived inputs of organic carbon and delivery of oxygen into an otherwise anaerobic soil, thereby decreasing the belowground carbon pool (Bernal et al., 2017;Jones et al., 2018;Mueller et al., 2015;Wolf et al., 2007).
Vegetation growth leads to a persistent oxygenated zone in wetland sediments (Marani et al., 2006;Boaga et al., 2014), and root induced priming has been shown to be a key factor in regulating the direction of change in terrestrial carbon stocks (Groenigen et al., 2014). However, the importance of priming relative to other drivers remains unexplored in coastal carbon pools.
Vegetation type may also be a strong driver of marsh carbon accumulation under SLR and eCO2.
Due to the relatively stressful conditions in coastal marshes, many plant species have evolved to use the C4 photosynthetic pathway, which unlike C3 plants, utilizes a CO2 concentrating mechanism that negates CO2 limitation. Thus, as opposed to studies of the effects of eCO2 on elevation in C3 marshes, C4 vegetation has been shown to exhibit little to no response to eCO2 (Bernal et al., 2017;Morris & Bowden, 1986;Mueller et al., 2015). Previous modeling neglects these differences in vegetation community response for simplicity and because of inherent difficulties separating the effects of eCO2 in mixed C3/C4 communities (Ratliff et al., 2015). Thus, a gap remains between short-term field experiments that show the importance of vegetation type, and long-term model experiments that suggest the importance of elevation dependent feedbacks. Here, we demonstrate with a novel soil-cohort model that eCO2 and SLR interact synergistically to increase soil carbon burial, driven by shifts in plant community composition, that facilitate an ever-expanding soil volume.
Our model is designed to simulate changes in the elevation and carbon content of a soil column at a single point on a marsh surface through time, and in response to environmental drivers such as SLR and eCO2. Following previous soil-cohort approaches (Kirwan & Mudd, 2012;Morris & Bowden, 1986), soil cohorts are built annually through the deposition of mineral and organic sediment on the marsh surface and each cohort expands and contracts through time according to organic matter production and decomposition within the soil column. As in previous approaches, mineral deposition and organic matter production vary with the depth and duration of tidal inundation of the marsh surface, and organic matter production and decomposition decrease exponentially with depth below the soil surface. These models demonstrate that marsh elevations and carbon stocks equilibrate to moderate rates of SLR, whereby rates of soil formation equal rates of SLR, but drown under higher rates of SLR (Kirwan & Mudd, 2012;Mudd et al., 2009).
Previous modeling efforts have focused on capturing only the most essential ecomorphodynamic interactions, neglecting many plant and microbial feedbacks that determine carbon preservation (Spivak et al., 2019). Here, we extend their utility by expanding the treatment of vegetation growth, belowground production, and decomposition to include nuanced feedbacks between vegetation type, eCO2, and organic matter priming. We consider two vegetation communities, a C4 marsh parameterized for Spartina patens and a C3 marsh parameterized for Schoenopletus americanus, both common tidal marsh species across North America. At elevations where species overlap, our model creates mixed communities with productivity and decomposition parameterizations weighted according to the relative species distribution. Organic matter decomposition rates increase as aboveground biomass increases (Jones et al., 2018;Mueller et al., 2015), reflecting priming of soil organic matter decomposition associated with root exudation and turnover. The model is used to explore the response of marsh soil carbon to interactions between SLR and eCO2 using observations from the Smithsonian Global Change Research Wetland (GCReW), an organic rich microtidal marsh on a tributary of the Chesapeake Bay (USA) that includes the longest running eCO2 experiment in the world (Drake, 2014).
We parameterized our model using empirical data and a model hindcast to represent conditions similar to a high marsh at GCReW. The entire GCReW site receives negligible mineral sediment, allowing us to simplify lateral gradients in sediment supply, and isolate the effects of dynamic organic matter cycling as drivers of marsh accretion and carbon accumulation. While other modeling experiments on the effects of eCO2 treated C3 and C4 species identically (Ratliff et al., 2015), long-term data at GCReW conclusively shows that eCO2 increases C3 biomass production with little effect on C4 biomass (Langley et al., 2009;Wolf et al., 2007). Our model therefore separates vegetation parameterizations for these fundamental plant functional types and additionally accounts for GCReW data that shows the CO2 fertilization effect on plant biomass is maximized at an intermediate inundation depth (Langley et al., 2013). Other empirical data used from the site includes tidal range, rooting depth profiles (Megonigal, 1999), species-specific relationships between aboveground biomass and elevation (Byrd et al., 2017), their responses to eCO2 (Drake, 2014;Groenigen et al., 2014;Langley et al., 2009), and the relationship between decomposition rate and aboveground biomass (Jones et al., 2018;Mueller et al., 2015;Supplementary Table 1).
In order to estimate the parameters for which we did not have data, chiefly the decomposition and turnover rates of belowground biomass, we performed a model hindcast and adjusted unknown parameters until the model produced stratigraphic characteristics consistent with GCReW. This consisted of a spinup period that created an organic rich soil profile in equilibrium with the local late-Holocene rate of relative SLR (~1 mm yr -1 ), and then a 150-year model run under the historic local rate of SLR (3.6 mm yr -1 ) observed in Annapolis, MD (NOAA, 2019). The model produced a final marsh elevation (0.34 m NAVD), accretion rate (~3.4 mm yr -1 ), and soil organic matter profile that are within the range of high marsh characteristics observed at GCReW today (Messerschmidt & Kirwan, 2020; Fig. 1). Also consistent with field observations at GCReW is the modeled loss of elevation relative to sea level (i.e. disequilibrium) (Fig. 1), as evidenced by the increase in flood tolerant C3 low marsh species encroaching into C4 high marsh habitat over the past three decades (Lu, in press). Taken together, these observations suggest that the model is generally capable of simulating soil-building processes at GCReW and other organic rich marshes.
To understand basic model behavior in a submerging marsh, we began our model experiments by subjecting a high elevation C4 marsh to a rate of SLR high enough to induce rapid drowning (25 mm yr - 1 ; Fig. 2). In these model runs, simulations begin with a marsh elevation (0.43 m NAVD) and organic rich soil profile created during the model spin-up period under a 1 mm yr -1 rate of SLR. Under ambient CO2 (aCO2), progressive inundation drives the conversion of a C4 marsh to a mixed community dominated by C4 vegetation (> 50%), then to a mixed community dominated by C3 vegetation, and finally to a C3 marsh that submerges (Fig. 2a). Accretion for the first 10 years is driven purely by organic inputs into a saturated high marsh that receives no tidal inundation. Although mineral accretion rates then increase through time, driven by longer and deeper flooding of the marsh surface, our parameterization (suspended sediment concentration= 5 mg L -1 ; tidal amplitude = 0.22 m) limits mineral sediment deposition and ensures that organic matter accumulation dominates marsh elevation change.
Under the modeled instantaneous increase in sea level rise rate (1 to 25 mm yr -1 ), organic accretion rates initially increase from ~ 1 mm yr -1 to 4 mm yr -1 (Fig. 2a) as production of C4 vegetation exceeds decomposition (Fig. 2b). As vegetation shifts to more flood tolerant C3 vegetation with slower parameterized root turnover, organic accretion rates decline, and the marsh eventually drowns (Fig. 2a).
Total carbon summed over the entire profile increases throughout the experiment, even as instantaneous carbon accumulation rates fluctuate during conversion to C3 vegetation (Fig. 2b).
Under eCO2 conditions, the model predicts qualitatively similar results (i.e. identical accretion rates for the C4 vegetation community, fluctuation in organic matter accumulation driven by vegetation type, and eventual submergence of the marsh platform). However, the positive effect of eCO2 on C3 vegetation growth allows the marsh to survive longer than under aCO2 (Fig. 2c). Elevated CO2 prolonged a state change from tidal marsh to open water by over a decade under the accelerated rate of SLR applied in this modeling exercise (25 mm yr -1 ), a response that would likely translate into several decades under most contemporary SLR scenarios. This behavior is driven by a bigger difference between enhanced production and decomposition than under aCO2 (Fig. 2d), and a more persistent mixed community where rapid C4 turnover accompanies C3 growth enhanced by eCO2, resulting in a synergistic enhancement of organic matter accumulation (Fig. 2c). Total organic matter and carbon summed across the soil profile are higher under eCO2 (Fig. 2d) than aCO2 (Fig. 2b), indicating that the priming effect of plant biomass under eCO2 does not completely offset the increase in organic production.
Next, we conducted three separate model runs at 25 mm yr -1 to explore how our novel parameterizations influence model behavior relative to approaches used in previous models, and in particular, the approaches used in the only other eCO2-informed tidal marsh model (Ratliff et al., 2015).
Our model differs from previous models in two key ways: 1) we parameterize C3 and C4 species separately to include a CO2 fertilization effect on C3 species only, as opposed to the assumption that eCO2 affects all vegetation (Ratliff et al., 2015), and 2) our model incorporates a previously unexplored priming effect where decomposition fluctuates with plant productivity due to the introduction of fresh carbon, and/or radial oxygen loss from roots and rhizomes (Jones et al., 2018;Mueller et al., 2015). This means our decay rate is dynamic, changing throughout a model simulation as opposed to a static value.
In these experiments, both species-specific eCO2 responses and organic matter priming lead to faster marsh drowning relative to the constraints of previous model (Fig. 3). However, we also find that priming has a much stronger effect than the species-specific eCO2 parameterization (Fig. 3). This surprising behavior is illustrated by the substantial difference in biomass produced between runs with and without species-specific effects, and yet little difference in the time until marsh submergence.
These strikingly similar results occur because the increase in biomass of the C4 community under eCO2 is counterbalanced by an increase in decomposition in our model (Fig. 3). This highlights an important aspect of our model, that production and decomposition are tightly coupled through time due to the plant-mediated priming effect and suggests that previous model results may overestimate the positive effects of eCO2 on marsh resilience (Ratliff et al., 2015).
To understand the interactive effect of eCO2 and SLR on marsh resilience and carbon accumulation rate, we subjected a high elevation C4 marsh to progressively faster rates of SLR and ran the model until it either equilibrated to the new rate of SLR or drowned. These simple model experiments illustrate that carbon accumulation rates (Fig. 4) and total carbon in the soil profile (Fig. 5) increase with increasing rates of SLR until the point of marsh drowning. Interestingly, increases in soil carbon occur even as soils became more mineral rich, driven by surficial sediment deposition that increased with tidal inundation (Fig. 5). Although this model behavior is consistent with faster burial and more efficient carbon preservation, we suggest it is more likely due to mineral deposition rates that increase more quickly than organic matter production rates, resulting in a relative decrease in percent organic matter (i.e. a decrease in carbon concentration). In contrast to terrestrial systems in which carbon accumulation is driven by changes in carbon concentration (Lu et al., 2019;Stewart et al., 2007), our results uniquely illustrate that marsh carbon pools are driven by changes in soil volume.
Our finding that carbon accumulation rates increase with the rate of SLR is consistent with previous modeling and stratigraphic observations that attribute accelerating rates of carbon accumulation to the reduction in carbon-sequestration saturation effects associated with an ever-expanding soil volume (Rogers et al., 2019). However, the explicit modeling of the interaction between eCO2 and SLR leads to new insights into the mechanisms responsible for increasing carbon accumulation. For example, based on inundation alone, organic matter production decreases as marshes transition from pure C4 to mixed communities (Fig. 2a), but under eCO2 C3 productivity more than compensates for the decline in productivity, and carbon accumulation rates increase (Fig. 2d). This finding is qualitatively similar to previous modeling demonstrating that SLR and temperature warming lead to an increase in carbon accumulation rates (Kirwan & Mudd, 2012), but in those model experiments warming did not enhance marsh persistence. Here, our results show that eCO2 extends marsh persistence where SLR drives the conversion of C4 vegetation to C3 vegetation that accumulates organic matter faster under eCO2 (Fig. 4). In these experiments, marshes drown when SLR rates exceed 4 and 11 mm yr -1 under aCO2 and eCO2 respectively. Although specific threshold rates of SLR and the quantitative effect of eCO2 on marsh accretion depend on model parameterizations, increased marsh resilience is generally consistent with empirical field experiments (Langley et al., 2009;Reef et al., 2017). Thus, the primary influence of eCO2 is to allow the marsh to survive faster rates of SLR, which in turn facilitates soil volume expansion, and faster carbon accumulation (Fig. 4). Results from a model experiment in which an organic marsh equilibrated to a sea level rise rate of 1 mm yr -1 was subjected to a sea level rise rate of 25 mm yr -1 in order to induce submergence under ambient and elevated CO2 conditions. Background panel colors represent vegetation types; a C4 (green), C4 dominant (>50%) mixed (C4MD; blue), C3 dominant (>50%) mixed (C3MD; yellow), and C3 (purple) community. a) Accretion rate (organic, mineral, and total) and water depth above marsh surface at mean high tide for a submerging marsh under ambient CO2 conditions, b) Organic matter dynamics (production, decomposition, and net accumulation), and total carbon (g m -2 ) in the marsh soil profile for a submerging marsh under ambient CO2 conditions, c) Accretion rate and water depth above marsh surface for a submerging marsh under elevated CO2 conditions, d) Organic matter dynamics and total carbon in the soil profile for a submerging marsh under elevated CO2 conditions. Note the differring scales of the x-axes, and that negative water depths indicate a supratidal position of the marsh relative to mean high tide. CO2 at sea level rise rates between 1 and 12 mm yr -1 . Marsh drowns when sea level rise exceeds 4 mm yr -1 under aCO2 (black line) and 11 mm yr -1 under eCO2 (red line). There is no effect of eCO2 on CAR at low sea level rise rates because the marsh equilibrates to elevations too high for C3 vegetation.
Our model is designed to simulate changes in the elevation and carbon content of a soil column at a single point on a marsh surface through time in response to sea level rise and elevated CO2. Following the approach taken in other marsh soil cohort models (Morris & Bowden, 1986;Mudd et al., 2009), our model considers the evolution of a soil column discretized into cohorts (Q(t)) that represent layers of soil of a given age (t). New cohorts are added annually to the surface of the soil column through the deposition of sediment from semidiurnal tides and advected lower in the soil column as they age, such that the oldest cohort (Q(1), deposited at t=1) is at the bottom and the most recently deposited at the top (Q(t+1). The thickness of each cohort evolves through time as live roots and decaying organic matter are deposited within the soil column. The vertical expansion and contraction of soil cohorts directly translates into changes in marsh elevation. Thus, much like previous models of saltmarsh vertical accretion (D'Alapos et al., 2007;Marani et al., 2007;Morris et al., 2002;Ratliff et al., 2015), marsh elevation change in our model occurs through the cumulative changes in mineral and organic inputs. At each time step, the model determines sediment deposition on the marsh surface, vegetation type (C3, C4, or mixed), and organic matter production and decomposition that together describe changes in marsh elevation through time. Here, we describe each of these processes in detail. Dimensions are denoted in square brackets of [M] for mass, [L] for length, and [T] for time.
Sediment is deposited on the marsh surface as a function of the height and duration of tidal flooding, and the availability of suspended sediment in the water column. Following previous approaches (D 'Alapos et al., 2007), the annual mass of sediment deposited (𝑞𝑞 𝑠𝑠 ) is calculated as:
where qs [M L -2 T -1 ] is sediment settling as defined by the product of the settling velocity (ws) [L T -1 ], and the instantaneous suspended sediment concentration Ct [M L -3 ] integrated over the tidal cycle dt [T]. Following previous approaches (Marani et al., 2010), Ct is constant on the rising flood tide, but declines throughout the ebbing phase as the difference between the depth integrated mass of sediment coming onto the marsh from tides with a fixed concentration of sediment (Co) and the cumulative mass of sediment deposited on the marsh surface at time t. In these simulations, we chose Co = 5 mg/L to represent organic rich marshes far from tidal channels, like those at the Smithsonian Global Change Research Wetland (GCReW). The mass of sediment deposited on the marsh surface depends on the height and duration of tidal inundation, which is calculated over one average tidal cycle and extrapolated to an annual time step by multiplying by the number of tidal cycles in a year. Other numerical models consider the influence of vegetation on sediment deposition, and lateral gradients in sediment supply (D'Alapos et al., 2007;Marani et al., 2010;Ratliff et al., 2015). In contrast, we intentionally focus our modeling efforts on conditions (low tide range, low Co) that lead to negligible mineral sediment deposition so that we can isolate the effects of dynamic organic matter cycling on carbon accumulation and elevation change.
The growth and type of vegetation is also related to tidal inundation. Following previous approaches (Morris et al., 2002;Mudd et al., 2009), the aboveground biomass of vegetation is defined by a parabolic relationship with marsh surface elevation relative to sea level (Supp. Fig. 1), which is a proxy for flooding duration and associated environmental factors. In our model, we use biomass-elevation relationships for two common species at GCReW (Byrd et al., 2017). These species represent plants with different carbon fixation strategies -Spartina patens (C4) found at the higher elevations with little flooding, and Schenoplectus americanus (C3) found at lower elevations where flooding is common. Marsh elevation relative to sea level in each time step is used to determine which vegetation community is present, or if a mixed community exists, as defined by the parabolic relationship between species biomass and elevation (Supp. Fig. 1). Aboveground biomass for each time step is given by:
where aboveground biomass (B) [M L -2 ] is a function of elevation relative to sea level (zsl) and 𝛼𝛼, 𝛽𝛽, 𝜃𝜃 are coefficients describing the relationship between elevation and biomass at GCReW (Suppl. Table 1; Supp. Fig. 1; Byrd et al., 2017). When species growth curves overlap, resulting in mixed communities, the percentage of each species is calculated as Bi / (Bi + Bj). These percentages are used to calculate a weighted average for all species-specific parameters related to productivity and decomposition in mixed communities. Under eCO2 conditions, aboveground biomass in the model increases in C3
species by a factor of ~ 1.3 and does not increase in C4 species, as has been observed at GCReW (Drake, 2014;Langley et al., 2009).
Belowground biomass is estimated from aboveground biomass, where rhizomes and roots are considered separately. Rhizome biomass (𝑅𝑅 ℎ ) [M L -2 ] is proportional to aboveground biomass:
where g = 1 for ambient conditions and g = 2 for eCO2. Root biomass (𝑅𝑅 𝑜𝑜 ) [M L -2 ] is estimated using a balanced growth model that describes a functional equilibrium in which the mass and uptake of carbon by leaf tissue (B) is balanced by the mass and uptake of nutrients (µ) by root tissue (Reynold & D'Antonio, 1996;Reynold & Thornley, 1982;Agren & Ingestad, 1987).
where 𝜏𝜏 is carbon uptake net respiration, 𝜎𝜎 is an optimal tissue C:N ratio, and µ is a nitrogen uptake rate. Although nutrient uptake by root tissues should be dynamic and species specific, there is not enough information on S. patens and S. americanus belowground processes to treat these properly in the model. Instead, we use Michaelis-Menton kinetics to calculate a temporally constant N uptake rate based largely on measurements for S. alterniflora (Bradley & Morris, 1990;Giurgevich & Dunn, 1981).
This approach necessarily leads to a constant root:shoot ratio for each vegetation type and CO2 scenario (R:S = 1.10 for C4 vegetation, 1.79 for C3 vegetation, and 2.15 for C3 vegetation under eCO2). Together, root and rhizome biomass calculations produce a total belowground:aboveground biomass ratio of 2.10 for C4 vegetation, 2.79 for C3 vegetation, and 4.15 for C3 vegetation under eCO2. Like estimates of turnover described below, these model parameters are highly uncertain. Nevertheless, the belowground:aboveground biomass ratios used in the model are consistent with speciesspecific field measurements from the GCReW site and other Chesapeake Bay brackish marshes. For example, measured root:shoot ratios in these marsh types are generally between 2-5, higher in C3 than C4 vegetation, do not change consistently with flooding (Kirwan & Gutenspergen, 2015), and increase with eCO2 in C3 vegetation (Langley et al., 2013). Following (Morris & Bowden, 1986;Mudd et al, 2009), the total biomass of roots and rhizomes is distributed through the soil profile to a maximum rooting depth below the soil surface as:
where 𝛾𝛾 [L] is the depth at which belowground biomass decreases by approximately onethird (Mudd et al, 2009), and ds is depth below the soil surface. Ro and Rh are calculated at the surface (Roo, Rho) and in each cohort following the scheme of Morris & Bowden, (1986). The masses of live roots and rhizomes (𝑅𝑅 𝑜𝑜 , 𝑅𝑅 ℎ ) [M L -2 ] are then summed to determine the total biomass produced in each cohort and each time step.
Organic matter accumulation within the soil column depends on root and rhizome turnover, the organic content of allochthonous sediment deposited on the marsh surface, and decomposition of organic matter. Organic accumulation within a soil cohort (qo) [M L -2 T -1 ] is given as:
where To and Th [T -1 ] are the annual turnover of root and rhizome biomass respectively.
Belowground biomass enters the soil profile and contributes to organic matter fluxes through turnover. We assume that roots turnover more quickly than rhizomes, and that C4
plants have faster turnover than C3 plants (Table 1). This assumption follows the general observation that root turnover increases with decreasing diameter class (Gill & Jackson, 2000), and measurements at GCReW that indicate S. americanus rhizomes are larger in diameter than S. patens rhizomes (Curtis et al., 1990). Nevertheless, belowground turnover is a poorly understood process, and these parameters were necessarily chosen based on model hindcasts rather than direct measurements. The parameter a represents a fraction of organic matter that is composed of non-decaying components of plant tissue (e.g. silica) that remain as ash during loss-on-ignition (LOI) analyses. In the model simulations presented here, we assign a=0.08 based on LOI derived soil organic fractions that almost never exceed 0.92 at GCReW (Messerschmidt & Kirwan, 2020) and in diverse wetlands across the United States (Morris et al., 2016). The coefficient j is the organic fraction of suspended sediment so that 𝑗𝑗𝑞𝑞 𝑠𝑠 [M L -2 T -1 ] represents allochthonous organic matter. For simplicity, we chose an arbitrary low organic fraction of suspended sediment (j = 0.05) so that organic matter accumulation in the model is driven by autocthonous rather than allochthonous processes. Finally, 𝐷𝐷 [M L -2 T -1 ] is the total amount of organic matter decomposition in a soil cohort as described below.
Following previous approaches (Morris & Bowden, 1986;Mudd et al., 2009;Kirwan and Mudd, 2012), organic matter produced through root and rhizome turnover
) is split into fast decaying (labile), and slow decaying (refractory) pools. The model experiments here follow the parameterization of Mudd et al., (2009),
where the labile organic fraction is 0.84 and the refractory component is 0.16. We additionally consider labile (0.1) and refractory (0.9) fractions of allochthonous organic matter deposition (jqs). This parameterization is based on the observation that allochthonous organic matter is typically old and recalcitrant (Hopkinson et al., 2018;Van de Broek et al., 2018) but is of limited importance in these simulations where parameters are chosen to minimize allochthonous organic sediment deposition.
The decomposition rate of organic matter in each soil cohort (D) follows a linear decay model,
where Qo [M L -2 ] represents the total amount of organic matter in each soil cohort (i.e. qo summed through time). Like previous models, labile and refractory pools have different decay coefficients (kl and kr) [T -1 ]. However, we also modify the decay coefficient in each cohort by its depth below the surface, and the amount of aboveground biomass to represent the effects of soil priming via radial oxygen loss from plant roots. As plant biomass increases, aerobic leakage into the rhizosphere leads to an increase in decomposition (Jones et al., 2018;Mueller et al., 2015;Wolf et al., 2007). For each organic matter pool, a reference decay coefficient (kref) at the surface of the soil profile is Live roots (Ro) and rhizomes (Rh) are additionally considered in this calculation, so that our organic fraction is consistent with sediment core analyses (i.e. LOI) that typically do not separate live belowground biomass. The fraction of organic matter is converted to the fraction of carbon (cf) [unitless] following Craft et al., (1991):
The total amount of carbon in a given cohort (Corg) [M L -2 ] is given as:
and the total amount of carbon within the soil column is given as the sum of all carbon in each cohort:
where n is the number of soil cohorts. Finally, the carbon accumulation rate [M L -2 T -1 ] for a given year represents the total change in the mass of carbon in the soil column, and is calculated as Model experiments begin with a spin-up period in which the model creates an organic rich soil profile that develops under a constant, slow rate of sea level rise (1 mm yr -1 ). The model spinup starts with an assigned 1m soil profile composed of 1,000 1 mm thick cohorts with no organic matter. The model then populates the initial mineral stratigraphy with organic matter that evolves dynamically as the balance between organic matter production and decomposition. The model spin up ends when marsh accretion rates equilibrate to the rate of sea level rise and organic matter accumulation rates