Investigation of Barrier Island Highway and Marsh Vulnerability to Bay-Side Flooding and Erosion

Abstract

Coastal highways along narrow barrier islands are vulnerable to flooding due to ocean and bay-side events, which create hazardous travel conditions and may restrict access to surrounding communities. This study investigates the vulnerability of a segment of highway passing through the Pea Island National Wildlife Refuge in the Outer Banks, North Carolina, USA. Publicly available data, computational modeling, and field observations of shoreline change are synthesized to develop fragility models for roadway flooding and marsh conditions. At 99% significance, peak daily water levels and significant wave heights at nearby monitoring stations are determined as significant predictors of roadway closure due to flooding. Computational investigations of bay-side storms identify peak water levels and the buffer distance between the estuarine shoreline and the roadway as significant predictors of roadway transect flooding. To assess the vulnerability of the marsh in the buffer area, a classification scheme is proposed and used to evaluate marsh conditions due to long-term and episodic (storm) stressors. Marsh vulnerability is found to be predicted by the long-term erosion rate and distance from the shoreline to the 5 m depth contour of the nearby flood tidal channel. The results indicate the importance of erosion mitigation and marsh conservation to enhance the resilience of coastal transportation infrastructure.

1. Introduction

Coastal highways and bridges are critical links in the transportation network for the movement of goods and people on a daily basis and especially for access by emergency services during post-disaster response. Along narrow barrier islands, this infrastructure is typically low-lying and near the water, restricting connectivity and posing continuous risks to operability, maintenance, and resilience. Recent studies predict that the flooding of coastal highways and bridges will continue to increase with rising sea levels and intensifying storm surges due to climate change [1].

Previous studies have analyzed the vulnerability of coastal transportation infrastructure, including bridges [2,3,4,5,6,7,8,9], near-coast roadways [10,11,12,13,14,15], and other infrastructure [16] subject to coastal storms or multi-hazard environmental conditions. Vulnerabilities are often presented in the form of probabilistic fragility functions relating the likelihood of damage or failure to a hazard intensity measure. For example, Kameshwar and Padgett [3] developed parameterized fragility functions assessing the probability of structural failure of highway bridges subject to earthquake and hurricane events. Another study developed a model for coastal bridges subject to sea level rise, landscape changes, and flooding; structural failure was estimated based on storm surge, waves, and inundation duration [4]. Beyond assessing bridge survival or failure, Padgett et al. [2] derived fragility curves considering bridge parameters and environmental conditions, predicting the likelihood of bridges exceeding a given damage state based on a four-point damage scale ranging from minor to complete. Storm surge and the number of spans were significant predictors of damage for the bridges considered in the study.

In addition to bridges, the effects of coastal flood hazards on highways and road networks have also been analyzed [12,17]. In an assessment of the cascading effects of hurricane waves and surges on physical (buildings and roadways) and social systems, Fereshtehnejad et al. [12] assessed roadway failure in Galveston Island, TX, USA, using a fragility model based on the distance to the roadway from the Gulf of Mexico and inundation duration. Modes of roadway failure including surface layer loss, cracking and potholes, and base failure were considered in the model. A failure model for roadways subject to surface runoff-driven flooding events was developed by Wang et al. [10]; the model distinguished between direct and indirect roadway failures based on the roadway’s flooded condition and connectivity to other roadway segments.

Few coastal highway vulnerability studies have also considered local morphological features that affect roadway fragility to flood impacts. For example, Nasrallah [11] used remotely sensed data and a morphological numerical model to forecast the storm impacts on coastal dunes that can lead to overwash and increased vulnerability of the North Carolina (NC) Highway 12. Three vulnerability indicators for coastal roadways in barrier islands based on island width, dune crest elevation above the roadway, and distance from the edge of the pavement to the ocean shoreline were developed by Velásquez-Montoya et al. [14]. Another study examined the effectiveness of distinct morphological indicators in predicting storm impacts on barrier island roadways, with the distance from the edge of the pavement to the dune toe identified as the most effective indicator of highway vulnerability [18]. However, these studies have mainly focused on ocean-side events and features (i.e., dunes), leading to a knowledge gap in the effects of bay-side storm events on coastal transportation corridors based on the extent and condition of the bay-side shoreline. Similar to ocean-side storms, bay-side storm events can cause significant issues for transportation facilities, including flooding, shoreline erosion, the deposition of sediment and vegetation, and wave damage in areas with a large fetch.

Previous studies have considered the vulnerability or response of barrier island and marsh shorelines to anthropogenic stressors [19,20,21], sea level rise [22,23,24], and episodic and long-term processes [25,26,27]. Many studies present marsh response parameters, such as the conversion to open water, erosion or accretion, or the change in elevation or shoreline position [22,23,24,25,27], based on aerial imagery [28], field observations [29,30], and/or computational modeling [31]. However, few studies have synthesized remotely sensed, field-based, and computationally modeled data to stochastically investigate the effects of erosive stressors on the marsh condition as healthy or eroded in varying degrees of severity.

This paper presents a vulnerability assessment for a coastal highway in North Carolina’s Outer Banks that focuses on bay-side impacts by considering publicly available roadway closure data and a synthesis of numerical model outputs and rapid-response field observations following two storm events. We focus on the roadway vulnerability to flooding, which disrupts access and exacerbates the degradation of transportation components. We further consider the interconnected performance of natural elements and transportation infrastructure by considering both the effect of the marsh buffer distance on the likelihood of roadway flooding and the vulnerability of the marsh itself to long-term and episodic erosion.

The remainder of this paper is structured as follows: Section 2 describes the study area, which comprises the section of a state highway passing through the Pea Island National Wildlife Refuge. This section of roadway is vulnerable to bay-side flooding due to the absence of dunes on the estuarine side, leaving the road unprotected during elevated bay-side water level conditions. Methods of (i) analysis of roadway closure based on publicly available traffic information data, (ii) computational modeling of bay-side storm scenarios, (iii) observations of shoreline change through assessments of aerial imagery and shoreline surveys conducted following two storm events, and (iv) fragility model development for characterizing roadway and marsh vulnerability are detailed in Section 3. Section 4 presents results of long-term and short-term shoreline change and significant predictors of roadway or marsh failure, showing fragility curves derived for roadway flooding and marsh erosion. Finally, Section 5 discusses the broader implications, considerations, and conclusions of this study.

2. Study Site

The NC 12 highway is the only coastal roadway connecting the Outer Banks of North Carolina from the communities of Corolla to Hatteras and providing direct access to the barrier islands from the mainland. The projected 2025 average daily traffic volume of NC 12 is estimated to be 9600 vehicles per day and 15,400 vehicles per day during the summer [32]. Given the location of this major roadway along narrow sections of barrier islands and its exposure to storms and high-water events, the road has been the subject of multiple vulnerability studies in the last two decades [13,33,34]. Such studies have identified several vulnerable hotspots where the roadway is subject to frequent flooding and sand burial due to overwash. Some of these vulnerable hotspots are located towards the northern end of Hatteras Island, where the Pea Island National Wildlife Refuge is located.

The specific stretch of shoreline and roadway analyzed here is located on the bay side of the northernmost end of Hatteras Island. The roadway section corresponds to the approach of the southern terminus of the Marc Basnight Bridge (Figure 1), located just south of Oregon Inlet. The bridge is located in Dare County, for which the United States Census Bureau [35] reports a 2021 population of 37,826; however, the daily population can increase to 225,000 to 300,000 during the summer tourist season from June to August [36]. The nearest communities to the study area are Nags Head, located 22 km to the north, and Rodanthe, located 21 km to the south; restricted access to this section of roadway cuts off access to towns further south, including Salvo, Avon, Buxton, and Hatteras Village. While these communities are exposed to coastal flood hazards due to hurricanes, nor’easters, or erosion, they are particularly impacted by damage or closure of this stretch of highway, which may delay or limit access by emergency response teams in the event of closure due to hazardous conditions or adversely affect tourism during summer months.

At this section of the barrier island, the ocean shoreline is accreting due to the presence of the terminal groin on the north end of the island [37]. However, the estuarine shoreline has been eroding at rates up to 3 to 4 m/year [38]. This shoreline erosion has been attributed to the morphological evolution of the adjacent Oregon Inlet and the rotation of the main channel that causes the southernmost flood channel of the inlet (location shown in Figure 1C) to encroach into the down drift back barrier [39]. This region was also identified as a potential location for barrier island breaching by the Federal Highway Administration (FHWA) and the North Carolina Department of Transportation (NCDOT) [32] due to its proximity to Oregon Inlet and the possibility of storm surge flows through adjacent shore-normal estuarine channels.

Although the study site is located on the bay-side of the barrier island system, its proximity to a tidal inlet makes it responsive to both ocean and bay conditions. Monitoring of such conditions at the Albemarle–Pamlico Sound (bay side) and ocean side is available via permanent observational stations. Four kilometers north of the study site, on the northern bay-side of Oregon Inlet, there is an NOAA Tides and Currents station (Oregon Inlet Marina, NC-Station ID: 8652587, herein referred to as the “marina tide gauge”) that has recorded hourly and six-minute water levels since April 1994 and January 1996, respectively, and two-minute scalar average wind speeds, two-minute vector average wind directions, and maximum five-second wind gusts since November 2007. Station 44095-Oregon Inlet, NC, herein referred to as the “waverider buoy”, is owned by the University of North Carolina Coastal Studies Institute and provides the closest wave measurements to the study site, including significant wave height, dominant wave period, and mean wave direction every 30 min. The station is about 18.5 km offshore of the study site at a water depth of 18.3 m and has been recording wave parameters intermittently since 2012 (Figure 1).

The top 10 highest total water levels recorded at the marina tide gauge have been above 1.02 m (referenced with respect to the North American Datum of 1988 (NAVD88)). The top three total water levels recorded, referenced with respect to NAVD88, were those corresponding to Hurricane Irene in 2011 (2.070 m), Hurricane Floyd in 1999 (1.525 m), and Hurricane Michael in 2018 (1.445 m). Numerical simulations of Hurricane Irene indicate that wind-generated surface waves and wind-driven storm surges are some of the most important contributors to extreme flooding along estuarine shorelines [40]. In addition to extreme seasonal events, long-term processes such as relative sea level rise contribute to the potential vulnerability of the stretch of shoreline. The linear relative sea level trend observed at the marina tide gauge is 5.32 ± 1.12 mm/year [41], with low (17th percentile) and high (83rd percentile) relative sea level rise contributions projected for 2050 ranging from 0.29 m to 0.45 m considering a 0.3 m global sea level rise and from 0.42 m to 0.71 m considering a 2.0 m global mean sea level rise [42].

3. Methods

3.1. NCDOT Traveler Information Management System Roadway Closure Analysis

NCDOT supports a Traveler Information Management System (TIMS), available to the public at DriveNC.gov. This system posts real-time notifications of closure or hazardous conditions along North Carolina roadways, enabling travelers to modify behavior accordingly [43]. The notifications and records of closures can also provide a record of occurrence frequency and the location of hazardous conditions. A record of the TIMS data was provided to the authors by NCDOT, comprising incidents on NC 12 along the Pea Island National Wildlife Refuge from September 2017 to November 2019. The dates of closure or hazardous conditions were recorded, and the corresponding daily maximum significant wave height and water level data were obtained from the waverider buoy and the marina tide gauge records, respectively.

3.2. Numerical Model and Description of Storm Scenarios

Field observations of water levels during 4 weeks in 2019 suggest that there are differences in the water level signals between the study site located south of Oregon Inlet and the marina tide gauge located north of the inlet. Such differences in water level signals may have consequences for predicting flooding [44]. The tidal amplitude at the study site is two to three times larger than at the marina tide gauge. The lowpass-filtered subtidal water levels (48 h cutoff period) are correlated with an r² = 0.42, suggesting differences in response to winds, waves, and surges. Given the spatial variability of the water level signals around the inlet, a numerical model was used to predict flooding scenarios for a range of conditions and to specifically investigate flooding by bay-side storm events.

A hydrodynamic, two-dimensional Delft 3D [45] model coupled with a wave model Simulating Waves Nearshore (SWAN) [46] is used in this study. The numerical model resolution varies from 470 m offshore to 15 m and 20 m within the inlet and the marshes behind NC 12. The computational domain includes a large portion of the Albemarle–Pamlico Sound and extends to the edge of mainland North Carolina and 35 m along the barrier island system. Wetting and drying thresholds were adjusted to better represent the flooding extent caused by Hurricane Irene (2011). Implementation of mean depth at the grid cell faces and a threshold depth of 0.01 m resulted in the best match of the simulated flooding extent and the bay-side wrack line as observed from aerial imagery taken after Hurricane Irene.

As part of the calibration and validation of the numerical model, simulations’ outputs were compared to in situ water levels and depth-averaged velocities for 30 days in the fall of 2019 and 15 days in the summer of 2020. The Willmott Skill Scores [47] for both periods for water levels ranged from 0.87 (very good) to 0.94 (excellent), providing confidence in the simulated hydrodynamics at the site. More details on the numerical model set up, calibration, and validation are presented in [14,48]. In addition to comparisons of water levels and currents, the simulated flooding extent was compared against 11 days of field measurements of high-tide flooding on the study site during October 2019 [44]. Peak measured water depths near station S3 on the marsh ranged from 0.3 m to 0.8 m. The root-mean-square (RMS) difference between observed and simulated peak water depths at each high tide for this location was 0.09 m, with the model tending to overestimate small peak water depths (<0.5 m) and show better agreement for higher peak water depths. The R² was 0.80, suggesting good agreement between the simulated and observed flood elevations on the marsh.

A total of 42 storm scenarios with varying water levels and wave conditions at the boundary (water levels at the bay boundaries ranging in 0.5 m intervals from 0.5 m to 3.5 m and significant wave heights at the ocean boundary ranging in 1 m intervals from 2 m to 7 m) were generated in the numerical model. The storm scenarios were set up to simulate bay-side water levels and ocean waves below, equal to, and above those created by Hurricane Irene (2011), which is the hurricane that has generated the largest bay-side surge since 1979 [49].

Monitoring stations were set up in the numerical model along ten cross-shore transects extending from the location of the NC 12 highway to the shoreline; the locations of each numerical model output station are shown in Figure 1C. Each numerical model output station recorded the occurrence (or not) of flooding at that station during a storm scenario. The placement of the stations along cross-shore transects allowed for the investigation of the marsh buffer distance’s impact on mitigating roadway flooding, with buffer distances determined as the perpendicular distance from the roadway station to the estuarine shoreline. For the roadway transects considered, existing marsh buffer distances ranged from 38 m to 563 m. Ten numerical model output stations were established along each the roadway, the marsh, and the estuarine shoreline (30 total) to assess the frequency and extent of marsh and roadway flooding during storm events (Figure 1C).

3.3. Shoreline Change Analysis

Bimonthly aerial images from NCDOT taken from 26 September 2003 to 16 April 2021 were digitized in order to identify historic positions of the estuarine shoreline. Shoreline change rates were calculated as a linear regression between the shoreline position and time using the Digital Shoreline Analysis System (DSAS) developed by the United States Geological Survey (USGS) [50].

In addition to aerial datasets, local estuarine shoreline surveys were conducted periodically from 2019 to 2021 to evaluate seasonal shoreline changes and the effects of storms, including a nor’easter in November 2019 and bay-side Hurricane Isaias in August 2020 (

Table 1. Dates of shoreline surveys performed.

Survey Dates	Purpose/Event	Max Water Level at Marina Tide Gauge between Surveys (m NAVD88)	Max Significant Wave Height Waverider Buoy between Surveys (m)
13 November 2019 (pre-storm) 25 November 2019 (post-storm)	Nor’easter 16–20 November 2019	0.6 (24 November 2019)	7.6 (16 November 2019)
22 June 2020	2020 Baseline	---	---
12 August 2020 (post-storm)	Hurricane Isaias 4 August 2020	0.9 (4 August 2020)	4.4 (4 August 2020)
8 March 2021	2021 Baseline	---	---
3 May 2021	2021 Final Shoreline	0.5 (29 April 2021)	5.2 (19 March 2021)

). Trimble R10 and R12 Global Navigation Satellite System (GNSS) Global Positioning Systems (GPS) with mobile connection to continuously operating reference stations (CORS) were used to survey the scarp edge of the marsh or the location of the dense vegetation for all but the May 2021 survey. These systems have a reported horizontal accuracy of 8 mm and vertical accuracy of 15 mm. The May 2021 survey was conducted with a hand-held Trimble R1 with maximum precision of 50 mm. Effort was made to perform all surveys as close to low tide as possible. It is estimated that there was up to 30 cm of uncertainty in visual identification/interpretation of the marsh edge. In some cases, areas were inaccessible due to hazardous conditions and were not surveyed. These areas and areas that were not consistently identified as either the marsh edge or the water line (i.e., sandy portions of the estuarine shoreline) were not considered in the analysis due to differences in the measurement location not related to erosion.

The data points surveyed with the GPS system were imported into ArcGIS and connected by lines using GIS tools. DSAS was used to create a series of transects with 5 m spacing along the shorelines [50]. The distance from the shoreline to the baseline was calculated for each date to estimate the marsh edge shoreline change.

3.4. Channel Bathymetry Comparison

Two bathymetric datasets were collected in the study area as part of a collaboration made possible via the During Nearshore Event eXperiment (DuNEX), a multi-institutional research program organized by the US Coastal Research Program [51]. An initial survey was conducted on 10 October 2019 by staff from the National Science Foundation (NSF) Natural Hazards Engineering Research Infrastructure (NHERI) Rapid Response Research (RAPID) program [52]. The survey was performed using the NHERI RAPID program’s Z-Boat 1800 with a single-beam echo sounder. The Z-Boat was remotely controlled by NHERI staff and performed soundings, which were located using an onboard Digital Global Positioning System (D-GPS). A second bathymetric survey was conducted on 20 April 2021 by Woods Hole Oceanographic Institution researchers (data release, [53]). The survey was completed using a single-beam echo sounder with an onboard GPS locator mounted on a remotely driven vessel [54]. For the initial survey, vertical elevations were adjusted to NAVD88 using VDatum and validated using water level data from the marina tide gauge. For the second survey, Post-Processed Kinematic (PPK) processing methodologies were used with the CORS reference station (NCBI) located 10 km from the site to measure the vessel’s position with 3 cm to 5 cm vertical and horizontal accuracy. These surveys were used to determine the position of the deepest part of the channel and to track channel migration along the study area via changes in the 5 m NAVD88 depth contour adjacent to the estuarine shoreline. This contour was chosen because it marks the boundary of the deeper portions of the channel, where velocities are higher.

3.5. Empirical Fragility Derivation

Three sets of empirical fragility curves were developed: two sets of curves were derived to predict roadway vulnerability to flooding based on either publicly available data or numerical model outputs, and one set of fragility curves was developed to identify marsh vulnerability to erosion. Fragility curves were derived by fitting available data to the Gaussian probability distribution, consistent with previous studies for engineering applications that fit damage data to normal or lognormal probability distributions [2,55,56,57]. Fragility curves derived using Method 1 were based on publicly available data from 2017 to 2019 of daily environmental conditions from the marina tide gauge and waverider buoy and roadway closure and hazard information for NC 12 from NCDOT TIMS data. Fragility curves created using Method 2 considered numerical outputs at the study site from bay-side storm scenarios. Finally, fragility curves were created to assess marsh vulnerability (Method 3) based on assessments of the marsh condition considering shoreline surveys following storm events, long-term erosion rates, and proximity of the nearby channel. For Methods 1 and 2, we define “failure” as the occurrence of flooding on a section of roadway, leading to either a traffic closure/hazard report in the TIMS data (Method 1) or a period of flooding at a roadway numerical model output station in the numerical model simulations (Method 2). Therefore, fragility curves present the probability of flood-based roadway closure due to a vector of environmental and location variables. For Method 3, failure was defined based on empirical classification of the marsh as either “healthy”, “eroded”, or “severely eroded” (

Table 2. Objective marsh classification scheme for categorization as healthy, eroded, or severely eroded.

System Rating	Shoreline Condition	Marsh Condition	Example
Healthy	Gentle slope; plant growth on or adjacent to shoreline; minimal to no exposed root mat.	Slope less than 1:30 *; all or majority of marsh above 0.61 m (2 ft) elevation; consistent plant growth throughout marsh, including juvenile plants; intact root mat.
Eroded	Scarp less than 0.30 m (1 ft); evidence of offshore (in-water) plant growth. May include: evidence of undercutting or cracks, chunks of marsh breaking off along shoreline.	Slope between 1:30 and 1:10 *; 50% of marsh above 0.61 m (2 ft) elevation; evidence of dead or otherwise removed plants; intact or exposed root mat. May include: signs of semi-regular flooding, evidence of channel incursion or paleo inlets.
Severely Eroded	Scarp greater than 0.30 m (1 ft); evidence of offshore (in-water) plant growth. May include: visible chunks of marsh sloughed off into water.	Slope greater than 1:10 *; less than 50% of marsh above 0.61 m (2 ft) elevation; significant evidence of dead plants or no plant growth; exposed root mat or no evidence of root mat present. May include: signs of regular flooding, evidence of channel incursion or paleo inlets.

* Slope parameter as defined in [58].

). These classifications were developed based on a literature review [30,58,59] and field observations at the study area. A detailed classification of the marsh shoreline was performed by the research team in May 2021 and is used in the fragility derivations.

The fragility curves based on publicly available data (Method 1) defined failure as a roadway designated as closed or hazardous in the TIMS data. Thus, curves assessed the vulnerability of the roadway to flooding or overwash causing hazardous travel conditions affecting roadway functionality. Independent variables measured at the marina tide gauge included the maximum daily water level WL_max, peak daily 5 s wind gust, V_wind, and corresponding wind direction θ. Water levels were referenced with respect to NAVD88. These variables were considered in addition to the maximum daily significant wave height H_s,max and corresponding dominant wave period T_pd measured at the waverider buoy. While TIMS data provide the county and often nearby cities of reported incidents, the precise locations of closures along NC 12 were not able to be determined, and all flood-driven closure events (from either the ocean or the bay) were considered. Therefore, landscape variables such as marsh or beach buffer distances to the roadway were not able to be disaggregated.

For Method 2, numerical model outputs from storm scenarios allowed for the investigation of environmental conditions specifically leading to bay-side flooding. Roadway section failure was determined based on a numerical model output of whether a station was shown as flooded (failure) or remained dry over the duration of a storm scenario. Flooding was the most likely indicator of roadway closure in this situation, as the limited dimensions of the infrastructure (two-lane roadway with minimal shoulder) provide minimal opportunity to maintain traffic during inundation events. Variables considered in the numerical model-based fragility curves included significant wave height at the boundary H_s,boundary, peak water level at the marina tide gauge WL_max,marina, significant wave height at the shoreline H_s,shoreline, and peak water level at the shoreline WL_{max,shoreline}. The resolution of the numerical model outputs further allowed for consideration of the marsh buffer distance X_marsh as a potential predictor of roadway vulnerability, where X_marsh is defined as the perpendicular distance between the numerical model output station and the estuarine shoreline.

For fragility curves based on marsh conditions (Method 3), two definitions of failure for marshes were considered: one considering failure when the marsh segment was classified as “severely eroded”, and one considering failure when the marsh segment was classified as “eroded” or “severely eroded” per Table 2. These conditions are important for both the persistence of the marsh and the performance of the vegetation in shoreline stabilization and infrastructure protection. Possible predictor variables influencing marsh failure included (i) the distance from the marsh shoreline to the shoreline-adjacent 5 m depth contour in the channel X_5mcontour based on bathymetric measurements taken in either October 2019 or April 2021 (Figure 1), (ii) the slope of the channel between the 2 m contour and 5 m contour m_2mto5m, (iii) the rapid response erosion rate RR determined using linear regression of shoreline positions between November 2019 and March 2021 (Table 1), (iv) the long term erosion rate LTR measured from satellite images of the shoreline at low tide taken between 2003 and 2021, and (v) the percentage of time T_BSS>0.2 that the modeled bed shear stress exceeded a critical threshold (0.2 N/m²) near the marsh shoreline. This threshold was defined based on the median grain size d₅₀ for the study site of 205 μm [44] and the 0.2 N/m² minimum critical shear stress for sands with d₅₀ = 200 μm, as reported from flume experiments [60,61].

Backward multiple regression was used to determine fragility models for each track within each method, and variables were assessed for importance based on their statistical p-value considering the 99% percent significance level [62]. Univariate regressions testing the significance of individual variables were also considered for roadway and marsh vulnerability. For the three sets of fragility curve derivations, multiple variable “tracks” were evaluated such that only one variable representative of an independent predictor was considered in any regression model. For example, marsh failure was modeled as a function of distance to the 5 m depth contour in the channel. While distances were available from both the 2019 and 2021 bathymetry data, only one distance (from either the 2019 or the 2021 bathymetry data) was included in a given regression analysis. The resulting fragilities were then compared to determine which variable (e.g., X_{5mcontour,2019}, the distance from the shoreline to the 2019 bathymetry 5 m depth contour or X_{5mcontour,2021}, the distance from the shoreline to the 2021 bathymetry 5 m depth contour) was the more significant predictor of marsh vulnerability. The goodness of fit was assessed based on each model’s deviance and R² value, which describes the proportion of the variance in the data that is explained by the predictor variables.

4. Results

4.1. Long-and Short-Term Shoreline Changes

Figure 2 shows the long-term shoreline change rates along the estuarine shoreline as determined from the shorelines derived from aerial imagery. The most severe erosion was observed in the northern portion of the study area, with rates of over 4 m/year of shoreline recession. Along the more southern portions of the study site, the shoreline was observed to be stable to slightly accreting (less than 0.6 m/year of accretion). The shoreline classification on 3 May 2021 according to Table 2 is also shown in the figure. Areas of ongoing long-term erosion were often classified as severely eroded or eroded, and sections observed to be stable were typically classified as healthy.

The results of the estuarine shoreline surveys illustrated that changes in shoreline position over the shorter term could vary significantly depending on the characteristics of the events between surveys (Table 1). Despite this variation, evidence of marsh erosion was observed during each of the rapid response shoreline surveys (following the nor’easter in November 2019 and Hurricane Isaias in 2020), including marsh platform cracking and scarps.

During the nor’easter in November 2019, ocean-side flooding and sand overwash led to the closure of the NC 12 Highway from the Basnight Bridge to Rodanthe, NC, from 17:00 November 16 to 10:00 November 20 [63,64]. During the rapid response shoreline survey on November 25, many of the areas that exhibited signs of ongoing undercutting in the pre-storm survey had eroded further, with escarpments observed throughout most of the shoreline except within the southern pocket beach. The primary mechanism of erosion appeared to be undercutting and slumping of the marsh platform, leading to sections of marsh breaking away from the shoreline. The average marsh edge erosion during the event, excluding the southern pocket beach region, was approximately 2.5 m along the surveyed area. The recession was generally between 1 m and 3 m, with a maximum of 5.3 m measured just north of the pocket beach.

Hurricane Isaias, with a track landward of the barrier island system (Figure 1), affected the study area on 4 August 2020, and the rapid response shoreline survey was conducted on 12 August 2020. Results of the shoreline survey indicated that the highest rates of erosion were in the northern to middle sections of the study area, consistent with historical analyses. In general, there was more erosion during the November 2019 nor’easter than during Hurricane Isaias. The average shoreline recession observed during Hurricane Isaias was approximately 0.8 m, less than half that observed during the 2019 nor’easter. It is theorized that this difference in shoreline response could be because of the longer duration that the nor’easter affected the shoreline (~5 days) compared to the duration that the hurricane affected the study area (~1 day). Additionally, differences in water levels during each storm may have affected erosion rates. Hurricane Isaias’s maximum water level was approximately 0.3 m higher than that of the 2019 nor’easter (Table 1). This higher water level may have reduced the erosion by inundating the marsh and decreasing the impact of waves on the marsh edge. Likely, a combination of these factors led to the reduced shoreline recession during Hurricane Isaias compared to the 2019 nor’easter.

The 2021 survey events were widely spaced and reflected ongoing difficulties with travel during the coronavirus pandemic. Between 8 March and 3 May 2021, there were several smaller nor’easter events with a maximum water level at the marina of 0.5 m on 29 April 2021, and a maximum significant wave height at the waverider buoy of 5.2 m on 19 March 2021 (Table 1). The marsh edge change between these surveys ranged from a maximum of approximately 2 m of recession in the northern area to 3 m of advance in the southern area due to seasonal marsh growth.

4.2. Marsh and Roadway Flooding due to Storm Scenarios

Numerical model output variables included the peak water level at the marina tide gauge, the peak significant wave height and water level at each estuarine shoreline numerical model output station, and the occurrence of roadway or marsh flooding at each of the respective stations. For the range of storm scenarios considered, peak water levels at the marina ranged from 0.46 m to 2.27 m, and resulting peak significant wave heights and peak water levels at the shoreline numerical model output stations ranged from 0.01 m to 0.11 m and 0.51 m to 1.90 m, respectively. Excluding the cases driven by the lowest surge level (0.5 m), maximum water levels at the marina station were roughly 60% of the magnitude of the surge level forced at the bay-side boundary. Given the extreme bay-surge conditions imposed in the simulations, waves at the waverider buoy had minor contributions to water level, with 0.08 m variability due to waves offshore of Oregon Inlet in the cases driven by the lowest surge level (0.5 m), and decreased in contribution as the surge level inside the sound increased (0.006 m variability in maximum water level at the marina tide gauge is due to offshore waves).

Of the ten numerical model output stations established along the length of the roadway (Figure 1), stations R3, R4, R6, R7, R8, and R9 did not have any flooding for any of the storm scenarios considered. Stations R0 and R5 flooded when the marina tide gauge water levels exceeded 1.91 m NAVD88. Stations R1 and R2 were also flooded when marina water levels exceeded 2.27 m NAVD88. Numerical model output stations along the landward marsh edge suggested that all marsh stations flooded when the marina water levels exceeded 0.46 m NAVD88, except for stations S0, S1, and S2, which flooded when marina water levels exceeded 1.18 m NAVD88, 1.91 m NAVD88, and 0.50 m NAVD88, respectively.

4.3. Fragility Curves

Results of the multivariate logistic regression for publicly available data (Method 1) indicate that the maximum daily water level at the marina tide gauge and maximum significant wave height at the waverider buoy are significant predictors (p < 0.01) of roadway closure due to either ocean or bay-side flooding (

Table 3. Summary of fragility model derivation methods, definition of failure, data sources, variables considered, significant variables based on p < 0.01, and model R² values.

Method	Failure Definition	Data Sources	Variables Considered	Significant Variables (p-Value)	R2
1	Roadway closure due to flooding	Publicly available from TIMS, marina tide gauge, waverider buoy	Hs,max, WLmax, Tpd Vwind, θ	Hs,max (1.14 × 10−12) WLmax (1.17 × 10−4)	0.26
2	Roadway transect flooding	Numerical model outputs	Hs,boundary, WLmax,marina, Hs,shoreline, WLmax,shoreline, Xmarsh,	WLmarina (8.2 × 10−9) Xmarsh (2.6 × 10−6)	0.48
3	Marsh condition as severely eroded	Bathymetry data (2019, 2021), aerial shoreline imagery (2003–2021), rapid response shoreline measurements	X5mcontour,2019, X5mcontour,2021, m2mto5m, RR, LTR, TBSS>0.2.	X5mcontour,2019 (2.48 × 10−4)	0.46
				LTR (1.39 × 10−4)	0.52
	Marsh condition as eroded or severely eroded			X5mcontour,2019 (3.11 × 10−4)	0.35
				LTR (1.39 × 10−7)	0.52

). Figure 3 provides an example of fragility curves for Method 1, showing the probability of roadway closure P(f)₁ for four different water levels as a function of significant wave height. Larger significant wave heights and larger water levels cause an increased probability of failure (i.e., roadway closure due to flooding). The ranges of water levels (referenced with respect to NAVD88) and significant wave heights considered in the fragility model are −0.20 m < WL_max < 1.43 m and 0.44 m < H_s,max < 7.59 m, respectively. The R² value for the fragility model is 0.26, indicating that 26% of the variance in the data is accounted for by maximum daily significant wave height and water level at the waverider buoy and marina tide gauge, respectively.

While a model including peak wind gusts, wind direction, water levels, and significant wave heights improves the R² value compared to the model considering only water levels and significant wave heights by 0.05 (from 0.26 to 0.31), both wind speed and wind direction are not significant predictors of roadway closure due to flooding (p > 0.01), although wind direction is statistically significant at the 0.05 level. Univariate models considering wind speed or wind direction individually indicate that wind speed is a significant predictor of roadway closure due to flooding, with p < 0.01. However, wind direction is not a significant predictor, potentially because bay-side and ocean-side flooding events are not distinguished in the fragility model based on TIMS roadway closure data. These results suggest that combinations of high wind speeds and directions may be correlated with water levels at the marina tide gauge, with high marina water levels associated with sustained periods of high westerly winds, while wind direction does not predict roadway flooding without the co-occurrence of high directional wind speeds.

The low R² of the final model considering peak water levels and significant wave heights (0.26) may partially be due to differences in local drivers of flooding (e.g., the water level at the shoreline of Pea Island [44]) from what was measured at publicly available gauges. In addition, other variables not able to be determined from publicly available data (e.g., buffer distance, the elevation of the flooded roadway) are likely to contribute to roadway vulnerability. Similarly, temporal considerations may also contribute to roadway inundation, such as previous flooding or rainfall events that saturate the soil, the duration of sustained directional winds, or the duration of flooding and elevated significant wave height conditions.

Considering the results of numerical model outputs for storm scenarios (Method 2), roadway vulnerability to bay-side flooding is dependent on the peak water level at the marina tide gauge and the buffer distance between the roadway transect and shoreline (p < 0.01) for the marina water levels ranging from 0.4 m to 2.2 m NAVD88 and buffer distances ranging from 38 m to 563 m. Significant wave height at the boundary and significant wave height at the shoreline are not significant indicators of bay-side flooding. The little significance of waves in this method results from simulations only accounting for bay-side storms, where the main flooding driver is a bay-side surge that counteracts the effects of waves as they try to propagate through the neighboring Oregon Inlet. While both the water level at the marina tide gauge and water level at the shoreline are found in separate models to be significant predictors of bay-side roadway flooding, both models exhibit similar performance (R² = 0.48), even though there may be local variability in water levels at the shoreline that is not captured consistently at the marina tide gauge. The water level at the marina tide gauge is selected as the predictor variable in Figure 4 due to it being readily obtained from publicly available data.

Figure 4A shows the probability of roadway transect flooding P(f)₂ as a function of water level at the marina tide gauge for a buffer distance of 50 m. The 95% confidence intervals are shown as dashed curves, and the solid fragility curve indicates that the probability of roadway flooding increases with increased water levels at the marina tide gauge. Figure 4B shows the effect of buffer distance on roadway transect flood vulnerability, depicting fragility curves for four buffer distances from the roadway based on the water level at the marina tide gauge. The probability of roadway transect flooding increases with a decreasing marsh buffer distance: for a water level at the marina tide gauge of 2.5 m (NAVD88), the probability of roadway transect flooding increases from 0.18 for a 200 m buffer to 0.95 for a 50 m buffer. These results highlight the importance of mitigating erosion to maintain large buffer distances between the bay-side shoreline and the roadway.

For fragility curves investigating marsh failure (Method 3), the horizontal distance from the marsh shoreline to the 5 m contour in the channel (ranging from 13.2 m to 66.5 m in 2019 data and 7.8 m to 49.5 m in 2021 data), the long-term erosion rate (ranging from −2.3 m/year (accretion) to 5.6 m/year), and the slope between the 2 m and 5 m contours in the channel offshore of the marsh segment (ranging from 0.12 to 0.58) obtained from the 2021 bathymetry data are identified in univariate regression as significant variables predicting marsh failure (p < 0.01). While the distance to the 5 m contour is identified as significant in both the 2021 and 2019 bathymetry data, the fragility model using the 2019 bathymetry data provides a better description of the proportion of the variance in the data (severely eroded model, R² = 0.46, severely eroded or eroded model, R² = 0.35) compared to the fragility model using the 2021 bathymetry data (severely eroded model, R² = 0.40, severely eroded or eroded model, R² = 0.22). The better performance of the fragility model derived based on 2019 bathymetry data compared to the fragility model based on the more recent survey suggests a lag between channel proximity (determined from bathymetric measurements) and shoreline erosion on a temporal scale of several months to years.

The 2 m contour was not continuous, resulting in insufficient data to construct a fragility model using the 2019 bathymetry data. However, sufficient data were available to calculate the slope between the 2 m and 5 m contours obtained from the 2021 data, which is identified as a significant predictor of the marsh condition. The performance of the fragility model based on channel slope (severely eroded model, R² = 0.33, severely eroded or eroded model, R² = 0.37) is similar to the fragility model based on the distance to the 2019 5 m contour, slightly under-performing in its prediction of marsh classification as severely eroded and slightly improving the prediction of marsh condition as either eroded or severely eroded.

The best overall predictor variable based on statistical significance and R² value is the long-term erosion rate, determined from overhead imagery of the marsh shoreline obtained at low tide between 2003 and 2021 (severely eroded model, R² = 0.52, severely eroded or eroded model, R² = 0.52). Notably, multivariate regression considering both long term erosion rate and distance from shoreline to the 5 m depth contour produces a fragility model with slightly better performance compared to the univariate model (severely eroded model, R² = 0.54, severely eroded or eroded model, R² = 0.53), but with non-significant p-values, indicating the collinearity of these two predictor variables. Fragility curves are shown in Figure 5 for marsh classification as severely eroded (black curve with markers) or eroded/severely eroded (black curve) as a function of (A) distance to the 2019 5 m contour X_{5mcontour,2019} and (B) the long-term erosion rate LTR. Shoreline data are shown as colored symbols, using a similar color scheme classification as in Figure 2B. As indicated in Figure 2, many of the areas classified as severely eroded or eroded (red squares and orange triangles in Figure 5, respectively) are associated with locations of high long-term erosion rates in the northern to central sections of the study area. These areas, similarly, are associated with closer proximity to the channel, as indicated by the 5 m depth contour (Figure 1).

5. Discussion and Conclusions

5.1. Implications for Planning

Empirical fragility curves derived from publicly available data and storm scenario simulations indicate the importance of measurements at nearby monitoring stations in predicting roadway inundation or closure due to flooding, particularly for bay-side events. Therefore, the results of this study may contribute to risk-management programs in the area: transportation planners may identify elevated water level conditions and take precautionary action to mitigate roadway flooding or prevent unsafe travel conditions, and coastal managers may identify adaptation alternatives to improve the resilience and robustness of transportation infrastructure. The fragility curves based on numerical simulations identified marsh buffer distance as a significant predictor of bay-side roadway flooding. An increased buffer distance of 150 m (from 50 m to 200 m) for a 2.5 m water level reduces the likelihood of roadway transect flooding by over 75% for the hydrodynamic conditions considered here, highlighting the importance of a healthy marsh buffer between the shoreline and the roadway for mitigating flooding impacts.

SLR may exacerbate the vulnerability of coastal transportation infrastructure by inundating marshes and reducing the buffer distance between the shoreline and the roadway. Interactions between marshes and developed near-shore infrastructure must also be considered. While no adverse effects of the roadway on marsh erosion were observed in this study area (i.e., marsh erosion was driven by proximity to the channel and long-term erosion rates related to channel velocity and/or sediment budget), the effects of coastal squeeze by near-shore infrastructure may limit the ability of vegetation to adapt to rising sea levels [65,66]. The marsh’s condition is an essential component of roadway vulnerability (or robustness) to bay-side flooding for the ranges of hydrodynamic conditions considered here, in addition to its significance for ecosystem services, such as a habitat for migratory birds and loggerhead turtles [67] and carbon sequestration [68], which are particularly important given the location of the study area in a National Wildlife Refuge. Therefore, the results indicate the importance of marsh monitoring and adaptive management through conservation, restoration, and erosion-mitigation measures.

As identified in the marsh vulnerability analysis, long-term erosion rates and proximity to a tidal inlet flood channel affect the marsh’s classification as severely eroded or eroded. Therefore, steps may be required to reduce long-term erosion rates by relocating the channel or providing some type of shoreline edge stabilization. The rapid-response erosion rate was not identified as a significant predictor of marsh condition in the fragility model, but as observed in pre-and post-storm surveys for the 2019 nor’easter and 2020 Hurricane Isaias, episodic events caused further erosion of already-eroding areas of the marsh shoreline. These observations suggest that the pre-storm condition of the marsh may affect the severity of impact from episodic events: an already-eroding marsh is vulnerable to more erosion during storms, while a healthier marsh (with more established vegetation and a gentle shoreline slope) may be more robust and resistant to storm damage. In other words, for this study area, while the main drivers of the shoreline condition are long-term erosion rates associated with the proximity of the inlet flood channel, episodic events can exacerbate existing erosion issues. A systems approach must be used to consider marsh vulnerability and its connection to the vulnerability of coastal transportation infrastructure.

5.2. Study Considerations

While the fragility models considered here identified significant variables influencing roadway flooding or marsh condition, several idealizations and assumptions were made in the fragility models’ derivations. First, the models were derived by fitting fragility data to a Gaussian distribution and assuming that data are normally distributed. Future work may consider fitting data to alternative distributions, such as the Weibull distribution or generalized extreme value distribution [3,12]. Similarly, variables not considered in fragility model derivation due to unavailability or insufficient data may be significant contributors to roadway or marsh fragility. The relatively low R² values for all fragility models considered indicate that other variables are required to explain the variance in the fragility models or that larger sample sizes are required to improve confidence in fragility model outputs. For example, publicly available data could not identify locations of roadway flooding, and therefore effects of buffer distance or flooded roadway elevation could not be included. While roadway elevations were available in the numerical model, marsh elevations were variable between the shoreline and roadway transect, and the range of elevations for roadway numerical model output stations was small (<1.0 m). Therefore, marsh buffer distance was selected as the landscape parameter mitigating roadway flooding, although greater elevations in the marsh or roadway are expected to further mitigate flooding impacts.

Similarly, sustained periods of high directional winds likely play an important role in bay-side roadway flooding, and wave direction may be important in driving flooding from the bay or ocean. While wind speed and direction are included in the fragility model based on publicly available data, multi-variate regression indicates a correlation between hydrodynamic conditions (water levels at the marina tide gauge and significant wave heights at the waverider buoy) and wind speed and/or direction. For example, high water levels at the marina tide gauge are associated with high sustained westerly wind events. The significance of wind speed and direction variables are likely reduced in the fragility model based on publicly available roadway closure data, which does not distinguish between bay-side and ocean-side flooding-related closures. Fragility model derivations based on numerical model outputs indirectly account for wind speed and direction, considering water levels at the marina tide gauge and shoreline and significant wave heights at the shoreline and boundary. Peak directional wind speeds must be sustained for durations sufficient to generate significant fetch-generated waves and water level setup. Therefore, future work may consider explicitly including threshold wind speeds, directions, and durations in fragility investigations.

Other temporal considerations likely play a role in both roadway vulnerability to flooding and marsh vulnerability to erosion, such as the duration of elevated water levels and/or wave heights and pre-storm marsh or roadway conditions based on the timing and frequency of previous rainfall or inundation events. As indicated in the analysis of short-term shoreline change following episodic events, longer-duration, lower-intensity storms (e.g., 2019 nor’easter) may have a more significant effect on event-driven shoreline erosion (and flooding) than shorter duration, higher peak-intensity storms (e.g., 2020 Hurricane Isaias). Processes occurring on longer temporal scales, such as climate change, long-term scour, and infrastructure deterioration, as well as the occurrence of multiple hazards, should also be considered [8,9]. Future work may thus consider the importance of time-dependent predictors in roadway or marsh fragility.

The models presented here are case-specific to the study area considered and only valid for the range of water levels, wave heights, and landscape conditions observed at nearby monitoring stations or tested in the numerical model. Future work may consider process-based relationships between sediment transport and hydrodynamic forcing (e.g., current, wave, and water level effects on bed shear stress or sediment transport). The relationships identified here as significant must be verified with additional modeling or field observations for a range of sites and conditions before these features can effectively be incorporated into design guidance [58,69].

While this study identified key relationships between roadway and marsh vulnerability, environmental conditions, and landscape features, other considerations and vulnerabilities may make climate change adaptation in the area more complex. For example, the subsidence [70] of the Outer Banks contributing to relative sea level rise may create long-term challenges for adaptation in the study area, particularly considering issues of a coastal squeeze if the marsh is not able to retreat due to the presence of the roadway. In addition, implications of interventions at regional scales must be considered for both updrift and downdrift locations. Critical next steps for future work include the integration of fragility curves presented here, coupling the marsh vulnerability to the resulting disruption in roadway functionality and, ultimately, the losses suffered by affected communities. Risk-based approaches may further help to understand hazard probabilities, community exposure, and vulnerability. Finally, the implications of mitigation actions (or inactions) on populations must be considered to ensure equitable shoreline management approaches. Communication among stakeholders, engineers, and decision-makers is essential to develop equitable and robust solutions for future climate change challenges.

5.3. Study Contributions

One contribution of this study is the proposed methodology for characterizing the marsh shoreline as “healthy,” “eroded,” or “severely eroded” (Table 2). While guidance exists for shoreline assessment that provides criteria for healthy or eroding marshes based on marsh slope or scarp height [58] and other studies have evaluated shoreline condition following extreme events [30,59,71], vulnerability to erosion [72,73,74,75,76,77], and recovery after storms [71,78], a standardized engineering methodology is required for assessing shoreline condition considering the status of erosion, vegetation persistence, and other landscape factors. Therefore, Table 2 may be refined, expanded, and generalized for evaluation of the condition of marsh shorelines or other nature-based shoreline-protection alternatives.

This study further provides a framework for assessing the vulnerability of coastal transportation infrastructure to chronic and extreme bay-side flooding events as a function of hydrodynamic characteristics and interconnection with existing landscape geomorphological processes. We present observations, modeling results, and empirical fragility curves showing (i) relationships between roadway flooding, water levels, wave heights, and/or marsh buffer distances and (ii) relationships between marsh condition, long-term erosion, and proximity to an encroaching flood channel. The results may be particularly useful for emergency and transportation planners; generally, roadways are closed when flooding or debris on the roadway creates unsafe driving conditions. The fragility curves presented here allow for the identification of water levels at which roadway flooding exceeds a threshold likelihood, providing infrastructure agencies with the ability to address potentially hazardous driving conditions or to mitigate potential flooding in advance of actual impacts (e.g., dune reinforcement). Methods for reducing roadway vulnerability (e.g., by increasing the marsh buffer distance between the roadway and the bay-side shoreline) are also identified; these results suggest the potential of nature-based alternatives, such as healthy wetland systems, to bolster the resilience of coastal transportation infrastructure. While the results are specific to the study area, the curves may be tested, calibrated, and validated for other locations and predictor variables to inform future planning and flood-risk-management efforts.

As coastal communities consider adaptation pathways to manage future coastal flood hazards in the face of climate change and sea level rise, the interconnectivity between the shoreline condition and infrastructure vulnerability may play a greater role in flood risk management. Understanding these processes may allow decision-makers to best leverage nature-based and conventional infrastructure to improve the resilience of coastal transportation infrastructure and surrounding communities.