Reviews and syntheses: ²¹⁰Pb-derived sediment and carbon accumulation rates in vegetated coastal ecosystems – setting the record straight

Vegetated coastal ecosystems, including tidal marshes, mangroves and seagrass meadows, are being increasingly assessed in terms of their potential for carbon dioxide sequestration worldwide. However, there is a paucity of studies that have effectively estimated the accumulation rates of sediment organic carbon (C_org), also termed blue carbon, beyond the mere quantification of C_org stocks. Here, we discuss the use of the ²¹⁰Pb dating technique to determine the rate of C_org accumulation in these habitats. We review the most widely used ²¹⁰Pb dating models to assess their limitations in these ecosystems, often composed of heterogeneous sediments with varying inputs of organic material, that are disturbed by natural and anthropogenic processes resulting in sediment mixing and changes in sedimentation rates or erosion. Through a range of simulations, we consider the most relevant processes that impact the ²¹⁰Pb records in vegetated coastal ecosystems and evaluate how anomalies in ²¹⁰Pb specific activity profiles affect sediment and C_org accumulation rates. Our results show that the discrepancy in sediment and derived C_org accumulation rates between anomalous and ideal ²¹⁰Pb profiles is within 20 % if the process causing such anomalies is well understood. While these discrepancies might be acceptable for the determination of mean sediment and C_org accumulation rates over the last century, they may not always provide a reliable geochronology or historical reconstruction. Reliable estimates of C_org accumulation rates might be difficult at sites with slow sedimentation, intense mixing and/or that are affected by multiple sedimentary processes. Additional tracers or geochemical, ecological or historical data need to be used to validate the ²¹⁰Pb-derived results. The framework provided in this study can be instrumental in reducing the uncertainties associated with estimates of C_org accumulation rates in vegetated coastal sediments.

1 Introduction

Recognition of the globally significant role of vegetated coastal habitats, including tidal marsh, mangrove and seagrass, as sinks of carbon dioxide (CO₂) (Duarte et al., 2013) has generated interest in the potential to mitigate CO₂ emissions through the management of these ecosystems, an approach termed the “blue carbon” strategy (Duarte et al., 2013; Mcleod et al., 2011; Nellemann et al., 2009). However, efforts to include vegetated coastal ecosystems in existing carbon mitigation strategies have met with an important limitation: there is a paucity of estimates of C_org sequestration rates, particularly in seagrass habitats (Johannessen and Macdonald, 2016, 2018; Macreadie et al., 2018).

Two interrelated measurements of importance to the rate of C_org sequestration are the sediment C_org content and the sediment accumulation rate. To date, most research has focused on the first term, which provides information about the C_org stock sequestered in sediments (Howard et al., 2014; Pendleton et al., 2012). However, C_org stocks alone cannot be used to fully assess the C_org storage capacity or to establish comparisons among sites. Measurements of C_org accumulation rates (CARs) address the question of how much C_org is sequestered in a specified time period and quantify the ongoing sink capacity. In general, CAR is obtained in one of two ways: by measuring the concentration of C_org in sediments and ascribing dates to either the entire profile of interest or to specific intervals; or by estimating sediment accumulation rates. Determination of mean CAR is partially dependent on the timescale of interest and the dating methods used. ²¹⁰Pb has been shown to be an ideal tracer for dating aquatic sediments deposited during the last ca. 100 years, providing a time frame compatible with management actions (Marland et al., 2001) and enabling the determination of CAR and its changes with time due to natural or human impacts. Due to the relatively long integration period (decades to a century), mean ²¹⁰Pb-derived CAR estimates are not affected by interannual variability, hence allowing for the assessment of shifts in C_org accumulation from the “baseline” condition (such as the C_org that naturally cycles through an ecosystem; Howard et al., 2017). Although several review papers have elaborated the applications of ²¹⁰Pb as a tracer in lacustrine and marine environments (Appleby, 2001; Baskaran et al., 2014; Du et al., 2012; Kirchner and Ehlers, 1998; Mabit et al., 2014; Sanchez-Cabeza and Ruiz-Fernández, 2012; Smith, 2001), little attention has been paid to the potential limitations of the ²¹⁰Pb dating method in vegetated coastal sediments. Experience shows that vegetated coastal environments often prove to be more challenging than lake or marine sediments (Saderne et al., 2018).

Table 1 Common values of the main parameters for vegetated coastal sediments (seagrass, mangrove and tidal marshes): average dry bulk density (DBD), range of sedimentation rates (mean values are provided in brackets) and organic matter (OM) content, median organic carbon (C_org) contents, and decay rate of buried C_org (from aboveground biomass to refractory sediment C_org).

^a Seagrass (Fourqurean et al., 2012); mangrove (Donato et al., 2011) and tidal marsh (Craft, 2007; Hatton et al., 1983). ^b Seagrass (Duarte et al., 2013), mangroves (Breithaupt et al., 2012) and tidal marsh (Ouyang and Lee, 2014). ^c Seagrass (Koch, 2001); mangrove (Breithaupt et al., 2012); tidal marsh (Cochran et al., 1998; Ember et al., 1987). ^d Seagrass (Fourqurean et al., 2012); mangrove (Breithaupt et al., 2012); tidal marsh (Chmura et al., 2003). ^e Seagrass, mangrove and tidal marsh (Lovelock et al., 2017b).

Vegetated coastal ecosystems may act as closed systems in which sediment accumulation is mainly associated with the buildup of autochthonous organic and inorganic material (McKee, 2011). In this situation, ²¹⁰Pb is deposited primarily from atmospheric fallout at steady state, with no postdepositional mobility except for physical or biological mixing of the sediments (e.g., Alongi et al., 2004; Cochran et al., 1998; Marbà et al., 2015). In some cases, however, the process responsible for incorporating ²¹⁰Pb into the sediments might be more complex. Vegetated coastal ecosystems may receive both autochthonous and allochthonous sediments from the upstream catchment, coastal erosion, from the offshore zone during storm events (Turner et al., 2007) or in response to land use change (Mabit et al., 2014; Ruiz-Fernández and Hillaire-Marcel, 2009). Their sediments might be reworked through the action of fauna (bioturbation), tides, currents and waves as well as through boat anchoring, dredging or fishing activities (e.g., Mazarrasa et al., 2017; Sanders et al., 2014; Serrano et al., 2016a; Smoak et al., 2013). Effects associated with climate change, such as sea level rise and extreme climatic events, may also have an impact on rates of production and decomposition of organic matter (OM) and on sediment and C_org accumulation (Alongi et al., 2008; Arias-Ortiz et al., 2018; Mudd et al., 2010). In such instances, sediment redistribution processes and complex accretion dynamics may violate some of the assumptions of the ²¹⁰Pb dating models, producing erratic ²¹⁰Pb specific activity profiles that are difficult to interpret.

Sediments of vegetated coastal ecosystems are known to be heterogeneous, consisting of coarse-grained deposits or bedrock covered by fine-grained sediments that settled as vegetation established (McGlathery et al., 2012; Olff et al., 1997). The percentage of living (e.g., roots) and recently formed organic material is greatest in the upper 10 cm and may be affected by varying inputs of detrital sediment within vegetated coastal ecosystems and by its relative rate of decomposition. While tidal marsh and mangrove sediments have relatively high organic matter content (on average 25 %) (Breithaupt et al., 2012; Cochran et al., 1998), mineral deposits account for the majority (>85 %) of the accumulated substrate in seagrass sediments (Koch, 2001; Mazarrasa et al., 2015) (Table 1). ²¹⁰Pb has a strong affinity for fine sediments (Chanton et al., 1983; Cundy and Croudace, 1995; He and Walling, 1996a) and organic matter (Wan et al., 2005), and thus any changes in these parameters due to sediment redistribution processes or to natural heterogeneity may also result in unique types of ²¹⁰Pb specific activity profiles in sediment cores of vegetated coastal ecosystems, adding complexity to the determination of sediment model age and sedimentation rates.

Here, we present how the processes of mixing and changes in the sedimentation rate, erosion, grain size heterogeneity and OM decay impact the depth distribution of ²¹⁰Pb in vegetated coastal sediments. We compare sediment and C_org accumulation rates derived from disturbed and ideal ²¹⁰Pb profiles and assess the discrepancy between them. Through a set of simulations based on examples from the literature and using various ²¹⁰Pb dating models, we assess the limitations that apply to the determination of last-century C_org accumulation rates in such ecosystems. Finally, we provide guidance on complementary analyses to accompany the ²¹⁰Pb dating technique that can improve sediment and derived C_org accumulation rate estimates.

Figure 1 ²¹⁰Pb cycle and idealized ²¹⁰Pb specific activity profile in sediments. Images of vegetated coastal habitats: Tracey Saxby, Integration and Application Network, University of Maryland Center for Environmental Science (http://ian.umces.edu/imagelibrary/, last access: 6 September 2018).

Table 2 Summary of the main ²¹⁰Pb-based models for sediment dating (adapted from Mabit et al., 2014).

(1) Non-postdepositional redistribution; (2) constant ²¹⁰Pb_xs fluxes at the SWI; (3) constant MAR. C_m: ²¹⁰Pb_xs activity concentration in sediments at mass depth m.
I: total inventory of ²¹⁰Pb_xs. I_m: excess inventory accumulated below depth m; k_m: effective mixing coefficient (Dρ²); m_a: mass thickness of top sediment zone.
Φ: Flux of ²¹⁰Pb_xs onto the sediment; g: fraction of ²¹⁰Pb_xs flux distributed within a certain mass depth. F_i: additional supply of ²¹⁰Pb_xs to layer i.

1.1 ²¹⁰Pb dating models

²¹⁰Pb is a naturally occurring radionuclide with a half-life of 22.3 years. Particles that sink and accumulate in the bottom of aquatic systems scavenge the ²¹⁰Pb present in the water column due to the decay of ²²²Rn in (1) the atmosphere and ulterior dry and wet deposition and (2) in the water. This is known as “excess ²¹⁰Pb” (²¹⁰Pb_xs) and is added to the “supported ²¹⁰Pb” (²¹⁰Pb_sup), which is continuously produced by the in situ decay of ²²⁶Ra in bottom sediments. The accumulation of sediments over time ideally generates a decreasing distribution of ²¹⁰Pb specific activity as a function of depth (or cumulative mass in g cm⁻²) governed by the decay of the ²¹⁰Pb_xs, as illustrated in Fig. 1. The ²¹⁰Pb dating models are based on the interpretation of the ²¹⁰Pb_xs rate of decline in a sediment profile. Under ideal conditions, ²¹⁰Pb can accurately date sediments back to about seven half-lives, i.e., about 150 years (the “dating horizon”), whereby the measurement uncertainty becomes too large to detect any ²¹⁰Pb_xs. However, chronologies reaching back that far might be rarely achievable in vegetated coastal sediments as these contain relatively low concentrations of ²¹⁰Pb. The basis of the distribution of ²¹⁰Pb_xs in sediments can then be described as (Koide et al., 1972)

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\rho }C}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\cdot \left(D_b\mathit{\rho }{\displaystyle \frac{\partial C}{\partial Z}}\right)-{\displaystyle \frac{r\partial \mathit{\rho }C}{\partial z}}-\mathit{\lambda }\mathit{\rho }C,\end{array}$$

where ρ is sediment bulk density (g cm⁻³), C is the specific activity of ²¹⁰Pb_xs (Bq g⁻¹), z is depth below the sediment–water interface (cm), D_b is the sediment mixing rate (cm² yr⁻¹), r is the sedimentation rate (cm yr⁻¹), λ is the ²¹⁰Pb decay constant (yr⁻¹) and t is time (years). Commonly, depth (z) is represented as cumulative mass (m) to correct for compaction. Cumulative mass (g cm⁻²) results from the multiplication of z and ρ, and sedimentation rates are expressed as mass accumulation rates (MAR) in g cm⁻² yr⁻¹, which can be described as $\mathrm{MAR}=\mathit{\rho }(v+q)$, where v and q are the accretion and compaction velocities, respectively (Abril, 2003) (Eq. 2).

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial C}{\partial t}}={\displaystyle \frac{\partial }{\partial m}}\left(k_{\mathrm{m}}{\displaystyle \frac{\partial C}{\partial m}}\right)-\mathrm{MAR}{\displaystyle \frac{\partial C}{\partial m}}-\mathit{\lambda }C,\end{array}$$

where k_m is an effective mixing coefficient (g² cm⁻⁴ yr⁻¹).

The ²¹⁰Pb technique was first applied by Koide et al. (1972) to date marine sediments. Since then, a family of dating models has been used to interpret the ²¹⁰Pb_xs depth distribution in marine and freshwater sediment cores, increasing in variety and complexity and involving a large diversity of postdepositional redistribution processes (Table 2). However, there are three models that are most widely used and described here: the constant flux–constant sedimentation (CF–CS) model (Krishnaswamy et al., 1971), the constant rate of supply (CRS) model (Appleby and Oldfield, 1978) and the constant initial concentration (CIC) model (Robbins, 1978). Although these three models each have specific assumptions, they share the following: (1) the deposition of ²¹⁰Pb_xs is at steady state, (2) there is no postdepositional mobility of ²¹⁰Pb and (3) the deposition of ²¹⁰Pb_xs is ideal; i.e., new radioactive inputs are deposited above the previously existing material. In the simplest of the cases, the CF–CS model assumes constant ²¹⁰Pb_xs depositional flux (Φ) and MAR. In this case the ²¹⁰Pb specific activity at the surface (C₀: Bq kg⁻¹) is constant and decreases exponentially with cumulative mass. The depth of burial m is related to the elapsed time since burial through the rate of sedimentation (MAR) (Table 2). If there is mixing at the surface of the core, the mean MAR can be calculated from the ²¹⁰Pb_xs specific activity profile below the surface mixed layer (SML). If the specific activity of ²¹⁰Pb_xs declines in sections, showing two or more exponentially decaying segments, then a mean MAR can be derived for each segment (Goldberg et al., 1977). In this way the model is to some degree able to cope with temporal variations in mass accumulation rates.

Variations in the accumulation rate may occur in response to natural processes or anthropogenic influences. Under these conditions, the CRS or CIC models could be suitable. The CRS model assumes a constant flux of ²¹⁰Pb (Φ) to the sediments over time. The initial specific activity is variable and inversely related to MAR (higher MAR leads to lower ²¹⁰Pb_xs specific activity and vice versa). The dating is based on the comparison of ²¹⁰Pb_xs inventories (I_m; Bq m⁻²) below a given depth (integration of ²¹⁰Pb_xs specific activity as a function of the cumulative mass) with the overall ²¹⁰Pb_xs inventory in the sediment core (I) (Table 2). The accurate determination of the ²¹⁰Pb_xs inventories is of critical importance and required for the application of the CRS model (Appleby, 2001).

The CIC model could be a better choice at locations where sediment focusing is a major factor, where event deposit layers are present, or if significant hydrologic changes have occurred or there are hiatuses in the sediment record caused by erosion events (Appleby, 2008). The CIC model assumes that the initial ²¹⁰Pb_xs specific activity at the sediment–water interface is constant with time irrespective of the sedimentation rate so that the ²¹⁰Pb_xs flux covaries with MAR. This model permits an estimation of the age (t) at any depth at which ²¹⁰Pb has been measured if the initial specific activity is known (Table 2). However, the CIC model requires a monotonic decrease in ²¹⁰Pb_xs specific activity down-core for age reversals to be avoided, which is rare in most vegetated coastal sediments. In that event, the calculation of mean accumulation rates alone using the CF–CS model would be a more reasonable approach, as it might be too ambitious to calculate a detailed stepwise chronology based on an often limited number of data points decreasing monotonically.

While the CIC and CF–CS models have typically been used in the marine environment, the CRS model is the most preferred in lake sediments and it is becoming widely applied in estuarine environments and vegetated coastal ecosystems (Andersen, 2017; Breithaupt et al., 2014). Some of the reasons could be that it suffers less from problems associated with non-monotonic features in the ²¹⁰Pb record and is relatively insensitive to mixing (Appleby et al., 1983; Appleby and Oldfield, 1992; Oldfield et al., 1978). The selection and use of a particular model should be based on the nature of the ²¹⁰Pb_xs specific activity and MAR. For further details on the main aspects relevant to the application of ²¹⁰Pb dating models in lake or estuarine environments we recommend two detailed and comprehensive book chapters by Appleby (2001) and Andersen (2017). Here, we focus specifically on analysis of the ²¹⁰Pb dating of sediments in vegetated coastal ecosystems.

2 Methods

We performed a literature review of studies on sediment accumulation in vegetated coastal ecosystems in the Web of Science^™ (http://webofknowledge.com, last access: 23 August 2018) with the following keywords: mangrove sediment; salt marsh, saltmarsh or tidal marsh sediment; seagrass sediment AND (Boolean operator) ²¹⁰Pb, Pb-210 or lead-210. The search produced 86, 223 and 27 results, respectively, all of them using one or more of the three models described above, probably due to their simplicity, with some exceptions such as Klubi et al. (2017) that additionally uses the TERESA model (Table 2). From the literature review, we identified the most common sedimentary processes that result in anomalous types of ²¹⁰Pb_xs specific activity profiles with depth. These could be summarized in five main processes: mixing, increasing sedimentation, erosion, changes in sediment grain size and decay of organic matter (OM). Then, we simulated the target processes on initial undisturbed seagrass, mangrove and tidal marsh sediments to determine the potential discrepancies in MAR between ideal and anomalous ²¹⁰Pb_xs profiles and analyze the limitations of the ²¹⁰Pb dating technique to derive CAR in these ecosystems. Throughout this work, we refer to “deviation” in MAR or CAR as the difference between the value computed and the value estimated from the ideal ²¹⁰Pb_xs profile.

Table 3 Summary description of the numerical simulations conducted to test for the effects of sedimentary processes on ²¹⁰Pb_xs specific activity profiles in seagrass and mangrove–tidal marsh sediments. MAR and CAR results derived from simulated profiles were compared with MAR and CAR estimates derived from the ideal ²¹⁰Pb_xs profiles reported here; k_s is the decay rate of the refractory sediment organic matter (OM) under anoxic conditions and k_ox is that in oxic conditions; k_lb is the decay constant of the labile OM derived from seagrass and mangrove–tidal marsh ecosystems (0.01 and 0.03 yr⁻¹, respectively).

2.1 Numerical simulations

All simulations started from an ideal ²¹⁰Pb_xs profile, complying with all assumptions, that was then manipulated to reflect the potential effect of each process. The ideal ²¹⁰Pb_xs profile was modeled considering the following: (1) a constant flux of ²¹⁰Pb_xs (Φ) of 120 Bq m⁻² yr⁻¹, i.e., the average global atmospheric flux reported by Preiss et al. (1996); (2) a MAR of 0.2 g cm⁻² yr⁻¹ and dry bulk density (DBD) of 1.03 g cm⁻³ to represent seagrass sediments; and (3) a MAR of 0.3 g cm⁻² yr⁻¹ and DBD of 0.4 g cm⁻³ to represent mangrove–tidal marsh sediments based on the typical values representative of these ecosystems (Duarte et al. 2013) (Table 1). Simulated surface activity per unit area of ²¹⁰Pb_xs (I₀; in Bq m⁻²) in ideal profiles was estimated through Eq. (3). Then Eq. (4) was applied to estimate ²¹⁰Pb_xs activities per unit area along the ideal profile (Table S1 in the Supplement).

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }{I}_{\mathrm{0}}={\displaystyle \frac{\mathrm{\Phi }}{\mathit{\lambda }}}\left(\mathrm{1}-e^{-\mathit{\lambda }{m}_{\mathrm{0}}/\mathrm{MAR}}\right)\text{(4)}& {\displaystyle }& {\displaystyle }{I}_m=I_{\mathrm{0}}\cdot e^{-\mathit{\lambda }m/\mathrm{MAR}}\end{array}$$

Activities of ²¹⁰Pb_xs per unit area (I_m) were then converted to activities per unit mass, C_m in Bq kg⁻¹, by dividing I_m by the cumulative mass (m) and converting grams to kilograms at each layer. Ideal profiles were then altered to simulate the following processes or scenarios: mixing (surface and deep mixing), increasing sedimentation (by 20 %, 50 %, 100 % and 200 %), erosion (recent and past), changes in sediment grain size (coarse and heterogeneous) and OM decay (under anoxic and oxic conditions and with labile OM contribution in sediments containing 16.5 % and 65 % OM) (Table 3; refer to Appendix A for a detailed description of the methodology used to conduct each simulation).

The CF–CS and CRS dating models were applied to the simulated ²¹⁰Pb_xs profiles to determine the average MAR for the last century (Table 2). The CIC model was excluded from the simulations presented in this study because in erratic ²¹⁰Pb_xs profiles (1) the CRS model would lead to more reasonable approaches when the flux of ²¹⁰Pb_xs is constant and (2) when that is not the case (e.g., simulations of erosion or heterogeneous grain size), the determination of mean accumulation rates alone by the CF–CS model would be a more reasonable approach. The models were applied in accordance with the simulated process. For instance, MAR was determined below the surface mixed layer in mixing simulations using the CF–CS model and piecewise in those with a change in average MAR (Appendix A). However, the models were also applied considering that (1) ²¹⁰Pb_xs profiles of mixing simulations were generated by increasing MAR and vice versa, and (2) erosion was not a factor in simulated scenarios H to J. This was done to test the potential disagreement in MAR and CAR relative to the ideal profile if the incorrect process was assumed and dating models were applied. Once the dating model was established, the C_org accumulation rate (CAR) was estimated through Eq. (5) assuming average sediment C_org contents of 2.5 % in seagrass and 8 % in mangrove–tidal marsh in both ideal and simulated sediment profiles. Under ideal conditions, CARs were 50 and 240 g C_org m⁻² yr⁻¹ in seagrass and mangrove–tidal marsh sediments, respectively. While this overall model structure was used in all simulated scenarios, MAR and CAR rates under ideal conditions varied from those reported above in increasing MAR and OM decay simulations to represent real increases in accumulation, changes in OM content and associated losses of sediment mass with depth (Table 3).

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathrm{CAR}={\displaystyle \frac{\sum _{n=i}^t(\mathit{\%}{C_{\mathrm{org}}}_i\cdot m_i)}{m_t}}\cdot {\mathrm{MAR}}_t,\end{array}$$

where (% C_orgi⋅m_i) is the mass per unit area of C_org at layer i (g C_org m⁻²), m_t is the cumulative mass over the period (t) (g m⁻²) and MAR_t is the mass accumulation rate of the period of interest (n−t) (g m⁻² yr⁻¹). When CAR is examined over the last 100 years, m_t is the cumulative mass down to the ²¹⁰Pb_xs horizon (i.e., the depth at which ²¹⁰Pb_xs activity approaches zero) and MAR_t is the mean mass accumulation rate.

Figure 2 Diagnostic features for seven distinct types of ²¹⁰Pb_xs profiles in vegetated coastal sediments (based on the literature and results from simulations in this study) and recommended actions to interpret each profile type and the sedimentary processes most likely involved. Characteristics of each profile type are explained in Sect. 3.1. The solid line in type VII represents the ²¹⁰Pb_xs specific activity profile at a reference undisturbed site.

3 Results and discussion

3.1 Types of ²¹⁰Pb_xs specific activity profiles

Seven distinct types of ²¹⁰Pb_xs specific activity profiles can be identified in vegetated coastal sediments based on examples from the literature (Fig. 2). Type I is produced by constant sediment accumulation in steady-state conditions (i.e., “ideal” profiles). The other six types of ²¹⁰Pb_xs profiles summarize the most common disturbances encountered in vegetated coastal sediments that are related to the presence of mixing (physical or bioturbation), increasing MAR, erosion or alteration by intrinsic features of sediments such as heterogeneous grain size distribution and decay of OM.

Type II (dark blue profile) illustrates a moderate decrease in the slope of ²¹⁰Pb_xs specific activities in the upper part of the sediment core, which is often attributed to higher MAR (Cearreta et al., 2002; Haslett et al., 2003; Swales and Bentley, 2015), but can also be related to a mixing process (Gardner et al., 1987).

Type III, showing constant ²¹⁰Pb_xs specific activities along the upper part of the core overlaying an exponential decaying trend, is usually interpreted as the outcome of mixing as a result of bioturbation or sediment resuspension, redeposition, and reworking (Jankowska et al., 2016; Sanders et al., 2010a; Serrano et al., 2016a; Sharma et al., 1987; Smoak and Patchineelam, 1999). This profile type has also been related to rapid accumulation of homogeneous sediment or to recent increases in MAR (Walsh and Nittrouer, 2004).

Type IV profiles show a reverse ²¹⁰Pb_xs pattern at the surface and have been attributed to a variety of factors. Similar to type III, this profile type can be caused by mixing processes in vegetated coastal ecosystems (Sanders et al., 2010b; Serrano et al., 2016a; Yeager et al., 2012) or by the deposition of allochthonous older material (Johannessen and Macdonald, 2018). It could also be produced by an acceleration of the sedimentation rate, as interpreted by Greiner et al. (2013), Smoak et al. (2013) and Bellucci et al. (2007) in seagrass, mangrove and tidal marsh, respectively, or by the decay of OM, as modeled by Chen and Twilley (1999) and Mudd et al. (2009) and observed by Church et al. (1981) in tidal marsh sediments containing >30 % OM in top layers. Additionally, type IV profiles could also be explained by nonideal deposition (i.e., a fraction of the new ²¹⁰Pb_xs input onto the sediment is not retained at the surface but penetrates to deeper layers), a process reported in peatlands and in sediments with very high porosities (>90 %) at the sediment–water interface (Abril and Gharbi, 2012; Olid et al., 2016).

Type V profiles show scattered ²¹⁰Pb_xs specific activities, which might reflect the periodic occurrence of processes that can cause type III or IV profiles and are often interpreted as evidence of repetitive reworking in the overall mixed sediment column (Alongi et al., 2001; Serrano et al., 2016a; Smoak and Patchineelam, 1999). However, this profile form has also been explained by the deposition of ²¹⁰Pb_xs outpacing its decay (λ=0.03111 yr⁻¹) (Alongi et al., 2005) or by a heterogeneous grain size sediment distribution with depth (Chanton et al., 1983; Kirchner and Ehlers, 1998; Sanders et al., 2010a), which could indicate varying ²¹⁰Pb_xs fluxes due to flood events, major land use changes or changes in vegetation cover (Appleby, 2001; Marbà et al., 2015).

Figure 3 Simulated ²¹⁰Pb_xs specific activity profiles of mixing (a), increase in sedimentation rates (b) and erosion processes (c) in vegetated coastal sediments. Several dry bulk density (DBD) and mass accumulation rates (MARs) are used to represent the effects of these processes in seagrass sediments (left panels: DBD 1.03 g cm⁻³; MAR = 0.2 g cm⁻² yr⁻¹; C_org=2.5 %) and in mangroves and tidal marsh sediments (right panels: DBD: 0.4 g cm⁻³; MAR = 0.3 g cm⁻² yr⁻¹; C_org=8 %). Bar charts illustrate the deficits in ²¹⁰Pb_xs inventories and C_org stocks after erosion events. See Appendix A for a detailed description of each scenario.

Types VI and VII represent low ²¹⁰Pb_xs activities with depth, apparently showing low, negligible modern net accumulation of sediments. Such profiles are usually related to an abundance of coarse sediments or to erosion processes, as shown in tidal marsh sediments (Ravens et al., 2009) and bare sediments that were previously vegetated with seagrass in Greiner et al. (2013), Marbà et al. (2015) and Serrano et al. (2016a).

These examples identified from the literature reveal that various sedimentary processes might produce similar types of ²¹⁰Pb_xs specific activity profiles. Any particular ²¹⁰Pb_xs profile can accommodate a range of mathematical modeling approaches (see below), which leads to the development of differing chronologies and MAR estimates. Hence, the identification of the process driving accumulation and causing variation in the ²¹⁰Pb_xs record aids in the determination of the C_org accumulation rates.

Figure 4 Simulated ²¹⁰Pb_xs specific activity profiles resulting from changes in sediment composition and organic matter decay. Coarse homogeneous grain size (a) and heterogeneous grain size with depth (b); triangles and dots represent a ²¹⁰Pb_xs profile in sediments consisting of fine grains (<63 µm) or sands (>125 µm) at surface layers, respectively. (c, d) Organic matter decay from a starting level of 16.5 % and 65 %, respectively (considering different scenarios described in Appendix A), in seagrass (left panels: DBD 1.03 g cm⁻³; MAR = 0.2 g cm⁻² yr⁻¹) and mangrove–tidal marsh sediments (right panels: DBD: 0.4 g cm⁻³; MAR = 0.3 g cm⁻² yr⁻¹).

Figure 5 Ratio of average 100-year C_org accumulation rates (CAR) between simulated and ideal ²¹⁰Pb profiles produced by various sedimentary processes in seagrass (a, b) and mangrove–tidal marsh habitats (c, b). (a, c) The correct process is assumed and models are applied. (b, d) An incorrect process is assumed and models are applied accordingly. Error bars represent the result of error propagation. Uncertainties for mean MAR were derived from SE of the regression and SE of the mean using the CF–CS and CRS models, respectively. Ratios of simulated to ideal sedimentation rates (MARs) are equal to those for CAR, determined from multiplying MAR by the fraction of C_org in sediments (Eq. 5), which was considered constant between ideal and simulated profiles. In simulations of increasing sedimentation and organic matter decay, new MAR and CAR were estimated for ideal ²¹⁰Pb profiles to represent real changes in accumulation, organic matter decay and associated changes in sediment mass with depth.

3.2 Simulated sediment and C_org accumulation rates (MAR and CAR)

We ran simulations for sedimentary processes (mixing, enhanced sedimentation, erosion) and heterogeneous sediment composition with depth (grain size distribution and OM decay). Results of the modeled ²¹⁰Pb_xs profiles are summarized in Figs. 3 and 4 and Tables S1–S7. We estimated mean 100-year MAR and CAR for the simulated profiles by applying the CF–CS and CRS models, and results were compared with those from their respective ideal non-disturbed ²¹⁰Pb_xs profiles. The estimated deviations in accumulation rates from those expected under ideal conditions are shown in Fig. 5 for seagrass and mangrove–tidal marsh ecosystems. These deviations are driven by variations in MAR caused by anomalies in ²¹⁰Pb profiles as the C_org fraction $\left(\frac{\sum _{n=i}^t(\mathit{\%}{C_{\mathrm{org}}}_i\cdot m_i)}{m_t}\right)$ was considered to be the same in both ideal and disturbed sediment profiles.

3.2.1 Mixing

Simulations of surface mixing (A and B in Fig. 3a) yielded ²¹⁰Pb_xs profiles similar to types II and III (Fig. 2), while deep mixing (scenario C) led to stepwise ²¹⁰Pb_xs profile forms similar to type V. Calculated MAR and CAR deviated between 4 % and 80 % from the ideal value in seagrass sediments, while such deviation was negligible (≤7 %) in mangrove–tidal marsh sediments due to the smaller proportion (5 %–10 %) of the ²¹⁰Pb_xs profile affected by mixing (Fig. 5a and c). In both cases, higher deviations from the expected rates were associated with deep mixing and with the use of the CF–CS model since this model interprets any divergence from the ideal exponential decrease in the ²¹⁰Pb_xs activity with depth to reflect random variation. In contrast, the CRS model is based on the ²¹⁰Pb_xs inventory (I) that is unaffected by vertical mixing.

Profiles of mixing in sediments could be equally explained by an increase in the sedimentation rate in recent years. If the incorrect process was assumed and dating models were applied accordingly, mean MAR and CAR were largely overestimated in seagrass sediments, by 20 %, 30 % and 95 %, using the CF–CS model in surface (scenario A, B) and deep mixing simulations, respectively (Fig. 5b). In mangrove–tidal marsh sediments, overestimation in mean MAR and CAR was substantial (30 %) when deep mixing was considered to be caused by an increase in MAR (Fig. 5d). A process mismatch between mixing and increased sedimentation in recent years did not cause large deviations (between 2 % and 5 %) in MAR and CAR derived by the CRS model. The CRS model outputs are similar if mixing or changes in accumulation rates are present, albeit ages within the mixed layer cannot be reported if mixing occurs.

3.2.2 Increasing sedimentation rates

Simulated increases in MAR from 20 % to 200 % (scenarios D to G, Fig. 3b) resulted in similar profile forms to those simulated with surface mixing. Increases in MAR were modeled over the last 30 years, a period over which more than a 100 % increase (2-fold) was needed to produce a reversal of ²¹⁰Pb_xs specific activities with depth under this simulation (type IV profiles). The influence of a change in MAR was better captured with the CRS model. The CF–CS model, in contrast, failed to account for rapid and large increases in MAR. Deviations from the ideal value ranged from 0 % to 15 % in scenarios D and E (20 % to 50 % increase in MAR) and were up to 60 % for a 100% increase in MAR (scenario F). Calculated MAR in scenario G (200% increase in MAR) was underestimated by a 30%, as piecewise dating is not applicable in profiles with constant or reversed activities of ²¹⁰Pb_xs with depth. In such situations, additional tracers or time markers are required to estimate MAR and CAR in the layer of constant ²¹⁰Pb_xs activities (see Sect. 4.2). Deviations from the expected value ranged from 0 % to 4 % when using the CRS model (Fig. 5a and c). Results were similar for both ecosystem types. If the recent increase in MAR was interpreted as mixing, the mean MAR and CAR were underestimated by between 10 % and 30 % in both habitat types using the CF–CS model (Fig. 5b and d). In contrast, the deviation from the ideal value was ≤5 % if the CRS model was applied.

3.2.3 Erosion

We ran three simulations (H, I and J) to represent recent (H) and past erosion events (I and J) (Fig. 3c). Simulations of erosion yielded lower ²¹⁰Pb_xs specific activities than those of the ideal reference profile (type VII, Fig. 2), and ²¹⁰Pb_xs dating horizons were found at shallower depths in these simulations (Fig. 3c). Consequently, ²¹⁰Pb_xs inventories (I) in eroded profiles were lower than expected (reference ideal profile I_ref: 3900 Bq m⁻²). Inventories of simulated seagrass sediments had a deficit of 2400 Bq m⁻² (60 %), 1250 Bq m⁻² (30 %) and 600 Bq m⁻² (15 %) in erosion scenarios H, I and J, respectively, while these deficits were of 900 Bq m⁻² (22 %), 700 Bq m⁻² (19 %) and 600 Bq m⁻² (15 %) in mangrove–tidal marsh sediments. Because seagrass ecosystems have lower sedimentation rates, a greater proportion of the ²¹⁰Pb_xs inventory was found in the top 10 cm of the sediment column and thus missing because of erosion. Simulations of past erosion events, which can be identified deeper in the profile, produced breaks in the slope of the ²¹⁰Pb_xs profiles similar to those of type II (Fig. 2). Simulated erosion scenarios did not result in a large impact in the CF–CS-derived MAR and CAR estimates under the conditions of this simulation (Fig. 5). Derived MAR and CAR were underestimated by 7 % and 2 % in seagrass and mangrove habitats, respectively. The CRS model cannot be applied to eroded ²¹⁰Pb_xs profiles unless the missing inventory is known and the total (I) and depth-specific (I_m) ²¹⁰Pb_xs inventories can be corrected. Assuming erosion was not a factor, the application of the CRS model to our simulated profiles underestimated MAR and CAR by up to 25% in seagrass and by 10 % in mangrove–tidal marsh sediments (Fig. 5b and d). Therefore, we caution against the use of the CRS model in profiles that show ²¹⁰Pb_xs inventories that are lower than those expected via atmospheric ²¹⁰Pb_xs deposition or relative to nearby undisturbed sites (Fig. 3c). The magnitude of erosion is better estimated by the deficit in inventories of ²¹⁰Pb_xs, rather than by sedimentation rates. The comparison between sediment records can provide information about the degree of erosion (see Sect. 4.3). In our simulations, the C_org stocks over the last 100 years were 20 % and 5 % lower in seagrass and in mangrove–tidal marsh sediments, respectively, compared to the corresponding ideal non-eroded profile. Part of this is likely related to the fact that the concentration of C_org is not changed, which in reality may actually change since fine sediments, in which C_org is more efficiently adsorbed, are more easily eroded and OM is remineralized when exposed to oxic conditions during resuspension (Burdige, 2007; Lovelock et al., 2017a) (see simulations in Sect. 3.2.4 and 3.2.5). Consequently, losses of sediment C_org could be significantly larger, as shown in some recent studies (Macreadie et al., 2013, 2015; Marbà et al., 2015; Serrano et al., 2016a).

3.2.4 Sediment grain size distribution

Coarse sediments are often unsuitable for ²¹⁰Pb dating as they may lead to very low ²¹⁰Pb_xs specific activities. We simulated ²¹⁰Pb_xs profiles in a sediment consisting of 70 % coarse sand (scenario K, Fig. 4a). This led to diluted ²¹⁰Pb_xs specific activities with depth similar to those produced by erosion processes. In contrast to erosion simulations, coarse but homogeneous grain size distribution with depth did not have any impact on MAR and CAR estimated by the CF–CS model, since the dilution effect did not cause any anomaly in the slope of the ²¹⁰Pb_xs profile. However, the CRS model underestimated the sedimentation rate by 10 % in both habitats (Fig. 5). The reduction of ²¹⁰Pb_xs specific activities may cause the limits of detection of ²¹⁰Pb_xs (0.35 Bq kg⁻¹ in our simulations) to be reached at shallower depths than in the ideal profile. In this simulation, the limits of detection were 4 and 6 cm shallower in seagrass and in mangrove–tidal marsh sediments, respectively (Table S5a and b in the Supplement). This was conducive to the overestimation of the sediment age at bottom layers by the CRS model and underestimated mean MAR due to the omission of a higher fraction of the integrated ²¹⁰Pb_xs activity per unit area from I_m and I at depths greater than those at which the limit of detection was reached (MacKenzie et al., 2011). This effect is known as the “old date error” of the CRS model and can be corrected as described in Binford (1990) and Appleby (2001). Because we have assumed the same C_org content in ideal as in simulated profiles, CAR estimates vary similarly to MAR. However, C_org content would likely covary with grain size, and we therefore expect lower C_org content with coarser sediments (Dahl et al., 2016; Sanders et al., 2012).

Simulations of varying grain size distribution with depth (scenarios L, M and N) led to stepwise ²¹⁰Pb_xs profile forms (Fig. 4b). A sharp increase in ²¹⁰Pb_xs specific activities in surface layers can be produced by the presence of finer sediments in which ²¹⁰Pb is preferentially associated (scenario L). As a result, sedimentation rates were 2 % to 20 % lower than those estimated for the ideal profile in both habitat types using the CF–CS and the CRS models, respectively (Fig. 5). The CRS model assumes that ²¹⁰Pb_xs specific activities are inversely related to the sedimentation rate, and thus higher ²¹⁰Pb_xs activities led to lower accumulation rates.

When coarser sediments dominated at the surface layers (scenario M), the simulated profiles obtained were similar to those with mixing and accelerated sediment accumulation in recent years (types II, III and IV). The dilution of ²¹⁰Pb specific activities caused by the deposition of coarse sediments in surface layers was interpreted by the CRS model as an increase in the sedimentation rate; however, this effect was compensated for in part by the old date error. With coarser sediments at surface layers, the CF–CS model applied piecewise overestimated average MAR and CAR by only 1 % in both habitat types, while the CRS model resulted in a 5 % overestimation (Fig. 5). If changes in grain size are considered throughout the entire ²¹⁰Pb_xs profile (scenario N), the deviations in accumulation rates were up to 10 % using both models in both habitat types. Indeed, the deposition of coarse sediments may indicate exceptional increases in sedimentation in the case of storm surge deposits or pulsed sediment deliveries. In these cases, the CF–CS model could be applied if the event layer can be identified and can be subtracted to produce a corrected depth profile from which to determine the CF–CS-derived ages and mean mass accumulation rate. However, the presence of coarse sediments is often related to a reduction in the deposition of fine particles or to the transport and erosion of these in high-energy environments, leading to a variation in the ²¹⁰Pb_xs flux onto the sediment surface that is considered constant through time by the two dating models. Where heterogeneous sediment layers are present, some corrections, such as the normalization of ²¹⁰Pb_xs specific activities, are required before the application of any of the ²¹⁰Pb dating models to obtain more accurate estimates of MAR and CAR (see Sect. 4.4).

3.2.5 Organic matter decay

Two different scenarios with low and high sediment organic matter (OM) content (16.5 % and 65 %, respectively) were modeled in relation to OM decay. In both scenarios simulated MAR and CAR were overestimated relative to those derived from ideal profiles. Variation in OM decay (from a starting level of 16.5 %) only slightly affected ²¹⁰Pb_xs specific activity profiles (Fig. 4c), causing a small overestimation of MAR and CAR of between 2 % and 5 % in both habitats and by both models under any of the rates of decay considered in this simulation (0.00005, 0.0005 and 0.01−0.03 d⁻¹) (Fig. 5a and c). OM decay in very rich organic sediments (65 % OM) caused increased ²¹⁰Pb_xs specific activities at surface (scenarios R and T) and subsurface sediments in which the decay of OM is greater, leading to a reversal of the ²¹⁰Pb_xs profile in simulated scenario S. Derived CARs were 20 %–30 % higher as estimated by the CF–CS model and 10 %–20 % using the CRS model in both habitat types (Fig. 5a and c).

Mass accumulation in vegetated coastal ecosystems is the result of the balance between material accretion (detritus and sediment) from autochthonous and allochthonous sources, decomposition, and erosion (e.g., Mateo et al., 1997). Assuming there is no erosion, the estimates of MAR and CAR by means of ²¹⁰Pb are the net result of mass accumulation with time and hence integrate both the burial and decomposition of organic matter over a centennial timescale. Because mean CAR rates are based on the C_org presently available and not the amount originally deposited, their determination will be dependent on the timescale over which they are calculated.

3.2.6 General remarks

Among the various ecosystems considered here, average last 100-year MAR and CAR derived from the CF–CS and the CRS models were less vulnerable to anomalies in mangrove–tidal marsh compared to seagrass sediments. Higher sedimentation rates lead to deeper ²¹⁰Pb_xs dating horizons and thus the fraction of the ²¹⁰Pb profile affected by anomalies was lower in mangrove–tidal marsh than in seagrass sediments. Anomalies caused by deep mixing or 2- to 3-fold acceleration in sedimentation had larger effects on the CF–CS-derived accumulation rates, while alterations caused by heterogeneous grain size composition or underestimation of ²¹⁰Pb_xs inventories primarily affected the CRS model results (Fig. 5). Care must be taken in these cases since the discrepancy in mean MAR relative to the ideal value could range between 20 % and 80 %. Our simulations showed that the decay of OM results in an overestimation of the accumulation rate, which was most severe in very rich organic sediments regardless of the model used (>50 % OM). However, this effect could reasonably be ignored in most cases since vegetated coastal ecosystems rarely contain OM concentrations >25 % (Table 1), for which the deviation in computed MAR was below 10 %. Overall, simulations showed that the variability in MAR and hence CAR due to sedimentary processes and differences in sediment composition was moderately low when appropriate dating models were applied and interpreted. Deviations in the determination of MAR and CAR, generally within 20 %, confirmed that the ²¹⁰Pb dating technique is secure (Fig. 5). However, failure to account for the correct process affecting ²¹⁰Pb specific activity profiles could lead to deviations in mean MAR and CAR exceeding 20 % (Fig. 5c, d).

MAR and CAR were most overestimated, from 20 % to 95 % in simulations with low accumulation rates, when acceleration was interpreted in mixed ²¹⁰Pb_xs profiles and the CF–CS model was applied piecewise. Deep mixing confounded with an increase in MAR generated the largest overestimation of mean CAR in both habitat types. In contrast, if mixing was assumed in ²¹⁰Pb_xs profiles showing a recent increase in MAR, mean accumulation rates were underestimated by up to 30 % using the CF–CS model below the surface mixed layer. Indeed, the CRS model was less sensitive to anomalies in ²¹⁰Pb_xs profiles; however, its application requires accurate determination of the ²¹⁰Pb_xs inventory at each depth (I_m) and in the entire record (I), which can be problematic, for instance, when all samples along a sediment core have not been analyzed or when sediment erosion has occurred at the core location. When the total ²¹⁰Pb_xs inventory is underestimated, be it through erosion, poor detection limits or insufficient core length, this generates erroneous dates and underestimation of average MAR and CAR. Underestimation of accumulation rates will depend largely on the proportion of the missing fraction of the ²¹⁰Pb_xs inventory from I_m and I. In our simulations, MAR and CAR were underestimated by 10 % to 25 %. While uncertainties within 20 % might be acceptable for the determination of mean MAR and CAR over a centennial timescale, they may not allow for the determination of a detailed geochronology, historical reconstruction, or ascertained rates of change and fluxes at specific times. In that event additional tracers or geochemical, ecological and historical data need to be used to validate the ²¹⁰Pb-derived results and reduce uncertainties caused by anomalies in ²¹⁰Pb_xs profiles in vegetated coastal sediments.

4 Approaches and guidelines

Retrieving reliable CAR depends on the correct determination of MAR and the diagnosis of the intervening sedimentary process. However, similarities in simulation outcomes and variations associated with anomalies in ²¹⁰Pb_xs profiles point to the need for additional sources of evidence to discriminate between alternative processes and constrain ²¹⁰Pb-derived estimates. ¹³⁷Cs or other independent radioactive tracers can be used to corroborate ²¹⁰Pb geochronologies. However, in its absence, geochemical information combined with knowledge of events related to land use and/or environmental changes (e.g., by means of aerial photographic evidence; Swales et al., 2015) can also be used as a tool to validate ²¹⁰Pb geochronologies and interpret ²¹⁰Pb_xs profiles. In Fig. 2 we have summarized the steps to characterize ²¹⁰Pb profiles and the sedimentary processes most likely involved and suggest several techniques to complement the ²¹⁰Pb dating method to obtain reliable MAR and CAR.

Prior to analysis, researchers can have control over some factors, such as coring, sampling or sample handling, that can create artifacts in ²¹⁰Pb profiles and therefore contribute to dating error. Guidelines for core sampling for the analysis of ²¹⁰Pb and other radionuclides have been described in detail, for example in Brenner and Kenney (2013) and in the technical report IAEA-TECDOC-1360 (2003). Some knowledge of the expected sedimentation rate is useful to decide how to section a sediment core for ²¹⁰Pb measurements, as is knowledge of the length that a core must have to reach the depth of the ²¹⁰Pb_xs horizon. Low sedimentation rates (∼1–2 mm yr⁻¹) and/or coarse sediments may imply that the ²¹⁰Pb-datable part of sediment cores is limited to the very top centimeters. In such situations, fine sectioning intervals (0.5 cm) would be required. Longer cores (of about 100 cm) should be collected if high sedimentation rates are expected (several mm yr⁻¹) so that the entire ²¹⁰Pb_xs inventory is captured and the CRS model can be applied. These can be sliced at thicker intervals without compromising the temporal resolution of the ²¹⁰Pb record. If the order of magnitude of the sedimentation rate is not known a priori, it is best to choose fine sampling intervals (e.g, at 0.5 cm along the upper 20 cm, at 1 cm from 20 to 50 cm, and at 2 cm below 50 cm) to ensure sufficient resolution.

After collection, a visual description (e.g., color, sediment texture, presence of roots, organisms or layers) of the sediments and measurement of parameters such as water content, OM and grain size are relatively low-cost actions that provide information to interpret ²¹⁰Pb distribution and the pattern of accumulation. Indeed, the type of sediment (e.g., fine vs. coarse, rich in carbonates, homogeneous or with organic debris embedded) is a factor that should be considered (IAEA-TECDOC-1360, 2003). Coarse particles or coarse-grained carbonates for which ²¹⁰Pb_xs is less preferentially adsorbed (Wan et al., 2005) may hinder the detection of any ²¹⁰Pb_xs in vegetated coastal sediments. In such situations, the analysis of ²¹⁰Pb in the smaller sediment fraction (i.e., <63 µm or <125 µm) is recommended to concentrate ²¹⁰Pb and reduce the dilution effect caused by coarse fractions. However, the application of the CRS model would then be limited to those cases in which the ²¹⁰Pb_xs is contained entirely in the sieved sediment fraction. Sieving combined with ²¹⁰Pb dating has been applied to mangrove sediments from arid regions where ²¹⁰Pb_xs flux is low (Almahasheer et al., 2017) and to carbonate-rich seagrass sediments in Florida Bay (Holmes et al., 2001). Similarly, large organic material such as roots and leaves should be removed from the sediment samples prior to ²¹⁰Pb analyses as these may contribute to the dilution of the ²¹⁰Pb_xs specific activity.

The analytical methods for ²¹⁰Pb measurements can also be chosen depending upon the amount of sample available and its expected specific activity. While the indirect determination of ²¹⁰Pb through alpha spectrometry of its granddaughter ²¹⁰Po requires a small amount of sample (150–300 mg) and will provide a significantly better limit of detection (<1 Bq kg⁻¹), the direct determination of ²¹⁰Pb by gamma spectrometry can simultaneously provide data for supported ²¹⁰Pb (²²⁶Ra) and relevant radionuclides, such as ¹³⁷Cs, ²²⁸Th, ⁷Be and ⁴⁰K, to validate the ²¹⁰Pb geochronologies. For a detailed description of the analytical methods and their advantages and disadvantages, see Corbett and Walsh (2015) and Goldstein and Stirling (2003).

4.1 General validation of ²¹⁰Pb models

4.1.1 Artificial radionuclides

Independent validation of the chronology is essential to ensure a high level of confidence in the results (Smith, 2001). Varved sediments used to validate chronologies in lakes do not occur in vegetated coastal sedimentary sequences, and thus transient signals such ¹³⁷Cs or ^239+240Pu become the most commonly used option to validate ²¹⁰Pb chronologies (Lynch et al., 1989; Sanders et al., 2010a). ¹³⁷Cs and ^239+240Pu were released to the environment through the testing of high-yield thermonuclear weapons in the 1950s to 1960s and can be used as chronometers in sediments either by assuming that the peak in activity corresponds to the fallout peak in 1963 or 1965 in the Northern and Southern Hemisphere, respectively, and/or that the depth of its first detection corresponds to the onset of fallout in the mid-1950s (Ribeiro Guevara and Arribére, 2002; Stupar et al., 2014). In addition, ¹³⁷Cs can also display a peak of elevated activity in sediment cores from Europe, corresponding to the emissions caused by the Chernobyl accident in 1986, which can also help to validate ²¹⁰Pb chronologies (Callaway et al., 1996).

However, the use of ¹³⁷Cs might have some limitations in vegetated coastal sediments. Two-thirds of the ¹³⁷Cs activity released due to the tests in the atmosphere decayed after 6 decades, rendering the identification of peaks and their correspondence to the mid-1950s and 1960s depths more difficult to determine. In addition, the detection of ¹³⁷Cs is more difficult in sediment cores from habitats located in the Southern Hemisphere and near the Equator. The low ¹³⁷Cs bomb-test fallout and Chernobyl inputs in these regions (Kelley et al., 1999; Ruiz-Fernández and Hillaire-Marcel, 2009), the greater solubility of ¹³⁷Cs in seawater (K_d=10² to 10³, Bruland, 1983), and the presence of sands and carbonates, particularly in seagrass sediments (Koch, 2001), are conditions that do not favor the adsorption of ¹³⁷Cs (He and Walling, 1996a) and may lead to its mobility (Davis et al., 1984). This effect could be intensified in the intertidal zone, which is not permanently submerged due to periodic changes in the water table. High contents of OM can also affect the distribution of ¹³⁷Cs in sediments as it is preferentially accumulated in leaf litter and may be absorbed by living roots (Olid et al., 2008; Staunton et al., 2002). In addition, decomposition of the organic phase in organic-rich sediments may cause the mobility of this radionuclide (Davis et al., 1984). These factors together may compromise the use of ¹³⁷Cs to validate ²¹⁰Pb geochronologies in vegetated coastal ecosystems. In contrast, Pu isotopes (²³⁹Pu half-life = 24 100 years and ²⁴⁰Pu half-life = 6500 years), although they are also dependent on the distribution of bomb-test fallout, would appear to offer several advantages over ¹³⁷Cs in these environments, since ^239+240Pu is relatively immobile under both freshwater and saltwater conditions (Crusius and Anderson, 1995). For instance, Sanders et al. (2016) determined sedimentation rates and ^239+240Pu penetration depths to study nutrient and CAR in intertidal mangrove mudflats of Moreton Bay, Australia. Nevertheless, because of the limitations in validating older ²¹⁰Pb dates near the base of the core and the low inventories of bomb-test fallout in sediments of Southern Hemisphere latitudes, alternative tracers might need to be used.

4.1.2 Geochemical information on sediments

Besides the irregular shape of ²¹⁰Pb_xs profiles, the absence of a secondary radioactive tracer to validate ²¹⁰Pb results can make interpretation even more complicated. However, geochemical information in sediments can provide the potential for an additional temporal frame and can also help to explain sedimentary processes that could be misinterpreted (e.g., mixing, increasing MAR, higher primary productivity or reduction of sediment supply). Analyses of additional proxies (pollen, diatom, nutrient concentrations, stable isotopes or trace metal records; López-Merino et al., 2017) that are based on well-described historical events at the study sites (e.g., pollution, crops and land clearance) could be used in the absence of secondary radioactive tracers to corroborate ²¹⁰Pb-derived dates and accumulation rates. For instance, stable Pb isotopes or total Pb concentrations in sediments are related to the history of the use of leaded gasoline in the area and can be used to identify age marks corresponding to peaks in its use or changes in lead sourcing. An example can be found in seagrass sediment cores from Florida Bay, USA (Holmes et al., 2001), and in Gehrels et al. (2005) that combines a marsh elevation reconstruction with a precise chronology derived from pollen analysis, stable isotopes (²⁰⁶Pb, ²⁰⁷Pb), ²¹⁰Pb and artificial radionuclides (¹³⁷Cs, ²⁴¹Am). Additionally, profiles of trace and heavy metals and of carbon δ¹³C and nitrogen δ¹⁵N isotopic composition of OM provide information about environmental changes for which historical information may be well known, i.e., human settlement, onset of tourism industry, temporal evolution of cropland areas or histories of variation in plant communities (Garcia-Orellana et al., 2011; Mazarrasa et al., 2017; Ruiz-Fernández and Hillaire-Marcel, 2009; Serrano et al., 2016b).

4.2 Mixing or rapid sedimentation

The methods described above for the general validation of ²¹⁰Pb models can also serve to discriminate between mixing or increasing MAR in recent years. ¹³⁷Cs and ^239+240Pu can also be used as tracers of bioturbation (Crusius et al., 2004) or acceleration of MAR during the past 50 years (Appleby, 1998; Cearreta et al., 2002; Lynch et al., 1989; Sharma et al., 1987). For instance, a demonstration of acceleration versus fast mixing could be supported when it is possible to find the distinct ¹³⁷Cs or ^239+240Pu peaks in the same zone where ²¹⁰Pb_xs activities are constant (Appleby, 2001). Changes along the profiles of geochemical elements consistent with shifts in ²¹⁰Pb_xs specific activities can often be associated with changes in sedimentation or erosion processes. For instance, instantaneous depositional event layers can be identified in the sedimentary record as isolated minima of ²¹⁰Pb_xs specific activities (Jaeger and Nittrouer, 2006; Smoak et al., 2013), but also as variations in grain size composition, OM, water content or dry bulk density (Smoak et al., 2013; Walsh and Nittrouer, 2004) (Sect. B1). Changes in sediment mineralogy can be discerned trough X-ray radiographs, X-ray fluorescence and CAT scans (described below), but also through other radionuclides, like ²²⁶Ra and ⁴⁰K, the profiles of which can be measured together with those of ²¹⁰Pb through gamma spectrometry. In particular, ⁴⁰K is also part of the mineral matrix and is often used as a surrogate for the lithogenic sediment fraction (Garcia-Orellana et al., 2006; Peterson, 2009; Xu et al., 2015).

4.2.1 Geophysical analyses

Prior to core sectioning and subsampling, nondestructive geophysical analyses such as X-ray radiographs, X-ray fluorescence (XRF), CAT scans (computerized axial tomography) or magnetic susceptibility can be conducted to identify changes in the composition of sediments with depth and in MAR or provide evidence of mixing. For instance, using X-ray radiographs many features and physical sedimentary structures may be visible (Sun et al., 2017) and, if preserved, could support the interpretation of a rapid increase in sedimentation (Walsh and Nittrouer, 2004). Pulsed sediment deliveries or erosion could be identified by discontinuous physical stratification, and sediment mixing by the presence of active burrows or the absence of sedimentary stratification (Chanton et al., 1983).

4.2.2 Short-lived radionuclides (²³⁴Th, ²²⁸Th, ⁷Be)

Radionuclides such as ²³⁴Th, ⁷Be and ²²⁸Th with properties such as particle reactivity and relatively short half-lives (24.1 days, 53.3 days and 1.9 years, respectively) are suitable to quantify sedimentation processes at scales from several months (²³⁴Th and ⁷Be) to a decade (²²⁸Th). Excess ²³⁴Th and ⁷Be are sensitive indicators of mixing in the zone of constant, scattered or reversed ²¹⁰Pb_xs specific activity profiles and are the most widely used (types II, III and IV, Fig. 2) (Cochran and Masqué, 2005; Sommerfield and Nittrouer, 1999). Demonstrating the presence of any of these short-lived radionuclides can give confidence that there is little material missing from the top of the sediment record and no recent erosion, which is essential for the application of the CRS model. An example is documented by Smoak and Patchineelam (1999) for a ²¹⁰Pb_xs profile affected by bioturbation in a mangrove ecosystem in Brazil (Sect. B2).

Recent increases in MAR could be estimated in vegetated coastal ecosystems from the slope of the best-fit line of the plot of excess ²²⁸Th against cumulative mass, as Alongi et al. (2005) showed in a mangrove ecosystem in Jiulongjiang Estuary, China (Sect. B3). However, the use of ²²⁸Th to derive recent increases in sedimentation is restricted to habitats with high accumulation rates (i.e., >4 mm yr⁻¹) with the last 10 years found in the upper centimeters. Its application is also limited due to the often very low excess activities in coastal waters and the constraint that sediments must contain a significant lithogenic-to-detrital fraction. In general, the use of short-lived radionuclides might be indicated to assess the magnitude of mixing or recent erosion in vegetated coastal sediments.

Mixing, either due to bioturbation or hydrodynamic energy, is the most common process affecting vegetated coastal sediment records. Although the presence of vegetation and anoxic sediments tends to reduce the depth of sediment mixing (Duarte et al., 2013), the mixed layer can extend to depths of 10–15 cm in marine sediments (Boudreau, 1994). If surface mixing occurs, valid estimates of sedimentation rates (within 5 % variability as shown in Sect. 3.2.1) can still be obtained using the dating models described above; however, this is only possible in sediments for which ²¹⁰Pb_xs is buried below the mixed layer prior to decay; i.e., the residence time of sediments in the mixed layer must be shorter than the effective dating timescale (∼100 years) (Crusius et al., 2004). In the example from Smoak and Patchineelam (1999) (Sect. B2) in which mixing extends to a depth of 11 cm, the sedimentation rate had to be higher than 1.1 mm yr⁻¹ in order for ²¹⁰Pb to be a useful chronometer (residence time in the mixed layer is 110 mm to 1.1 mm yr⁻¹ = 100 years, which is within the effective dating timescale of ²¹⁰Pb).

4.2.3 Maximum penetration depth of ²¹⁰Pb_xs

A chronology cannot be estimated if mixing affects the whole or the vast majority of the sediment record, such as in the deep mixing simulation in seagrass sediments in this work. However, information such as the total historical inventory of elements, like nutrients accumulated at a site, and the maximum conservative sedimentation rate can still be estimated. The penetration–depth method (Goodbred and Kuehl, 1998; Jaeger et al., 2009) uses the maximum penetration depth of ²¹⁰Pb_xs (depth of disappearance) as a marker horizon for sediments that are ∼100 years old. Low surface ²¹⁰Pb_xs specific activities can greatly restrict the age of the ²¹⁰Pb dating horizon; therefore, this is an issue that should be considered when establishing the age of the ²¹⁰Pb_xs horizon. For surficial specific activities less than ∼100 Bq kg⁻¹ this could be as few as three to four ²¹⁰Pb half-lives, i.e., 65–90 years. By locating the dating horizon independently of the alteration of sedimentary processes and of assumptions of the CF–CS or CRS models, an upper estimate of the average MAR can be derived. Note that by using this method, the rates of change or fluxes cannot be estimated and these types of ²¹⁰Pb_xs profiles may be of little use in establishing chrono-stratigraphies since they are unlikely to have good records of other environmental parameters.

4.3 Erosion: ²¹⁰Pb_xs inventories (I)

Assessing the extent of erosion requires the comparison of the ²¹⁰Pb_xs inventories between reference, i.e., undisturbed locations (I_ref), and eroded sites (I). Because ²¹⁰Pb_xs is particle reactive, once deposited in sediments, its subsequent lateral redistribution is primarily controlled by resuspension and transport processes, and thus a deficit in ²¹⁰Pb_xs inventories relative to undisturbed sediments may indicate loss or mobilization of sediment particles. This approach has been used in terrestrial soils (Martz and Jong, 1991; Walling et al., 2003) and more recently to assess erosion of seagrass sediments (Greiner et al., 2013; Marbà et al., 2015; Serrano et al., 2016a) (Sect. B4). Because the ²¹⁰Pb_xs inventories at a reference undisturbed location may be spatially variable, we recommend the use of a reference inventory value based on several cores (i.e., mean ± 2 SE). The consistency of the resulting reference inventory value can then be assessed by comparing it with the ²¹⁰Pb_xs inventory measured in a nearby undisturbed soil characterized by minimal slope or with that expected from the local atmospheric flux of ²¹⁰Pb_xs (Φ). See Preiss et al. (1996) for global and regional ranges of atmospheric fluxes. The expected inventory (I_ref) can then be derived as $I_{\mathrm{ref}}=\mathrm{\Phi }/\mathit{\lambda }$, where λ is the decay constant of ²¹⁰Pb (0.03111 yr⁻¹).

4.4 Heterogeneous sediment composition

4.4.1 Normalization of ²¹⁰Pb_xs specific activity

Dating models assume rapid and nondiscriminatory removal of radionuclides from the water column regardless of major changes in grain size or OM content along a sediment record. Radionuclide adsorption onto sediments is strongly governed by the binding capacity of the settling particles (Cremers et al., 1988; Loring, 1991), and thus its scavenging is increased by fine-grained texture (He and Walling, 1996a) and OM particles (Yeager and Santschi, 2003). Variations in the influx of these particles into vegetated coastal sediments may proportionally affect the influx of particle-bound ²¹⁰Pb_xs (as long as it is still available), thus violating the assumption of constant flux of the CRS and CF–CS models and leading to subsections and irregularities in ²¹⁰Pb_xs profiles. Constant or reversed patterns in ²¹⁰Pb_xs activity profiles, which could be easily mistaken for reworked deposition, could be caused, for instance, by vertical fluctuations of grain size due to seasonal variations in sediment discharge or reoccurring tidal currents. Sediment studies often attempt to minimize these effects by normalizing radionuclide specific activities to granulometric or geochemical parameters that reduce the influence of preferential adsorption by fine sediments and OM (Álvarez-Iglesias et al., 2007; Loring, 1991; Wan et al., 2005), allowing researchers to obtain ²¹⁰Pb_xs profiles showing an exponential decreasing trend with depth (Kirchner and Ehlers, 1998; Sun et al., 2017). Radiometric applications in coastal sediments have traditionally opted for grain size normalizers such as the <4, <63 µm fraction or Al content (Álvarez-Iglesias et al., 2007; Sanders et al., 2010b; Sun et al., 2017; Walsh and Nittrouer, 2004), while in dynamic, sandy-rich coastal systems where the mud fraction is small, normalization by OM content has been shown to also be effective (Van Eaton et al., 2010). Equation (6) can be used to normalize ²¹⁰Pb_xs specific activities (²¹⁰Pb_xs-NORM in Bq kg⁻¹) by grain size fractions, OM content or other geochemical parameters that control the variation of the input of ²¹⁰Pb_xs.

$$\begin{array}{}\text{(6)}& {\displaystyle }{}^{\mathrm{210}}{\mathrm{Pb}}_{\mathrm{xs}\text{-}\mathrm{NORM}}={}^{\mathrm{210}}{\mathrm{Pb}}_{\mathrm{xs}\text{-}\mathrm{MEAS}}({\mathrm{NP}}_{\mathrm{AVG}}/{\mathrm{NP}}_m),\end{array}$$

where ²¹⁰Pb_xs-MEAS is the measured specific activity of the bulk sample at depth m, and (NP_AVG∕NP_m) is the ratio between the core average normalizing parameter and its content at depth m. For instance, multiplication by this ratio corrects measured ²¹⁰Pb activities for variations in OM with respect to an average core value.

4.4.2 ²²⁶Ra specific activity profiles

²¹⁰Pb_xs specific activity is determined by subtracting supported ²¹⁰Pb, assuming it is in equilibrium with ²²⁶Ra, from total ²¹⁰Pb specific activity. This is straightforward when gamma spectrometry is employed since the total ²¹⁰Pb and ²¹⁰Pb_sup (i.e., ²²⁶Ra) can be quantified simultaneously. On some occasions, particularly when ²¹⁰Pb is determined by alpha spectrometry, ²²⁶Ra is not measured, and supported ²¹⁰Pb is most often determined from the region of constant and low ²¹⁰Pb specific activities at depth or, alternatively, from a number of determinations of ²²⁶Ra via gamma spectrometry or liquid scintillation counting (LSC) along the core. This method assumes that ²²⁶Ra and ²¹⁰Pb_sup are constant throughout the sediment core (Binford, 1990). However, this might not be always the case, especially in heterogeneous profiles consisting of a variety of sediment types (Aalto and Nittrouer, 2012; Armentano and Woodwell, 1975; Boyd and Sommerfield, 2016) or in records containing episodes of rapid sedimentation (Chanton et al., 1983). In addition, equilibrium of ²¹⁰Pb_sup with ²²⁶Ra might be compromised in surface sediments for which ²²²Rn is deficient (Appleby, 2001). Although variations in ²²⁶Ra specific activities with depth are small in most cases, accurate determination of ²²⁶Ra might be crucial in sediments with low total ²¹⁰Pb (e.g., due to the presence of coarse sediments), for which slight variations in the ²¹⁰Pb_sup fraction may result in significant errors in the estimation of ²¹⁰Pb_xs and derived accumulation rates (Diemer et al., 2011). Therefore, we recommend measuring ²²⁶Ra specific activity profiles or, at least, using depth-specific ²²⁶Ra values at several depths along a sediment profile to estimate ²¹⁰Pb_xs.

5 Conclusions

The ²¹⁰Pb dating technique provides crucial information for the study of carbon sequestration in vegetated coastal sediments and can also provide accurate geochronologies for the reconstruction of environmental processes. However, ²¹⁰Pb-based geochronologies may be difficult to conduct in mangrove, tidal marsh and seagrass ecosystems where unaltered sedimentary records are rare.

Shallow vegetated coastal sediments are often affected by a number of processes such as mixing and bioturbation, accelerated sedimentation, or erosion and might be composed of heterogeneous sediments. These factors may lead to anomalies in the ²¹⁰Pb_xs specific activity profiles and thus produce erroneous geochronologies and biased mean last-century MAR and CAR. Discrepancies in mean MAR and CAR between irregular and ideal ²¹⁰Pb profiles simulated in this study are within 20 % if the intervening sedimentary process is well diagnosed. Otherwise, these deviations may range between 20 % and 100 %, with higher errors associated with the application of (1) the CF–CS model in records showing intense mixing or large increases in MAR and (2) the CRS model in incomplete ²¹⁰Pb_xs sediment profiles. Additional tracers or geochemical, ecological or historical data can be used to identify the process causing anomalies in ²¹⁰Pb_xs profiles and reduce uncertainties in derived accumulation rates. Using the procedures in Sect. 4, researchers have been able to obtain reliable chronologies and CAR in vegetated coastal sediments. Special caution, however, should be applied for sites at which sediments might be altered by multiple processes (leading to profile types V or VI shown in this study) and when other chronological tools or time markers are not available (e.g., ¹³⁷Cs). Sites that have slow accumulation rates and/or intense mixing may not be datable and derived CAR estimates may be largely overestimated. Mistakes would include assigning discrete ages in mixed sediments or extrapolating an age–depth model for a core that should be considered undatable to depths down the core or to nearby sites. While attention should be paid to the limitations of ²¹⁰Pb-derived results in vegetated coastal ecosystems, the guidelines provided here should help interpret complex ²¹⁰Pb profiles obtained from vegetated coastal sediments and develop a strategy to strengthen the evaluation of MAR and CAR.

Appendix A: Simulation methods

A1 Mixing

To simulate surface mixing (scenarios A and B), we estimated the accumulated ²¹⁰Pb_xs activity per unit area over the top 5 cm of the ideal ²¹⁰Pb_xs profile (I_5 cm: 2126 Bq m⁻² in seagrass and 723 Bq m⁻² in mangrove–tidal marsh sediments) (Table S1a and b in the Supplement). We split this inventory within the upper 5 cm using a random function, the outputs of which fell within the standard deviation (±SD) of the mean of the ²¹⁰Pb_xs activities in the upper 5 cm (±107 Bq m⁻² in seagrass; ±9 Bq m⁻² in mangrove–tidal marsh sediments). To simulate deep mixing (scenario C), we followed the same methodology but we randomly split the ²¹⁰Pb_xs inventory within the upper 15 cm, which is a depth reported as deep mixing in seagrass (Serrano et al., 2016a), mangroves and tidal marshes (Nittrouer et al., 1979; Smoak and Patchineelam, 1999) and is characteristic for marine sediments globally (Boudreau, 1994). We ran the simulation several times until we obtained three scenarios (A, B, C) of mixing encompassing a range of surface mixed layers (SMLs) (Table S2a and b). Mixing A (k_m: ∞ g² cm⁻⁴ yr⁻¹) consisted of constant ²¹⁰Pb_xs specific activities with depth in surface layers; mixing B (k_m: 20–23 g² cm⁻⁴ yr⁻¹) was characterized by a decrease in the slope of ²¹⁰Pb_xs in top layers; and mixing C represented deep mixing from the sediment surface down to 15 cm (k_m: 6–25 g² cm⁻⁴ yr⁻¹). ²¹⁰Pb_xs activities per unit area (A) were converted to ²¹⁰Pb_xs specific activities (C) in Bq kg⁻¹, which we averaged every two layers to represent smooth transitions. Sedimentation and derived CAR were estimated from the modeled profiles using the CF–CS and the CRS models. The CF–CS model was applied below the depth of the visually apparent SML (3 cm) in scenarios A and B to avoid overestimation of MAR. The CF–CS model was applied to the entire profile in deep mixing scenario C in seagrass sediments and below the apparent mixed layer (13 cm) in mangrove sediments. Deep mixing affected 10 % of the entire ²¹⁰Pb_xs profile of mangrove–tidal marsh sediments and 45 % of seagrass sediments. To account for the deviations in mean MAR and CAR associated with a process mismatch (i.e., as if considering the actual mixing to be caused by an increase in MAR), we applied the CF–CS model piecewise (scenarios B and C) and to the entire profile (scenario A). In the case of the CRS model, ages were determined at each layer and average centennial MAR was estimated by dividing the mass of sediment accumulated (g cm⁻²) down to 100-year depth by its age (i.e., 100 years) in all cases.

A2 Increasing sedimentation

We simulated an enhancement of MAR that could result, for instance, from increased sediment runoff due to coastal development by increasing the basal MAR (0.2 and 0.3 g cm⁻² yr⁻¹ in seagrass and mangrove–tidal marsh, respectively) by different magnitudes (20 %, 50 %, 100 % and 200 %). Increases in MAR were simulated over the top 6 cm and 23 cm of the idealized ²¹⁰Pb_xs specific activity profiles, which represent the last 30 years of accumulation in seagrass and mangrove–tidal marsh sediments, respectively. Last-century mass accumulation rates expected for ideal profiles were estimated by dividing the accumulated mass down to a 100-year depth (derived from gradual increases in MAR) by its age (Table S3a and b). ²¹⁰Pb_xs specific activities (C_m) as a result of increased MAR were estimated through Eq. (A1) for each layer. Simulations of increasing MAR generated four profiles per habitat type (scenarios D, E, F and G) (Fig. 3b). Average MAR and CAR were estimated from the modeled profiles using the CF–CS and CRS models. The CF–CS model was applied piecewise in scenarios D, E and F and below the layer of constant ²¹⁰Pb_xs in scenario G.

$$\begin{array}{}\text{(A1)}& {\displaystyle }{C}_m={\displaystyle \frac{\mathit{\lambda }\cdot I_m}{\mathrm{MAR}\cdot \mathrm{10}}},\end{array}$$

where λ is the decay constant of ²¹⁰Pb (0.0311 yr⁻¹) and I_m is the ²¹⁰Pb_xs inventory accumulated at layer m. The value 10 allows for the unit conversion to Bq kg⁻¹.

We also estimated mean MAR and CAR assuming that the process causing scenarios D, E, F and G was mixing. For this, we applied the CF–CS model below the surface mixed layer (6 and 23 cm in seagrass and mangrove–tidal marsh sediments, respectively). The shift in the slope of the ²¹⁰Pb_xs profile in scenario D in seagrass sediments was minimal, and hence we applied the CF–CS model to the entire ²¹⁰Pb_xs profile, as this would likely be the method applied by most researchers in a real case. The CRS model was run similarly if mixing or changes in accumulation rates are expected, ages were determined at each layer and average centennial MAR was estimated by dividing the mass of sediment accumulated (g cm⁻²) down to 100-year depth by its age (i.e., 100 years). If mixing is expected, ages within the mixed layer cannot be reported.

A3 Erosion

Erosion in vegetated coastal sediments can occur due to high-energy events (Short and Wyllie-Echeverria, 1996), vegetation loss and subsequent destabilization of sediments (Marbà et al., 2015) or mechanical disturbances (e.g., Serrano et al., 2016c). We ran three simulations to represent recent (H) and past erosion events (I and J) (Fig. 3c). We started with an ideal ²¹⁰Pb_xs profile with a total initial ²¹⁰Pb_xs inventory of 3900 Bq m⁻². To simulate erosion, we removed the ²¹⁰Pb_xs inventory accumulated in the top 0–5 cm (H), middle 5–10 cm (I) and 10–15 cm sections (J) in sediments from both habitat types (mangrove–tidal marsh and seagrass). The resulting ²¹⁰Pb_xs activity per unit area (Bq m⁻²) was converted to ²¹⁰Pb_xs specific activities (Bq kg⁻¹) by dividing by the corresponding mass depth (g cm⁻²) at each section after correcting the latter for the loss of sediment layers (Table S4a and b). ²¹⁰Pb specific activities were averaged every two layers to simulate smooth transitions rather than a sharp discontinuity after an erosion event. We estimated the resulting average MAR and CAR using the CF–CS model (applied piecewise in erosion scenarios I and J). The CRS model should not be applied in simulated erosion scenarios since the overall core inventories (I) are incomplete. However, we ran the CRS model to test the errors associated with its application in eroded sediments assuming that erosion is not a factor.

A4 Changes in sediment grain size

We simulated various ²¹⁰Pb_xs profiles with changes in sediment grain size distribution using the approach described by He and Walling (1996), in which the specific surface area of particles exerts a primary control on the ²¹⁰Pb_xs adsorbed:

$$\begin{array}{}\text{(A2)}& {\displaystyle }C\left(S_{\mathrm{sp}}\right)=\mathit{\mu }\cdot S_{\mathrm{sp}}^{\mathrm{0.67}},\end{array}$$

where C is ²¹⁰Pb_xs specific activity (mBq g⁻¹), S_sp is the specific surface area of the sediment particles (m² g⁻¹) and μ is a constant scaling factor depending upon the initial ²¹⁰Pb_xs activity per unit area (mBq m⁻²). The ²¹⁰Pb_xs specific activity in bulk sediments can also be represented by Eq. (A2), replacing S_sp by the mean specific surface area S_mean (m² g⁻¹) of the bulk sample. In this work, we estimated μ at each layer of an ideal ²¹⁰Pb_xs profile in seagrass and mangrove–tidal marsh sediments if S_sp throughout the core is 0.07 m² g⁻¹, corresponding to a mean particle size of 63 µm. The surface area can be estimated as (Jury and Horton, 2004)

$$\begin{array}{}\text{(A3)}& {\displaystyle }{S}_{\mathrm{sp}}={\displaystyle \frac{\mathrm{3}}{\mathit{\rho }\cdot r}},\end{array}$$

where ρ is the density of the sediment particles and r is the mean radius of sediment particles, which are considered spherical.

We estimated the weighted-mean specific surface area of a very coarse sediment composed of 70 % coarse sand (500–1000 µm), 20 % medium sand (250–500 µm) and 10 % silt (4–63 µm) (S_mean=0.0153 m² g⁻¹) through Eq. (A3) (size scale: Wentworth, 1922). The bulk density (ρ) of sediment fractions was considered: 1.03 g cm⁻³ for silt, 1.6 g cm⁻³ for medium sand and 1.8 g cm⁻³ for coarse sand. First, we simulated ²¹⁰Pb_xs profiles as a function of the specific surface area by applying Eq. (A2) to an ideal ²¹⁰Pb_xs specific activity profile (scenario K) (Table S5a and b). Second, we simulated a shift to sandy and clayey sediments in surface layers, which could result after the restoration or loss of vegetated coastal ecosystems. The percentages of sands and clay along the core were changed using a random function (from 60±20 % in surface to 15±5 % in bottom layers; scenarios L and M) (Table S6a and b). The shift was simulated at the same age depth (30 years before collection) in all scenarios and habitat types. Finally, we simulated a heterogeneous grain size distribution along the entire sediment profile intercalating sand and clay layers randomly with depth (scenario N) (Table S6a and b). The mass depth term was corrected in each case for changes in grain size, which lead to variations in DBD with depth. The bulk density (ρ) of sediment fractions was considered: 0.4 g cm⁻³ for clays and 1.6 g cm⁻³ for medium sands. In addition, the value of μ was readjusted at each sediment depth of the ideal profile to represent non-monotonic variations in cumulative dry mass. ²¹⁰Pb_xs specific activity profiles were estimated as a function of the specific surface area that was estimated at each layer according to the various proportions of clay and sand. The average MAR was estimated using the CF–CS and CRS models. The CF–CS model was applied piecewise in simulated scenarios L and M.

A5 Organic matter decay

²¹⁰Pb_xs in vegetated coastal sediments is deposited in association with mineral particles but also with organic particulates (Krishnaswamy et al., 1971; Yeager and Santschi, 2003). Once buried, sediment organic matter (OM) content usually decays with sediment depth and aging due to remineralization of labile fractions, leading to an enrichment of ²¹⁰Pb_xs specific activities. We simulated the resultant ²¹⁰Pb_xs profiles derived from this process in two sediments with different OM contents (16.5 % and 65 %). The first value (16.5 % OM) is within the usual range of tidal marsh–mangrove and in the high range for seagrass sediments (Fourqurean et al., 2012) (Table 1). The second value (65 % OM) represents an extreme scenario based on existing studies in seagrass and mangrove ecosystems (Callaway et al., 1997; Serrano et al., 2012). The simulations were run under three OM decay constants assuming the following: (1) the whole pool of OM is refractory under anoxic conditions, decaying at a rate of 0.00005 d⁻¹ in seagrass and in mangrove–tidal marsh sediments (Lovelock et al., 2017); (2) 50 % of the refractory pool is exposed to oxic conditions, decaying at a rate of 0.0005 d⁻¹ in mangrove–tidal marsh sediments; and (3) 50 % of the OM pool is labile and decaying fast, although exposed to anoxic conditions, at 0.01 and 0.03 d⁻¹ in seagrass and mangrove–tidal marsh sediments, respectively (Lovelock et al., 2017).

The ²¹⁰Pb enrichment factor (η) can be determined for a given time after deposition as

$$\begin{array}{}\text{(A4)}& {\displaystyle }\mathit{\eta }\left(t\right)={\displaystyle \frac{{\mathit{\chi }}_{\mathrm{s}}+{\mathit{\chi }}_{\mathrm{org}}\cdot e^{-k_{\mathrm{org}}\cdot t}}{{\mathit{\chi }}_{\mathrm{s}}+{\mathit{\chi }}_{\mathrm{org}}}},\end{array}$$

where χ_s is the mineral fraction of sediments, χ_org is the organic fraction of sediments at time 0, k_org is the decay constant of the OM in sediments and t is time and can be estimated as m∕MAR. As time (t) increases the exponential term tends to zero, and hence the OM stored in the sediment reaches a constant value at which it is no longer decomposed. We assume that the remineralized OM leaves the sediment as CO₂, but in fact a fraction (f) would transform to mineral matter as ${\mathit{\chi }}_{\mathrm{s}}\left(t\right)={\mathit{\chi }}_{\mathrm{s}\left(\mathrm{0}\right)}+f\cdot {\mathit{\chi }}_{\mathrm{org}\left(\mathrm{0}\right)}\cdot \left(\mathrm{1}-e^{-k_{\mathrm{org}}\cdot t}\right)$. In our simulations f=0 was assumed.

Then, the ²¹⁰Pb_xs specific activity of a sample of age t with initial specific activity C₀ is

$$\begin{array}{}\text{(A5)}& {\displaystyle }{C}_t={\displaystyle \frac{C_{\mathrm{0}}\cdot e^{-\mathit{\lambda }t}}{\mathit{\eta }\left(t\right)}},\end{array}$$

and the total mass accumulated with depth (M) above a layer of age t is

$$\begin{array}{}\text{(A6)}& {\displaystyle }M=\mathrm{MAR}\cdot {\mathit{\chi }}_S\cdot t+\mathrm{MAR}\cdot {\mathit{\chi }}_{\mathrm{org}}\cdot e^{-k_{\mathrm{org}}\cdot t}\cdot t.\end{array}$$

MAR was estimated using the CF–CS and CRS models. The CF–CS model was applied below the ²¹⁰Pb_xs reversed profile in scenario S. CAR was estimated through Eq. (5). Organic matter (%OM) in mangrove–tidal marsh sediments was transformed to %C_org using Eq. (A7) (Kauffman and Donato, 2012). In seagrass sediments we applied the relationship reported by Fourqurean et al. (2012) (Eq. A8) (Table S7a and b).

$$\begin{array}{}\text{(A7)}& {\displaystyle }& {\displaystyle }\mathit{\%}{C}_{\mathrm{org}}=\mathrm{0.415}\mathit{\%}\mathrm{OM}+\mathrm{2.89}\text{(A8)}& {\displaystyle }& {\displaystyle }\mathit{\%}{C}_{\mathrm{org}}=\mathrm{0.43}\mathit{\%}\mathrm{OM}-\mathrm{0.33}\end{array}$$

For this simulation new MAR and CAR were estimated as derived from ideal ²¹⁰Pb profiles to represent changes in organic matter content due to decay and associated losses of sediment mass with depth. This resulted in lower ideal MAR in seagrass and mangrove–tidal marsh sediments (seagrass: 0.17 and 0.07 g cm⁻² yr⁻¹; mangrove–tidal marsh: 0.25 and 0.10 g cm⁻² yr⁻¹ in OM decay simulations starting at 16.5 % and 65 % OM, respectively) (Tables 3 and S7).

Appendix B

B1 Case study of a sedimentation event

Hurricanes and cyclones can lead to the sudden delivery of large amounts of sediments and nutrients to mangroves and tidal marshes, which in turn can result in enhanced production (Castañeda-Moya et al., 2010; Lovelock et al., 2011). Smoak et al. (2013) obtained a ²¹⁰Pb_xs specific activity profile consistent with a large pulse of sediment delivered to fringing mangroves in the Everglades, Florida (Fig. B1). The specific activity of ²¹⁰Pb_xs was vastly different (several times lower) in sediments accumulated during the event. The sediment accumulation rate estimated by the CRS model for the upper part of the sediment record was 6 times that of background levels, resulting in a doubled accretion rate due to the high bulk density of the delivered sediments (Castañeda-Moya et al., 2010). C_org concentrations in the abruptly accumulated sediments were lower (5 %) than those of the sediments beneath the event layer (20 %–25 %). In fact, event deposits could consist of coarse sediments (for instance, sand and shell sediment layers deposited during storm events characteristic of offshore environments; Swindles et al., 2018), but also of fine sediments that could present lower ²¹⁰Pb_xs specific activity compared to surrounding layers (e.g., siltation events due to clearing of the catchment area; Cambridge et al., 2002; Serrano et al., 2016d). Indeed, if the initial ²¹⁰Pb_xs specific activity (C₀) is known, the CIC model could be useful to constrain dating when it is difficult to precisely define the thickness of such deposits. Otherwise the CF–CS model could be applied if the event layer is identified (e.g., using XRF, ²²⁶Ra or granulometry) and can be subtracted to produce a corrected depth profile from which to determine derived CF–CS ages and mean mass accumulation rates.

Figure B1 ²¹⁰Pb_xs, mass sedimentation rates (MAR), dry bulk density and C_org content in a mangrove sediment core at the Everglades, Florida. The gridded area represents the period 2002–2009, when Hurricane Wilma (2005) delivered a large pulse of sediment (Adapted from Smoak et al., 2013).

B2 Case study of mixing

An example of bioturbation processes is documented by Smoak and Patchineelam (1999) in which they showed a mixed ²¹⁰Pb_xs profile down to 11 cm of depth in a mangrove ecosystem in Brazil evidenced from the ²¹⁰Pb, ²³⁴Th and ⁷Be specific activity profiles (Fig. B2). The ²¹⁰Pb_xs activities decreased exponentially below the surface mixed layer, resulting in an estimated accumulation rate of 1.8 mm yr⁻¹. In the upper layers the ²¹⁰Pb_xs followed a complex pattern, with alternate relative maxima and minima, which could be representative of varying conditions of fluxes and sediment accumulation rates, the presence of coarse sediments, or physical or biological mixing. However, ⁷Be penetrated down to 4 cm of depth and excess ²³⁴Th was detected only in the surface layer. Sediments that are buried for a period of more than 6 months will have undetectable ⁷Be, and hence its presence at 4 cm of depth indicated that the activity of benthic communities had remobilized it downwards to a much greater degree than sedimentation.

Figure B2 ²¹⁰Pb_xs specific activity profile affected by bioturbation. Short-lived ⁷Be and excess ²³⁴Th specific activity profiles are indicators of mixing in the zone of constant ²¹⁰Pb_xs (0–5 cm) (Adapted from Smoak and Patchineelam, 1999).

B3 Case study of rapid sedimentation rates

Alongi et al. (2005) studied the rates of sediment accumulation at three mangrove forests spanning the intertidal zone along the south coastline of the heavily urbanized Jiulongjiang Estuary (China). Mass accumulation rates (MARs) were rapid and one of the ²¹⁰Pb_xs specific activity profiles showed scattered concentrations with depth. This could be related to either a very high MAR during the last decades or an intense mixing down-core. However, the excess ²²⁸Th specific activity profile, determined from the difference between the total ²²⁸Th and ²²⁸Ra activities in the sediment, showed a clearly decaying trend down to 15 cm (Fig. B3). The exponential decay curve fitted to the excess ²²⁸Th profile yielded an accumulation rate of 10 cm yr⁻¹, which was consistent with the ²¹⁰Pb specific activity profile. Therefore, the evidence provided by excess ²²⁸Th indicated that a very high MAR was the most plausible process responsible for the sediment record.

Figure B3 Vertical specific activity profiles of excess ²¹⁰Pb and ²²⁸Th in core 3564 from Alongi et al. (2005), produced by a rapid mass accumulation rate.

B4 Case study of erosion

Incomplete inventories of ²¹⁰Pb_xs indicative of erosion can be illustrated by the measured ²¹⁰Pb specific activity profiles in sediments from Oyster Harbor (Albany, Western Australia), some of which have been devoid of seagrass vegetation since the 1980s due to eutrophication (Marbà et al., 2015). The measured ²¹⁰Pb_xs specific activities in unvegetated sediments were relatively low, and the horizon of ²¹⁰Pb_xs was detected at a shallower sediment depth than in neighboring sediments, in which seagrass meadows persisted (Fig. B4). The inventory of ²¹⁰Pb_xs in the unvegetated sediment exhibited a deficit of 722 Bq m⁻² compared to that in the vegetated site. This deficit results not only from the lack of accumulation of ²¹⁰Pb_xs while sediments were unvegetated (30 years; atmospheric flux of 25 Bq m⁻² yr⁻¹), but also from the subsequent sediment erosion. These results, combined with C_org analyses, showed that unvegetated sediments had an average deficit in accumulated C_org stocks of 2.3 kg C_org m⁻² compared to vegetated sediments over the last ca. 100 years. This deficit has been produced due to seagrass loss in the 1980s but is equivalent to a loss of approximately 90 years of C_org accumulation.

Figure B4 Comparison of ²¹⁰Pb specific activity profiles and inventories of ²¹⁰Pb_xs and organic carbon (C_org) between a vegetated and unvegetated site. The grey area indicates supported ²¹⁰Pb specific activity (Adapted from Marbà et al., 2015).