Setting a reference for wetland carbon: the importance of accounting for hydrology, topography, and natural variability

Our study took place on the Delmarva Peninsula (Maryland, USA) within the mid-Atlantic Coastal Plain province (figure 1(a)). The Peninsula contains thousands of freshwater depressional wetlands ('Delmarva bays' [56]), which exist naturally as either forested or emergent (vegetated) systems. Forested bays are small (∼1 ha) closed-canopy pools with little to no herbaceous vegetation, while emergent wetlands are larger with diverse meadow-like vegetation in distinct zonal communities [51, 57]. Soils come from sandy or silty parent materials and are dominated by mineral components, although organic soils (i.e. histosols) develop in some natural bays [58]. Delmarva bays interact with shallow groundwater, and annual hydrologic patterns are characterized by seasonal (winter–spring) flooding followed by evapotranspiration-driven drawdown in the summer and fall [59]. Owing primarily to agriculture, more than two thirds of Delmarva bays have been at least partially destroyed, and those remaining have likely been hydrologically altered [56]. Recently, however, the Peninsula has seen net wetland gain [60], due in part to large-scale restoration efforts (e.g. the US Department of Agriculture's Conservation Reserve Program) [61]. Common techniques used to hydrologically restore Delmarva bays include ditch plugging, soil surface excavation, and installation of water control structures [45, 61].

Zoom In Zoom Out Reset image size

Our study includes data from two sampling sets to capture the range of variability in Delmarva bays. We initially chose three restored and three nearby least-disturbed ('natural') wetlands to assess spatial variability in SOC. In this set, hereafter called the 'emergent wetlands', we sampled from three vegetation zones: woody, herbaceous, and open water. After finding high variability in natural wetland SOC, we sampled an additional 19 natural wetlands (hereafter 'forested wetlands') at two sites with similar vegetation, edaphic and topographic conditions ('BC' and 'JL'; figure 1). Although emergent wetlands result from most wetland restoration/creation projects, forested wetlands are more common in the study region [52, 62] and were our focus for further baseline data.

In forested wetlands, duplicate soil cores (⩾1 m depth) were taken from the basin center and separated into five depth intervals (0–10, 10–30, 30–50, 50–75, and 75–100 cm) for analysis. To determine SOC, we first measured organic matter content in all samples by mass loss on ignition (LOI) for 16 h at 400 °C [63]. We then chose a 90 sample subset spanning the range of LOI values to measure C content using dry combustion (CHN analyzer, LECO Corp, St. Joseph, MI). We fit a linear regression between LOI and C content in the sample subset and used the fit equation (R² = 0.99) to estimate C content from LOI in all samples. Bulk density was determined using the core method [64], with linear correction for compaction [65].

Stocks were calculated by multiplying C content, bulk density, and sample length. For each core, SOC stocks were summed over all depth increments to 1 m, and wetland values were taken as the average of duplicate cores. The mean difference between duplicates was 3.6 kg C m⁻² (11.6% relative difference).

Because our initial objective was to examine soil spatial variability, our sampling approach differed in emergent wetlands. For each wetland and vegetation zone, we selected one random location and took three soil cores spaced 10 m apart (total n = 9 cores per wetland). Cores were taken to a minimum depth of 50 cm (mean = 67 cm) and separated by naturally-occurring horizons. Organic matter was determined by LOI, and C content was estimated using the linear regression fit with forested wetland samples. To compare with forested wetlands, we used a data harmonization approach: SOC density was modeled as a function of LOI using forested wetland data (sensu [66, 67]), and 1 m stocks were estimated by integrating depth functions fit with horizon values [68]. Details of soil data harmonization are in supplementary methods.

Continuous water level data was collected from 2018 to 2019 in a subset of study wetlands. Shallow wells were established in 12 forested wetlands and instrumented with pressure transducers (HOBO water level loggers; Onset Computer Corporation, Bourne, MA). Water levels were aggregated to daily averages for each wetland. Where necessary, data gaps were filled using nearby wells not included in this study.

We chose six metrics to capture different aspects of the hydrologic regime with potential importance in C cycling [70]. From 2019, which had average rainfall (1114 mm; figure 4(a)), we calculated: median, maximum, interquartile range, inundation duration (i.e. hydroperiod), and exposure date (i.e. the first day water level dropped below the surface). From 2018, which was wetter than average (1555 mm) but did not require gap filling, we calculated recession rate as the median daily water level decline during the growing season (May–October). We included only days when water levels declined in all wetlands and excluded days following rainstorms. A summary of metrics can be found in table 1(a).

We analyzed wetland relative topography using a 1 m lidar-derived digital elevation model (DEM) of the study area [71]. Wetland basins were delineated with a Stochastic Depression Analysis tool [72], which has been used to map wetlands in the region [73]. Using the delineated wetlands and DEM, we calculated six terrain metrics (table 1(b)) associated with wetland presence [44]. These included standard descriptors (depth, area, and shape index) as well as other metrics that, based on literature reports, may contain unique information about wetland hydrology (relative topographic position, relative catchment area, and basin profile coefficient as indicators of landscape position; potential runoff inputs; wetland bathymetry; table 1(b)). Due to an extreme value, relative catchment area was transformed by $- 1/x$ for analysis. Geospatial data was processed and analyzed in QGIS [74] using tools from Whitebox [72] and SAGA [75], and in R [v4.2.0; 76] using the sf [77] and terra [78] packages. An extended description of geocomputation is in supplementary methods.

Relationships between forested wetland SOC stocks and hydrologic/terrain metrics were quantified using multiple linear regression. Our analysis focused on two models: one with hydrologic metrics ('hydrology model'; n = 12 wetlands) and one with terrain metrics ('terrain model'; n = 19 wetlands). We included a site term in each model to account for unknown differences between our sampling sites. All statistical analysis was done in R using the car [79] and rdacca.hp [80] packages.

We expected multicollinearity, so we examined pairwise correlations and variance inflation factors (VIFs) to select minimum subsets of hydrologic/terrain metrics that characterized differences among wetlands [81]. Since our analysis was exploratory, we did not do further variable selection [80, 82], and we present regression coefficients descriptively while recognizing their uncertainty due to small sample sizes [83, 84]. Fitted models were validated by visually inspecting residuals (figures S1 and S2).

The relative importance of metrics in regression models was evaluated using variance partitioning methods [85]. Specifically, we determined the average contribution of each metric to the total variance explained in the model over all subsets of explanatory variables (i.e. hierarchical partitioning [86]). Results are expressed as partial R² values for each metric, which sum to the multiple R² of the model, and as % R², which is the partial R² divided by the multiple R².

Commonality analysis was used to compare the variability in SOC stocks explained by hydrologic metrics, terrain metrics, and site differences [87]. First, we re-fit both models without site and chose the three metrics with the highest partial R² in each. We then fit a third model including the selected groups of hydrologic and terrain metrics, and site (n = 12 wetlands), and partitioned variance using commonality analysis. This technique quantified each group's unique contribution to the model as well as the contributions shared among groups [87]. As in hierarchical partitioning, contributions are expressed as proportions of the model's R². Pairwise relationships among all metrics were assessed with Pearson correlation coefficients.

In forested wetlands, SOC stocks ranged from 19.7 to 51.0 kg m⁻² (mean = 32.2 kg m⁻² SD = 8.40, CV = 26.1%; figure 2). Natural emergent wetland estimates spanned a similar range (14.0–49.4 kg m⁻²), and two of the three wetlands contained soils dominated by organic materials. Restored emergent wetlands had a narrower range of estimated SOC stocks (10.9–28.1 kg m⁻²), although several point estimates fell within the ranges of both types of natural wetlands (figure 2).

Zoom In Zoom Out Reset image size

There was high spatial variability in emergent wetland SOC. Mean estimated stocks were highest in open water (36.7 kg m⁻²) for natural wetlands and in herbaceous vegetation (24.1 kg m⁻²) for restored wetlands (figure 2). Natural wetlands had higher among-zone variability (i.e. mean of CVs for each wetland; 36.3%) and among-wetland variability (i.e. mean of CVs for each zone; 40.5%) compared to restored wetlands (among-zone CV = 18.9%, among-wetland CV = 27.3%). Variability among natural wetlands was driven more by open water (CV = 46.7%) and herbaceous (CV = 64.4%) than by woody vegetation (CV = 10.6%; figure 2). For the herbaceous zone, all restored SOC stock estimates fell within the natural range.

As expected, there was multicollinearity and strong correlations among hydrologic metrics (figure 3). Since water levels in all wetlands followed a similar seasonal pattern (figure 4(b)), we assumed correlated metrics reflected the same hydrologic dynamics, and proceeded with a subset containing duration, max, and recession. The strongest correlation in this subset was between duration and max (r = 0.44).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The hydrology model explained variability in SOC stocks (R² = 0.75, R² _adj.= 0.61), but individual metrics differed in their contributions (table 2). Recession was the most important, with a higher contribution (37.7%) than the other two combined (20.8%), while duration was relatively unimportant. Both recession and max were positively associated with SOC stocks (β = 0.52 ± 0.24 and 0.53 ± 0.23, respectively); in other words, wetlands with deeper water and slower drying rates tended to have higher SOC stocks. The negligible r² value of max, however, shows its contribution to the model (% R² = 13.4) was likely due to suppression of residuals rather than by direct association with SOC stocks. Recession was the only metric that explained variability in SOC stocks individually (r² = 0.34).

Site was an important component in the hydrology model (% R² = 41.6). SOC stocks were lower in JL than in BC (β = − 1.62 ± 0.62), and this difference was not explained by hydrologic metrics. Although we did not include interaction terms in the model, pairwise relationships between hydrologic metrics and SOC stocks appeared to depend on site (figures 5(a)–(c)). Within BC, all three metrics had higher slopes and explained more variability (mean r² = 0.69) than within JL (mean r² = 0.14).

Zoom In Zoom Out Reset image size

There were strong correlations between terrain metrics and hydrologic metrics (figures 4 and S3). The strongest correlation was between depth and max (r = 0.87), which was expected given that depth is effectively an estimate of max based on relative topography. Duration and recession were most strongly correlated with RTP and profile, respectively (r = −0.70 and 0.66). Among terrain metrics, catchment and depth were correlated (r = −0.79) so we excluded the former from further analysis.

As with hydrology, individual metrics contributed unevenly to the terrain model (table 3). Profile was the most important (% R² = 58.0) and was positively associated with SOC stocks (β = 0.74 ± 0.34); that is, wetlands that were more cylindrical (i.e. had flatter bottoms) tended to have higher SOC stocks (figure 5(d)). Despite correlating with hydrologic metrics, depth and RTP had minor contributions to the terrain model (% R² = 12.5 and 9.5). Site was a minor component (% R² = 7.6) and relationships between terrain metrics and SOC stocks did not appear to consistently differ between sites (figures 5(d)–(f) and S4).

Commonality analysis revealed that hydrologic and terrain metrics both explained a common portion of variability in wetland SOC stocks (figure 6). Hydrologic metrics were represented by duration, max, and recession, and terrain metrics by depth, profile, and shape. The largest component was that shared by the hydrologic and terrain groups, which contributed 42.9% of the model's R² value. Both groups had similar unique components (hydrologic = 12.8%, terrain = 12.0%). The unique component of site was approximately zero (0.2%); nearly all its contribution was shared by terrain or (to a lesser extent) hydrologic metrics. Variability shared among all groups was a small negative value (−5.4%), which was likely due to sampling error and can be considered zero [88].

Zoom In Zoom Out Reset image size

We measured a wide range of SOC stocks among wetlands that would be considered replicates for most purposes. Wetlands are often classified based on vegetation and hydrogeomorphology, and wetland types are widely used to approximate functions in lieu of empirical data [89]. However, despite their shared vegetation type and hydrogeomorphic class (depressional), as well as their proximity to one another, the 19 forested wetlands we sampled had SOC stocks spanning values characteristic of mineral wetlands (∼20 kg C m⁻²) to those of organic wetlands (>50 kg C m⁻²; table 4). Mean stocks in our forested wetlands were 32.2 ± 8.4 kg m⁻²; in comparison, a recent analysis of 529 inland wetlands, which spanned vegetation types and ecoregions throughout the conterminous US, estimated mean SOC stocks at 19.7 ± 10.6 kg m⁻² [90]. We compiled a list of relevant studies encompassing a range of wetland types and found that high within-type variability is common (table 4). This suggests that caution should be broadly applied when assuming that wetland types, as defined by existing classification systems, reflect variability in ecosystem processes [27, 91–93].

We also found that SOC varied within spatially heterogeneous wetlands. Accounting for spatial variability is essential for accurately assessing wetland functions and can help to identify drivers of function in reference systems [94, 95]. In our natural emergent wetlands, estimated SOC stocks varied among vegetation zones, but in an inconsistent pattern among wetlands. Spatial variability was lower in restored emergent wetlands, which aligns with studies indicating that soil properties are more spatially homogeneous in restored wetlands [96, 97]. While vegetation communities are thought to influence wetland SOC storage [7, 98, 99], our results suggest that vegetation type may not be a reliable indicator of within-wetland SOC variability. We did not directly measure productivity, but we expect its variability to be dominated by vegetation type, as reported elsewhere (e.g. [100]). Instead, other factors (e.g. flooding patterns and soil properties that affect oxygen availability and decomposition [101]) likely contribute to spatial variability within these wetlands.

Our observations of high natural variability support a revised strategy for defining reference conditions (e.g. [102]). In small datasets like ours, high variability limits the ability to evaluate restored wetlands relative to references. Incorporating attributes like hydrology can help interpret reference information (as we address in the following section), but significantly increasing precision in reference conditions requires extensive data and a thorough understanding of landscape variability [103]. Constraints on time and funding often put these conditions out of reach for practitioners. With outcomes unpredictable, restoration assessment would benefit from a more rigorous treatment of uncertainty and a greater focus on the natural range of variability [104, 105]. For example, strategies like setting targets for groups of wetlands (rather than for individual wetlands) and reducing the precision of targets acknowledge ecosystem dynamism and promote heterogeneity [106, 107]. Importantly, our results also suggest that estimates of restored wetland 'C storage potential' should be accompanied by uncertainty analysis based on landscape-scale variability.

Water level recession rate had the strongest relationship with SOC stocks, while standard hydrologic descriptors, notably inundation duration (i.e. hydroperiod), explained relatively little variability among wetlands. Measures of duration and average water level are commonly used in wetland restoration to set performance standards and evaluate success [116, 117]. Although these metrics have been linked to wetland processes [118, 119], our results indicate that the most relevant aspects of the hydrologic regime may depend on the system and process of interest [54].

The importance of recession rate suggests a key role of drying processes in wetland SOC storage. Recession rate represents wetland response to negative water fluxes, mainly evapotranspiration and groundwater outflow [38]. Human management at the catchment scale has been shown to affect wetland recession rates by altering upland forest evapotranspiration [120]. Wetlands that dry faster may have soils with increased oxygen diffusion, increased temperature sensitivity of respiration, and more frequent wet–dry cycles; all of which can promote C loss through mineralization [121–123]. Recent reports have linked recession rates to other ecological and biogeochemical variables [124–126], which implies a broader importance in wetland function that should be examined in future studies.

Interestingly, site-scale differences appeared to modulate hydrology–SOC relationships. Such differences could include, for example, proximity to drainage ditches, absolute elevation, forest age, and adjacent land use [120, 127, 128]. While studies have established that large-scale factors like climate and geological setting can be important controls of SOC (e.g. [129]), these did not vary within our study area. Furthermore, since our forested wetlands lacked emergent vegetation and most had the same soil type (i.e. map unit; table S2), we assumed in situ production and soil properties were relatively uniform, but fine-scale spatial variability in these factors could affect SOC storage [40, 118, 130, 131]. For example, soil texture can vary among neighboring wetlands [132, 133], including those within the same soil map unit [134], and we did note that deeper layers at BC had finer textures than at JL (data not shown). Studies have shown that soil physical properties may interact with water level to influence wetland C processes [135, 136], so it is plausible that soil texture contributed to the site differences we observed. Similarly, though in situ production was likely a minor source of variability, terrestrial inputs are also important sources of C, and may have interactive effects with hydrology on wetland C budgets [137–139]. Regardless of its cause, spatial variability in wetland hydrology–SOC relationships may occur within limited areas and should be considered when setting expectations for restored wetlands.

Our study joins a growing list that have found complex or unclear relationships between wetland hydrology and SOC (e.g. [6, 7, 109, 140]). In addition to other sources of variability, these results suggest that hydrology is a multifaceted and indirect driver of wetland SOC [141, 142]. Hydrology influences wetland C processes through numerous mechanisms such as dissolved organic matter flows, temperature, and, predominantly, soil oxygen diffusion [139, 141, 143]. Despite the negative effect of water level on oxygen content (and thus decomposition [143]), persistent water saturation can maintain anoxia long after inundation recedes [144]. Further, inundation does not always slow decomposition [142]. Decomposition rates also depend on factors such as organic matter chemistry, physicochemical stabilization, microbial communities, and the availability of alternative electron acceptors [145–147]. Overall, though hydrology is a 'master variable' in wetland function, its breadth of biogeochemical influence makes it difficult to quantify its role in specific processes.

We found that hydrology and relative topography explained much of the same variability in wetland SOC stocks. Relative topography provides the physical basis for wetland hydrology [43, 148], and therefore may help characterize wetland ecosystem processes that respond to the hydrologic regime [149]. Accordingly, our results support evidence that cylindrical (i.e. flat-bottomed) wetlands have slower recession rates [150–152] and indicate that these wetlands may also have higher SOC stocks. Relative topography can also reflect other forms of landscape-scale variability [38, 153]; for example, the basin-forming processes that determine wetland bathymetry may also affect sedimentation and soil texture [69, 154]. Developing a mechanistic understanding of the relationships between relative topography, hydrology, and wetland processes is imperative to inform how restoration interventions can support functions like SOC storage.

Hydrology did not offer a clear advantage over relative topography in explaining SOC stock variability among our wetlands. This finding is significant in that relative topography is a far more practical tool for large-scale analyses. Terrain metrics can be rapidly computed over large areas with publicly-available remote sensing data, in contrast to water level data, which is more expensive and time-consuming to obtain [49]. It should be noted that human-made structures (e.g. ditches, tile drainage) that modify wetland hydrology may not be detected by remotely-sensed terrain metrics [128]; thus, care should be taken when scaling to broader areas.

Relative topography may even be a preferable indicator when only short-term hydrologic data is available. Although 2019 represented average water conditions, wetland hydrologic tendencies can shift over time with changes in catchment land cover, regional water withdrawals, and climate-driven precipitation patterns [119, 155, 156]. Hydrologic non-stationarity is especially likely in heavily altered areas, where contemporary conditions may uncouple from 'slow' variables like SOC that develop over decades to centuries [157, 158]. Relative topography is less responsive to nonlocal disturbance and could therefore more robustly represent reference wetland conditions at these time scales.