Evaluating the effects of land-use change and future climate change on vulnerability of coastal landscapes to saltwater intrusion

Global sea level is rising, but the consequences of sea level rise for freshwater dependent ecosystems depends on topography, hydrology and land use among other factors (IPCC, 2013). In coastal plain regions, which are characterized by low topographic relief, freshwater-dependent ecosystems such as swamp forests and freshwater marshes are particularly susceptible to saltwater exposure through the landward movement of marine salts via surface or ground water connections (Herbert et al. 2015, White and Kaplan 2017). Such exposure is one category of saltwater intrusion, a term often limited to studies of groundwater salinization in coastal regions. Here, however, we expand the term saltwater intrusion to include salinization of surface drainage ways, including natural stream channels as well as artificial canals and ditches (as per Herbert et al., 2015).

Where wetlands have been drained and tilled, elevation loss through soil oxidation and enhanced hydrologic connectivity via ditch and canal systems may exacerbate both the frequency and the extent of saltwater intrusion by facilitating the movement of water between interior coastal landscapes and the sea (Hackney and Yelverton, 1990; Turner, 1997; Doyle et al., 2007; Manda et al., 2014; White and Kaplan, 2017). The projected increases in drought severity and long-term sea level rise will act in concert with these land-use changes to extend the upstream and upland reach of saltwater with deleterious effects on vegetation, soil biogeochemistry and water quality (e.g., Anderson and Lockaby, 2012; Ardon et al., 2013; Helton et al., 2014).

Vulnerability of coastal landscapes to saltwater intrusion is a frequent topic of scholarly discussion (see Herbert et al., 2015); however, robust tools have not yet been developed to forecast the potential inland extent of saltwater intrusion under various scenarios of climate change. Sea level rise models simulate changes in water elevation and inundation (so-called “bathtub models”, Titus and Richman, 2001), but hydrologic mechanisms that allow marine salts to penetrate many kilometers inland are missing from most modeling efforts. Extending sea level rise predictions to include the much larger extent of saltwater intrusion vulnerability requires incorporating knowledge about hydrologic networks and flows. In low relief landscapes, it is primarily freshwater flow that limits the movement of saline water inland, against the intuitive land-to-sea gradient. During periods of extended drought, a lack of flow together with strong landward winds off the coast can push saline water inland (e.g., Day et al., 2007; Anderson and Lockaby, 2012; Ardon et al., 2013).

Here we develop a tool, the Saltwater Intrusion Vulnerability Index (SIVI), which combines the effects of sea level rise with freshwater flow routing and accumulation in order to assess the vulnerability of coastal ecosystems to saltwater intrusion via surficial drainage ways. In this study, we apply SIVI to a portion of the Albemarle-Pamlico Peninsula (APP), a coastal landscape in North Carolina where most of the landscape is below 1.5 m mean sea level in elevation with respect to National Geodetic Vertical Datum (Moorhead and Brinson, 1995; Titus and Richman, 2001). Recent projections from a range of emissions scenarios suggest that North Carolina will likely experience between 24 and 132 cm of relative sea level rise (SLR) by the year 2100 (Kopp et al., 2015), putting a significant fraction of the APP at risk of saltwater intrusion within the next century. Because of its low elevation and topographic relief, much of the APP landscape was historically occupied by wetlands (estimated wetland percent cover: 53%, wetland types below 1.5 m elevation: palustrine: 87%, estuarine: 23%) (Moorhead and Brinson, 1995). Extensive land use change in the region began in the 1920’s when private agriculture and forest companies built roads and canals to drain and transform the APP’s extensive wetlands for crop cultivation and plantation forestry (Carter, 1975; Richardson, 1983). Previous studies of ecosystem change in the APP have focused on changes in flooding patterns (Poulter and Halpin, 2008), vertical accretion of wetlands (Moorhead and Brinson, 1995), salinity effects on vegetation (Powell et al., 2016), coastal erosion (Phillips et al., 1993; Wang and Allen, 2008), and the biogeochemical impacts of saltwater intrusion (e.g., Ardon et al., 2013; Helton et al., 2014; Ardon et al., 2016, 2017).

In this study, we examine the extent to which artificial drainage infrastructure has altered the potential for saltwater intrusion in the APP. We hypothesize that current drainage infrastructure supporting agriculture in this region has increased the vulnerability to saltwater intrusion in surface waters in this region by altering patterns of flow routing and flow accumulation. We use a high resolution digital elevation model (DEM) and related geospatial datasets to identify and remove artificial drainages and other human infrastructure. We perform geospatial analyses on the original DEM with artificial drainages and the modified DEM without artificial drainage infrastructure. We use these analyses to estimate the effect of highly connected artificial drainages on hydrological flow in a low gradient coastal landscape, and we discuss implications for major land uses within the region.

We used a LiDAR-derived DEM raster dataset with pixel size of 6.1 meters covering the entire study area (approximately 652 km²; Figure 1a) in Tyrrell County, North Carolina (North Carolina Floodplain Mapping Program, http://www.ncfloodmaps.com). This DEM is part of a statewide dataset collected in the wake of major hurricane-related flooding in the 1990s. The vertical LiDAR accuracy of the DEM was +/– 0.13 m. 99% of elevations within the DEM range from 2.03 m above mean sea level (a.s.m.l.) to –0.26 m a.m.s.l. We acquired an artificial drainage network from the U.S. Geological Survey (USGS) National Hydrography Dataset (NHD), representing both natural and artificial water bodies in the region (US Geological Survey, 2013). Feature attributes in the NHD dataset identify artificial canals and drains.

The original, LiDAR-derived DEM represents the post-development DEM. The original DEM was used to simulate a landscape with no artificial drainages i.e., representing a landscape before development started in this region. To identify the pixels that contain artificial drainages, the original DEM was processed through a range filter, which generated a raster where each output pixels contained the difference between maximum elevation and minimum elevation for a 3-by-3 pixel neighborhood around the corresponding pixel in the original DEM. Range filter allows the identification of the pixels that contain artificial drainages. We then applied a threshold of the 95% percentile of the pixel elevations (approximately 0.5 m) to select the artificial drainages on the original DEM. We overlaid the artificial drainages from the NHD network and there was 85% agreement between the NHD network and the network of artificial drainages identified using the range filter with a 95% percentile criterion. We then interpolated the elevation of those pixels identified as artificial drainages using a median filter across a 30-by-30 pixel moving window. The resulting interpolated raster provides a reasonable estimate of the landscape without artificial drainages (Figure 1b) in the form of a DEM which represents the pre-development conditions on the APP. The raster had the same spatial extent and pixel resolution as the original DEM. MATLAB R2014a (Mathworks, 2014, Natick MA) was used for these processes and a workflow is included in Figure S1.

Flow accumulation in the pre- and post-development landscapes were simulated using the ArcHydro extension in ArcGIS 10.5.1 (ESRI, 2010, Redlands CA). ArcHydro tools have widespread applications in hydrologic modeling and analysis. Both pre- and post-development DEMs were filled to remove sinks with undefined flow direction. We used the artificial drainage network shown in Figure 1a to burn artificial drainages (Hellweger, 1997) into the post-development DEM to force flow routing in a landscape with artificial drainages. In addition to burning the stream or drain pixels, we used linear interpolation to create a smooth slope of five pixels on each side of the burned stream pixels. This step was necessary to ensure that flow direction and flow accumulation rasters (described below) reflected actual hydraulic gradients observed in artificially drained coastal environments (O’Doherty, 2013).

For both pre- and post-development DEMs, the flow direction for each pixel was computed. The flow direction raster was used to compute a flow accumulation raster with D8 connectivity (a pixel is connected to its cardinal and diagonal neighbors), which represents the number of upstream pixels that flow into each downslope pixel. Streams were then delineated from the flow accumulation raster by applying an area threshold of 0.37 km². This threshold was determined through visual inspection of natural drainages found on the pre-development DEM and consistent with drainage area of first order watersheds in low-gradient coastal plain in eastern North Carolina (e.g., Chescheir et al., 2003) and in Maryland (e.g., Angier et al., 2005) and thus represents a realistic average stream generation threshold for the area (Figure 1b). The artificial and natural drainages from Figure 1a and 1b were used to calculate the drainage density (drainage length per unit land area) of the landscape.

We used information on elevation and flow accumulation to create a spatially explicit saltwater intrusion vulnerability index (SIVI) for the study area. The index is determined by the ratio of saltwater intrusion risk, defined as the reciprocal of depression-filled elevation above sea level (z^-1), to freshwater subsidy potential derived from the log transformed flow accumulated area.

We based the numerator of SIVI on the depression-filled elevation to account for the fact that some low-lying pixels are topographically isolated from the coastline and, therefore, protected from inundation until water level reaches a higher threshold than the actual elevation of the pixel. The reciprocal of this value acknowledges that higher depression-filled elevations are less vulnerable to saltwater exposure.

A log-transformed flow accumulation (A) raster describes the topographic concentration of surface flow on the landscape (Seibert and McGlynn, 2007). This raster was modified by a series of three low pass filters to better capture the relatively diffuse nature of flow accumulation in this low-relief environment. The resulting raster (F_a) represented the ability of freshwater flows to counteract the encroachment of saltwater into the interior of the landscape (Warner et al., 2005). The complete SIVI was thus computed as,

The SIVI units of m^–1 do not have a strict physical interpretation but represent the low-lying pixel’s elevation and the accumulation of upstream pixels flowing into that same pixel that can counteract saltwater intrusion. We use the numerical values of SIVI to compare scenarios – higher SIVI values indicate greater vulnerability to saltwater intrusion, and lower values indicate less vulnerability.

We selected four, square areas within the study area (2.28 km² each; 249 × 247 pixels) representing a gradient in drainage density (natural to artificial). These areas represented the range of drainage densities observed, and for each area we evaluated the effects of artificial drainage infrastructure on elevation (z), flow accumulation (A), and SIVI. We used the land cover data from the NOAA’s Coastal Change Analysis Program (C-CAP) for the four major land-covers across the entire study area. The four classes and their areal coverages for the study area included wetlands (59%), agriculture (25%), shrub/scrub (8%) and forest (4%). Of the 59% wetland cover, palustrine wetland had an areal coverage of 99.3% over the study area. The estuarine wetlands had less than 1% coverage for the study area. The shrub/scrub and forest land-cover classes represent the upland systems and do not include palustrine and estuarine forested, shrub wetlands which are included in the wetlands class. Also, the deciduous, evergreen, forest land cover types were combined for the forest land-cover class.

To test the difference between the pre-development and post-development elevation, flow accumulation, and SIVI, we performed a Two-Sample Kolmogorov-Smirnov test with a null hypothesis that the hydrological variables in the pre-development scenario were from the same continuous distribution as the post development scenario. Similarly, we performed a Wilcoxon Rank Sum Test to determine differences in the median of SIVI between land cover classes.

Traditional approaches to model validation are difficult to apply in the case of this saltwater intrusion vulnerability model. We do compare our SIVI estimates to the results of a large-scale survey of surface waters and soil solution throughout the peninsula (Figure S2). Although we make these comparisons, we approach such an exercise with limited expectations. There are two major limitations on the predictive capacity of our SIVI estimates. First, we lack key information about drainage control structures. Our model identifies areas of the coastal plain that are susceptible to saltwater intrusion because of their elevation and surface water connectivity to the estuary but we lack empirical data on interruptions to these surface water connections. Control structures such as flood gates and pump stations are regular features of this landscape, many of the agricultural fields in this region require active drainage to maintain their suitability. Our SIVI estimates suggest what is possible in the absence of these control structures. The second limitation of our model is that saltwater intrusion, unlike sea level rise, does not represent a permanent state change. Instead, saltwater intrusion may occur episodically as a result of storm surges or seasonally as a result of drought (Ardon et al., 2013; Manda et al., 2018). The salinity signals that can be detected in the small number of permanent monitoring locations in the region during these events are subsequently lost as rainfall flushes these marine salts back out to sea. The extent to which a residual chemical signal is left behind by saltwater intrusion will depend upon drainage and precipitation – with the very features that make an area susceptible (hydrologic connectivity) tending to reduce the longevity of elevated salinity.

For the post-development scenario, the drainage density was 2020 m per km² whereas for the natural pre-development scenario, the drainage density was 650 m per km² (Table 1). The probability density function (PDF) of elevation for the pre- and post-development scenarios are shown in Figure 2a. The post-development elevation does not reflect the ditch burning processing step carried out during flow accumulation. Most of the artificial drainages were located on land within 1 m of mean sea level. Similarly, flow accumulation (A) of DEM with or without artificial drainages is given in Figure 2b. The flow accumulation values for the post-development scenario are abundant and less than those for the pre-development scenario. The elevation and flow accumulation in the pre-development scenario are significantly different at p < 0.01 (Table 2).

The SIVI values derived from the pre- and post-development scenarios in Figure 1 are shown in Figure 3a (SIVI in the presence of artificial drainages) and Figure 3b (SIVI in the absence of artificial drainages). In both the pre- and post-development scenarios, land closer to the coastline tends to be more vulnerable to saltwater intrusion than the inland areas. Key characteristics of SIVI from the pre- and post-development scenario, including the maximum, 75^th percentile, median, 25^th percentile and minimum SIVI values are shown in Table 1. The SIVI for post-development DEM tends to be higher (all except minimum value) than SIVI for the pre-development DEM (Table 1). Figure 3c and 3d show a ΔSIVI raster (pre-development DEM subtracted from post-development DEM) for two sub-sections, where the positive values (post development SIVI > pre-development SIVI) are shown in warmer colors and median ΔSIVI value was 0.3. Overall, 66% of the study area had positive ΔSIVI values, meaning that vulnerability to saltwater intrusion increased for most of the study area following development.

We selected four 2.28 km² sub-sections within the study area based on different combinations of natural and artificial drainage densities (Figure 4a). The sub-sections, ranked in order of their artificial drainage densities were North (6155 m per km²), South (3985 m per km²), West (1352 m per km²) and East (0 m per km²). Their natural drainage densities were as follows, West (1704 m per km²), East (595 m per km²), South (476 m per km²), and North (240 m per km²). The North sample had the greatest number of artificial drainages but the lowest natural drainage density (Δ drainage density = artificial – natural; North (5915 m per km²), South (3509 m per km²), West (–352 m per km²) and East (–595 m per km²). The East sample on the contrary had no artificial drainages and was the only sample within our study area comprising mostly of unditched land. Elevation, flow accumulation and SIVI for all samples are included in Figure 4b. The elevation, flow accumulation, SIVI for all sub-sections with artificial drainages were significantly different than without artificial drainages (Table 2).

We plotted the median values of z^–1 and F_a with lines showing ±1 standard deviation for four land-cover classes in the present-day, post-development scenario (Figure 5). Conceptually, the x-axis of this figure (F_a) represents freshwater subsidies that keep saltwater at bay, and increasing values along this axis represent increasing ability of the landscape to withstand droughts. The y-axis (z^–1) represents vulnerability to saltwater intrusion due to low elevation, and greater y values are associated with increased vulnerability to sea level rise. The slope between any point and the origin (z^–1/F_a) is SIVI (equation 1), and four reference values for SIVI are shown as lines on Figure 5. The highest reference value of 5 for SIVI corresponds, approximately, to the most vulnerable quartile of pixels in our post-development scenario. The median SIVI for the wetland land-cover class is greater (p < 0.01) than the median for other agriculture, forest or scrub/shrub land cover classes. The median SIVI for palustrine wetland is also greater (p < 0.01) than the median SIVI of estuarine wetland (Figure 5 inset).

Through the construction of canals and ditches, humans have increased drainage density by approximately three-fold across the study area. Our results from the post-development scenario (2 km per km²) together with other estimates of artificial drainage density on the APP of between 3 km per km² and 12.5 km per km² (Heath, 1975; Poulter and Halpin, 2008; Poulter et al., 2008) show that the region’s flow network has been highly altered by human activity. Earlier work has found that high artificial drainage densities of coastal wetlands are associated with increased runoff and reduced evapotranspiration in wetlands (Richardson, 1983; Richardson and McCarthy, 1994; Holden et al., 2004). This hydrological behavior is quite different from unaltered coastal wetlands, which tend to have relatively high rates of evapotranspiration, low runoff ratios, and slow rainfall-runoff responses that help maintain freshwater flows in coastal streams during inter-storm periods and droughts (Sun et al., 2002, 2006; Lu et al., 2003). These freshwater flows represent F_A in our vulnerability calculations and altering their status by increasing drainage densities has the potential to impact vulnerability of these landscapes to saltwater intrusion. The altered drainage densities, particularly in wetlands, could have serious water quality implications given the important role of wetlands in retaining sediments and nutrients generated in surrounding landscapes (Mitsch and Gosselink, 2008).

Although our analysis does not account for the loss of elevation due to draining and decomposition of peatland soils, this result is expected, in part, as the elevations of artificial drainages tend to be lower than the lands they replace. However, the draining of the peatland soils to support agriculture and forestry in the coastal wetlands can also lead to rapid oxidation of formerly saturated soils, which causes subsidence and further declines in elevation (Richardson, 1983; Holden et al., 2004). The loss in elevation in artificial drainages with concomitant subsidence owing to land-use change then could lead to further loss of elevation. The loss of elevation therefore increases the threat of saltwater intrusion for these already vulnerable areas.

Although changes to elevation across the region may appear subtle, flow accumulation is highly altered by artificial drainage (Figure 2b). These changes have major implications for water flow in the landscape. Because of the extensive network of artificial drainages (Figure 1a), the natural flowpaths (i.e., route of surface water flow), that previously traversed un-channelized portions of the landscape, now terminate in nearby canals or ditches. These changes short-circuit flows that previously moved slowly and diffusely across un-channelized landscapes are now short circuited, reducing the network distance between the coast and inland portions of the landscape.

Conceptually, the short-circuiting of the hydrological flow occurs when unditched areas accumulate very little flow en route to artificial drainages. The larger flow accumulation values of the post-development scenario are concentrated in fewer artificial drainages compared to the pre-development scenario. The average flow accumulation is one order of magnitude greater for the pre-development scenario (8.44 × 10⁴ m²) than for the post development scenario (3.84 × 10³ m²). It is important to note that while it was commonly understood that draining of coastal wetlands lowered the elevation by oxidation of soil, our results show that dissecting the landscape with artificial drainages can have a much larger impact on vulnerability to saltwater intrusion by short-circuiting hydrological flows.

We observed that areas closer to the coastline are more vulnerable in both pre- and post-development scenarios, because they tend to be lower in elevation (Figure 3). These low elevation areas closer to the coast in Northern Tyrrell County are also susceptible to high rates of shoreline erosion (Eulie et al., 2013). Although our vulnerability analysis does not consider erosion susceptibility, we acknowledge that this is a widely reported and well-understood issue for coastal ecosystems and human communities. Given that our analyses do not include coastal erosion, our estimates of vulnerability are likely conservative. Further inland, elevations tend to be greater and flow accumulations larger, hence the vulnerability of the landscape is relatively lower. Though the spatial pattern of vulnerability appears similar in both scenarios (Figure 3), the SIVI values are generally higher in the post-development DEM (median SIVI: 1.62) than the pre-development scenario (median SIVI: 1.28, p < 0.01, Table 2) because the flow accumulation is short-circuited (Table 1). This increase in vulnerability is represented by ΔSIVI (Figure 3c and 3d) and we found that areas in the post-development DEM have higher SIVI than the pre-development scenario (66% of the study area had positive ΔSIVI). With a positive median ΔSIVI value of 0.3, the effect of artificial drainages in increasing vulnerability to saltwater intrusion in this low-lying coastal system is therefore documented. To contextualize SIVI with the ecological effects of saline water on freshwater dependent landscapes, we compare post-development SIVI estimates with ionic concentration of saline water (chloride, sulfate and divalent cations) in soil solution (Figure S2). We observe that SIVI is positively correlated with Cl^– and SO₄^2– concentrations. Despite the caveats we mentioned earlier regarding the stochastic nature of saltwater intrusion events and a lack of complete assessment of human water infrastructure on the peninsula, the linear regression models suggest that SIVI explains more of the variability (higher r²) in saltwater ions than elevation (Table S1).

An important caveat to this finding is that our study does not account for the role of water control structures such as pumping stations and one way flap gates already in place to limit saltwater intrusion. These control structures are difficult to map comprehensively for large regions such as the APP. Such water control structures are one strategy already in place in the APP and elsewhere to reduce vulnerability of coastal landscapes to saltwater intrusion. However, there are significant logistical and financial constraints to install and manage water control structures at every entrance point on the APP to limit saltwater intrusion. For example, Poulter et al. (2008) found that out of 87 potential entrance points on Dare County (east of our study area), six locations control the flow inland. We also note that the same structures, which can serve as barriers to saltwater intrusion during normal weather conditions, may become overtopped or otherwise compromised during severe coastal storms and may then artificially impede the natural drainage of saltwater from the landscape following the event.

Distributions of elevation and flow accumulation within our four high resolution sampling areas (spanning a range of natural and artificial drainages, Figure 4b) are similar to the elevation and flow accumulation PDFs for the entire study area (Figure 2). Though there are no artificial drainages within the East sample, the flow accumulation of East sample with no artificial drainages is similar to the flow accumulation of other samples with artificial drainages because there are artificial drainages on the edge of the sample area. This result shows that even in the absence of local artificial drainages, local flow accumulation can be affected by artificial drainages elsewhere due to the regional nature of flowpaths and flow accumulation. In other words, the effects of artificial drainages are not localized, but instead manifest at larger spatial scales. For all four of the samples, we found that median SIVI values are lower (p < 0.01) for the pre-development DEM (North pre: 0.42, North post: 0.60; East pre: 1.84, East post: 2.25; West pre: 0.62, West post: 0.75; South pre: 0.46, South post: 0.65), elevation, flow accumulation and SIVI are significantly different without artificial drainages (Table 2) and that the difference in median SIVI values is highest in East sample.

As elevation and freshwater subsidies decrease for specific landscape positions, saltwater intrusion vulnerability increases (Figure 5). Relative to other land cover types, these results show that wetlands on the APP are more vulnerable to saltwater intrusion (p < 0.01). Because palustrine wetlands represent ~99% of the wetlands land-cover class, the median SIVI of palustrine wetlands is equal to the median SIVI of wetlands. The palustrine wetlands however are vulnerable compared to the estuarine wetlands (Figure 5 inset) as the median SIVI of palustrine wetland is higher (p < 0.01) than the median SIVI of estuarine wetland. Because sea level rise reduces the relative elevation of all landscapes and increases the elevation-related component of SIVI, all land cover types should become more vulnerable over time (Figure 5) in the absence of any acceleration of vertical sediment accretion.

If the pace of historic development of artificial drainage on the APP continues for agriculture and timber industry, it is likely that all land cover types will continue to have reduced freshwater subsidies over time, and thus be increasingly vulnerable during drought periods, moving them toward the left side of Figure 5. It is important to note that with SIVI, we cannot forecast how far to the left the shifts will occur for each land-cover type. Droughts are projected to become more common in the southeastern United States during the 21^st century (Carter et al., 2014) and incidence of storm surge flooding in coastal North Carolina is expected to rise with sea level (Kopp et al., 2015). While we do not know the exact timing of droughts and storm surge flooding, there is evidence that droughts have already exacerbated the effects of the sea level rise by pushing saline water inland (Williams et al., 1999; Doyle et al., 2007; Ardon et al., 2013). A recent study on the APP showed that wind tides pushed saltwater in the inland canals and wetlands which then persisted in the canals for up to 10 days (Manda et al., 2014). Besides wind tides, the connectivity of artificial drainage network to open water serve as conduits to bring in saltwater and increase salinity in coastal plain (up to of 6 parts per thousand) during periods of droughts (e.g., Ardon et al., 2013). Because SIVI considers the role of freshwater subsidies in mitigating saltwater intrusion, it could be a useful tool for land managers and decision makers to identify the vulnerable portion of the landscape. Identifying the most vulnerable locations with SIVI is the first important step in developing a monitoring plan to minimize saltwater intrusion inland by installing and managing water control structures at entrance points that control the flow.

Increasing vulnerability of wetlands, forests and other land cover types to saltwater intrusion will have negative ecological effects. For example, the submergence of wetlands due to rising sea levels (Moorhead and Brinson, 1995), conversion of freshwater wetlands to brackish marshes because of saltwater intrusion and SLR (Craft et al., 2009) and forest die-off because of enhanced exposure to salinity (Williams et al., 1999; Day et al., 2007) have implications for carbon cycling and ecosystem services. Most notably, these landscapes serve an important function for maintaining water quality in a region that is surrounded by shallow estuaries that are sensitive to nutrient loading, eutrophication, hypoxia, and harmful algal blooms (Mitsch and Gosselink, 2008; Herbert et al., 2015). These estuaries are important economic engines for the region, supporting important commercial and recreational fishing industries.

The societal impacts of saltwater intrusion and SLR, however, extend beyond the fishing economy. The region has already experienced abandonment of agricultural fields due to saltwater intrusion and reduced soil fertility (Moorhead and Brinson, 1995; Hackney and Yelverton, 1990), and a recent study by Hauer et al. (2016) estimated that 94% of the Tyrrell County residents will be exposed to catastrophic flooding by year 2100. Thus SLR, increased artificial drainage, and drought means that all land cover types will continue to become more vulnerable to saltwater intrusion, and their vulnerability will increase at a pace that is greater than if either SLR or artificial drainage occurred separately. These conclusions support other work showing that the incidence and severity of saltwater intrusion in coastal regions will increase in the coming decades as a result of climate change and human alteration to the hydrological cycle (Neubauer and Craft, 2009; Herbert et al., 2015; Tessler et al., 2015).

In this study, we showed that artificial drainages not only increased the drainage density of the landscape but also altered flowpaths in ways that reduce the ability of freshwater flows to counter the effects of saltwater intrusion. Even in areas that lacked local artificial drainages, hydrologic flows were affected due to the large-scale nature and connectivity of flow networks. Together, these alterations to the landscape increased the overall vulnerability of our study region to saltwater intrusion, which we quantified using the spatially explicit saltwater intrusion vulnerability index. Vulnerability to saltwater intrusion following artificial drainage increased more for palustrine and estuarine wetlands than for agriculture, forest and shrub land cover types. We show how sea level rise, continued infrastructure development and drought all have the ability to stress the coastal landscape by increasing the vulnerability of the region to saltwater intrusion.

The data presented in this work are available through UC Press Dash Digital Repository at https://doi.org/10.15146/R3539J.

We would like to thank TL Jass for his assistance with field work and laboratory analyses. We are grateful to Ellen R. Herbert and an anonymous reviewer for their comments.

The funding of this work was provided by National Science Foundation (EF-1427188 and EAR 1462169).

The authors have no competing interests to declare.

• Contributed to conception and design: REE, ESB, MA, TKB, JPW

• Contributed to acquisition of data: REE, SMA, MGS, EAU

• Contributed to analysis and interpretation of data: AB, REE, ESB, MA, SMA, MGS, EAU, TKB

• Drafted and/or revised the article: AB, REE, ESB, MA, JPW, TKB

• Approved the submitted version for publication: AB, REE, ESB, MA, SMA, MGS, EAU, TKB, JPW