Millennial-Scale Carbon Storage in Natural Pine Forests of the North Carolina Lower Coastal Plain: Effects of Artificial Drainage in a Time of Rapid Sea Level Rise

Abstract

Coastal forested wetlands provide important ecosystem services along the southeastern region of the United States, but are threatened by anthropogenic and natural disturbances. Here, we examined the species composition, mortality, aboveground biomass, and carbon content of vegetation and soils in natural pine forests of the lower coastal plain in eastern North Carolina, USA. We compared a forest clearly in decline (termed “ghost forest”) adjacent to a roadside canal that had been installed as drainage for a road next to an adjacent forest subject to “natural” hydrology, unaltered by human modification (termed “healthy forest”). We also assessed how soil organic carbon (SOC) accumulation changed over time using ¹⁴C radiocarbon dating of wood sampled at different depths within the peat profile. Our results showed that the ghost forest had a higher tree density at 687 trees ha⁻¹, and was dominated by swamp bays (Persea palustric), compared to the healthy forest, which had 265 trees ha⁻¹ dominated by pond pine (Pinus serotina Michx). Overstory tree mortality of the ghost forest was nearly ten times greater than the healthy forest (p < 0.05), which actually contributed to higher total aboveground biomass (55.9 ± 12.6 Mg C ha⁻¹ vs. 27.9 ± 8.7 Mg ha⁻¹ in healthy forest), as the dead standing tree biomass (snags) added to that of an encroaching woody shrub layer during ecosystem transition. Therefore, the total aboveground C content of the ghost forest, 33.98 ± 14.8 Mg C ha⁻¹, was higher than the healthy forest, 24.7 ± 5.2 Mg C ha⁻¹ (p < 0.05). The total SOC stock down to a 2.3 m depth in the ghost forest was 824.1 ± 46.2 Mg C ha⁻¹, while that of the healthy forest was 749.0 ± 170.5 Mg C ha⁻¹ (p > 0.05). Carbon dating of organic sediments indicated that, as the sample age approaches modern times (surface layer year 2015), the organic soil accumulation rate (1.11 to 1.13 mm year⁻¹) is unable to keep pace with the estimated rate of recent sea level rise (2.1 to 2.4 mm year⁻¹), suggesting a causative relationship with the ecosystem transition occurring at the site. Increasing hydrologic stress over recent decades appears to have been a major driver of ecosystem transition, that is, ghost forest formation and woody shrub encroachment, as indicated by the far higher overstory tree mortality adjacent to the drainage ditch, which allows the inland propagation of hydrologic/salinity forcing due to SLR and extreme storms. Our study documents C accumulation in a coastal wetland over the past two millennia, which is now threatened due to the recent increase in the rate of SLR exceeding the natural peat accumulation rate, causing an ecosystem transition with unknown consequences for the stored C; however, much of it will eventually be returned to the atmosphere. More studies are needed to determine the causes and consequences of coastal ecosystem transition to inform the modeling of future coastal wetland responses to environmental change and the estimation of regional terrestrial C stocks and flux.

1. Introduction

Wetlands occupy 4–6% of Earth’s land area and store C in the range of 202–535 Gt [1,2]. The importance of wetlands in providing habitats for diverse plants and animals, and a range of other ecosystem services, is well-recognized, but a sufficient understanding of the aboveground and belowground carbon cycling of wetlands is not well-understood [3,4]. Wetlands are a significant type of ecosystem of the southeastern USA coastal plain. However, natural and anthropogenic disturbances threaten the ecological integrity and C storage function of coastal wetland ecosystems.

Climate change threatens wetland ecosystems by exacerbating environmental stressors, which could have effects on the biodiversity, wetland productivity, and resilience to other stressors [5,6,7,8]. How climate change affects such wetland ecosystem functions has been given relatively little attention [9,10]. Climate change is expected to accelerate sea level rise (SLR) to between 0.4 and 1.2 m by 2100 [6], and alter the frequency and intensity of coastal storms [11]. Recently, it was recognized that the formation of “ghost forests”, created by the submergence of low-lying land, is one of the most obvious indicators of climate change [12]. SLR also increases saltwater intrusion into freshwater ecosystems [13,14,15] and alters the distribution and quantity of carbon within coastal wetlands, causing species shifts and landward migration, contributing to the direct loss of wetland area [10]. Shifts in species composition in response to SLR may cause lower marsh species to replace upper marsh species as wetlands move landward and as the continental margin adjusts [16].

In addition to these climatic impacts, the draining of coastal wetlands across the southeastern USA persisted well into the 1970s and has only slowed in the past 30 years [12]. During that time, the US government decided to drain the wetlands to decrease the groundwater level for flood control and more viable crop and timber production, construction of access roads, and other land uses [17]. Eastern North Carolina was extensively ditched over many years, and one area formerly managed for timber is now the Alligator River National Wildlife Refuge (ARNWR) in Dare County, administered by the US Fish and Wildlife Service [18]. Typical lower coastal plain road construction consisted of digging a ditch connected to a major tributary, sound or other waterway to provide drainage, and using the resulting “borrow” material as the bed for the adjacent road. Such roads run throughout ARNWR and across the southeastern lower coastal plain [18].

The effect of ditching on groundwater hydrology and drainage for a given site varies depending on the location relative to the ditch. Visual observations, whether from the ground or remotely-sensed (Figure 1), indicate that in eastern NC, vegetation on the side of the road with a ditch is often dead or dying, while on the other side the vegetation is much healthier/abundant/productive. Ditches were implicated in allowing saltwater to move inland into ecosystems that are not salt-tolerant [12], yet disentangling the hydrologic effects from the salinity effects of coastal forest decline is still a developing science.

Wetland productivity is mediated by complex bio-geomorphic feedbacks, whereby plants accumulate organic matter and trap inorganic sediments to maintain their elevation in relation to local groundwater or adjacent water bodies [19,20]. The soil accretion capacity of coastal wetlands is a key factor in the overall resilience to environmental change [20,21,22]. Quantifying primary biomass production and soil organic carbon (SOC) storage in wetland forests is therefore critical for predicting the fate of coastal wetlands and the associated ecosystem services as the continental margin adjusts to SLR.

Wetlands contain 20–30% of the estimated 1500 Pg of global soil carbon [4] making them a significant global storehouse for C. The anoxic condition of wetland soils slows down the decomposition process, resulting in the accumulation of organic matter. Thus, wetlands can store a large amount of carbon, making them an important sink of atmospheric carbon. Wetland C stocks vary in response to factors, including salinity, regional climate, water chemistry, soil type, lateral and vertical variability in soil properties, and vegetation, which complicates quantification of process rates [23,24]. These complexities limit our understanding of the quantity and distribution of stored carbon in wetlands [2,4]. Therefore, a more accurate approach to wetland soil carbon accounting must be used to develop models relating the C dynamics to temporal trends in wetlands. Techniques such as radiocarbon dating markers to determine the geochronology of a series of plant macrofossils are useful tools to estimate carbon accumulation rates, and were used to integrate the ages of fossils from 2000–3000 years ago [25]. The radiocarbon dating of plant material found at depth in peat profiles near ARNWR was also used to recreate historical patterns of past vegetation dynamics/marsh accretion as a means of sea level reconstruction over the recent geologic past [19,25].

In this study, we quantified the vegetation structure and C storage of two adjacent natural wetland pine forests as affected by millennial-scale SLR and century-scale artificial drainage at ARNWR, in eastern North Carolina, USA. We estimated the biomass of vegetation at different levels of the overstory canopy and evaluated tree mortality and species composition. We also estimated the total carbon stock in vegetation and soils and how these changed as a function of distance from local water bodies/open water. Further, we assessed the rate of SOC accumulation using ¹⁴C radiocarbon dating on wood samples collected at different depths within the peat profile down to the underlying Pleistocene mineral sedimentary deposits, spanning approximately 1800 years BP, and related findings to local sea level reconstructions. We hypothesized that the ghost forest, which has a visibly high mortality (e.g., Figure 1), would have lower total C stocks (in vegetation and soils) compared to adjacent healthy forests. We also hypothesized that ecosystem C storage would be positively correlated with the distance from open water, due to protection from salinity and/or hydrologic stress. Finally, we hypothesized that the recent decades-long visible increase in ghost forests would be associated with recent decadal-scale change in the rate of local SLR as revealed by the ¹⁴C dating estimation of marsh accretion over the past two millennia.

2. Methods

2.1. Study Site

This study was conducted in the northern reaches of the Alligator River National Wildlife Refuge (ARNWR) in Dare County, North Carolina (Figure 2), which spans over 61,600 hectares [18]. The wetlands physiography of ARNWR is a product of post-Wisconsinian glaciation, as coastal rivers cut channels through the landscape at a time of lower sea level that consequently filled with marine sediments overlain by accumulating organic peat [26,27]. The soil at this site is classified as a Belhaven series histosol (loamy, mixed, dystic, thermic Terric Haplosaprists) with a highly decomposed organic matter layer underlain by loamy marine sediments [28]. This part of the refuge, adjacent to the water bodies of Albemarle Sound and East Lake, was historically dominated by extensive pond pine forest due to this species’ tolerance of hydric soil conditions and dependence on fire for regeneration [18,29]. An abandoned fire road with adjacent drainage ditch, built in 1965 (Ed Sawyer Road) (Figure 1 and Figure 2), provided an ideal location to observe the effects of altered hydrology and the distance from surrounding water bodies on pond pine forest carbon storage and ecosystem transition. Common overstory species at the study site include swamp bay (Persea palustris), pond pine (Pinus serotina), red maple (Acer rubrum) and sweetgum (Liquidambar styraciflua), in addition to numerous understory shrubs and herbaceous plants. The elevation is approximately 30 cm above mean sea level, average annual precipitation is 1163 ± 49 mm (1981–2018) and average annual temperature is 15.7 °C ± 0.4 °C (2009–2018).

2.2. Overstory Tree Biomass and C Content

Vegetation sampling occurred throughout 2014 and 2015. Paired, circular (15 m radius) vegetation sampling plots were installed in healthy forest and ghost forest along both sides of the fire road, starting near the bank of East Lake and progressing inland (Figure 2). Overstory tree species were identified and live/dead status was determined. All tree heights (living and dead) were measured with a Nikon laser hypsometer, and diameters were measured at breast height (1.3 m above ground level) with a D-tape. Dead trees were assigned to a decay class using USDA Forest Service Forest Inventory and Analysis (FIA) protocols [30] based on visible characteristics (

Table 1. Allometric equations and decay classes used to estimate live biomass and C content of overstory tree species at Alligator River National Wildlife Refuge in eastern, North Carolina, USA [30,32].

Live Biomass Regressions
Species	Allometric Equation	Notes
Pond pine (P. serotina)	y = exp(2.0171 + 2.3373 × ln(dbh))	Equation for P. rigida
Red maple (A. rubrum)	y = 2.52363 × (dbh2)1.9648
Sweetgum (L.sty raciflua)	y = 1.822108 × (dbh2)1.2635
Bay (Persea spp.)	y = −13.388 + 6.82 × (dbh2)
FIA decay class	Qualification	C remaining (%)
0	Alive	100
1	Dead, small branches attached, knife does not penetrate	97
2	Dead, no small branches, has large branches, little penetration	97
3	Dead, no large branches, height intact, more penetration	86
4	Dead, top broken, easy penetration, highly decayed	53

). Overstory standing tree biomass was estimated using species-specific allometric relationships [31] for all overstory species (Table 1). Allometric regressions for estimating biomass of pond pine do not exist in the literature, representing a pressing research need, so equations were used for the taxonomically closely-related pitch pine (Pinus resinosus L.) [32]. Live tree C content was calculated as 50% of biomass [32], and that of dead trees as the calculated live C content adjusted by the % remaining by decay class.

2.3. Mid-Story and Understory Biomass and C Content

Mid-story biomass was sampled using 5 m radius circular plots nested within each overstory plot. Diameter of woody plants up to 10 cm diameter was measured at 10 cm above the ground. On-site, species-specific, allometric equations were developed to estimate mid-story, woody plant biomass. Fifty mid-story plants ranging in size and species were destructively harvested, and heights, diameters and total fresh weight were measured directly. Wood samples were cut from each plant to determine fresh-to-dry weight ratios, applied to the green biomass estimates. For small mid-story plants, biomass was estimated using the exponential equation y = 0.0944exp^(0.7102x) (R² = 0.89), where x is the diameter-at-breast-height (dbh). For large mid-story plants, the linear equation y = 1.3229x − 2.3518 (R² = 0.78) was used. Understory, primarily herbaceous, plant biomass was estimated by destructively harvesting 1 × 1 m clip plots in the center of the mid-story plots, dried to constant mass at 65 °C, and weighed. Carbon content of all mid-story and understory plants was assumed to be 50% of live biomass.

2.4. Soil Sampling

Soil samples were taken with a McCauley auger at 10 cm increments to a depth of 1 m at one location in each overstory plot. Samples were sieved (<2 mm) and dried to constant mass at 65 °C, and weighed. Bulk density was determined by dividing mass by volume for each depth increment. Each soil increment was analyzed separately for carbon concentration, content and stable C isotopes using a Picarro G2201-isotopic CO₂/CH₄ analyzer with combustion module at the Tree Physiology and Ecosystem Science Laboratory, North Carolina State University, USA. Percent soil carbon and bulk density were then used to estimate the total soil carbon content. The soil bulk density and carbon content were determined using the following formula:

$$Soil bulk density \left({\mathrm{g} \mathrm{cm}}^{-3}\right)=\frac{Ovendry sample mass \left(\mathrm{g}\right)}{Sample volume \left({\mathrm{cm}}^3\right)}$$

$$Soil carbon content\left({\mathrm{g} \mathrm{C} \mathrm{cm}}^{-2}\right)=bulk density \left({\mathrm{g} \mathrm{cm}}^{-3}\right)\times soil depth \left(\mathrm{cm}\right)\times \% organic C$$

The total soil carbon pool was determined by summing the carbon mass of each of the sampled soil depths extrapolated into per hectare basis (Mg C ha⁻¹).

2.5. Belowground Wood Sampling

Wood fragments were recovered from soil samples at various depths within the 1 m peat profile sampled with the McCauley auger, and down to 2.3 m depth using a bucket auger. This was the maximum peat depth encountered at the site, furthest inland from the shoreline, below which lay Pleistocene mineral sedimentary deposits. The deepest wood sample collected was embedded in the mineral sediments, in what appeared to be a buried soil A-horizon, suggesting presence of an upland forest at some point in the past. We used the radiocarbon dating method to determine the age of the buried wood samples [33]. Sixteen of these samples were sent off for ¹⁴C dating at the National Ocean Sciences Accelerator Mass Spectrometry Facility at the Woods Hole Oceanographic Institution.

2.6. Data Analysis

We used ANOVA and Tukey’s HSD Test for comparison between paired plots, testing the level of significant differences at a 95% confidence level. Smoothed-curve fittings were carried out with locally weighted logarithmic regression in ggplot2 package [34]. Contour plots were used to present a 2-dimensional surface by plotting the aboveground C and SOC as contours in x-axis (latitude) and y-axis (longitude). Contour plots were generated using plotly [35], tidyverse [36], and reshape2 [37] packages. All analyses were processed in R version 4.1.1 [38].

3. Results

3.1. Overstory Species Composition, Tree Density, Mortality and Diameter

Overstory species composition, tree density, and mortality differed significantly between the ghost forest and the healthy forest (

Table 2. Overstory species composition, tree density, and decay class of dead trees in healthy forest and ghost forest in the northern part of the Alligator River National Wildlife Refuge in eastern North Carolina. Different letters after a total value denote significant differences at p < 0.05 between ghost and healthy forest.

Location	General Species Composition	No. of Trees Per ha−1			* Decay Class of Standing Dead Trees
		Healthy/No Decay	Standing Dead/Dying	Total (Healthy + Standing Dead)
Ghost forest	Swamp bay	176	210	388	2
	Maple	12	24	36	1
	Pitch Pine	85	178	263	2
	Total	273 a	412 a	687 a
Healthy forest	Maple	63	2	65	1
	Pitch Pine	148	42	190	1
	Sweetgum	10	0	10
	Total	220 b	44 b	265 b

* Decay class 0—alive; 1—dead but small branches attached; 2—dead with only large branches left, 3—large branches gone; 4—top broken off, very decayed.

). The ghost forest had a tree density of 687 trees ha⁻¹ (Table 2), characterized by 56% swamp bay, 38% pond pine, and the remaining 5% was red maple. Sixty percent of the overstory (or 412 trees ha⁻¹) in the ghost forest was standing dead or dying, half of which were swamp bay, 43% were pond pine, and the remainder were red maple. The most prominent decay class of swamp bay and pond pine in the ghost forest was class 2, characterized as dead with only large branches left. In contrast, the healthy forest had a tree density of 265 trees ha⁻¹ (Table 2), with 72% pond pine, 24% red maple and 4% sweetgum. The amount of standing dead trees in the healthy forest was only 17% (42 trees ha⁻¹). Most dead trees were pond pine of decay class 1 (dead, small branches attached). Interestingly, there were no swamp bay trees in the healthy forest, a species that dominated the overstory of the ghost forest (e.g., species replacement during ecosystem transition). While the stand density of the ghost forest was over 2.5 times higher than in the healthy forest, mortality was nearly ten times greater (Table 2). The ghost forest was also characterized by having a significantly smaller diameter and height of trees than the healthy forest (p < 0.05; Figure 3).

3.2. Overstory, Mid-Story and Understory Biomass

Ghost forest and healthy forest differed significantly in total overstory and midstory biomass, but not in understory biomass (

Table 3. Summary of biomass in each vegetation strata by species in healthy forest and ghost forest in the northern part of the Alligator National Wildlife Refuge in eastern North Carolina. Different letters after a biomass value denote significant differences at p < 0.05 between healthy forest and ghost forest.

	Species	Overstory Biomass (Mg ha−1)			Mid-Story Biomass (Mg ha−1)	Under-Story Biomass (Mg ha−1)
		Alive	Standing Dead	Total BIOMASS (Alive + Dead)
Ghost forest	Swamp bay	16.35	10.08	26.43
	Maple	2.69	1.42	4.11
	Pine	11.87	13.57	25.43
	Total	30.91 a ± 6.96	25.06 a ± 6.26	55.97 a ± 12.61		3.65 a ± 2.70	8.32 a ± 4.1
Healthy forest	Maple	8.45	0.23	8.68
	Pine	15.42	2.96	18.38
	Sweetgum	0.93	0.00	0.93
	Total	24.81 b ± 7.25	3.19 b ± 1.65	27.99 b ± 8.75		14.06 b ± 7.30	7.41 a ± 1.28

). On a per tree basis, the healthy forest had a higher average overstory tree biomass (0.11 Mg ha⁻¹), which is 18% higher than in the ghost forest (0.09 Mg ha⁻¹) (p < 0.05). However, given the higher tree density in the ghost forest (Table 1), the overall total overstory biomass in the ghost forest was 55.97 ± 12.61 Mg ha⁻¹ compared to 27.99 ± 8.75 Mg ha⁻¹ in the healthy forest (Table 3). Almost 44% of the total overstory biomass in the ghost forest was in standing dead trees, however, with a total dead biomass of 25.06 ± 6.26 Mg ha⁻¹. In contrast, of the 27.99 Mg ha⁻¹ total overstory biomass in the healthy forest, a total of 88% was in live trees (24.81 ± 7.25 Mg ha⁻¹). The remaining 11% was in standing dead trees, primarily pine (Table 3). Interestingly, the mid-story biomass (3.65 ± 2.70 Mg ha⁻¹) of the ghost forest was only 25% of the amount in the healthy forest (14.06 ± 7.30 Mg ha⁻¹), possibly due to the dense shrub canopy hindering the development of a mid-story in the ghost forest. Finally, understory biomass in the ghost forest averaged 8.3 Mg ha⁻¹ compared to 7.4 Mg ha⁻¹ in the healthy forest, but the difference was not significant due to the relatively large variance.

The ecosystem transition from the healthy forest to the ghost forest resulted in a major shift in overstory species composition and forest structure, from dominance/co-dominance by pond pine/red maple to that of woody bay shrubs/(declining) pond pine, and resulting shifts in biomass distribution (Table 3). In the healthy forest, the total live overstory biomass averaged 24.8 Mg ha⁻¹, with pond pine and red maple accounting for 62% and 34%, respectively. In ghost forest, the total live overstory biomass of 30.9 Mg ha⁻¹ was comprised of 53% bay shrubs and 38% pond pine, although the pond pine was in rapid decline (Figure 1). With this transition, the structure of the forest changed from median diameter/height of 24 cm/11.5 m in the healthy forest to 20 cm/8.5 m in the ghost forest (Figure 3), in addition to the increased density. Finally, the proportion of dead overstory biomass to total overstory biomass shifted from 0.16 in the healthy forest to 0.60 in the ghost forest, signifying impending large losses of ecosystem C as the dead material decomposed in the coming years.

3.3. Ecosystem Carbon Stocks in Vegetation and Soils

As with biomass, the total ecosystem C stocks differed significantly between the healthy forest and ghost forest, as did the distribution between different compartments (Figure 4). The total aboveground C content of ghost forest, 33.9 Mg C ha⁻¹, was 9.2 Mg C ha⁻¹ greater than that of the healthy forest at 24.7 Mg C ha⁻¹. The C content of dead trees in the ghost forest contributed to almost 44% of the total aboveground C, while that contribution in the healthy forest was only 11%. The percentage contribution of overstory live tree C to the total aboveground C was similar in both forests, at 45–50%, but the amount was double in the ghost forest at 27.9 Mg C ha⁻¹, compared to healthy forest at 13.9 Mg C ha⁻¹. The percentage contribution of the mid-story to total aboveground C differed greatly between the healthy forest (28%) and ghost forest (5%), as did the amounts, at 7.0 and 1.8 Mg C ha⁻¹, respectively (Figure 4). Finally, the understory C content of both forests was similar at 3.7 Mg C ha⁻¹ in the healthy forest and 4.2 Mg C ha⁻¹ in the ghost forest, contributing 15% and 12% to total aboveground C, respectively.

3.4. Change in Overstory C Stocks with Distance from Open Water

The visible ecosystem transition occurring at our study site (Figure 1 and Figure 2), was accompanied by a progression of increasing vegetation biomass C and soil organic C with a distance from open water. Averaged over all sampling plots, the total C in overstory trees averaged 13.6 Mg C ha⁻¹ in the ghost forest, compared to 25.0 Mg C ha⁻¹ in healthy forests (Figure 5a,c). There was a correlation (R² = 0.52) between the increasing live overstory C and increasing distance from open water in the healthy forest but not with the ghost forest (R² < 0.01; Figure 6). There was a strong relationship between the increasing soil organic C (SOC) with increasing distance from open water in the ghost forest (R² = 0.60) compared to the healthy forest (R² = 0.30; Figure 5c,d and Figure 6).

3.5. Soil Carbon Stocks

The mean soil organic carbon stocks across the plots in the ghost forest was 412 Mg C ha⁻¹, which was 37.5 Mg C ha⁻¹ greater than the healthy forest: 374.5 Mg C ha⁻¹ (Figure 7). In both forests, the soil C content in the top 50 cm was less than the lower 50 cm. Across the entire vertical 1 m soil profile, the ghost forest had a higher soil C than the healthy forest. The mean soil C in the top 50 cm in the ghost forest was 367.2 Mg C ha⁻¹ and, in the bottom 50 cm, 456.9 Mg C ha⁻¹, whereas the healthy forest had 324.7 Mg C ha⁻¹ and 424.3 Mg C ha⁻¹, for top 50 cm and below 50 cm depths, respectively.

3.6. Wood Sample Carbon Dating and Soil Accumulation

The radio-isotope of C (¹⁴C) in wood samples recovered from the peat during soil sampling was used to age the soil profile and calculate peat accumulation rates (

Table 4. Depth of soil where 16 wood samples were collected, the estimated age (approximate calendar year AD) and calculated soil accumulation rate using the method of McTigue et al. (2019) [26].

Soil Depth (cm)	Age (Apprx. Calendar Year BP 2015)	Soil Accumulation Rate (mm/Year)
−85	1790	3.78
−90	1380	1.42
−90	1835	5.00
−92	1830	4.97
−95	1485	1.79
−95	1870	6.55
−160	865	1.39
−168	865	1.46
−169	805	1.40
−175	1415	2.92
−178	1155	2.07
−179	875	1.57
−189	1120	2.11
−190	1775	7.92
−195	265	1.11
−234	235	1.31

). Wood samples ranged in age from ca. AD 235 to ca. AD 1870, with a fairly strong relationship between sample age and peat soil depth in soil (Figure 8a). It is notable that the deepest (−234 cm) and oldest (AD 235) wood sample was recovered from the A-horizon of a buried mineral soil profile, while all the rest were recovered from the overlying peat. As the sample age approached modern times, we found a decreasing soil accumulation rate (y-axis near 0, Table 4; Figure 8b). The rate of soil accumulation over time ranged from 2.59 to 7.92 mm year⁻¹. To obtain the annual soil carbon accumulation rate, we divided the average carbon stock between the ghost forest and healthy forest (393 Mg C ha⁻¹) by the age of the wood samples following the methods of [25]. The resulting average annual carbon accumulation rate was roughly 0.21–1.67 Mg C ha⁻¹ year⁻¹.

4. Discussion

4.1. Species Composition and Forest Mortality

Historical ditching and drainage along lower coastal plain roads and other infrastructure alter site hydrology, and potentially saltwater exposure regimes, in ways that affect contemporary ecosystem function, with implications for responses to future environmental change. The ditch, which provides drainage for the road, may also allow the inland intrusion of fresh- or saltwater during storm surges further than would normally occur. In contrast, the roadbed can act as a dam, blocking floodwater from inundating ecosystems on the other side. In the current study, this altered hydrology due to a roadbed–ditch system resulted in the conversion of what we termed healthy pond pine forest into actively transitioning ghost forest. The strong differences in species composition between the ghost forest and healthy forest was the result of this disturbance to the hydrologic and, potentially, the salinity regimes. The composition of wetland vegetation is strongly controlled by the wetland water level and hydroperiod [39,40]. Physical and environmental stressors may exclude several plant species, and the duration and intensity of flooding serves as environmental determinants of plant species composition [9]. Swamp bay, pond pine, red maple, and sweetgum are typical overstory tree species found within the Alligator River National Wildlife Refuge [41]. The healthy forest studied here was dominated by a pond pine overstory with the other species as associated sub-canopy co-dominants typical of much of the refuge. However, the adjacent ghost forest, just on the other side of the roadbed–ditch, was characterized by the almost complete mortality of the pine overstory, and the replacement of the understory by dense thickets of swamp bays and other woody shrubs. This type of ecosystem transition (replacement of forest by shrub-dominated systems) was also typical over much of the refuge, and was due to the shrubs having a higher flood tolerance, even when compared to the moderately flooded tolerant pond pine, as was reported for other similar systems [3].

According to personnel from the refuge and the ARNWR Comprehensive Management Plan [18], the roadbed–ditch system was constructed at the study site in 1965. However, it is interesting to note that the remotely-sensed signs of forest decline only began to appear after the mid-2000s (Figure 1). This suggests the environmental stress factors leading to large-scale ecosystem transition became severe only relatively recently. Although pond pine is considered an invader species in wetland coastal plain settings, it is unusual for this species to mature in waterlogged areas [42]. This supports our hypothesis that the 67% overstory mortality in the ghost forest could, at least in part, be due to the increasing inundation at the site in recent years. The decreasing canopy cover of the dying pond pine results in high understory light levels, allowing for very flood-tolerant swamp bays to flourish, changing the species composition. The continuing dominance of a pond pine overstory in the healthy forest implicates the roadbed–ditch system as a causative factor driving this ecosystem transition by slowing site drainage, thus lengthening the inundation period (hydroperiod), leading to tree mortality. However, the 54% mortality of swamp bays and other woody shrubs in the ghost forest, indicates that the hydrological stress is exceeding the tolerance limits of even these very flood-tolerant plants as part of the ecosystem transition to ghost forest. This suggests that the current shrub-dominance is only a transitory state, and that further ecosystem transition occurs with continuing increases in the severity of the driving environmental stress factors.

In addition to hydrologic stress, saltwater intrusion was implicated as a factor driving the formation of ghost forests along the southeast Atlantic and Gulf coasts [12]. In the early stages of saltwater intrusion, live trees may suffer reduced annual growth, and then forest distress becomes more visible. Young trees die conspicuously and new tree recruitment ceases [43,44]. This cessation of recruitment skews tree age distributions towards older age classes until these older remaining trees eventually die, forming relict trees that will become ghost forests. Although we did not measure salinity exposure regimes in the current study, accordingly, trees are most vulnerable to salinity stress during germination and as seedlings [45]. Additionally, even the most flood-tolerant species require periodic soil aeration for successful seed germination [3]. Therefore, the mortality occurring in our study may have been driven by the synergistic impacts of salinity and inundation, a similar condition reported by other studies [45,46,47]. High mortality rates and the subsequent large amount of forest carbon losses in coastal forests was attributed to sea level rise, drought [48], hurricane/storms [10,13], and hydrological connectivity from dams, ditching, and canals [49]. The increasingly apparent forest mortality and ecosystem transition, of which the current ghost forests appear to be at a seral stage, suggest that the environmental drivers recently surpassed the threshold of tolerance of coastal pond pine in eastern North Carolina, accelerated by the effects of human infrastructure on hydrology and, potentially, salinity regimes.

4.2. Aboveground Biomass and Carbon Content of Vegetation

The proliferation of swamp bays and the decline of overstory pines with ecosystem transition to ghost forests has implications for aboveground biomass and the carbon content of these coastal ecosystems. The capacity of forest ecosystems to store C in biomass varies depending on the species composition, age, and population density of each layer of the plant community, N:P ratio, total nutrient supply, morphological and physiological traits of plants, the amplitude of water level, water pH and electrical conductivity [39,50]. Despite the fact that the healthy forest had a higher average aboveground biomass on a per tree basis, the ghost forest, with a higher number trees of smaller diameter, had a higher total aboveground biomass than that of the road side. The intraspecific competition of trees in higher population densities usually led to a higher total aboveground carbon content [51].

We did not measure the water table depth in this study; however, a carbon flux monitoring station within ARNWR recorded a prolonged hydroperiod (2009–present) during non-growing seasons [52,53,54]. This monitoring station is located inland of the ARNWR [44]. Since our study plots were located in the outer part of ARNWR and were closer to a body of water (East Lake), we supposed our study sites would have experienced more severe, prolonged flooding than the carbon monitoring site. Hydrological conditions significantly influence wetland functioning [55]. Water level affects oxygen availability, decomposition rates, and nutrient mobilization from the sediment [56,57] and thus affects biomass production [58,59]. The water level also influences the ability of seeds to germinate and establish [60].

A healthy forest has a multi-layered structure that includes overstory, mid-story, and understory vegetation. Forest structure is a very critical component affecting natural tree regeneration and forest stability [61,62]. However, our results showed an imbalance in forest structure in both the ghost and healthy forest. The overstory trees have a far greater density than the mid-story and understory stratum. In the ghost forest, only a small percentage of regeneration species live past the mid-story stage, as shown by low biomass in the mid-story level (Table 3). An increase in overstory structure induces a reduction in mid-story or understory vegetation due to competition for resources such as light, nutrients, and space. Maintaining enough understory vegetation biomass should be considered in the management goals of any forest ecosystem.

A higher overstory canopy density and biomass in both the ghost and healthy forests in our study ultimately has a negative effect on the understory layer. This negative impact significantly limits the restoration of proper stand structure in this wetland forest. However, there are also positive effects with rising overstory tree density as this increases the tree number, allocation fraction of photosynthetic products in stems, and the aboveground net primary productivity [63]. As a result of the interaction between the positive and negative effects of tree density, the aboveground biomass of trees increases with a higher tree density but at a declining rate after a period of time [51].

4.3. Soil Carbon

Our results confirmed the relevance of freshwater wetlands within ARNWR as a C sink. Total carbon (C) stocks in the study sites ranged from 749.0 ± 170.5 Mg ha⁻¹ in the healthy forest to 824.1 ± 46.2 Mg ha⁻¹ in the ghost forest. These values are higher than those published for other eastern Mountains and the Upper Midwest, averaging 478 ± 58 Mg ha⁻¹ to 539 ± 47 Mg ha⁻¹ in the top 100 cm soil depth [4]; 533Mg C ha⁻¹ in organic wetlands of southeastern Atlantic coastal plains [64]; and 586 Mg C ha⁻¹, 100 cm soil depth [65] in forested wetlands in New England, USA. Our total soil C was also higher than (195 ± 25 Mg ha⁻¹) that reported in the Interior Plains of the Dakota Praire [66,67] and 198 ± 21 and 216 ± 30 Mg ha⁻¹ in the Coastal Plains and West [4]; 340 ± 94 Mg ha⁻¹ in the Gulf of Mexico measured to soil depths of 100 and 200 cm [16]; 109 Mg C ha⁻¹ in flats and levees; and 193 Mg C ha−1 in mineral wetlands of the southeastern Atlantic coastal plains [64].

The depth of peat horizons is mainly determined by the activity of soil organisms, the repeated interruption of humus accumulation by water dynamics, soil micro-relief, the distance to bodies of water, production and decomposition dynamics, organic and mineral depositional history [2,16], which might cause the low soil C stocks found in top 0–50 cm soil horizons in both the ghost and healthy forests. By contrast, the higher soil carbon below (50–100 cm depth) in our study indicates the episodic burial and preservation of soil C due to anaerobic soil conditions that slow the mineralization of the buried C [64]. However, a study reported that soil layers below 30 cm deep contain substantial cumulative reservoirs of carbon, with 65% of the total wetland soil carbon stored between 30 and 120 cm [4]. It must be noted that we only determined C stocks of up to a 100 cm depth and that total stocks may be even higher at deeper layers. Soil sampling in deeper layers is likely needed for the more effective evaluation and modeling of wetland soil carbon dynamics.

4.4. Soil Accumulation over Time

The organic sediment of the soils we analyzed was deposited over 1870 years, which indicated the soil’s capacity to sequester C in the long-term (soil-carbon-memory effect) [2]. Fixed C may be released as CO₂ by microorganisms during decomposition or stored belowground, where anaerobic decomposition occurs slowly, and recalcitrant C may remain for millennia [16]. The ongoing delivery of sediment due to flooding leads to sediment deposition, the burial of organic matter, and the vertical accretion of marsh surfaces, thus allowing carbon to accumulate over long periods [4]. Increasing rates of sea level rise may also contribute to soil accretion by increasing the duration of inundation and increasing sediment deposition on marsh surfaces [68]. Anthropogenic impacts such as the presence of ditches and dikes also have a strong association with affecting soil carbon deposition [4]. A slightly lower soil carbon stored in the healthy forest might be due to the presence of ditches increasing stream and base flows, thereby increasing the annual discharge, which leads to lower (that is, drier) groundwater levels. This lower groundwater level over time could increase soil carbon oxidation and affect soil carbon stores [4]. More stored carbon in the ghost forest might be due to the occurrence of constantly saturated conditions, allowing for more efficient C mineralization at the soil surface and limiting microbial respiration at depth [64], and some contribution of dying overstory trees to soil C inputs.

Our mean soil accumulation rate was 2.59 mm year^−1, yet we found a decreasing trend as our site aged closer to modern times some 235 to 265 years ago (1.11 mm year⁻¹ to 1.31 mm year⁻¹) (Table 4). This decline in trend may occur because the vertical accretion rates overlap with the range of sea level rise rates. Numerical models predict that the maximum vertical accretions rates 2100 years ago were generally 5–30 mm year⁻¹ [69]. Once these maximum limits are achieved, wetland drowning can be observed [12,70]. When these threshold rates of sea level rise are exceeded, wetlands must migrate laterally into submerging uplands to survive [12].

Vegetation at our site persisted for at least 1870 years, resisting drowning to sea level rise for many centuries because natural ecogeomorphic feedbacks allowed coastal habitats to vertically accrete and keep pace with the sea level [71]. However, a study reported that sea levels have risen at an average rate of 2.1 mm year⁻¹, representing the steepest century-scale increase of the past two millennia [72], and this rate was initiated between AD 1865 and 1892. Another study also reported a rise in sea level from 2.2 mm year⁻¹ to 2.4 mm year⁻¹ around the 18th century [73]. If we follow the result of these studies with a sea level rise rate of 2.1 mm year⁻¹ to 2.4 mm year⁻¹, this means that the sediment accumulation in our study sites of less than 1.5 mm year⁻¹, during the turn of 18th century, was less than the estimated rate of sea level rise (2.1–2.4 mm year⁻¹). These estimates suggest that our site may have already been slightly submerged long before and has transitioned from a depositional to an erosive environment and experienced a drowning effect [12], slowly migrating laterally towards uplands to survive as it approaches modern times. Recent global modeling suggests wetland migration into submerging uplands is the single most significant factor influencing wetland area over time, and that global wetland area could increase by up to 60% by 2100 for a 1.1 m sea level rise [74]. Our study site also appeared to be undergoing a transition from organic to mineral accumulation based on the decreased organic matter and sediment accumulation heading towards the open water. This stratigraphic sequence was indicative of submerging marsh [16].

Our estimated carbon accumulation rate (CAR) of 0.21–1.67 Mg ha⁻¹ year⁻¹ aged 235–1870 years old at both the ghost and healthy forests (Table 4 aligns with other studies conducted elsewhere in the United States. Johnson et al., 2007 [75] collected and aged a New England salt marsh core to 3700 cal BP and reported a CAR of 0.4 Mg ha⁻¹ year⁻¹. Brevik and Homburg, 2004 [76] collected a core in Southern California that formed over the course of 5000 years and measured a CAR of 0.3 Mg ha⁻¹ year⁻¹. Drexler, 2011 [77] measured CAR in cores >6000 years old and found CARs from 0.38 to 0.79 Mg ha⁻¹ year⁻¹. A study also found a CAR of 0.39 Mg ha⁻¹ year⁻¹ when integrated over ~2400 years at 2.2 m depth of salt marsh [25]. These results suggest that carbon storage in marshes can be enhanced by SLR-driven OM burial. We caution that our sampling sites only covered a few plots with a small sampling size, and our estimates came with some degree of uncertainty. Telford et al., 2004 [78] compared different age–depth modeling techniques for radiocarbon-dated sediments and found that uncertainty is often underestimated, and the errors were surprisingly large, especially when there were few dated levels.

5. Conclusions

Our study quantifies the impact of roadside ditching and drainage on ecosystem transition from pond pine to shrub pocosin along the lower coastal plain of North Carolina. Ditches and the roadbeds protect forest on the roadside but accelerates the ecosystem transition to a ghost forest on the ditch side, which has implications for the aboveground C storage and cycling. Our study helped us to better understand how coastal wetland soil carbon accumulation changed over time and the response mechanism of vegetation to these changing hydrological perturbations. Our investigation contributes to providing information on the current status and past trends in coastal wetland soil dynamics that may provide valuable insights into the future responses of the wetland ecosystems to sea level rise and climate change. The mismanagement of the wetland ecosystem may convert the present soil C pools into net C sources in the future. Therefore, conservation strategies for the wetland ecosystem should, to a greater extent, incorporate the functions of disturbed and undisturbed wetland forests as C sinks. Flood control structures prevented the migration of outer bank marshes into inland wetlands. The management of wetland ecosystems must therefore consider the cost of constructing more conventional flood control structures to the cost of damages associated with flooding and the danger of inland forested wetlands becoming ghost forests. Choices between defending the coast from sea level rise and facilitating ecosystem transition play a critical role in determining the fate and function of coastal wetland ecosystems. Accurate carbon accounting in wetlands is vital to reduce the risk of climate change contributions by identifying and protecting wetland landscapes that hold disproportionately large soil C stocks.