Uncertainty in United States coastal wetland greenhouse gas inventorying

As in the previous NGGI, we mapped a subset of palustrine wetlands categorized as coastal lands because their tidal elevation was lower than MHHWS. However, uncertainties in digital elevation model (DEM) elevations and in mapping tidal height were not previously included in the NGGI uncertainty analysis (EPA 2017). We enhanced the inventory by creating a probabilistic coastal lands map (supplemental information: section 2.1).

For wetland surface elevation data we used DEMs that were created using Light Detection and Ranging (LiDAR) and were aggregated by the National Oceanic and Atmospheric Association (NOAA) for their Sea-Level Rise Viewer (supplemental table 1). DEMs were created to Federal Emergency Management Administration accuracy standards (ASPRS 2004, Coveney 2013). DEMs have a nominal root mean square error (RMSE) of 0.185 m for low-relief areas and assume no bias (supplemental table 1). However, wetland vegetation and soil introduce system-specific bias and random error (Chassereau et al 2011) not captured by the nominal accuracy reporting. We corrected for a mean error of 0.173 m and estimate a RMSE of 0.205 m for wetland surfaces based on a weighted average of results from multiple US-based studies (supplemental table 2). We created a map of MHHWS heights using empirical Bayesian kriging to interpolate between NOAA tide gauges. We also created a corresponding uncertainty map incorporating random error in LiDAR mapping, datum transformations (Schmid et al 2013, Leon et al 2014), and distance between tide gauges. We combined the DEMs, the MHHWS map, and the associated uncertainty surfaces into a single spatial layer representing the probability of elevation being below MHHWS (figures 1, 2).

For palustrine wetlands, we treated mapped area as a random variable. For each of 111 palustrine wetland categories we extracted pixel counts by probability class for the coastal lands map intersecting the C-CAP class and represented mapped area as a normal distribution approximated from the multiple binomial distributions (supplemental information: section 2.1). The means and standard deviations for all 111 palustrine wetland stable classes and palustrine wetland change events are reported in supplemental table 2.

We calculated an estimated area from mapped area (Olofsson et al 2014, Byrd et al 2018) by combining accuracy assessment matrices (McCombs et al 2016) with area data from C-CAP (supplemental table 1) (supplemental tables 4–5). C-CAP did not assess classification accuracy for all individual land cover change events between 2006–2011. Instead there is an overall accuracy assessment for 2011 classification and one for the 2006–2011 generalized 'change' or 'no change' categories.

The accuracy assessment matrix records counts for all instances of mapped classes—what a datapoint was mapped as—and reference classes—what it actually was (supplemental tables 4–5). We converted the accuracy assessment matrix from counts to proportional areas (Olofsson et al 2014, Byrd et al 2018), and calculated the estimated proportional area for each class as the reference class' column sum in the proportional area matrix. We used estimated and mapped area at the full map scale to calculate an estimated to mapped area ratio (r). For each 2006–2011 C-CAP class, we used the appropriate r to scale mapped area by the 2011 class. We then used a second r value from the 'change' and 'no change' matrix to scale again based on change detection. Additional detail on how we calculated proportional area accuracy assessment matrices and class-specific scaling factors are available in the supplemental information: section 2.2).

We represented uncertainty in estimated to mapped area ratio by representing each mapped class in the accuracy assessment count matrix as a multiple multinomial distributions, a distribution that describes counts falling into two or more categories as a random variable (supplemental information: section 2.2).

We estimated soil carbon stock change in wetlands remaining wetlands and lands converted to wetlands as annual carbon burial rate from a literature review of lead-210 (²¹⁰Pb) dated cores (supplemental information; section 2.3.1). ²¹⁰Pb-based measurements typically integrate carbon burial over a century, compared to cesium-137 (¹³⁷Cs)- and artificial plot-based measurements, which integrate carbon burial over multi-decadal to annual time scales; therefore we assumed ²¹⁰Pb-based rates are more representative of long-term storage rates. We described soil carbon burial using a lognormal distribution because observed removals can not be negative when strictly relying on dated sediment profiles, observed values were always greater than zero, and the data show a positively skewed distribution (figure 3; table 1).

Zoom In Zoom Out Reset image size

For soil carbon stock change associated with wetland loss, we used average soil carbon density values reported by Holmquist et al (2018) to characterize the CO₂ emission rate (table 1). Holmquist et al (2018) determined that soil carbon density did not vary significantly by depth, and that the probability distribution of soil carbon density was described well by a normal distribution, truncated so that values could not be less than zero. They also determined utilizing a single average value for all wetlands was more parsimonious and precise than stock estimates based on available maps of soil carbon.

The previous NGGI (EPA 2017) made two assumptions about carbon changes during wetland conversion events that were not considered in the error propagation. First, the depth of soil lost to conversion was based on a range of values reported for aquaculture and salt production pond construction (0.5–2.5 m; IPCC 2014) but was fixed to 1 m. In the NGGI, this value was applied to wetland areas that converted to open water as indicated by C-CAP. Because wetland to open water conversion events were dominant in our accounting and the IPCC depth intervals for degradation were largely not applicable, we represented uncertainty regarding this assumption by using a uniform distribution ranging from 0.5–1.5 m (table 1) to represent a wide distribution centered on 1 m. This uncertainty reflected a consensus from our coauthor group and reflected an expert assumption rather than data, as we could not readily locate or ingest any relevant data. The NGGI also assumed that 100% of the carbon released by conversion from coastal wetlands to open water is lost to the atmosphere. However (Lovelock et al 2017) reviewed available studies and estimated 25%–50% of terrestrial carbon delivered to the marine environment was buried in ocean sediments (Baldock et al 2004, Cai 2011, Blair and Aller 2012). Therefore we represented the fraction lost back to the atmosphere as a uniform distribution ranging from 50%–75% (table 1).

We utilized biomass data from (Byrd et al 2018) to generate emissions factors for emergent wetlands. We accounted for forested wetland biomass using a synthesis of tree diameter at breast height (DBH) for mangrove and tidal freshwater forested plots, then converting DBH to above ground biomass using allometric equations cited within the data source, or originating from a similar representative study (supplemental information: section 2.3.2). We represented scrub/shrub data using a subset of the Byrd et al (2018) biomass data, plots that were dominated by the shrub Iva frutescens, and a subset of the forested biomass dataset, plots in which average tree heights were lower than 5 m. We converted biomass to organic carbon using a conversion factor of 0.441 (Byrd et al 2018). We represented above-ground biomass with lognormal distributions because the data exhibited skewed positive distributions (table 1; supplemental figure 2).

For CH₄ fluxes, we utilized a synthesis of annual CH₄ fluxes compiled by (Poffenbarger et al 2011) and further developed as part of the 2nd State of the Carbon Cycle Report (Windham-Myers and Cai in Revision) (supplemental information: section 2.3.3). Although IPCC guidance recommends separating CH₄ emissions by salinity class using an 18 ppt threshold (IPCC 2014), C-CAP's two salinity categories are not optimized for this purpose. We instead had to represent CH₄ emissions with separate estuarine and palustrine emissions factors based on a 5 ppt salinity threshold (Dobson et al 1995) (figure 4).

Zoom In Zoom Out Reset image size

We represented CH₄ fluxes using a normal distribution for estuarine wetlands because while the vast majority of sites indicated a net emissions scenario, one oligohaline site in New Jersey displayed net-uptake of CH₄ for much of the two years reported (Weston et al 2014) (figure 3). We represented palustrine CH₄ emissions using a lognormal distribution because flux values had a skewed positive distribution and there were no instances of net-uptake of CH₄ (figure 3; table 1). We estimated the GWP of CH₄ as 25 CO₂e ${{{\rm{CH}}}_{4}}^{-1}$ for consistency with the NGGI (IPCC 1997, EPA 2017) even though IPCC 5th Assessment Report recommends updated conversions (28 CO₂e ${{{\rm{CH}}}_{4}}^{-1}$ or 34 CO₂e ${{{\rm{CH}}}_{4}}^{-1}$ with feedbacks; table 1) (Pachauri et al 2014).

We propagated uncertainty using a Monte Carlo analysis (Ogle et al 2003, IPCC 2006, Metsaranta et al 2017). We calculated the inventory (equation (2)) 10 000 times, simulating the underlying data using random draws from the probability distributions for 145 random variables (supplemental information: section 2.3.4): including normal distributions for the mapped area for each of 111 possible palustrine stable and change classes (supplemental table 2) and multinomial distributions used to randomly draw accuracy assessment matrices for twenty-three 2011 C-CAP land cover classifications (supplemental table 4), and 2006–2011 change and no change categories (supplemental table 5).

We also propagated uncertainty for nine emissions factors or emission factor components (table 2). For normally distributed variables we randomly drew the same number of datapoints from literature review from the probability distribution then represented the emissions factor or component as the mean of the randomly drawn data. For uniform distributions, we randomly drew a single value. For emissions factors that were lognormally distributed we randomly redrew the underlying data as in normal distributions but represented the central tendency of using the exponentiated logmean. This choice is consistent with IPCC Wetlands Supplement guidance, however arithmetic means are often used for lognormally distributed emissions factors (Levy et al 2017). Because the goal of this paper is to quantify the effect of assumptions on the inventory, we repeated the uncertainty analysis using the arithmetic mean of lognormally distributed values (supplemental information: section 3.2; supplemental figure 4).

We performed a one-at-a-time sensitivity analysis (Metsaranta et al 2017), meaning we categorized sensitivity of the US scale emissions and removals to assumptions, datasets, and mapping accuracies by manipulating one input at a time and recording the effect. For each random variable we re-calculated the coastal wetland total GHG emissions and removals using the 0.025 quantile and 0.975 quantile values from Monte Carlo analysis, while fixing all others at their median value. We reported sensitivity of the inventory to each input as the difference in the total flux between using the input's minimal and maximal settings.

The sensitivity analysis also helped test the effect of some of the fundamental assumptions. For example, CH₄ fluxes need to be converted to CO₂e, and there is controversy about whether to use the GWP (25 CO₂e ${{{\rm{CH}}}_{4}}^{-1}$) (IPCC 2014) or the Sustained Global Warming Potential (SGWP; 45 CO₂e ${{{\rm{CH}}}_{4}}^{-1}$) and Sustained Global Cooling Potentials (SGCP; 203 CO₂e g ${{{\rm{CH}}}_{4}}^{-1}$) which more effectively represent the system (Neubauer and Megonigal 2015). We quantified the effect of that choice by calculating the inventory using a GWP and median values for all other inputs and then recalculated changing only the GWP to SGW/CP (Neubauer and Megonigal 2015). Also, we tested the assumption of relying on the coastal lands definition for determining how much palustrine wetland area to include in the inventory compared to a tidal wetlands definition from the National Wetlands Inventory (NWI) (Hinson et al 2017, Holmquist et al 2018, Najjar et al 2018). For this alternative analysis, we included all C-CAP palustrine wetlands intersecting an NWI-based tidal wetlands map (Holmquist et al 2018) and treated all palustrine mapped areas as fixed. In the sensitivity analysis we calculated the difference in total inventory between the default settings and the NWI based mapping strategy. We also we repeated the sensitivity analysis using the arithmetic mean of lognormally distributed values, and discuss the results further in the supplemental information (section 3.2; supplemental figure 5).

The uncertainty and sensitivity analysis presented herein suggest that uncertainty could be reduced at the scale of the contiguous US primarily by improving data availability and process-based models for CH₄ emissions, CH₄ radiative forcing, and carbon burial rates. Net-wetland CH₄ emission combines CH₄ production by methanogenic archaea under anaerobic conditions, CH₄ oxidation and consumption by methanotrophic bacteria mainly under aerobic conditions, and CH₄ transport to the atmosphere (Conrad 1989, Whalen 2005). Major controls of these processes include: water table position; soil temperature; sulfate supply and potential production of hydrogen sulfide, a methanogen toxin, for which salinity is a proxy for; vegetation, including both biomass and species composition, which may facilitate CH₄ transport from soil production sites to the atmosphere; and primary production of vegetation, since new photosynthate may be a substrate for methanogenesis (Wang et al 1996, Walter and Heimann 2000). Large discrepancies have also been noted between chamber and eddy covariance measurements of CH₄ fluxes (Hendriks et al 2010, Krauss et al 2016), suggesting the need for additional comparisons between these two methods.

The use of GWPs serves an important policy need because GWPs are transparent and tractable. However, GWPs are an oversimplification because modeling CO₂e in power units (W m⁻²) that relate directly to radiative forcing is several steps removed from actual climate impacts such as changes in temperature, precipitation, and sea level. The SGW/CP model is equally transparent and tractable, but more closely represents reality by acknowledging that changes in GHG emissions persist over several years (Neubauer and Megonigal 2015). Therefore we recommend that SCW/CP's should be considered for adoption by the IPCC. When considering the consequences of GHG inventory data beyond the IPCC context, ecosystem scientists and policy analysts should discuss metrics that are independent of time frames, such as switchover time, as they are more informative of the long-term impacts (Frolking and Roulet 2006). Our uncertainty analysis is focused on variables that are inputs to GWP and SGW/CP models, but there is an ongoing need to address the uncertainty introduced by using these models to underpin climate policy.

Currently, IPCC guidance recommends applying separate carbon burial rates to different wetland types and ecoregions to increase accuracy. However, multiple studies suggest other relevant geographic and methodological factors need to be considered in the US inventory. In some locations, accelerating sea-level rise is expanding the area conducive to carbon burial, potentially increasing carbon burial rates (Kirwan and Mudd 2012, Hill and Anisfeld 2015). A sensitivity analysis of the marsh equilibrium model highlighted relative sea-level rise, plant productivity and relative tidal elevation as dominant drivers of carbon sequestration in stable wetlands (Morris and Callaway in Press). Elevation/inundation gradients were correlated with sediment accretion dynamics in San Francisco Bay (Callaway et al 2012). Finally, there are many ways to measure carbon burial that integrate different time scales: decades- ¹³⁷Cs, centuries- ²¹⁰Pb, or millennia- ¹⁴C (Turetsky et al 2004). We recommend that future studies rectify the complex interactions between regional variability in relative sea-level rise, plant productivity, local elevation/inundation dynamics, and the potential effects of measuring this carbon burial using differing methods.

Some inputs in the inventory could be improved by targeted studies and additional data collection, including soil depth affected by conversion to open water, and percent carbon returned to the atmosphere upon loss.

Available data on elevation loss due to the diking of wetlands for agriculture (Drexler et al 2009), and the mass lost to 50 cm depth following vegetation die off (Lane et al 2016), are not suitable proxies for the vast majority of losses occurring from 2006 to 2011, estuarine emergent to open water conversions resulting from hurricane impacts and erosion in the Gulf of Mexico (Couvillion et al 2011 ). Although average carbon mass at depth in wetland soils is well constrained for coastal wetlands (Holmquist et al 2018, Sanderman et al 2018), the sensitivity of this carbon stock to different disturbances across regions, relative elevations, and time is not well known.

Uncertainty in assumptions about carbon loss is not unique to this study and was discussed explicitly in a recent global analysis of soil and biomass loss from mangrove conversions (Sanderman et al 2018), which report that the rate and forms of carbon loss may depend on soil type and depth (Donato et al 2011). Because assumptions about loss events vary from study to study, and because of the fact that these assumptions are dominant contributors to uncertainty (figure 7), future research should prioritize empirical and modeling studies that constrain depth and percent carbon loss due to wetland conversion events.

The wetlands supplement of the IPCC report provides two CH₄ emissions factors for wetlands, one for fresh to brackish conditions and another for higher salinity (18 ppt threshold) (Poffenbarger et al 2011, Bridgham et al 2013). However, C-CAP salinity categories do not match these categories, instead mapping estuarine and palustrine (5 ppt threshold; figure 3). This inconsistency limits our ability to confidently assess the true GHG balance for saline wetlands at the national scale. We propose developing maps and data to support at least three categories of salinity—saline (>18 ppt), brackish (0.5–18 ppt), and fresh (<0.5 ppt)—in order to reduce uncertainty in landscape scale CH₄ emissions from coastal wetlands (figure 4).

Existing remote sensing approaches for vegetation and inundation dynamics could improve mapping both CH₄ emissions and carbon burial rates. Recent strides in mapping coastal wetland vegetation biomass (Byrd et al 2018), vegetation species classification (Immitzer et al 2016) and seasonal dynamics (Mo et al 2015) could provide more detailed vegetation descriptions that would be a proxy for salinity zones. For inundation/elevation regimes, extensive coastal DEMs are available, but lack the accuracy to adequately map tidal flooding depth and inundation time at relevant scales and could be improved by integrating additional remote sensing and modeling (Hladik et al 2013, Parrish et al 2014, Buffington et al 2016). Future studies should quantify the precision needed for DEMs in the tidal zone. Currently soil emissions factors are calculated using tabular data, however improvements in mapping should be leveraged to support spatially-explicit approaches in future versions of the inventory incorporating trends in productivity and seasonality (Knox et al 2017), variation in carbon mineralization rates (Mueller et al 2018), edaphic factors and geomorphology (Rovai et al 2018). Many improvements may be forward-looking and hindcasting may not be appropriate (Byrd et al 2018), and spatially-explicit approaches should only be utilized only if they actually do improve precision and accuracy of inventorying compared to simpler approaches (Holmquist et al 2018).

Biomass changes were not a top contributor to uncertainty, but changes in forested and scrub/shrub biomass were the ninth and fifteenth contributors to uncertainty respectively. This study quantified the effect of uncertainty by upscaling means and uncertainties from multiple field studies, however remote sensing approaches using LiDAR, RADAR, object based image detection, and optical remote sensing, can all be used to characterize biomass changes on local to regional scales (Byrd et al 2018). Future studies could expand the uncertainty and sensitivity analysis to capture the effect that uncertainties in genus-specific assessments of wood density (Jenkins et al 2003), biomass carbon content (Byrd et al 2018), and the contributions and decay rates of downed wood (Krauss et al 2018).

C-CAP's accuracy was not a dominant contributor to the overall uncertainty in the inventory, but we were only able to quantify this from 2006– 2011. C-CAP is available for the entire CONUS coastal zone from 1996–2011, and trends were extrapolated out back to 1990 and forward to 2015 for the NGGI inventory. Future studies are needed to assess accuracy for earlier time steps.