Cohort Marsh Equilibrium Model (CMEM): History, Mathematics, and Implementation

Marsh accretion models predict the resiliency of coastal wetlands and their ability to store carbon in the face of accelerating sea level rise. Most existing marsh accretion models are derived from two parent models: the Marsh Equilibrium Model, which formalizes the biophysical relationships between sea level rise, dominant macrophyte growth, and elevation change; and the Cohort Theory Model, which formalizes how carbon mass pools belowground contribute to soil volume expansion over time. While there are several existing marsh accretion models, the application of these models by a broader base of researchers and practitioners is hindered because of (a) limited descriptions of how empirically derived ecological mechanism informed the development of these models, (b) limitations in the ability to apply models to geographies with variable tidal regimes, and (c) a lack of open‐source code to apply models. Here, we provide for the first time an explicit description of a mathematical version of the Cohort Theory Model and a numerical version of a combined model: the Cohort Marsh Equilibrium Model (CMEM) with an accompanying open‐source R package, rCMEM. We show that, through this “depth‐aware” model, we can capture how tidal variation impacts broad patterns of marsh accretion and carbon sequestration across the United States. The application of this model will likely be imperative in predicting the fate and state of coastal wetlands and the ecosystem services they provide in an era of rapid environmental change.


Introduction
Coastal marsh ecosystem dynamics are highly dependent on sea level, building elevation relative to the rate of sea level rise (i.e., accretion) or becoming overwhelmed by the inundation of the sea leading to coastal retreats.The resiliency of coastal marshes in the face of accelerating sea level rise, driven by anthropogenic greenhouse gas emissions, is a pressing societal concern (Fagherazzi et al., 2012;Kirwan & Megonigal, 2013;Kirwan et al., 2016).In an era of rapid environmental change and wetland loss, predicting the state and fate of coastal marshes using mechanistic models is of great importance.
The mechanism of marsh accretion is composed of two main, entwined factors: surface mineral deposition from tidal inundation and organic carbon content of the subsurface soil governed by plant inputs and decomposition (Kirwan & Megonigal, 2013;Figures 1a and 1b, Figure S1 in Supporting Information S1).Sediment suspended in the coastal ocean is carried over the marsh surface during tidal cycles and, as the water slows due to friction with the marsh surface and plants, settles onto the surface of the marsh providing fresh mineral deposition (Mudd et al., 2010).Flooding of the marsh surface is also a major driver of plant growth in marsh ecosystems, wherein the plant community is adapted to a window of optimal inundation such that plant productivity falls off if the marsh is too dry or too wet (Figure 1a).Plant productivity provides inputs to the soil subsurface via root growth and turnover.Organic matter inputs into the soil are reduced via decomposition by microbial organisms, and the balance of these inputs and outputs governs the total organic content of the soil (Figure 1b).Organic material is less dense than the mineral sediment input at the top of the soil (Morris et al., 2016), which changes the volume of the soil profile and results in a correlated change in the surface elevation of the marsh.Changes in marsh elevation then feedback to plant growth rates, as surface elevation and sea level rise jointly impact the amount of flooding plants experience.(a) Peak aboveground biomass is a function of the depth below mean high tide (i.e., flooding), which is both a function of the elevation of the marsh surface and tidal dynamics.The solid line represents where the marsh is stable and the dashed lines represent instability (i.e., plants can no longer produce enough biomass keep pace with sea level rise).(b) Soil cohorts are composed of four mass pools: live belowground biomass, fast decaying organic matter, slow decaying organic matter, and mineral sediments.Cohorts are tracked by their age (in years).At each annual timestep (t) a new layer of mineral sediment is added to the surface of the marsh.Then, belowground biomass is turned over within the rooting zone and is portioned into fast and slow organic matter pools, and is decayed.The proportion of fast organic matter that is decayed is then respired into greenhouse gases.Then, the cohort "ages" and the change in volume of the cohort due to biomass allocation and decay adjusts the depth of the top and bottom of the cohort in the subsequent timestep t + 1. (c) Timeline of journal articles and other model documentation of the building and fusion of the Marsh Equilibrium Model and the Cohort Theory Model from Morris et al. leading up to the development of the R package rCMEM, presented here.Events that were particularly influential in updating the model are shaded in blue.Small diagrams represent major foci of the papers (from left to right): development of the Cohort Theory Model, relationships between precipitation and macrophyte productivity, development of the Marsh Equilibrium Model, multispecies version of the Marsh Equilibrium Model, online interface of the Marsh Equilibrium Model, coupling of aboveground and belowground accretion processes, and the characterization of mineral and organic matter packing densities.
There are myriad marsh accretion models that characterize biophysical relationships between flooding, mineral deposition, soil biogeochemistry, and plant growth (e.g., Fagherazzi et al., 2012;Kirwan et al., 2010;Rietl et al., 2021;Swanson et al., 2014) and are used to predict marsh elevation gain through time and quantify the carbon sequestration potential of wetlands (Holmquist et al., 2018;Windham-Myers et al., 2018).Marsh accretion models are often composed of plant growth and organic matter components derived from two classic models: (a) the functional relationship between flooding and plant aboveground biomass (i.e., Marsh Equilibrium Model, Morris et al., 2002;Figure 1a) and (b) the tracking of belowground carbon pools via annualized "cohorts" (i.e., Cohort Theory Model, Morris & Bowden, 1986; Figure 1b).For more than three decades, empirical and modeling work has been built onto the original modeling scaffolds of the Marsh Equilibrium Model and the Cohort Theory Model, reflecting relevant advances in the fields of wetland biogeochemistry and geomorphology as new mechanisms emerge over time (Figure 1c; Morris, 2006Morris, , 2007Morris, , 2010;;Morris & Bradley, 1999;Morris & Bowden, 1986;Morris & Haskin, 1990;Morris et al., 2002Morris et al., , 2012Morris et al., , 2016;;Windham-Myers et al., 2018).Together, the Marsh Equilibrium Model and the Cohort Theory Model are critical for mapping changes in wetland area across space and/or time (e.g., Schile et al., 2014), predicting global carbon flux (e.g., Mack et al., 2023), and for understanding the role of organismal response to global change on emergent ecosystem processes (e.g., Vahsen, Blum, et al., 2023).
Despite considerable advances in marsh accretion models as detailed below, there are several aspects of the development of these models that limit their utility to a wider audience.First, to date, there has not been a full description of the mathematical dynamics of the combined Marsh Equilibrium Model and Cohort Theory Model despite its application in several research endeavors.Second, the application of the model across different geographies has not been particularly intuitive given the impact of locally driven tidal dynamics and their projections into the future given climate change.Finally, a focus on numerical simulations of the model for prediction and projection has limited mathematical advances to the model that would allow for the approximation of parameters that are difficult to estimate empirically.
Here, we detail the development of the Cohort Marsh Equilibrium Model (CMEM), a fusion of the Marsh Equilibrium Model and Cohort Theory Model, including both a mathematical formulation and release of the R package, rCMEM, as an open-source model (Holmquist et al., 2022).This effort represents the first articulation of the Cohort Marsh Equilibrium Model which has been previously used in a variety of studies to predict marsh accretion and carbon accumulation in several geographies across the United States (e.g., Schile et al., 2014;Vahsen, Blum, et al., 2023).We introduce the historical development of CMEM to provide context to its structure and mechanisms by tracing its origin from the original Marsh Equilibrium Model and Cohort Theory Model, identifying advances in the biogeosciences literature that align with major updates in the model (Figure 1c).Then, we present a novel mathematical version of the Cohort Theory Model, followed by a numerical model that informs the R package rCMEM, highlighting pertinent ecological drivers like tidal regime.Our goal is to provide modular and open-source tools for hindcasting and forecasting the linked processes of marsh resiliency and carbon sequestration given projected increases in the rates of sea level rise.Open-source tools like rCMEM are imperative in an era of unprecedented environmental change, wherein collaboration from multiple research groups is needed to improve prediction.

Background of the Cohort Marsh Equilibrium Model Parent Models: The Cohort Theory Model and the Marsh Equilibrium Model
The Cohort Marsh Equilibrium Model presented here can be seen as the culmination of research endeavors over the past several decades (Figure 1c).Morris and Bowden (1986) established the Cohort Theory Model, which breaks down marsh soils into mass cohorts that are tracked through the soil column by their depth and age.Each year, a "cohort" of mineral sediment settles on the surface of the soil column at age zero and depth zero to some depth, defining a single annual cohort.This mineral-defined cohort then gains organic inputs from the roots and rhizomes within the rooting zone and loses organic carbon via decomposition (Figure 1b).After the initial deposition, the cohort is isolated from the surface via new mineral deposition and sees no more mineral inputs.However, the cohort continues to receive root inputs and decomposition losses if it is within the rooting zone (Figure 1b).Due to this addition and removal of organic matter, the volume of the cohort expands and shrinks over time as it "moves" down the soil column (Figure 1b).The key feature of the Cohort Theory Model is that all matter and bulk density is distributed downcore via individual soil cohorts and thus, the depth of a given soil cohort is dependent on the dynamics (and resulting bulk density) of younger cohorts deposited above.
Much of Morris and colleagues' relevant work in the 1990s shifted to understanding what drives variation in aboveground biomass density, which influences organic matter inputs driving accretion rates (Figure 1c).Longterm experiments were critical to parameterizing interannual shifts in Spartina alterniflora primary production and subsequently concluded that these shifts were primarily linked to inundation rather than to shifts in salinity as previously proposed (Figure 1a; Morris & Bradley, 1999;Morris & Haskin, 1990;Morris et al., 2002).Morris et al. (2002) described the relationship between sea level and plant production as a hump-shaped parabola with an optimal growth condition (i.e., the peak of the parabola in Figure 1a) bound by a minimum and maximum depth below high tide within which vegetation can persist (i.e., x-intercepts of the parabola in Figure 1a), which was formalized into the first version of the Marsh Equilibrium Model.An important finding of this modeling exercise was that the flooding depth at which aboveground biomass is maximized is not likely optimal for the long-term trajectory of the marsh: the stabilizing feedback between the marsh surface elevation, biomass production, and sea level can be broken at certain flooding depths (Figure 1a, dashed line; Morris et al., 2002).Within the next few years, Morris enhanced the Marsh Equilibrium Model by extending the biomass function into a multi-species scenario (Morris, 2006(Morris, , 2007)).
While the Cohort Theory Model explicitly modeled belowground biomass as a major contributor to vertical accretion and carbon accumulation, the Marsh Equilibrium Model recognized how aboveground biomass mediated accretion through the sediment trapping from tidal flows.Subsequent work by Morris et al. integrated these two perspectives (Figure 1c).For example, in Morris et al. (2012) the amount of mass that enters into the marsh soil over time was partitioned into separate aboveground and belowground components.Further, Morris et al. (2016) developed a novel approach for calculating accretion rates using the separate, bulk self-packing densities of organic and inorganic material as well as loss-on-ignition (a proxy for organic matter fraction), which allowed for better predictions across different sea level rise scenarios and across different concentrations of suspended sediment in the water column.

Mathematical Model Formulation
Here we present a novel mathematical model to generalize the dynamics of the mass pools belowground.A key development of this mathematical model is the introduction of an age for depth substitution where instead of the depth of a specific cohort, we shift perspective to the age of the mineral component of a cohort.This mirrors the original Cohort Theory Model implementation and then builds on a discrete formulation by allowing the age of that cohort to shrink to an arbitrarily small size (i.e., 1 year to <1 day) to create a continuous age axis over which we can integrate.Given a long enough consistent mineral and root input (i.e., constant marsh elevation relative to sea level), we can then derive the "steady state" or stable profile of the surface cohorts providing expected downcore soil organic matter distribution and mineral age.Under these constant conditions we can also calculate an accretion rate (i.e., change in total soil volume over time) that would imply a constant sea level rise for this hypothetical marsh.
We derive the mathematical model from the following set of assumptions: 1.All mineral and organic deposits comprising the cohort of age a at time t are located at the same depth x(t, a) in the soil column.Cohorts of distinct ages do not mix.Each soil cohort in this model is initiated at age a = 0 when the associated mineral mass enters the system as a "fresh" laminae of minerals that are being laid on the soil surface (Figure 1b).The particular components of each cohort are the rapidly decaying (C f ) and slowly decaying (C s ) carbon deposits, the mineral deposits (M), and the roots and rhizomes (R). 2. The organic carbon enters the system through the influx of roots and rhizomes at the corresponding depth of the soil cohort.In this formulation, the "age" of the organic carbon is defined as the age of the mineral deposits located at the same depth x in the soil column.Thus, the concept of the "cohort age" is closer to the inorganic isotopic date (e.g., 210 Pb or 137 Cs) rather than the organic isotopic date (e.g., 14 C). 3. The organic and mineral components at depth x in the soil column have densities ρ c (x) and ρ m (x), respectively.4. The volumetric fraction R = R(x) of roots within the soil column is a function of depth x only.Specifically, the total volume occupied by the roots V R in a vertical soil column of fixed cross-sectional area A and depths z ∈ [0, x] is given by 5. The mass-age density of the mineral deposits of age a at time t is given by M(t, a), so that the total mass of mineral deposits of all ages up to age a (equivalently, located at depths z ∈ [0, x(t, a)]) is given by ∫ a 0 M(t,s) ds, where s is a dummy variable for integration, and this mineral mass occupies the respective volume V M , where ds.
6. Similarly, the mass-age density of the organic (carbon) deposits of age a at time t is given by C f (t, a) + C s (t, a), so that the total mass of carbon deposits of all ages up to age a (equivalently, located at depths z ∈ [0, x(t, a)]) is given by ∫ a 0 {C f (t, s) + C s (t, s)} ds, and this organic mass occupies the respective volume V C , where 7. The conservation equation represents the volume of a vertical soil column of fixed cross-sectional area A and depths z ∈ [0, x(t, a)] as a sum of the volumes occupied by carbon deposits, mineral deposits, and the roots, respectively: (1) In the subsequent analysis, we will scale the model so that A = 1 and neglect the effects of compression so that the densities of all deposits remain constant throughout the entire soil column.With these simplifications, the conservation equation can be expressed in the integral form as follows: (2) 8. We further assume that the volumetric fraction of root density decays linearly with depth: so that the roots are present in the soil column only up to depth x max .
For mathematical convenience, we introduce a function which represents the total volume occupied by carbon and mineral components at depths z ∈ [0, x].Our assumption that 0 < r m < 1 implies that ϕ(x) is a monotonic increasing function of x with ϕ(0) = 0 and . Furthermore, the inverse function can be explicitly expressed as where α = 1 r m and β = r m 2x max .
To obtain a differential form of the conservation equation, we first make a substitution u(t, a) = ϕ(x(t, a)), so that and then differentiate on both sides to arrive at The time-age dynamics of the mineral pool is modeled by the equation where ζ(t) corresponds to the surface mineral loading rate at time t.This equation can be solved directly to yield the solution This formulation assumes that there is no loss of mineral matter (e.g., from wave action) once sediment is deposited on the marsh surface.The time-age dynamics of the fast and slow carbon pools is modeled by the equations where, after simplifications, we have Here, we also assume that the fast carbon component decays at an exponential rate k f , while the decay of the slow carbon component in the timescale of our model is neglected altogether.We argue that the assumption of negligible decay from the slow carbon pool is reasonable because the realistic timescale of marsh development/ degradation is on the order of decades to centuries which is much faster than the typical timescale of slow carbon decay (i.e., millennia; Kirwan & Mudd, 2012).
Combining the above equations, we obtain the differential form of the cohort model: Journal of Geophysical Research: Biogeosciences where we can track the depth of the cohort of age a via x(t, a) = ϕ 1 (u(t, a)).
We have purposefully neglected addressing units explicitly in the above formulation, consistent with common mathematical presentations, to allow the reader to focus on the kinetics and formulation of the model.In practice, the model proceeds on an annual timestep (a and t are measured in years), mass pools are measured in grams, and depth (and thus area and volume) is measured in centimeters.

Stationary Age Distribution With Temporally Constant Mineral and Root Inputs
It is clear that the system (Equations 4-7) does not allow a true equilibrium solution because the mineral content in the soil column is being top-loaded at a time-dependent rate and as a consequence, the total mineral content is steadily increasing.In addition, the "slow" carbon component is being constantly supplied by the root structure within the soil column and its decay rate is simply neglected in our model.As a consequence, the total "slow" carbon content is also steadily increasing, so "slow" carbon in our model is rather a permanent carbon content.
Nonetheless, if we assume that the loading rate is constant M(t, 0) = ζ > 0 (a situation which we refer to as "the steady-state conditions"), then we can find the corresponding solution to Equations 4-7 where the agedistribution (and thus also depth-distribution) of the cohorts is time-invariant (i.e., stationary).
Assuming such steady-state conditions allows to us uncouple the age-dynamics of the organic components and the roots from the mineral dynamics and reduces Equations 4-7 to a system of ordinary differential equations: Due to the nonlinear nature of the system (Equations 8-10), we cannot obtain closed-form expressions for the components of the time-invariant age distribution.However, in the next section we discuss some relevant properties of such distribution for data estimation.

Parameter Estimation Informed by the Mathematical Model
In this section, we show how Equations 8-10 from the mathematical model can be used to infer some of the parameters from the data without solving these nonlinear equations explicitly.For instance, it follows directly from Equations 8-10 that all three functions u(a), C f (a) and C s (a) are increasing initially when a > 0 is small, and thus all three functions are positive on the interval (0, +∞).Furthermore, since u′(a , that is, u(a) is increasing at least as fast as a linear function, and there exists a unique critical value a max such that u(a max ) = u max and u(a) > u max for all a > a max .This, in turn, implies that no roots are present in cohorts older than a max , that is, R(a) = 0 for a > a max , where we use R(a) in place of R ϕ 1 (a))) for notational brevity.Consequently, for all ages a > a max , C f (a) decays exponentially to zero while C s (a) remains constant: We let C(a) = C f (a) + C s (a) be the cumulative content of carbon matter in age cohort a.We also let a x be the age at which the organic matter is maximized, that is, The critical value a max naturally corresponds the maximal age of the root biomass, thus 0 < a x < a max .
Due to Equation 11we have C f (∞) = 0 and C s (∞) = C s (a max ), and Hence, we have that C(∞) = C s (∞), which allows us to express Now we can estimate the rate of organic decay k f by evaluating or equivalently, To evaluate the fast carbon root allocation fraction f f , we integrate Equation 10 from 0 to ∞ to find that Substituting C s (∞) = C(∞), we obtain the estimate Finally, we integrate Equation 10 from 0 to a x , and then use Equation 13 to express and then use the equation to estimate the rate of root turnover as Equations 11-14 demonstrate that the downcore patterns of soil organic carbon density and mineral age can be used to determine key parameters in the Cohort Marsh Equilibrium Model, with the critical assumption that the Journal of Geophysical Research: Biogeosciences 10.1029/2023JG007823 marsh has kept up with sea level rise consistently over the last several decades (or as long as it takes for the mineral cohort to travel past the rooting zone and lose the associated fast carbon).In practice, this key assumption may not hold for a location of interest and classic parameterizations via incubation studies may be more appropriate.

Forcings and Feedbacks
A key component of Cohort Marsh Equilibrium Model is a feedback between the sea level and marsh platform elevation via changes in biomass and mineral inputs, which is derived from the Marsh Equilibrium Model.We note that this feedback is not explicitly represented in the mathematical formulation presented above (although the authors hope to see this extended in future works).Thus, the purpose of the mathematical model described above is to better characterize mass pools, while the purpose of the numerical model in this section is to make predictions of changes in mass pools over time, and in turn, changes in marsh surface elevation and carbon accumulation, while accounting for biophysical feedbacks between plant production, sea level rise, and marsh surface elevation.
The numerical model also represents what is featured in the R package rCMEM, wherein there are two main dynamic modules of the model-the sediment and biomass modules-both of which are forced by sea level rise.
We include a supplemental vignette in Supporting Information S2 that walks through a full workflow of using CMEM to aid users of the R package for the current version.In this manuscript, we do not include the function names and calls specific to the rCMEM package as this is a living codebase that may be edited and versioned in the future.Instead, we walk through the mathematical formulation of the sediment and biomass input modules as a function of tidal inundation.We note that the units for mass (grams) and length (centimeters) are kept consistent throughout the inputs, internal processing, and outputs of CMEM and rCMEM.

Dynamic Feedbacks: Sediment Module
The sediment module calculates the amount of inorganic sediment (i.e., minerals) captured on the marsh surface due to tidal flows.The amount of water above the marsh surface for a set area controls the amount of sediment delivered and is governed by variability across daily and monthly tidal cycles (Figures 2a and 2b), and trends across decadal and centennial timescales influence annual tidal predictions (Figures 2c and 2d).Below, we walk through classical formulations for tidal inundations at sub-annual, annual, and multi-year scales.Then we discuss the sediment input at the marsh surface as a function of inundation time, water column depth, and capture rate.

Sub-Annual Tidal Inundation
A single tidal flood time for a marsh at elevation Z(t) is broken down into three cases, where ω is the fraction per tide that the marsh spends below the water at year t given mean high Z h (t) and low Z l (t) water.Alternatively, a trigonometric function can be used to characterize fractional flooding time (sensu Hickey, 2019), which is a new functionality within rCMEM.
In the trigonometric case, flooding time is calculated as a function of the absolute value of the rising time ω r minus ϕ, which is the time of one half of a tidal cycle, and the falling time ω l for tidal cycle class (sensu Hickey, 2019), Journal of Geophysical Research: Biogeosciences 10.1029/2023JG007823 VAHSEN ET AL.H where H 1 and H 2 are intermediate functions used to translate Z, Z h , and Z l data to rising and falling times (ω r and ω l , respectively).When ϕ is set to 0.5 (i.e., half of a tidal cycle), the value of ω(t) is equal to the fractional flooding time.
This trigonometric option allows sedimentation processes to be more sensitive in high marshes by smoothing transitions at the outer boundaries of the flooding extent, but does not drastically influence calculated fractional flooding times.Fractional flooding time is capped at 1 if elevation is lower than lowest flood elevation (100% flooded), and 0 (0% flooded) if the elevation is greater than the highest water level.

Annual Number of Tidal Cycles
The number of tides per year n t also regulates the amount of flooding on the marsh surface.Estimating n t is specific to geography.Semi-diurnal and mixed semi-diurnal tidal regimes are the dominant tidal type in the U.S. Southeast-wherein the Marsh Equilibrium Model was first calibrated (Morris et al., 2002)-and have two floods per lunar day.In contrast, diurnal tides, which are dominant along much of the U.S. Gulf Coast, have only one flood per lunar day.Estimating the flooding time for these regions would thus be inaccurate using the same number of tidal floods as for semi-or mixed semi-diurnal tides.
This version of CMEM and rCMEM allows for n t to be specified to better align with observed tidal scenarios across geographic locations by introducing three alternative tidal cycle scenarios governing annual tidal inundation patterns: mean high tides (352.657 per year), mean higher high tides (327.937 per year), and mean higher high spring tides (24.720 per year) (Figure 2b).The fractional flooding time ω is then calculated and summed across scenarios resulting in a total fractional flooding time.This new adaptation within rCMEM especially improves the realism of the model when the starting elevation of the marsh surface is above the high water line (i.e., a supratidal marsh) because it allows for the marsh surface to flood periodically, allowing some mineral sedimentation and adjusting the functional response of vegetation to that flooding.

Multi-Year Tidal Amplitude Change
High tide flood dynamics can also vary on an 18.61-year cycle known as the lunar nodal cycle (Figure 2c) following where Z h (t) is the high tide level, Ω μ is the mean tidal amplitude, λ is the amplitude of the 18.61 years cycle, θ is the phase of the 18.61 years cycle, and Z μ (t) is mean sea level, all at time t.

Mineral Inputs
The amount of available mineral sediment inputs m a at an elevation in a given year to be delivered to the marsh surface is calculated as where the available suspended sediment input is the fractional flooding time ω multiplied by the suspended sediment concentration m and the capture rate of the marsh τ, if the fractional flood time ω (see Equation 15) is smaller then the inverse of the number of times the sediment settles out of the water column for each tidal cycle (i.e., capture rate, τ).Alternatively, if the inverse of capture rate τ is less than the fractional flooding time, then the full available sediment is deposited on the surface.In practice, capture rate (τ) serves as a tuning parameter and while it has a basis in physical processes, it will only be meaningful if constrained empirically.
The sediment that is delivered to the surface of the marsh for any given year is then where M is the mineral input, m a is the available suspended sediment input (see Equation 22), n t is the number of tides per year, and the height of the average water column over the marsh is the elevation of high tide Z h (t) minus the elevation of the marsh Z(t) divided by 2 to account for cyclic oscillations.

Dynamic Feedbacks: Biomass Module
Biomass production is a function of flooding and is optimized at an intermediate level of flooding-higher in the tidal frame plant growth is stunted and lower in the tidal frame productivity is reduced and mortality increases (Figure 1a).Because the elevation of the surface, the elevation of sea level, and the tidal range vary from year to year and across geographies, we calculate the elevation of the biomass parabola as dimensionless (Z*), relative to the tidal range.Relative tidal elevation is calculated as where Z* is relative tidal elevation, Z is the surface elevation, Z μ is mean sea level, and Z h is the mean high water level.Z* values are convenient for comparisons across geographies: Z* values between 0 and 1 are between mean sea level and the elevation at which the surface receives twice daily tidal floods.Z* > 1 indicates elevations above the mean high water line, and Z* < 1 denotes elevations that are approximately below low tide.
In the original Marsh Equilibrium Model, the parabolic relationship between Z* and aboveground biomass depended on the elevation limits of vegetation Z min and Z max and the maximum amount of biomass that can be produced across elevations B max (Figure 3, black curve).Under a scenario of increased sea level rise (as is generally predicted), Z min will characterize the survival potential of the marsh because marshes that do not produce vegetative biomass will transition into a mudflat.In this version of the Cohort Marsh Equilibrium Model, we introduce an additional fourth parameter Z peak, or the elevation at which biomass production is optimal (i.e., where B = B max ) to allow for an asymmetric "parabola shape" (Figure 3, gray curves).We define Z range , Z sum , and Z prod to simplify the calculation of aboveground biomass via a piecewise function. .
Aboveground biomass is then converted to belowground biomass via a root-to-shoot ratio (φ; technically a ratio of root and rhizome biomass to aboveground biomass), such that Total root and rhizome biomass (r tot ) is the integrated sum of the root density described in Equation 3, implying that r m = 2 r tot (Z * ) Ax max .An alternative implementation of an exponential decay for this root distribution is included in the R package, rCMEM.
To account for species composition in overall ecosystem biomass production we consider species specific parameters: Z max , Z min , Z peak , B max , φ, and x max (Equations 25-28).We then constructed a "competition" function, where the total aboveground biomass is equal to the predicted aboveground biomass of the most dominant species at the current level of flooding allowing for multiple species to dominate within their niche flooding level (Figure 4; sensu Morris, 2006).

Boundary Conditions and Numerical Discretization
In prior sections we dealt primarily with range change for the various pools in the models and how they interrelate.We now consider three key boundary conditions for the Cohort Marsh Equilibrium Model-initial soil organic matter and mineral pools downcore, sea level rise, and the fixed bottom assumption.
To generate the initial distribution of soil organic matter and mineral pools downcore, we first assume that the surface elevation of the marsh is below the high tide and then run the simulation forward under a prescribed sea level rise scenario.During year one, we assume there is a single mineral cohort based on annual mineral inputs calculated from the initial elevation.Cohorts are then added to the profile until three conditions are met: (a) the profile meets a minimum age of 50 years, (b) the profile is deep enough to accommodate the maximum rooting depth, and (c) the fast and slow organic matter pools at the bottom of the profile have stabilized.The initial elevation can also be outside the range of the biomass parabola (Figure 3), however only a mineral sediment lamina will be added to the column until the relative tidal elevation crosses back into the range of the biomass parabola.An alternative initialization could leverage the 'stable' formulations given in Equations 8-10, however the above scenario allows for historical sea level patterns to be directly incorporated.
Sea level rise is then projected as a quadratic increase over time (Sweet et al., 2017).Mean sea level at time Z μ (t) is calculated as where Z μ (0) is initial mean sea level, r 0 is the initial rate of sea level rise, R is the total amount of sea level rise, and T is the total number of years in the scenario.Finally, to calculate the change in marsh elevation we assume that the bottom cohort is stationary and that any change in the total core volume results in a change in the surface of the marsh.

Example Scenario Simulation: Geographic Specific Parameterization
To examine the impact of realistic differences in sea level rise on marsh accretion and soil carbon storage, we simulated the numerical version of the Cohort Marsh Equilibrium model across coastal cities in the United States informed by site-specific tidal histories, while keeping other model parameters constant.Specifically, we extracted verified tidal data from 1983 to 2020 to inform the tidal parameters of seven locations from the contiguous United States which represented a variety of coastal geographies: Seattle, Washington; San Francisco, California; Port Isabel, Texas; Pensacola, Florida; Charleston, South Carolina; Annapolis, Maryland; and Portland, Maine (Figure 5a, Table 1; Stephenson, 2016).
We fit harmonic constituents to verified tidal data from 1983 to 2020 using the R package TideHarmonics (Stephenson, 2016) and used a moving window analysis to identify tidal data for each site: average mean tidal level (which we treated as equivalent to mean sea level), mean high water, mean higher high water, and mean higher high water spring tide.For the two gauges with diurnal tidal regimes (Port Isabel, TX and Pensacola, FL), we only reported mean sea level and mean high water data (Table 1).Across sites, we used mean tidal level at 2000 as the initial mean sea level.We estimated initial sea level rise for each gauge as the slope of a fitted linear relationship between median decimal year and mean tidal level from 1983 to 2000 and we calculated the annual mean tidal amplitude as the difference between annual mean tidal level and the mean high water data.
We applied a 100-year sea level rise scenario with the localized scenario at each location, using median probability realized concentration pathway (RCP) 4.5 (Kopp et al., 2014).To test whether the lunar nodal cycle was an important driver by location-and thus determine if it should be parameterized for the model simulations-we fit a sinusoidal function to the relationship between mean tidal amplitude and year using the non-linear solver nls function in R (version 4.3.2;R Core Team, 2023).If lunar nodal phase and amplitude terms were significantly different than zero (p < 0.05) then we used the estimated coefficient values as inputs, otherwise we set lunar nodal amplitude as zero (Table 1).
We derived the relationship between relative tidal elevation (Z*) and aboveground biomass for the simulations from a Spartina alterniflora mesocosm experiment at the North Inlet (Morris et al., 2013).We fit an asymmetric parabolic relationship (Equations 25-27) with a Bayesian hierarchical model to the biomass-elevation data from the mesocosm experiment via rjags (Plummer, 2021) S1).We used mean high water at 2020 and Z* of 0.83 (the median peak elevation from the biomass analysis Z * peak ; Table 2) to calculate a unique starting elevation for each site.While we use data from the Morris et al. (2013) Spartina alterniflora mesocosm experiment to inform the parameter values of the biomass-elevation relationship, we note that previous applications of the model have used data from other mesocosm experiments with different species (Vahsen, Blum, et al., 2023) or from observational data that characterize a parabolic relationship with tidal elevation and aboveground biomass (Schile et al., 2014).All other parameters and inputs were kept constant for the simulations across different geographic locations (Table 2).
Finally, to visualize the underlying model dynamics that inform the emergent accretion and soil carbon storage rates, we visualized the four mass pools tracked by the Cohort Marsh Equilibrium Model-belowground biomass,  1 and 2).(a) Map of locations for which the model was simulated given geographically specified tidal parameters.(b) Resulting relative tidal elevation (dimensionless) for each location (calculated following Equation 24) across the simulation period.(c) Carbon flux (g m 2 yr 1 ; i.e., the change in the amount of carbon in the system across an annual timestep) for each location across the simulation period.
fast organic, slow organic, and mineral-for each cohort by depth for Annapolis, MD and Seattle, WA: two sites with different resilience to sea level rise based on the localized sea level rise scenarios applied here.

Numerical Model Simulations
The model simulations using data from seven NOAA tide gauges show that the variability in future marsh resiliency is strongly influenced by geographically explicit inundation drivers (Figure 5).Qualitatively, the simulated marshes in Annapolis, MD, Port Isabel, TX, and Pensacola, FL showed an eventual submergence of the marsh surface (Figure 5b; Relative Tidal Elevation <0) and loss of carbon sequestration capacity, as depicted by negative carbon flux peaks in Figure 5c, wherein plant biomass is fully depleted.The Charleston, SC and San Francisco, CA simulated marshes became more flooded over the course of the simulation and contributed progressively less carbon each year due to reductions in biomass production under increased flooding, however, these marshes did not become submerged within the 100-year simulation.The Seattle, WA and Portland, ME simulated marshes were the most resilient to projected localized sea level rise, maintaining their surface elevation nearly constant relative to sea level and maintaining carbon sequestration capacity.Finally, the importance of significant lunar nodal cycle oscillations on relative tidal elevation and carbon flux can be seen in Figures 5b and  5c for all locations except Annapolis, MD and San Francisco, CA, wherein lunar nodal phase and amplitude were not found to be significantly different than zero.Differences in relative tidal elevation and carbon flux coincided with differences in mass pools belowground, elucidating the cohort mechanism within the Cohort Marsh Equilibrium Model.In Figure 6, we compared the composition of cohorts downcore at the end of the simulation for Annapolis, MD and Seattle, WA-two marshes with different accretion and carbon storage trajectories as shown in Figure 5.At the end of the simulation, the simulated marsh in Annapolis submerged and thus was no longer producing any belowground biomass (Figure 6a) and all fast organic matter was decomposed (Figure 6b).Alternatively, the simulated marsh in Seattle maintained belowground biomass and fast organic matter pools that generally declined in mass moving downcore until the end of the rooting zone at 30 cm below the marsh surface.
Once the Annapolis marsh submerged and no longer supported vegetative production, new cohorts were only composed of mineral matter, which increased with increasing sea level rise as the marsh surface continued to decline in relative tidal elevation (Figure 6d).The lack of biomass production also meant that new cohorts in the Annapolis marsh following submergence were not contributing to slow organic matter pools, but due to its recalcitrant nature, slow organic matter was conserved within cohorts deeper in the marsh soil (Figure 6c).Alternatively, maintained vegetation in the Seattle marsh led to sustained contributions to the slow organic matter pool (Figure 6c) and relatively less mineral mass in the top cohorts due to being higher in relative tidal elevation (Figure 6d).

Discussion and Conclusion
The Cohort Marsh Equilibrium Model described here is the first open-source fusion of the Marsh Equilibrium Model (Morris et al., 2002) and Cohort Theory Model (Morris & Bowden, 1986) which revolutionized our understanding of the link between soil biogeochemistry and surface elevation in tidal ecosystems.These foundational models are merged through the development of a "depth-for-time" substitution, linking key feedbacks between soil cohort biogeochemistry, root inputs, and relative sea level.While versions of this joint model have been used to predict marsh resilience in the face of sea level rise (e.g., Schile et al., 2014), we hope that providing an open-source codebase accompanied with a more detailed discussion of the mechanisms and equations motivating the model dynamics will improve the usability of the model across the research community.
By introducing a novel mathematical version of the cohort theory model, we provide unique opportunities for data integration to estimate parameters in the Cohort Marsh Equilibrium Model that are usually difficult to estimate.Particularly, the decay rate of the fast pool κ f , turnover time of the root mass ψ, and allocation of the dead roots to the fast pool f f could be informed by the organic carbon by age slope across different sections of the soil core.Characterizing κ f , ψ, and f f is important because these parameters are typically informed by single studies or expert judgment, yet are integral in accurately estimating carbon stocks.For example, Rietl et al. (2021) showed through model simulations that characterizing species-specific decomposition rates was critical in capturing carbon storage dynamics in a coastal marsh in the Chesapeake Bay.Although more work is needed in developing the methodologies for estimating parameters from core data, we feel this is a promising line of research that could provide hyper-local parameters for soil kinetics.
We recognize that there are several fundamental assumptions that are made in this mathematical formulation that could be revisited in future work including the assumption of a negligible decay rate of the slow organic matter pool, a constant organic and mineral packing density downcore (i.e., ignoring the effects of compaction), and constant downcore decay rates (i.e., ignoring the effects of changing redox dynamics).Additionally, the mathematical formulation in its current state does not account for biophysical feedbacks between plant production, sea level rise, and marsh surface elevation, although this is a near-term goal for our research team.Despite these limitations, this Cohort Marsh Equilibrium Model formalizes past patterns in the parent models and allows us to integrate insights from both the Cohort Theory Model and Marsh Equilibrium Model under a single framework.
Our geographic simulations informing the Cohort Marsh Equilibrium Model with localized inundation patterns confirm previous findings that wetland resilience is particularly sensitive to localized sea level rise scenarios (Figure 5; Holmquist et al., 2021).Kirwan et al. (2016) similarly suggested that parameterizing dynamical models, like the Cohort Marsh Equilibrium Model, with localized sea level rise and tidal ranges is needed to correctly assess marsh resilience.Together, these findings call for site-specific management plans and tailored models to inform them.
We note that there are indeed still limitations to our numerical modeling framework in characterizing the local dynamics needed to accurately predict regional marsh resilience.First, this model is aspatial, that is, the Cohort Marsh Equilibrium Model does not account for horizontal transgression, which has been shown to be a critical mechanisms for predicting how marshes will expand or contract in area in response to sea level rise ( Holmquist et al., 2021;Miller et al., 2021).Second, the Cohort Marsh Equilibrium Model lacks some nuance in vegetation parameters.It does not explicitly account for sediment trapping by vegetation aboveground, a feature that is represented in similarly minded models (Mudd et al., 2010).Thus, this version of the Cohort Marsh Equilibrium Model could be underestimating the capacity for marshes to accrete via mineral sedimentation on the marsh surface.Additionally, the root-to-shoot ratio in the model is currently set to be static across inundation regimes, however there is evidence from mesocosm experiments that root-to-shoot ratios for some marsh species decline with increasing inundation (Vahsen, Kleiner, et al., 2023) and that variation in root-to-shoot ratios can lead to uncertainty in model predictions (Vahsen, Blum, et al., 2023).We hope that future work assessing the importance of and accommodating these mechanisms will be made easier via an open-source code base.
Overall, open-source code bases such as rCMEM allow users to modify existing model parameters and underlying code to suit local conditions, creating opportunities for relevant simulations on a more regional level.Open-source projects can support vibrant developer communities that intersect policy and academic research spheres.Finally, we see future improvements of model forecasts by adopting a full data-model integration approach (e.g., ecological forecasting; sensu Dietze, 2017), which would account for sensitivies and uncertainties of parameters and is made more feasible with an open-source version of Cohort Marsh Equlibrium Model.We are excited by this new potential for open science.

Figure 1 .
Figure 1.Conceptual diagrams of foundational components of the Marsh Equilibrium Model (a) and the Cohort Theory Model (b) and the timeline of model development and merging (c).(a) Peak aboveground biomass is a function of the depth below mean high tide (i.e., flooding), which is both a function of the elevation of the marsh surface and tidal dynamics.The solid line represents where the marsh is stable and the dashed lines represent instability (i.e., plants can no longer produce enough biomass keep pace with sea level rise).(b) Soil cohorts are composed of four mass pools: live belowground biomass, fast decaying organic matter, slow decaying organic matter, and mineral sediments.Cohorts are tracked by their age (in years).At each annual timestep (t) a new layer of mineral sediment is added to the surface of the marsh.Then, belowground biomass is turned over within the rooting zone and is portioned into fast and slow organic matter pools, and is decayed.The proportion of fast organic matter that is decayed is then respired into greenhouse gases.Then, the cohort "ages" and the change in volume of the cohort due to biomass allocation and decay adjusts the depth of the top and bottom of the cohort in the subsequent timestep t + 1. (c) Timeline of journal articles and other model documentation of the building and fusion of the Marsh Equilibrium Model and the Cohort Theory Model from Morris et al. leading up to the development of the R package rCMEM, presented here.Events that were particularly influential in updating the model are shaded in blue.Small diagrams represent major foci of the papers (from left to right): development of the Cohort Theory Model, relationships between precipitation and macrophyte productivity, development of the Marsh Equilibrium Model, multispecies version of the Marsh Equilibrium Model, online interface of the Marsh Equilibrium Model, coupling of aboveground and belowground accretion processes, and the characterization of mineral and organic matter packing densities.

Figure 2 .
Figure 2. Drivers of variation in water level across four temporal scales using tidal gauge data from NOAA Station 8665530 in Charleston, South Carolina.(a) Water level (m) derived from a data logger across a lunar day (24.83hr).(b) Water level (m) derived from a data logger across a year.Horizontal lines represent different tidal data: mean sea level (MSL), mean high water (MHW), mean higher high water (MHHW), and spring mean higher high water (MHHWS).(c) Tidal amplitude (m) across an 18.6 years lunar cycle.(d) Predicted change in mean sea level within a century.The black bolded line represents mean predictions and the shaded band represents a 95% confidence interval.Fluctuations in the bands represent the 18.6 years lunar cycle as highlighted in (c).

Figure 3 .
Figure 3. Predicted aboveground biomass g cm 2 ) as a function of relative tidal elevation (Z*) for three hypothetical biomasselevation functions.The classic Marsh Equilibrium Model parabolic relationship where Z * peak = Z * max +Z * min

Figure 4 .
Figure 4. Predicted aboveground biomass g cm 2 ) as a function of elevation (cm NAVD88) for two species (left; green lines) and the total aboveground biomass used to parameterize the model given the "competition function" (right; black line) as specified in the Cohort Marsh Equilibrium Model.

Figure 5 .
Figure 5. Resulting trajectories of Cohort Marsh Equilibrium Model simulations from 2000 to 2100 given the same initial elevation normalized to tidal amplitude (Relative Tidal Elevation [Z*] = 0.83) and locally specific tidal amplitudes, sea level rise scenarios, and lunar nodal cycle parameters (Tables1 and 2).(a) Map of locations for which the model was simulated given geographically specified tidal parameters.(b) Resulting relative tidal elevation (dimensionless) for each location (calculated following Equation24) across the simulation period.(c) Carbon flux (g m 2 yr 1 ; i.e., the change in the amount of carbon in the system across an annual timestep) for each location across the simulation period.

Figure 6 .
Figure 6.Predicted mass for each of the four mass pools: (a) belowground biomass, (b) fast organic, (c) slow organic, and (d) mineral, across cohorts in 2100 for Annapolis, MD and Seattle, WA.Points represent individual cohorts and are arranged on the y-axis by the depth at the top of each cohort.In 2100, the Annapolis simulated marsh became submerged and lost all vegetation, and the Seattle simulated marsh was still keeping pace with sea level rise and retained vegetation.Note that the scales of the x-axes vary across panels (a-d).
to inform parameter estimates of B max , Z * max , Z * min , and Z * peak following methods from LeBauer et al. (2013) (see Supporting Information S1: Text S1 and Table

Table 1
Tidal Driver Parameters and Inputs for Geographic Simulation of the Cohort Marsh Equilibrium ModelNote.r 0 = initial rate of sea level rise (cm yr 1 ), R = total amount of sea level rise (cm), Z 0 = initial elevation (cm), Z μ = initial mean sea level, Z μ,datum = mean sea level over the last recorded tidal datum period, Z h,datum = mean high water over the last recorded tidal datum period, Z hh,datum = mean higher high water over the last recorded tidal datum period, Z hhs,datum = mean higher high spring tide over the last recorded tidal datum period, λ = lunar nodal amplitude, θ = lunar nodal phase, n t = number of tides per year, τ = capture rate.All marsh surface and tidal elevation values are relative to the NAVD88 datum.Pensacola, FL and Port Isabel, TX do not have values for Z hh,datum and Z hhs,datum because they have diurnal tidal regimes.

Table 2
Non-Tidal Parameters for Geography Simulation Note.These parameters were held constant across the different geographies.