Importance of Pre-Storm Morphological Factors in Determination of Coastal Highway Vulnerability

Abstract

This work considers a database of pre-storm morphological factors and documented impacts along a coastal roadway. Impacts from seven storms, including sand overwash and pavement damage, were documented via aerial photography. Pre-storm topography was examined to parameterize the pre-storm morphological factors likely to control whether stormwater levels and waves impact the road. Two machine learning techniques, K-nearest neighbors (KNN) and ensemble of decision trees (EDT), were employed to identify the most critical pre-storm morphological factors in determining the road vulnerability, expressed as a binary variable to impact storms. Pre-processing analysis was conducted with a correlation analysis of the predictors’ data set and feature selection subroutine for the KNN classifier. The EDTs were built directly from the data set, and feature importance estimates were reported for all storm events. Both classifiers report the distances from roadway edge-of-pavement to the dune toe and ocean as the most important predictors of most storms. For storms approaching from the bayside, the width of the barrier island was the second most important factor. Other factors of importance included elevation of the dune toe, distance from the edge of pavement to the ocean shoreline, shoreline orientation (relative to predominant wave angle), and beach slope. Compared to previously reported optimization techniques, both machine learning methods improved using pre-storm morphological data to classify highway vulnerability based on storm impacts.

1. Introduction

Coastal road infrastructure represents one of the main components of cities located at or near the ocean. This infrastructure suffers inevitable damage due to natural hazards that pose not only life-threatening events to population, but also a burden to the satisfactory operation of infrastructure [1]. In the United States, there is an average general cost of 20.5 billion USD per hurricane each year at places exclusively located close to shorelines [2]. Particularly, the state of North Carolina, which has the second largest shoreline (484 km) of the US Atlantic coast, reported a value of 181 million USD, on average, over the years of 2016 and 2020, due to the presence of inclement weather events, such as hurricanes [3]. In recent years, local authorities have conducted prevention programs, including the regular measurement of the physical characteristics of coastal areas. These data-collecting efforts have guided researchers and transportation agencies to perform prevention-related activities, such as the development of vulnerability indicators for coastal roadways. The analysis of these indicators may provide a better understanding of the drivers of coastal roadway vulnerability and assist transportation departments with preventive actions in response to natural hazards.

Pre-storm coastal morphology, in partnership with storm properties (waves and storm surge), has long been assumed to control impacts on beaches, dunes, barrier islands, and the infrastructure present on those islands. Coastal morphology features, including island width, dune crest elevation, and distance from the road edge to the ocean shoreline, are indicators to assess whether road infrastructure is vulnerable to damage or not [4]. Vulnerability analysis reflects the susceptibility of the infrastructure to reduce its functionality, ranging from debris removal without interruption to temporary road closures [5]. Processes such as storm surges and sea-level rise are the main drivers of the overall damage to coastal infrastructure. Previous research works and agency reports agree that the increment in the number of storms and rising sea levels will likely affect coastal communities at different severity levels [6,7]. Furthermore, as mentioned by [8], there is a need for local coastal road systems to handle dynamic conditions. This work aims to contribute to the literature by developing a classifier model coupling past storm information with current and historical local morphology factors to inform stakeholders of the vulnerability of coastal roadways to hurricanes and tropical storms.

The monitoring of areas close to coastal roadways provides critical information about the past and current conditions of the existing road infrastructure. Especially in barrier islands, where bridges and roadways are likely the first infrastructure systems facing the ocean, monitoring indicators have provided managers valuable insights for beach restoration projects and natural habitat conservation programs [9,10]. Moreover, state agencies use field-measured data coupled with visualization from aerial photographs to monitor the morphological conditions of coastal zones, such as barrier islands, over time, especially after natural hazards events. Therefore, modeling and analyzing these dynamic data sets are the backbone of the current and future action plans to protect coastal infrastructure, such as roadways.

Recently, machine learning techniques have been explored to inform land cover classifications and storm responses [10]. A review of machine learning techniques to classify land cover is provided by [11], with the random forest algorithm having the highest accuracy level for this application. An examination of machine learning in land use and cover change detection and modeling was provided by [12]. These authors stated that machine learning has the strong potential to advance the modeling of land use and cover changes by identifying and incorporating new exploratory variables.

An application of machine learning to storm response and recovery, considering power utility, was detailed by [13]. The authors developed a data-driven predictive model to aid the utility in its emergency planning efforts. Another power utility application was presented in [14], with a focus on quantification of uncertainty in machine-learning-based prediction modeling. Other works have focused on predicting physical and social storm impacts; ref. [15] developed a convolutional-neural-network-based model to rapidly predict storm surge across an extensive coastal region using a storm track. Other studies have examined storm impacts and responses via social media data [16,17,18].

In parallel to the machine learning work, much work regarding the control of variable storm impacts has focused on damage to protective dunes [19,20]. These studies linked controls on hydrodynamic forcing to erosion of the dune. Recently, ref. [5] developed a data set of seven storms and the corresponding impacts on a coastal highway, along with a robust collection of varying morphological factors that could be used to predict impact. They used deterministic optimization techniques (i.e., mixed integer linear programming) and found that the distance from edge-of-pavement to dune toe, volume above mean high water between edge-of-pavement and ocean shoreline, distance from edge-of-pavement to ocean shoreline, and dune crest height above the road were the most skilled individual predictors of highway vulnerability. A multi-indicator function of dune toe elevation and distance from edge-of-pavement to dune toe was found to be more skilled than any of the individual indicators.

This paper aims to apply two machine learning techniques to the storm impact data set presented in previous studies [5], in order to classify roadway vulnerability and identify the indicators that best infer storm impacts. We applied two classifiers, K-nearest neighbor (KNN) and ensemble of decision trees (EDT), to classify whether different longitudinal segments or transects of a coastal roadway are considered vulnerable (class of interest) or not. Our models were applied to a 13-mile (21-km) expanse of state highway North Carolina 12 (NC 12), along Hatteras Island, NC (Figure 1). NC 12 is an essential, 148.0-mile-long (238.2 km) coastal highway that connects Corolla to Cedar Island on the northern side of the Outer Banks in the state of North Carolina, USA [4].

2. Materials and Methods

2.1. Dataset

This section describes the data set and methods used in the analysis. The data sets comprise seven groups, corresponding to seven storm events registered by the North Carolina Department of Transportation (NCDOT) between 2011 and 2019 [4]. The storms featured in this paper are summarized in

Table 1. Summary of storms and their impacts on the study area. Peak wave heights and water levels were obtained from the stations shown in Figure 1.

Storm #	Storm	Date and Duration	Storm Approach	Peak Wave Height (ft)	Max Water Level (ft NAVD88)	Storm Impacts to Study Area
1	Hurricane Irene	August 2011 50 h	Sound-side	21.4	3.0	Island breaches, erosive damage to the road, and overwash on road
2	Hurricane Sandy	October 2012 88 h	Ocean-side	18.8	4.5	Wide-spread overwash on road
3	Hurricane Arthur	July 2014 14 h	Sound-side	18.5	1.3	Erosive damage to pavement, overwash on road
4	Nor’easter	February 2016 55 h	N/A	15.7	3.5	Overwash on road
5	Hurricane Matthew	October 2016 153 h	Ocean-side	17.2	3.3	Overwash on road
6	Nor’easter	March 2018 148 h	N/A	17.5	3.7	Overwash on road
7	Hurricane Dorian	September 2019 47 h	Ocean-side	21.9	4.6	Overwash on road

. Of the seven storms in the data set, five were hurricanes, and their tracks, in relation to the study area, are shown in Figure 1. Two of the storms in the data set were winter extra-tropical storms, or “nor’easters”, which are not tracked in the same way as hurricanes. Locations of storm impacts to the study area’s highway, NC 12, were identified for each of the seven storms. Aerial imagery from the National Oceanic and Atmospheric Administration (NOAA) national geodetic survey and imagery, collected under the NCDOT Coastal Monitoring Program [4], were manually examined at the 1:1250 scale to locate road impacts along NC 12 (Table 1). These data collection efforts were described in detail in [5] and summarized in this section. A section of roadway was determined to be impacted if more than 50 percent of the roadway was covered with sand or visibly damaged. In general, post-storm imagery was available on the order of days post-storm, which precluded the inclusion of visible flooding as a storm impact, as it had already receded.

The indicators used to determine an impacted road transect are numerical values of the geomorphologic characteristics extracted from the coastal features of the area; they are listed in

Table 2. Geomorphological characteristics and units used as indicators for vulnerability classification.

Geomorphological Characteristic

(1) Island width (ft) from the estuary shoreline to the ocean shoreline (island width) (2) Dune crest elevation (ft) relative to NAVD 88 (dune crest elevation) (3) Road elevation (ft) relative to NAVD 88 (road elevation) (4) Difference in elevation (ft) between the dune crest and the road (dune crest height above road) (5) Distance from road edge-of-pavement (EOP) to ocean shoreline (ft) (EOP to ocean shore) (6) Volume above MHW between EOP and ocean shoreline (ft3/ft) (volume per length) (7) Distance from road EOP to estuarine shoreline (ft) (EOP to estuary shore)

(8) Distance from the dune toe to the ocean shoreline (ft) (beach width) (9) Beach slope found according to the USGS 2-point method [24] (beach slope) (10) Dune toe elevation (ft) relative to NAVD 88 (dune toe elevation) (11) Distance from road EOP to dune toe (ft) (EOP to dune toe) (12) The angular difference between shore-normal orientation and the weighted mean wave direction (degrees)

and illustrated in the schematic in Figure 2 [5]. Most of the geomorphological measurements provided in Table 2 are straightforward and visible in Figure 2. It is noted that the shoreline was identified as the visible wet-dry line on aerial photographs and is approximately the mean high-water line [21]. In this area, the mean high-water elevation is approximately 1.1 ft (0.34 m) above the North American Vertical Datum of 1988 (NAVD88) [22]. The “angular difference between shore-normal orientation and weighted mean wave direction”, however, is a bit complicated to illustrate. The hypothesis was that the shore’s orientation can be used as a proxy for exposure to average storm waves. Transects along the study area that are more exposed to direct storm waves are more likely to see significant erosion and overwash, which can affect the highway. To calculate the weighted mean wave direction (similar to the weighted mean wind direction used by [23]), wave data from the Oregon Inlet waverider buoy (NDBC station 8652587, Figure 1), between 2012–2020, were filtered to only include significant wave heights greater than 2 m, and wave direction was weighted by the square of significant wave height to better represent the directional distribution of wave energy. Circle statistics functions were used to ensure that angular directions were preserved. The mean wave direction was assumed to be constant across the relatively small (13-mile, 21-km) study area. The orientation of a line normal to the shore at each transect was computed at each transect along the study area. The angular difference between the wave height weighted mean wave direction and shore-normal orientation was computed on a transect basis. The measurement could vary from 0 to 90 degrees, where 0 degrees would suggest that the shoreline is most exposed to the average storm wave direction, and 90 degrees would be most sheltered from the mean storm wave impact.

The binary outputs of each data set assign 0 if the road transect was not impacted by a storm and 1 otherwise. These data sets have an imbalance problem, as the class of interest (impacted road transect or 1) happens less frequently than the not impacted road class (0). Therefore, the evaluation of the proposed methodology includes metrics of performance that reduce this implicit bias from the results, as explained below.

2.2. Correlation Analysis

A correlation analysis was performed to identify linear dependence among predictors prior to building one classification model, as explained below. The correlation matrix, where rows represent observations and columns the variables, shows the Pearson correlation coefficient calculated as shown in Equation (1):

$$\rho \left(A,B\right)=\frac{1}{N-1}{\displaystyle \sum }_{i=1}^N\left(\frac{A_i-{\mu }_A}{{\sigma }_A}\right)\left(\frac{B_i-{\mu }_B}{{\sigma }_B}\right)$$

(1)

where A and B are two random variables with N observations; µ and σ are the mean and standard deviations of A and B, respectively.

The preprocessing step filters predictors based on the correlation ($\rho$) and implements a feature selection subroutine. The Chi-square test reports low p-values when the corresponding predictor is dependent on the response variable and considered an important feature. If $\rho$ between two variables is equal to or higher than 0.95 (e.g., considered perfectly correlated), the most important feature reported by the Chi-square test remains in the data set.

2.3. Binary Classifiers

We apply two widely known classifiers methods for binary classification as the K-nearest neighbors (KNN) [25] and the ensemble of decision trees (EDT) [26]. These two classifiers were selected because KNN and EDT have been proven as classifiers that effectively handle imbalanced data sets [27,28,29]. For the KNN classifier, we filter perfectly correlated variables based on their rank correlation and select the most important predictors based on a feature selection subroutine that uses the Chi-square test [30]. This step avoids inflating the distance from correlated predictors. Decision trees do not assume relationships between features, but split data into subsamples that boost their classification performance. To overcome overfitting issues, related to creating a single decision tree, we use the EDT model that improves its performance, albeit decreasing interpretability. Additionally, these classifiers are not considered black boxes with random initial parameters, which makes them feasible to reproduce for similar potential problems.

2.3.1. K-Nearest Neighbor (KNN)

KNN is a simple supervised-learning classifier that assigns a new variable to the class with the most values from the closest k neighbors located on a search space, determined by the number of predictors [25]. The two inputs required to train and predict a KNN model are the number of neighbors (k) and type of distance used to determine the closest neighbors in a jth dimension, where j corresponds to the number of predictors (Figure 3). In this work, we performed an optimization subroutine to find the best hyperparameters of the KNN model for each storm. A Bayesian optimization algorithm [31] was applied to improve the performance of the KNN classifier model for each storm. The Bayesian optimization subroutine decision variables include the number of neighbors, distance metric, distance weight, and whether the data are standardized.

Table 3. KNN optimizable hyperparameters and search range.

Hyperparameter	Search Range
Number of neighbors	$\max \left(1,NumObservations/2\right)$
Distance metric	City block, Chebychev, Correlation, Cosine, Euclidean, Hamming, Jaccard, Mahalanobis, Minkowski (cubic), Spearman
Distance weight	Equal
	Inverse
	Squared inverse
Standardize data	True
	False

shows the search range used for the hyperparameter optimization of the KNN model [32].

KNN takes, as response, a binary variable called impact, where 1 represents the road is vulnerable to storms impacts and 0, otherwise. The model includes the most important features retrieved from the correlation analysis and feature selection subroutine as predictors. These features, as shown in Table 2, include (1) island width (ft), (2) dune crest elevation (ft), (3) road elevation (ft), (4) dune crest height above road (ft), (5) EOP to ocean shore (ft), (6) volume per length (ft³/ft), (7) EOP to estuary shore (ft), (8) beach width (ft), (9) beach slope, (10) dune toe elevation (ft), (11) EOP to dune toe (ft), and (12) shore orientation (degrees).

2.3.2. Ensemble of Decision Trees (EDT)

Decision trees are part of machine learning methods widely applied in the water quality, risk assessment, and forecasting domains to classify the risk of events, including pipe failure and water quality problems [27,28,29,33]. Decision trees graphically represent a classification or regression problem. The components of a single decision tree include a root node, internal nodes, and terminal nodes (see Figure 4). For classification models, a decision tree starts by assigning a class label to each leaf node. The non-terminal nodes, which include the root and other internal nodes, contain feature test conditions to separate records that have distinct characteristics. The root node has no incoming edges and can have multiple outgoing edges. Each internal node has exactly one incoming edge and two or more outgoing edges, and each leaf node has exactly one incoming edge and no outgoing edges [27,28,33,34]. In terms of bias and variability, a decision tree has low bias and high variance; therefore, averaging the result of many decision trees reduces the variance, while maintaining low bias. The application of ensembles enhances the performance of a single tree. Ensemble of decision trees (EDTs) group individual trees as one model to predict or fit numerical and categorical data (Figure 4). EDTs are built by comparing the “out-of-bag” (OOB) error value from all the individual decision trees and a voting combination of their results. EDTs have proved to be more robust when dealing with bias and variance, thus leading to generally a better prediction performance. However, EDTs are complex models to analyze, compared to individual decision trees [35].

EDTs start by identifying the root node. We used the well-known Gini index (G_idx) (Equation (2)) of a node as the split criterion to identify the root node of each tree in the ensemble. The G_idx with the lowest value is selected as the root node, and each tree grows until no more features remain for splitting [36]. EDTs are built using mainly two types of methods: bootstrap aggregation (bagging) and boosting. Bagging creates data set replicas using bootstrapping to incorporate single decision trees. A random selection of the observations with replacement is applied to create the bootstrap replicas. The most common application of bagging is the random forest method [36], where every tree of the ensemble randomly selects predictor variables for each split along the decision tree. Boosting methods, on the other hand, build the ensemble sequentially by using the output of the previous decision tree as input for the next one. The adaptive logistic regression (LogiBoost) algorithm, which is a variation of the adaptive boosting methods, can improve model performance when dealing with binary classification problems [37]. LogiBoost is a type of the widely applied adaptive boosting algorithm, where the objective function is to minimize the binomial deviance rather than the exponential loss, as shown in Equation (3):

$$G_{idx}=1-{\displaystyle \sum }_{i}pr{\left(i\right)}^2$$

(2)

where pr is the probability of an observation i classified into a particular class.

$${\displaystyle \sum }_{n=1}^N{w}_{n}log\left(1+\mathit{exp}\left(-2y_{n}f\left(x_n\right)\right)\right),$$

(3)

where:

• $y_n\in \left\{-1,+1\right\}$ is the true class label;
• $w_n$ are normalized weights;
• $f\left(x_n\right)\in \left(-\infty ,+\infty \right)$ is the predicted classification score calculated, as shown in Equation (4):

$$f\left(x\right)={\displaystyle \sum }_{t=1}^T{a}_t{h}_t\left(x\right),$$

(4)

where:

• $a_t=\frac{1}{2}\mathit{log}\frac{1-{\epsilon }_t}{{\epsilon }_t}$ are weights of the weak hypotheses ($h$) in the ensemble that are used to determine weighted error ($\epsilon$), as shown in Equation (5):

$${\epsilon }_t={\displaystyle \sum }_{n=1}^N{d}_n^{\left(t\right)}∎(y_n\ne h_t\left(x\right)),$$

(5)

where:

• $x_n$ is a vector of predictor values for observation n;
• $y_n$ is the true class label;
• $h_t$ is the prediction of learner with index t;
• $∎$ is the indicator function;
• $d_n^{\left(t\right)}$ is the weight of observation n at time step t.

In this study, we implemented a binary classification based on the analysis of the impacted roads with data collected from a set of storms that struck the North Carolina coast between 2011 and 2018.

2.3.3. EDT Settings

The number of trees or learners, their depth, and the overall learning rate are needed to set up an EDT model. Accuracy improves with the number of trees but also increases EDT’s complexity and computation time. Based on previous applications in water-related problems [27,28,29,33] and a series of preliminary analyses to determine the number of trees, a range between 10 and 100 trees was evaluated. A suitable option was found at 27 trees for the ensemble, as the performance did not improve with more trees. The maximum number of splits determines how deep a tree can be expanded. We analyzed a range between 10 and 100 observations, where 34 splits per tree produced the best results. Similarly, the learning rate that defined the step size of the error minimization was analyzed between 0.01 and 1.0, and a value of 0.65 was defined as suitable for the application of this model.

2.3.4. Model Training and Testing

The data sets for each storm were randomly partitioned into training (70%), validation (10%), and testing (20%) sets. Based on previous applications [33], each model was trained ten times to account for the stochasticity of the sampling method, and the average metrics of performance associated with each storm are reported as the results of the analysis. For the application of each classifier method, we optimize its hyperparameters, using the validation subset of each storm. Additionally, we compare the performance of these two classifiers using metrics that include the F1 score and precision-recall area under the curve (PR AUC) values, which are suitable for dealing with imbalanced data sets [38]. Finally, feature selection is presented to give a comprehensive analysis of the predictors with different effects on the model performance.

2.4. Metrics of Performance

The performance of classification models uses a set of commonly known metrics that are calculated using the confusion matrix, which shows the performance of a model in predicting samples within four classes (

Table 4. Confusion matrix of vulnerability. Adapted from [5].

	Classified Vulnerable	Classified Not Vulnerable
Observed vulnerable	True positives	False negatives
Observed not vulnerable	False positives	True negatives

). True positives (TP) and true negatives (TN) denote the number of positive and negative events that are correctly identified, respectively. In this research, positive events correspond to events with a value of 1, which means that the predicted and observed outputs are vulnerable, and negative events correspond to events with a value of 0, which means that predicted and observed outputs are not vulnerable. False positives (FP) denote the number of not vulnerable events incorrectly identified as vulnerable, and false negatives (FN) indicate the number of vulnerable events incorrectly identified as not vulnerable. Performance metrics, including accuracy (Equation (6)), precision (Equation (7)), recall (Equation (8)), and the F-1 score (Equation (9)), were calculated using the confusion matrix, as follows.

$$Accuracy=\frac{TP+TN}{TP+TN+FP+FN}$$

(6)

$$Precision=\frac{TP}{TP+FP}$$

(7)

$$Recall=\frac{TP}{TP+FN}$$

(8)

$$F-1Score=\frac{2\ast Precision\ast Recall}{Precision+Recall}$$

(9)

3. Results

3.1. Correlation Analysis and Feature Selection

This section reports on the correlation analysis and feature selection subroutine applied prior to the KNN implementation. Then, the performance metrics of the two machine learning methods used for the binary classification problem are discussed. The correlation analysis (Equation (1)) of the predictors shows that, for all storm events, Pearson’s correlation coefficient indicates perfect correlation among several variables. As expected, “(1) Island Width” is perfectly correlated to “(7) EOP to Estuary Shore”; “(2) Dune Crest Elevation” correlated to “(4) Dune Crest Height above Road”; and “(5) EOP correlated to Ocean Shore” to “(11) EOP to Dune Toe” (Figure 5).

Feature selection is used to determine which perfectly correlated variable remains in the predictors’ data set. Using the Chi-square test, the method reports that “Island Width” is a slightly better estimate than “EOP to Estuary Shore“ for storm events 1 to 5 (Figure 6a–e), and the opposite occurs for storms 6 and 7 (Figure 6f–g). Similarly, “Dune Crest Height above Road” is more important than “Dune Crest Elevation” for all storms, except for storm 5. Finally, “EOP to dune toe” is a better predictor than “EOP to Ocean” for all storms, except for storm 4 (Figure 6).

3.2. KNN and EDT Performances

The classifier performances (Figure 7) show that the F-1 score reported by the KNN classifier looks similar to the EDT classifier’s F-1 score across the seven storm events. However, the KNN F-1 score had a slightly higher value than EDT in all storms, except storm number 4. We performed a t-test to examine the statistical relevance of the difference between the reported F-1 scores. The following hypothesis was evaluated:

• H₀ = The pairwise difference between F-1 scores from KNN and EDT has a mean equal to zero at the 5% significance level. ${\mu }_{F-1KNN}-{\mu }_{F-1EDT}=0$
• H_A = The pairwise difference between F-1 scores from KNN and EDT has a mean not equal to zero at the 5% significance level. ${\mu }_{F-1KNN}-{\mu }_{F-1EDT}\ne 0$

The two-tailed t-rest resulted in a p-value of 0.04, indicating that the test rejects the null hypothesis (H₀) of having identical mean values, in favor of the alternate hypothesis (H_A) that the mean KNN F-1 score is significantly different than that of the EDT F-1 score at the 5% significance level.

The KNN and EDTs models’ performance, in terms of the F-1 score, shows that the models correctly classified a coastal road as vulnerable, based on the morphological indicators, at least 90% of the time. In the remaining 10%, the methods cannot classify an observation into the correct class. These high F-1 score values demonstrate the robustness of the classifiers to handle an imbalanced data set, where the class of interest (i.e., the road is vulnerable) is less frequent than the majority class (i.e., the road is not vulnerable to impact).

In terms of the area under the curve (AUC) reported by the classifiers, the EDT outperformed KNN classifier across the seven storms events (Figure 8). Both classifiers reported their maximum value during the storm event number 6, where KNN AUC represented 97% of EDT AUC. This difference increased for the lowest performance of both classifiers that occurs during storm event number 4, where KNN AUC only represented 86% of EDT AUC. Similar to the F-1 score analysis, the following hypothesis of the AUC values was evaluated:

• H₀ = The pairwise difference between AUC values from KNN and EDT has a mean equal to zero at the 5% significance level. ${\mu }_{AUCKNN}-{\mu }_{AUCEDT}=0$
• H_A = The pairwise difference between AUC values from KNN and EDT has a mean not equal to zero at the 5% significance level. ${\mu }_{AUCKNN}-{\mu }_{AUCEDT}\ne 0$

The two-tailed t-test resulted in a p-value of 0.0063, indicating that the test rejects the null hypothesis (H₀) of having identical mean values, in favor of the alternate hypothesis (H_A) that the mean KNN AUC value is significantly different than that of the EDT AUC value at the 5% significance level.

For further interpretation of the performance metrics’ values, each storm’s receiver operating characteristic (ROC) curve is presented (Figure 9). The predictability of the classifiers reports a recall or sensitivity higher than 90% at a very low FPR. Furthermore, the EDT model reports high F-1 scores and AUC values for all storms and samples selected, as explained in the “Model Training and Testing” subsection (Figure 10).

A comparison with a recent related application to the analysis of vulnerability assessment [5] was conducted to highlight the performance of the proposed two classifiers (

Table 5. Comparison of the F-1 score values of classifying highway vulnerability.

Storm	Parameter	[5]	KNN	EDT
1	F-1 Score	0.85	0.93	0.92
2		0.88	0.92	0.91
3		0.83	0.90	0.89
4		0.57	0.96	0.96
5		0.76	0.94	0.93
6		0.89	0.97	0.97
7		0.93	0.95	0.95

). KNN and EDT classifiers show an improvement in the classification of highway vulnerability, with respect to the previous application across all the storm events. For the F-1 score, KNN that reported the best performance, when compared to EDT, had an average F-1 score of 0.94, with a minimum of 0.90, corresponding to storm event 3. The average F-1 score reported by [5] was 0.82, with a minimum of 0.57, corresponding to storm event 4. Similarly, the proposed classifiers showed an improvement in the AUC values (

Table 6. Comparison of the AUC values of classifying highway vulnerability.

Storm	Parameter	[5]	KNN	EDT
1	AUC values	0.91	0.92	0.97
2		0.85	0.91	0.97
3		0.89	0.89	0.96
4		0.53	0.81	0.94
5		0.74	0.91	0.95
6		0.94	0.96	0.99
7		0.94	0.94	0.97

). For the AUC value, EDT that reported the best performance, when compared to EDT, had an average AUC value of 0.96, with a minimum of 0.94, corresponding to storm event 4. The average AUC value reported by [5] was 0.83, with a minimum value of 0.53, corresponding to storm event 4. The classifying capabilities of both the KNN and EDT methods outperformed the previous implementation [5] by an average of 10% across the storm events.

3.3. EDT Feature Importance Analysis

Different attributes contribute to the model with different magnitudes. For the EDT classifier, a feature importance analysis was conducted to determine how attributes affected the model performance for each of the seven storm events (Figure 11). The most important predictor across all the events is the “(11) Distance from road EOP to dune toe (ft)” (see Table 2). The second-best predictor depends on the evaluated storm. While “(5) Distance from road edge-of-pavement (EOP) to ocean shoreline (ft)” ranks second in storms 2 and 7, and the “(10) Dune toe elevation (ft) relative to NAVD 88” ranks second in storms 1 and 5. The rest of predictors that are part of the top three vary across the storm events and include “The angular difference between shore-normal orientation and the weighted mean wave direction (degrees)”, “(6) Volume above MHW between EOP and ocean shoreline (ft³/ft)”, “(1) Island Width (ft)”, “(8) Beach Width (ft)”, and “ (9) Beach Slope”.

4. Discussion

The results of coastal highway vulnerability studies inform decisions regarding when and where to implement projects. Examples of projects that could be informed by highway vulnerability studies include highway realignment, artificial dune construction, shoreline hardening, bridge construction, and regular maintenance. Maximizing the accuracy of vulnerability assessment is crucial for assuring the effectiveness of these projects and making the most out of limited transportation budgets.

A previous study [5] sought to evaluate geomorphological indicators and determine whether a weighted linear function of multiple indicators could improve the prediction of highway vulnerability, beyond that of individual indicators. Ref. [5] (2022) found that combining multiple indicators, namely the distance from road EOP to dune toe and dune toe elevation, did improve performance. However, it was not clear whether the attained accuracy (F1 score and PR AUC) was limited by the assumed function form, and there were open questions regarding whether more complex function form or machine learning methods could improve the performance.

This present study supported the hypothesis that machine learning methods would improve the classification of highway vulnerability as impacted or not. This study also affirmed the finding of [5], i.e., that distance from EOP to dune toe is a powerful and reliable predictor of highway vulnerability. Distance from EOP to dune toe was found, in this study, to be the most important predictor among most of seven storms of varying intensity and approach, relative to the study area. Some of the results seemed counterintuitive, such as the finding that dune crest elevation is not important for certain storms (Figure 6 and Figure 10). However, this area is highly managed, as discussed in [4], with frequent dune rebuilding. This may lead to dunes with relatively high crest elevations being quite close to the shoreline, where the dune is essentially a high, narrow barrier between storm waves and the roadway. In this situation, a high dune crest does not substantially lessen vulnerability.

Other factors, such as the angular difference between the shore-normal orientation and weighted mean wave direction, are shown to be important for some storms but not others. With regard to this angular difference, the relative angle between waves and structures has been shown to lead to complex effects by [39,40,41,42,43]. The research performed by [39] examined low-crested structures under an oblique wave attack. It was found that, for smooth (asphalt) structures, there was a strong dependency between the transmission coefficient and angle of wave attack, with the transmission decreasing with increasing incident wave angles. The stability of placed block revetments under oblique wave attack was investigated by [40], with the results showing that waves induce a pressure difference across the cover layer that is more complicated than that of perpendicular wave attack. Ref. [41] examined an oblique wave attack on open granular (sand) filters underneath armored slopes via physical model tests, showing that, for larger wave angles (where 0 degrees is perpendicular), the amount of filter erosion decreases. The effect of wave structure angles for tsunami loads on bridges was examined by [42,43], with the findings including the result that complex three-dimensional effects are generated by the interaction of the waves and structure. Ref. [44] examined the wave-structure interactions on bridge decks with varying geometries, also finding that loading is complex. Although these studies were not focused on ground-based roadways, similar complexity could be contributing to the results, showing the angular difference to be important for specific storms, while not as influential for others, and presenting an opportunity for further study.

While this research investigated 12 geomorphological factors that were applicable to the study area, there are other parameters cited in the literature that may be useful indicators for vulnerability assessments of roadways in other locations. A previous work [14] found that the width of elevated coastal decks affects the hydrodynamic loads imparted by waves on the structure, as well as that the ratio of “wavelength-to-deck width” is a critical factor. In the region of the present study, the roadway width did not significantly vary along the study area; therefore, this factor could not be examined. Similarly, dune vegetation was not included in this study. Past work [45,46,47,48] has shown that vegetation impacts wave-induced pressures, loads, and erosion on dunes and other structures. The fact that vegetation on dunes enhances their resiliency and dunes generally serve as protection for roadways on their landward side suggests that vegetation could have a significant effect on the vulnerability of coastal roadways. This factor was not included in the vulnerability classification models of this current study because it was difficult to adequately quantify vegetation density along a transect cross-section. In is noted that, for this specific study area, the sections that had more dense vegetation between the road and shoreline, which also had a large distance from roadway EOP to dune toe and/or roadway EOP to shoreline—this would lead to a high correlation of these variables. These two factors (roadway width and dune vegetation) could be crucial indicators in classifying roadway vulnerability for other study areas, where there may be substantial spatial variation not captured by other indicators. The study area extent was a limitation for the current research, which hindered our ability to more comprehensively assess additional vulnerability indicators. Future research could add data sets from other barrier islands, where additional variables, such as vegetation and roadway properties, could be included in the model to potentially improve the vulnerability assessment.

It is noted that the most important feature for both classifiers used in this study, distance from EOP to dune toe, can also be measured without digital computation, using aerial photographs. The development of the presented study allows researchers and practitioners to classify road vulnerability as impacted or not, based on local geomorphological observations. The models’ performance outperforms previous works significantly and can be applied to different locations where these features are available. Distance from EOP to ocean shoreline has been previously employed to evaluate highway vulnerability in this study area [4], as have long-term erosion rates and deterministic numerical modeling [49]. However, the results of this study indicate that coastal managers may be able to employ data-driven models to evaluate variability in coastal highway vulnerability on barrier islands. As remote sensing technology improves, and frequently updated topography becomes available [50,51], having a data-driven model that will classify portions of the roadway as vulnerable, without the high computational cost of deterministic numerical modeling, will be beneficial. The success of EDT in this study implies that this machine-learning model may be a useful tool for highway vulnerability studies. The model could be run with the most up-to-date geomorphology data available, in order to provide predictions regarding barrier island highway locations that are most susceptible to storm impacts. The application of the models used in this study are limited to the vulnerability classification of land-lying roadways located on wave-dominated barrier islands with sandy coastlines. Additionally, the models may not perform well for study regions where vegetated dunes are not well-correlated to distance from EOP to dune toe because vegetation was not included as a distinct variable within the model. The long-term effects of sea-level rise may bias the performance of the models, unless elevation datums are adjusted over time.