Toward an Integrated Probabilistic Coastal Vulnerability Assessment: A Novel Copula‐Based Vulnerability Index

Estimating the exposure of the coastal systems to natural hazards using coastal vulnerability models, which benefits from index‐based approaches and utilize information about the characteristics of the system, has become extensively adopted in the past few decades in coastal management and planning. However, the explanatory power of index‐based approaches and subjective selection of vulnerability factors are still in dispute. This study aims to introduce a stochastic coastal vulnerability model and assess its skill in characterizing and preserving simultaneous information about various comprising factors. Two common coastal vulnerability indices, Additive Coastal Vulnerability Index and Multiplicative Coastal Vulnerability Index are formed, and then their performances are compared to the proposed Probabilistic Coastal Vulnerability Index (PCVI) for the coastal counties of South Carolina. PCVI is developed based on the joint‐probability analysis of vulnerability factors using copula functions, which makes it capable of preserving the importance of multivariate information, and in turn, forms a more informative index. The performance of indices is benchmarked against post‐hazard flood maps and the cost of fatalities from Hurricane Florence (2018) and Hurricane Matthew (2016). The PCVI revealed more accurate results in terms of explaining the importance of vulnerability associated with biophysical and socio‐economic factors. The capability of PCVI to preserve multivariate vulnerability information offers a more pragmatic approach to reflect the exposure and adaptive capacity of coastal communities facing coastal hazards.

2 of 25 identify vulnerable communities that are susceptible and vulnerable to coastal hazards (Bevacqua et al., 2018). The predominant relationship between vulnerability caused by natural hazard exposures and the socio-economic conditions of the coastal communities, pressing the need for an integrated coastal vulnerability assessment. Such assessment needs a broader focus beyond just biophysical factors that can incorporate the spatial-temporal dynamics of socio-demographic, socio-ecological conditions (Sajjad & Chan, 2019), natural habitats of the coastal area (Sajjad et al., 2018) alongside the bio-geophysical exposures (Zhang et al., 2021) to understand the nexus between people and the coastal environment. The biophysical aspect can be described by the geology, geographical features, climate, and ecology of a coastal region, while socio-economic activities, such as land use change, population growth, socio-demographic status, and economic activities, can determine the sensitivity to hazards (Goharian et al., 2017;Wheater & Evans, 2009). Throughout ancient and even modern urbanization, the adaptive capacity of coasts has been undermined (Goh, 2019). The adaptive capacity of the coastal systems primarily depends on the socio-economic characteristics of the communities that reside in these regions (Berkes, 2007). To render the adaptive capacity as a vulnerability component, a multi-dimensional and integrated vulnerability framework that combines a set of biophysical, environmental, and socio-economic factors is required, beyond the conventional approach (Furlan et al., 2021). Appropriate consideration must be given to the biophysical and socio-economic aspects of vulnerability mapping to set up adaptation policies (IPCC, 1991).

Coastal Vulnerability Index
Index-based approaches for assessing coastal vulnerability, such as Coastal Vulnerability Index (CVI), provide vulnerability measures in a quantitative manner. Alternatively, there have been other approaches, such as Geographic Information System (GIS)-based decision support or dynamic model outcomes that include flood hazard models and continuous remote sensing observation to represent the coastal regions' vulnerability characteristics (Rajasree & Deo, 2020). The index-based vulnerability assessment approach continues to be widely used, as its application and visualization are easy, and it considers a variety of vulnerability-related indicators for the final assessment. Moreover, a CVI index can be assessed based on hierarchical aggregation or by developing an aggregated index construction. There are two approaches to form a CVI, deterministic or probabilistic. In the deterministic approach, two common CVI calculation methods are the additive vulnerability index calculation (Mclaughlin & Cooper, 2010), which is formed based on the arithmetic mean of consisting factors, and the multiplicative vulnerability index (El-Zein & Tonmoy, 2015), which is developed based on the geometric mean method. In the deterministic approach, various combinations of biophysical and socio-economic vulnerability indices, based on allotted weights, can result in CVIs (Choi, 2019). This result may not be an appropriate outcome to capture the multivariate nature of the vulnerability, while the probabilistic approach is able to capture the uncertainty or randomness associated with weight calculation (Engström et al., 2020) or by using joint probability distributions of factors (Hao & AghaKouchak, 2013). Thus, basic steps for combining the socio-economic and biophysical factors to determine the CVI for a holistic coastal vulnerability assessment involves the identification of the most accountable factors to indicate vulnerability, qualitative and quantitative contextualization of vulnerability indicators, weighting, factor normalization, and the final aggregation of indicators into one unique index (Choi, 2019). A variety of studies attempted to feature different factors, develop normalization techniques to assemble the vulnerability factors, and present the final index for CVI calculation (e.g., Bukvic et al., 2020;Kantamaneni et al., 2019;Nguyen et al., 2016). These methods are often different in their data processing, selection of representative grid resolution, quantifying factors, representing spatial and temporal (dynamic or static) scales, weighting methods, and aggregation framework.
Recent bilateral discussions have been shaped around how the nexus of socio-economic and biophysical factors should be managed in coastal systems, as described in , Smith and Rodríguez-Labajos (2021), and Almenar et al. (2021). Among the attempts pursued by Furlan et al. (2021) and Mclaughlin and Cooper (2010) conceptualized a CVI by the transformation of different dimensions of coastal vulnerability as sub-index (Mclaughlin & Cooper, 2010), for example, Social Vulnerability Sub-index, Biophysical Vulnerability Sub-index, and so on, and the CVI is developed by integration of resulting sub-indices to support the diverse characteristics of vulnerability as a multi-dimensional index. In practice, however, a more meaningful assessment measure is needed, which can reflect the combined and integrated effects of multiple dimensions of coastal vulnerability. Accordingly, this research aimed to address this gap by developing CVI for integrated coastal vulnerability analysis using a probabilistic CVI. A stochastic CVI can measure the simultaneous importance of the degree of vulnerability from multiple dimensions of vulnerability, and then combines them with CVI to inform managers (Yan et al., 2017).

Probabilistic Coastal Vulnerability Index
The interrelationship between coastal vulnerability factors can be formulated using the joint probability concept of copula functions. The copula concept is widely adopted for multivariate data analysis to model the joint probability between two or more random variables. Copula approach uses the marginal distribution of multiple variables to obtain their multivariate and joint probability distribution. The resultant multivariate index from the joint probability distribution will then be formed based on bivariate relationships. The joint probability concept has been applied for different purposes in the field of water resources engineering, such as drought monitoring (Hao & AghaKouchak, 2013), risk analysis of compound flooding (Wahl et al., 2015), vulnerability assessment of water supply systems (Goharian et al., 2018), hurricane hazard (Song et al., 2020), flood hazard (Didier et al., 2019), natural hazards (Mazdiyasni et al., 2017), sediment transport (Shojaeezadeh et al., 2018), and so on.
It is prominent that the growing body of literature on coastal vulnerability is looking more at the coping capacity of communities or the population's ability to mobilize assets or entitlements to deal with loss or to prevent future damage from hazards. This coping capacity of a community can be incorporated into the socio-economic factors, which reflects the adaptive capacity of local communities. Current methods have improved the identi fication technique of vulnerable coastal places by applying GIS, while the socio-demographic factors related to coastal vulnerability are increasingly getting importance as coastal vulnerability factors (Cutter & Emrich, 2017;Gupta et al., 2020). To prioritize the coping capacity of a coastal community, vulnerability can be determined as a multi-dimensional index using probabilistic and deterministic CVI considering biophysical vulnerability factors and socio-economic vulnerability factors. The CVI formulation using a deterministic approach continues to have inherent challenges in spatial mapping because any unaccounted socio-economic or biophysical characteristics may change the CVI outcome. On the other hand, the stochastic CVI can more reliably deal with uncertainty in factor quantification, that can be useful for long-term coastal planning and management. Therefore, existing multivariate vulnerability approaches using additive CVI (e.g., Gayen & Villalta, 2022;Mclaughlin & Cooper, 2010;Mullick et al., 2019) and multiplicative CVI (e.g., Gornitz, 1991;Mohd et al., 2019) require a comparative evaluation with probabilistic CVI.
Yet, few studies (Goharian et al., 2018;Hao & AghaKouchak, 2013) have explored the need for stochastic indices for system vulnerability analysis, the knowledge base needed to develop the stochastic CVIs to deal with coastal vulnerability analysis, integrating multiple vulnerability dimensions of socio-environmental systems in coastal vulnerability assessment still lacks a suitable framework (Kantamaneni et al., 2019;Nguyen et al., 2016). The ability of CVI to jointly represent the vulnerable coastal populations and biophysically vulnerable locations is rarely viewed in the context of both spatial contexts and coastal problems. In a recent study by , they improved the factors quantification approach for the integrated coastal vulnerability assessment in the South Carolina (SC) coast using a deterministic CVI which has a limitation on preserving the importance of simultaneous information associated with multivariate vulnerability factors. The main hypothesis of the current study is that the stochastic method-driven copula-based CVI can be more enhanced and informative for integrated coastal vulnerability analysis. Accordingly, two main research questions are explored in the current study to test the hypothesis: 1. How does the significance of the vulnerability measured by an index vary from deterministic CVI to stochastic CVI, and which metric is best suited for integrated coastal vulnerability analysis? 2. Between the deterministic and stochastic CVI, which one can better preserve the importance of multivariate vulnerability information?
Accordingly, there are two main limitations in the current state of knowledge on coastal vulnerability mapping that will be addressed in this research. The first limitation is determining the appropriate factor aggregation technique to represent the nexus between biophysical drivers and socio-economic components. A deterministic and probabilistic CVI framework is developed to address the problem. Another limitation is that the explanatory power of CVIs for interpreting disaster outcomes remains disputed, often in the context of complex socio-environmental systems. In this regard, the explanatory power of CVI and the dependence structure of biophysical and socio-economic factors can be better explored by comparing the Probabilistic Coastal Vulnerability index (PCVI) with the Additive Coastal Vulnerability Index (ACVI) and Multiplicative Coastal Vulnerability Index (MCVI).
To test the explanatory power of CVIs to see if these indices are informative and capable of keeping the importance of information from concurrent factors, two proxy data sets-the historical flood inventory and the cost of fatalities data associated with the coastal hazard-are used for evaluation purposes.

Methodology
Vulnerability in the current study refers to the degree to which a socio-environmental system is susceptible to natural hazards and the lack of coping capacity with coastal hazards or the adaptive capacities of coastal communities (Bevacqua et al., 2018;Sahin & Mohamed, 2014;Zhang et al., 2021). The coastal vulnerability is defined in the present study as a function of biophysical exposure, the sensitivity of the coastal systems, and coastal communities' lack of adaptive capacities. The vulnerability assessment follows an integrated vulnerability concept where vulnerability arises from a set of conditions and processes resulting from both biophysical and socio-economic factors. The coastal hazards, including hurricanes, coastal floods, storm surges, high winds, and SLR, are the focus of the biophysical exposures in the vulnerability assessment. The proposed vulnerability analysis was conducted for the South Carolina (SC) coast of the USA. There are two main groups of coastal vulnerability factors-biophysical and socio-economic vulnerability factors. The associated coastal vulnerability sub-indices for these groups are the Bio-Physical Coastal Vulnerability sub-index (BPVI) and the Socio-Economic Vulnerability sub-index (SEVI) (Mclaughlin & Cooper, 2010). BPVI looks for hydroclimate characteristics of coastal hazards and physical attributes of coastal areas that incite loss due to a coastal hazard. SEVI seeks to identify the populations susceptible to coastal hazards depending on their socio-demographic and economic values.
The selection of the vulnerability factors associated with biophysical and socio-economic vulnerability groups is based on the coastal vulnerability study of the SC coast done by .  developed a comprehensive and multi-dimensional CVI, combining physical, hydroclimate, socio-economic, ecological, and shoreline characteristics of SC, which consisted of approximately 21 vulnerability indicators. In the present study, to estimate the vulnerability of SC's coast, the biophysical vulnerability factors were chosen from the hydroclimate and physical vulnerability indicators. The factors associated with the biophysical vulnerability group are the number of historical coastal hazard events, hurricane track density, surge height, rainfall intensity, SLR, land use and land cover, available water storage in soil, and elevation ( Table 1). The socio-economic vulnerability factors consist of Social Vulnerability Index (SVI), the number of Historical and Archeological Structures (NHAS), and the cost of fatalities (Table 1). The socio-demographic characteristics of the SC coast are incorporated using the factor, SVI. In the current study, the SVI is a composite socio-demographic factor that is developed using the Probabilistic Principal Component Analysis (PPCA) method as described in . PPCA is a multivariate data analysis method that combines a probability model and a classical linearized projection for performing the principal component analysis (Khajehei et al., 2020;Tipping & Bishop, 1999). To derive the SVI, the socio-demographic factors provided by Cutter and Emrich (2017) for the SC coast are engaged. The socio-demographic data sets of Cutter and Emrich (2017) contained 29 socio-demographic variables attributed at the census block level which are originally derived from the US census data.
The economic dimension is further represented in the socio-economic vulnerability group adding two other factors, NHAS and the historical cost of fatalities. The cost of fatalities is chosen because damage to property and infrastructure, as well as the disruption of transportation systems, can have far-reaching economic consequences, both inside and outside the flood zone. Besides that, the SC coast is rich in historical and archeological sites. The SC coast has a total of approximately 8,749 NHASs. These archeological and historical sites provide significant socio-economic contributions to the SC economy since they receive lots of visitors throughout the world. The presence of cultural heritage works as an asset and contributes to a high level of socio-economic benefits. However, the SC coast has many NHAS that are susceptible to a high degree of coastal hazard vulnerability.
For developing CVI, a methodological framework is developed that consists of four major stages: (a) factors quantification, (b) factors aggregation to obtain sub-index, (c) CVI development, and (d) CVI evaluation. The methodological framework is shown in Figure 1. The stagewise CVI model set up is described as follows: 1. Factors quantification: The chosen factors (Table 1) are quantified by using a spatiotemporal data assimilation framework (Figure 1). In this study, the spatial data sets of the considered factors are processed using ArcGIS for the coastal counties of SC. Depending on an individual factor's contribution to overall vulnerability, relevant spatial analysis is carried out to quantify a factor. The spatial data sets are acquired for processing the spatial factors. Three different types of spatiotemporal data features are available for data processing, that is, the point, polygon, and raster. Since the vulnerability is quantified at the census block level, the polygon features are set as the final spatial factor quantification unit. Therefore, the data aggregation from raster and point features is attributed to the polygons to make the factor quantification process consistent. The number of coastal hazard events and cost of fatalities retrieved from NOAA storm event databases is retrieved as a point feature. Later, the spatial data sets are aggregated into a polygon (census block level) to quantify the spatial attribute by summing up all points in the corresponding polygon as total counts of coastal hazard events. Following this quantification approach, the polygon feature represents historical total counts of coastal hazard events. Further, the socio-demographic variables are available at the census block level, which is remained as it is because the factors quantification is done at the census block level. However, in cases of processing the raster data set, inconsistency is present while adjusting the spatial resolution to a polygon. For example, the gridded raster data set of rainfall intensity has inconsistency in the spatial resolution of raster data sets with the polygon. The inconsistency in temporal resolution of the raster data sets is quantified by 7-year of historical records of rainfall intensity. In such cases, the raster is handled using a spatial analysis tool in ArcGIS called "zonal statistics." The tool estimates the spatial and temporal average of the rainfall intensity at the census block level. To get an idea of the spatial distribution of rainfall the magnitudes of rainfall intensity are temporally averaged over the whole period. The gridded raster data sets, land use, and available water storage of the soil also have spatial inconsistency with the polygon. For quantifying the land use factor, a weighted Curve The spatial resolution is given according to raw input data sets. b Representative unit to quantify the degree of vulnerability of the census block. Number (CN) on each polygon is applied based on land use types (Table 1). The available water storage of the soil in each census block is quantified using the "zonal statistic" tool by summing up the gridded raster values of each soil grid to the corresponding census block. A detailed description of factor quantification is described in Table 1. 2. Factors aggregation to sub-index: After quantifying all spatial factors of a group the next stage is to aggregate all spatial factors using a spatial model to obtain the sub-index of each vulnerability group ( Figure 1). The main component in the spatial model developed in this stage consists with three main steps: (a) factor normalization using fuzzy functions, (b) determination of the spatial importance of each factor using objective weighting techniques as described in Section 2.1, and (c) aggregation of all fuzzy normalized spatial factors to a sub-index (BPVI and SEVI) using the objective weighting technique. The main advantages of using the fuzzy logic-based factor normalization approach are to determine the degree of vulnerability using a non-linear data transformation function, the ability to reduce the uncertainty of discrete scale factor normalization, and precise indexing during the factor aggregation by employing the vulnerability functions of both deterministic and stochastic . Details of the fuzzy functions are provided in . Fuzzy functions transform a vulnerability factor in a normalized range between 0 and 1 at a continuous scale to quantify the degree of coastal vulnerability. Fuzzy logic reflects the probability of a spatial entity being vulnerable assigning a membership function that results in a continuous scale, 0 to 1 rather than 0 or 1 of Boolean logic (Zadeh, 1975). Three fuzzy membership functions: fuzzy-large, fuzzy-small, and fuzzy-linear, are applied to the quantified factors for fuzzy normalization. Then, the normalized factors are aggregated using the objective weighting method which is adopted from Shannon's entropy concept (Shannon, 1948). When the higher input values of the quantified factors are more likely to be a member of the set, the fuzzy-large transformation function is employed, and when the lower input values of the quantified factors are more likely to be a member of the set, the fuzzy-small transformation function is used (Bahrani et al., 2016;Mullick et al., 2019). When a linear relationship between the quantified factors and coastal hazard vulnerability exists, the fuzzy-linear function is applied to perform a linear transformation of the quantified factors. 3. Index development: After aggregating all spatial factors associated with each vulnerability group to two different sub-indices: BPVI and SEVI, the next task is CVI development. The sub-indices obtained in the previous step are normalized using min-max standardizations into 0 to 1. Then, three vulnerability functions: additive, multiplicative, and copula are used to combine the sub-indices into PCVI, ACVI, and MCVI, respectively. The stochastic copula function is described in Section 2.3 and the other deterministic vulnerability functions are described in Section 2.2. After getting all CVIs, an intercomparison among the maps using the developed indices is illustrated. The relative degree of vulnerability in the CVI map was classified into five classes: low, moderate low, moderate, moderate high, and high. In the data specific classification case such as SVI, this study applied the natural break methods to find the break value of the vulnerability classes. For the intercomparison of the CVIs, an equal interval method is used. To accomplish this, first, the PCVI, ACVI, and MCVI are normalized following min-max standardization, then to compare the spatial pattern of different types of CVI, based on the histogram an equal interval of 0.2 (bin width). 4. CVI evaluation: The disaster outcomes in terms of flood maps and cost of fatalities of two hurricane events, respectively Hurricane Florence (2018) and Hurricane Matthew (2016) are used as a proxy to evaluate the CVI. For the evaluation of CVI, these event specific data sets were separated in the factor quantification stage for validation purposes. An accuracy assessment using a confusion matrix was used to evaluate the CVI. The development of the confusion matrix is described in Section 2.4.

Objective Weighting Method
Subjective weighting is frequently connected with biases coming from human perception and inadequate support of facts and evidence when selecting, prioritizing, and evaluating criteria., objective weight determination approaches are likewise reliable in the absence of specific preferences and expert or decision-maker opinions because they are based on mathematical models. Following the entropy weighting technique of the objective weighting method, initially, the normalized performance index ( ) is calculated, using the fuzzy normalized membership probability ( ) at each census block using Equation 1: where r is the fuzzy normalized membership probability of a census block in a range between 0 and 1 for the census block i (1 to m) and the spatial factor j (1 to n); is the normalized performance index of a census block i and the spatial factor j. The entropy value, , of the normalized performance matrix is calculated using Equation 2 as follows: Then, the entropic weight of each spatial factor j in a vulnerability group is calculated using Equation 3: (3) where = 1 − and ∑ =1 = 1 . The term dj denotes the vector of diversification degrees which reflects the difference of jth index value. The larger the dj value, the larger the weight of the j spatial factor in a vulnerability group. Finally, the sub-index of the corresponding vulnerability group is obtained using Equation 4 and following a weighted sum function.
where is the entropic weight of a spatial factor j; and is the fuzzy normalized membership probability of a census block and vulnerability factor j; and is the resulting sub-index obtained for a census block and vulnerability group k. In the current study, the resulting sub-index in the objective weighting method for the biophysical and socio-economic vulnerability groups are BPVI and SEVI, respectively.

Additive and Multivariate Factor Aggregation
The additive utility function (Equation 5) is the arithmetic mean, and the multiplicative utility function (Equation 6) is the geometric mean of the consisting sub-indices. The sub-index was obtained from Equation 4, in a factor aggregation from the objective weighting method. In this study, two sub-indices from the vulnerability groups, that is, BPVI and SEVI, were combined following these two methods.
where ACVI and MCVI are the Additive and Multiplicative Coastal Vulnerability Index, respectively.
is the weight (0-1) assigned to the corresponding vulnerability group, and k is the number of vulnerability groups. is the sub-index obtained for the k vulnerability group at census block i.
The sensitivity of the assigned weight is also tested for the CVI evaluation. The weights assigned in Equations 5 and 6 are changed in a range between 0 and 1 with an interval of 0.1 to perform the sensitivity analysis. The multiplicative vulnerability function (Equation 6) is not sensitive to changes in sub-index weights. In a multiplicative vulnerability function (Equation 6) MCVI is proportional to a constant term, "mth root" of the products of all weights, √ 1 × 2 × 3. . . . . . . . . . −1 × . When the weights of the sub-index, of the vulnerability group is changed, the resulting index values also changed proportionally to the constant term and after normalization, the associated function returns to the same index magnitude. This can be inferred that regardless of how the weights of the sub-index: ACVI or BPVI are changing in the multiplicative vulnerability function; it shows no sensitivity to the changing weights on resulting MCVI values. Therefore, the sensitivity analysis is just performed for the ACVI and is not required for the MCVI. The advantage of the PCVI is that sub-indices weighting is not required in the copula-based stochastic method.

Sklar's Theorem and Copula Functions
Sklar's theorem (Nelsen, 2006) shows that a function C, represented here by a copula function ( , for a pair of two random variables exists where their marginal cumulative distribution functions, and for all ( ) 2 , express the joint probability function, F BS , as follows: where, the B and S variables in Equation 8 represent biophysical and socio-economic vulnerability indices, respectively. The marginal cumulative distribution is found using the copula function, which can be developed as Equation 8: where G links the bivariate distribution to marginal distributions of variables B and S, and and are the marginals of those variables.
1. The BPVI and SEVI were calculated for each census block. A pairwise data set (BPVI, SEVI) for the copula analysis has been developed using 226 pairs of BPVI and SEVI where each pair preserves corresponding BPVI and SEVI information of each census block. It should also be mentioned that the study area contains a total of 226 census blocks. 2. The best marginal probability distribution fit for BPVI and SEVI was found. The MvCAT tool estimates the parameters of the 17 different continuous marginal distributions using the maximum likelihood algorithm. Akaike information criterion (AIC) and Bayesian information criterion (BIC) criteria were used to find the best-fitted distributions. 3. The best-fit distribution for each variable was then checked with the Anderson-Darling goodness of fit test to determine the statistical significance of the fitted distribution within a 5% significance level. 4. MvCAT is used to test the performance of 25 different copula families (Gaussian, t, Clayton, Frank, Gumbel, Independence, Ali-Mikhail-Haq, Joe, Farlie-Gumbel-Morgenstern, Plackett, Cuadras-Auge, Raftery, Shih-Louis, Linear-Spearman, Cubic, Burr, Nelsen, Galambos, Marshall-Olkin, Fischer-Hinzmann, Roch-Alegre, Fisher-Kock, BB1, BB5, and Tawn) to find the appropriate function that described the correlation structure. The performance of the various copula functions was evaluated in MvCAT using multiple goodness of fit measures, such as AIC, BIC, and Root Mean Square Error (RMSE). 5. Based on the statistical measures for the goodness of fit, the best copula function was chosen. As a result, the Roch-Alegre copula was determined to be the best fit among the other copula families (Table S1 in Supporting Information S1). The mathematical formulation of the probability density function (PDF) for the Roch-Alegre copula is presented in Equation 9.
where 1 ∈ (0, ∞) and 2 ∈ (1, ∞) are the copula parameters, and u and v are the optimally fit bivariate marginal distributions. 6. The optimal sets of best fit copula model parameters were derived using MvCAT tool, where the optimization problems can be solved by local optimization and Bayesian inference-based Markov Chain Monte Carlo methods. At this point, we used the local optimization method to find the copula model parameters in the Roch-Alegre copula function (Equation 9). The local optimization method utilizes a gradient-based, interior-point search method (Waltz et al., 2006) to solve the optimization problem by repeating the search for 30 iterations from different random initial points to reduce the chance of being trapped in a local minimum (Sadegh et al., 2017). 7. When the copula model parameters were set, the PDF and CDF of Roch-Alegre copula functions were used to simulate 100,000 sets of random points in a range between [0,1] to represent the dependence structures of BPVI and SEVI metrics.

Probabilistic Coastal Vulnerability Index
Finally, following Goharian et al. (2018), the PCVI was formulated with Equation 1 using the CDF. This function was used for the joint CDF, and to estimate the vulnerability surface, the PCVI was applied using the inverse function of a standard normal distribution ( −1 ). Probabilistic index has been successfully used before when conducting water supply system vulnerability assessments (Goharian et al., 2018), agricultural drought monitoring (Hao & AghaKouchak, 2013;Liu et al., 2016), meteorological drought assessment (Masud et al., 2015), power production risk and water stress (Ganguli et al., 2017), and compound drought and hot extremes (Hao et al., 2020).
where the PCVI denotes the Probabilistic Coastal Vulnerability Index, ( ) represents the CDF of BPVI and SEVI.

CVI Evaluation
To evaluate and validate all CVIs, the reference data for index validation was derived from flood maps and the cost of fatalities considering two hurricanes, that is, Hurricane Florence (2018) and Hurricane Matthew (2016). For this purpose, NASA-Moderate Resolution Imaging Spectroradiometer products from the Dartmouth flood observatory were retrieved for processing post disaster flood images. The images for Hurricane Florence (2018) and Hurricane Matthew (2016), were processed at 500 m resolution. These satellite images helped to validate the coastal hazard maps and the biophysical aspects of the CVI. ArcGIS was used to project retrieved inundation maps and to transform them into flood intensity maps. The vulnerability of each census block was defined based on the percentage of flood exposures per unit area per event in the census block. Larger inundated areas in each census block resulted in higher degrees of vulnerability in that block and vice versa.
Since the validation of the indices is more targeted to disaster outcomes, here the coastal flood magnitude and frequency are engaged to investigate the CVI inform ability. The footprint of flood extent is related to the biophysical nature of vulnerability. Besides that, the evaluation of an integrated CVI in the sense of a multivariate vulnerability index is incomplete in the absence of socio-economic validation aspect. Therefore, the cost of fatalities data sets during Hurricane Florence (2018) and Hurricane Matthew (2016) are used to quantify the socio-economic aspects of the CVI. The cost of fatalities is obtained from the NOAA Storm event database.
Natural Breaks methods based on Jenk's optimization (Chen et al., 2013) are applied to the reference data sets and CVIs (PCVI, ACVI, and MCVI) to get a binary classified output: vulnerable and non-vulnerable classes for each of these cases. The Natural Breaks algorithm minimizes the deviation within classes by maximizing the standard deviation between the classes (Chen et al., 2013). The binary classified indices are pairwise compared with the binary classified reference data sets for an instance of flood data sets and another instance of cost of fatalities. The confusion matrices for each index are estimated to evaluate the accuracy of the indices using the values of precision (Equation 11), recall (Equation 12), F1-score (Equation 13), and accuracy (Equation 14) (Han et al., 2011). The confusion matrix depends on four performance measures: true positive (TP), true negative (TN), false positive (FP), and false negative (FN), which were obtained previously from a pairwise comparison between the binary classified reference data sets and the indices. When a census block is recognized as vulnerable by both reference data sets and indices, then the census block is identified as a TP sample. When a census block is identified as non-vulnerable in both reference data sets and indices, then the census block is discerned as a TN sample. The FP sample is obtained in an instance when a census block is predicted as vulnerable by the indices, although it is identified in the non-vulnerable class according to the reference data sets. When a census block is predicted as non-vulnerable by the classified indices but according to the reference data sets it is identified in a vulnerable class, then the sample is called an FN sample. By summing up the total number of TP, TN, FP, and FN samples, the confusion matrix for an index can be calculated. The precision was used to determine the positive predictive value, which is the proportion of the total correctly identified by the CVI map and matched with the reference flood hazard or the cost of fatalities. Recall, also known as sensitivity, is the proportion of actual positive predictions (TP) that are identified correctly in the vulnerability maps. The F1 score is the weighted average of precision and recall, and it can be obtained using Equation 13.
where TP, FP, TN, and FN denote true positive, false positive, true negative, and false negative samples, respectively, in the validation data set. The values of precision, recall, and F1 score ranged from 0 to 1. A precision value of 1 indicates that the CVI predicted vulnerability is correct with respect to the reference hazard vulnerability.
Low precision values indicate the prediction of vulnerability. A recall value of 1 indicates that the CVI map correctly predicted the biophysical and socio-economic vulnerability, while lower values indicate an underprediction of vulnerability. The perfect precision, F1, and recall scores occur when they have magnitudes of 1.

Study Area
The study area (Figure 2) covered six coastal counties including Jasper, Beaufort, Colleton, Charleston, Georgetown, and Horry, along with 226 census blocks of SC, which have a combined area of 15,113 km 2 . The coastal region is gently sloping toward the Northern Atlantic Ocean. The coastal areas have numerous riverine estuaries, lagoons, salt marshes, and wildlife refuges. The majority of these coastal features are bordered by barrier islands and have been shaped over the last 15,000 years as a result of SLR. The barrier islands and salt marshes on the SC coast play a crucial role in maintaining a favorable environment for 291 ecological species . The Winyah Bay is the largest estuary in Northern SC, fed by Black Rivers, blackwater Waccamaw, and the PeeDee River (Mallin et al., 2000). The North and South Santee River discharges freshwater into the Winyah Bay from three coastal watersheds, Santee, Saluda, and Broad watersheds (Figure 2). Most port and harbor activities in SC occur in the Charleston Harbor area. The Ashley, Cooper, Waccamaw, and Santee are some of the major rivers in SC coasts ( Figure 2). The Edisto River, another major river in SC, passes through the southwestern part of Charleston County and is one of the longest, free-flowing blackwater River in US (Mallin et al., 2000). Further south, the St. Helena Sound and the Broad River Estuary enrich the SC coast with numerous tidal creeks, salt marshes, and islands.
Approximately 28% of the 5.02 million residents in South Carolina live in eight coastal counties, and this population has increased by 127% from 1970 to 2010, which placed SC's coast as the third highest populated area among the 31 coastal and Great Lakes states nationwide (Crossett et al., 2013). According to the US Census Bureau (2018), an increase of 0.62 million in population occurred between 1989 and 2018. Meanwhile, the SC coast is experiencing an accelerated SLR in recent decades, 0.25-2.0 m of SLR by 2100 on the SC coast (Daniels, 1992). In 2016 and 2017, the tidal flood days that occurred in the city of Charleston were 50 and 46 days, respectively (Hollman, 2020). Long-term, historic tide gauge records at Charleston Harbor, SC shows that the Charleston peninsula has an SLR of 3.32 mm/yr and the recorded SLR for the City of Beaufort is 2.75 mm/yr at the Fort Pulaski tide gauge (Morris & Renken, 2020;NOAA, 2019). It is expected that the annual days of nuisance (high tide) flooding, as well as the coastal inundation area, will increase in the future with SLR, which, in turn, emphasizes the need for a more precise and comprehensive coastal vulnerability assessment for SC's coast.

South Carolina's Coastal Vulnerability Factors
The factors considered in the BPVI, and SEVI are shown in Table 1 with the respective sources of data and their descriptions. The assigned fuzzy functions for each factor are also described in  . The cost of fatalities represents economic damage to property, crops, and casualties during disaster events. It is estimated from NOAA storm event databases using the available disaster records between 1996 and 2020.

Socio-Economic and Biophysical Vulnerabilities
About 13 principal components are retained in the PPCA analysis ( Figure S1 in Supporting Information S1). About 76% cumulative explained variance is achieved from the principal components ( Figure S1 in Supporting Information S1). The SVI map developed using the PPCA analysis is shown in Figure 3. A histogram (bin width = 0.2) is also provided in Figure 3 of SVI that shows the total counts of census blocks under different ranges. The SVI map shows the vulnerable census blocks' location is mostly in Beaufort, Charleston, and Horry County (Figure 3). The census blocks adjacent to the ocean demonstrate more social vulnerability than those inland. Hence, coastal hazards have a significant influence on the social vulnerability of the coastal areas.
The results from the estimated SEVI (Figure 4a) show that most of the populated areas in the study area, such as the City of Beaufort, the City of Charleston, the City of Georgetown, and the City of Myrtle Beach (Figure 4a), have low to moderate socio-economic vulnerability. Apparently, technological and structural integrity that exists in these communities, including early flood warning systems, evacuation facilities, transportation routes, emergency health services during coastal hazards, and long-term disaster management planning, have increased the adaptive capacity to coastal hazards in those coastal urban areas. However, between the metropolitan and the rural areas of SC, a large gap remains in the adaptive ability of coastal hazards. Besides the socio-economic factors, including the level of education, public health, per capita income, socioeconomic status, household composition, disability, minority status, people aged older than 65 or below 17, percentage of people below poverty estimate, unemployment rate, the proportion of disable population in a community, percentage of mobile homes, and high percentage of people with no or minimum education, are among the factors that contribute to the high SEVI in SC's coastal counties. The SEVI can be used as a measure of the SC coastal system's adaptive capacity against coastal hazards. According to SEVI estimates, 1.8%, 69.2%, 24.55%, 1.78%, and 2.67% of the SC coast are categorized as low, moderately low, moderate, moderately high, and highly vulnerable, respectively. The SEVI map also identified many disadvantaged communities which are moderate to highly vulnerable.
The BPVI map in Figure 4b shows the biophysical vulnerability of the coastal system. A total of 4,532 coastal hazard events during 1996-2019 have been reported by NWS storm event databases for the study area. Among these hazard events, 69.9% flooded, 22.34% were hurricanes, and 4.25% were high tide events. Looking at the past, Hurricane Florence (2018) itself caused $2.2 million of economic damage in this region. Currently, 84% of businesses in Charleston are exposed to areas prone to storm surges, along with approximately 71% of business facilities that are in the coastal and inland floodplains (CRC, 2017). Charleston county is mostly classified on the BPVI map as moderate to moderately highly vulnerable areas (Figure 4b). About 70% of houses in the City of Charleston are in the 100-year return period floodplains and 87% of houses overlap with storm surge-prone areas (CRC, 2017). The frequency of coastal floods in the City of Beaufort is increasing, with a 0.5-m SLR, resulting in an estimated 15% land loss in the city (CRC, 2017). The city has experienced 13 coastal flood events since 1980, with estimated surge height records of 0.8-1.52 m above Mean Higher High Water. In the City of Beaufort, 12%-14% of businesses are exposed to coastal flooding (Knapp et al., 2019). Prior analysis indicated that the total business sales volume in this city can face about $155.7 million in economic damages with 2 m of SLR (Knapp et al., 2019). However, it is not simply businesses and people living in urban areas who are in danger. In general, the coastal environment plays a significant role in the ecological balance of SC's coast.  Charleston County, which hosts the largest wildlife sanctuary on the SC coast, Cape Romain National Wildlife Refuge (NWR), plays a significant role in preserving the ecological balance. By 2100, an SLR of approximately 0.53 m will have submerged 51.4% of Cape Romain NWR, leaving the habitat of 52 ecological species in a highly vulnerable condition and four species on the brink of extinction (Daniels, 1992). Moreover, 80% of historical sites of the state are in Charleston County, which indicates that future climate change effects should be taken into serious consideration and appropriate actions should be taken to increase the adaptive capacity of the historical sites. Beaufort County (average BPVI = 0.79) has the highest biophysical vulnerability followed by Charleston County (average BPVI = 0.75). The BPVI map shows that the area adjacent to the shoreline is associated with moderately high to highly vulnerable due to wide exposure to the hurricane, SLR, and their proximity to the ocean as well. Particularly, the northeastern parts of Charleston County are associated with high vulnerability (Figure 4b). According to hurricane records 16 deadly hurricanes from 1851 to 2017 made landfall in Charleston County and 9 out of 16 hurricanes made landfall in the northeastern part of Charleston County. The BPVI map also shows that the urban areas on the SC coast, such as the City of Charleston, the City of Beaufort, and the City of Myrtle beach are recognized as highly vulnerable. The City of Charleston is associated with moderate to high vulnerability based on BPVI, while the whole peninsula is depicted as highly vulnerable (Figure 4b). The U.S. Army Corps of Engineers has already proposed a $1.8 billion coastal management plan for the Charleston peninsula to increase the adaptive capacity of the City by developing new structural measures, such as storm surge barriers, breakwater or wave attenuation structures, deployable floodwall, and levees, along with elevating the roads, bridges, and buildings (USACE, 2020). Colleton County (average BPVI = 0.36) appears to be the least vulnerable area based on BPVI; however, three Category 4 hurricanes that is, Hazel (1954), Gracie (1959), and Hugo (1989) made landfall in Horry County at Myrtle Beach, SC.

Joint Probability
The BPVI and SEVI are both fitted to generalized extreme value distributions (Figures 5a and 6a). Figures 5a  and 6a show the fitted probability distributions for BPVI and SEVI. The quantile-quantile (q-q) plot ( Figures 5  and 6b) compares the empirical probability with the fitted probability. The fitted marginal distribution of BPVI in q-q plots is perfectly matched with its empirical probability. However, the q-q plot shows a slight deviation of the empirical SEVI from its fitted probability function. The result of the goodness of fit test measures, AIC, BIC, and RMSE (see Table S1 in Supporting Information S1), demonstrates that among the 25 copula families investigated in this study, the Roch-Alegre copula performed the best.

PCVI
The correlation between BPVI and SEVI is determined using Spearman's rank-order and Pearson correlation, which are obtained, respectively, −0.33 and −0.34. These results suggest that a significant correlation between BPVI and SEVI exists to derive PCVI. The CDF of the Roch-Alegre copula is provided in a color-coded contour, as shown in Figure 7 also illustrates the relationship between the BPVI and SEVI to obtain the PCVI. The vulnerability surface plot in Figure 8 shows how non-normalized PCVI magnitudes vary in a bivariate relationship with BPVI and SEVI. The PCVI-based vulnerability interpretation is rather non-linear than the conventional additive or multiplicative vulnerability index estimations. The PCVI concept is developed by paying attention to the co-existence of vulnerability associated with the BPVI and SEVI, where their simultaneous information is incorporated using a copula function. One major advantage of adopting the copula approach is to get an idea of the relationship between biophysical and socio-economic vulnerability using the CDF (Figure 7) and the vulnerability surface of PCVI (Figure 8). Lower BPVI up to marginal probability 0.2 and higher SEVI above marginal probability 0.8, in combinations result in higher CDF and vice versa. For example, the marginal probability values of SEVI range between 0.8 and 1 and for BPVI between 0.1 and 0.2, which result in higher CDF (Figure 7). Marginal probability values of BPVI are estimated as 0.8-1 and for SEVI where it ranges between 0.12 and 2, also results in higher CDF (normalized value 0.9-1) (Figure 7). The CDF in Figure 7 approaches toward higher values with marginal values of BPVI greater than 0.7 and marginal values of SEVI greater than 0.65. This result suggests a combination of higher BPVI values (marginal probability above 0.7) and higher SEVI values (marginal probability above 0.65) form a higher PCVI. Figure 8 illustrates the inverse normal of CDF into PCVI. The transformation of the PCVI surface, transformed from a flatter surface (BPVI-SEVI marginal probability combination in moderate-moderate range) to a steeper surface (BPVI-SEVI marginal probability combination in high-high range). The lower probability of SEVI and BPVI produced low vulnerability scores and high-high combinations of these variables result in high vulnerability scores. The question of which approach is the more accurate vulnerability index may arise at this point; however, this question will be resolved in Section 4.5, where we examine the validity of all aggregation methods in this study. For instance, a coastal region may experience severe damage due to hurricane hazards, and if there is no socio-economic exposure in this area, the system would not be considered as vulnerable based on the PCVI. However, such an instance using ACVI and MCVI would result in moderate or high vulnerability. Therefore, identifying vulnerable areas using the PCVI index depends on the co-existence of socio-economic vulnerability and coastal hazards, along with the concurrent occurrences of a hazard in a coastal region. Using joint probability distributions, if either bivariate factor (BPVI or SEVI) has zero probability, regardless of the other factor's value, the system may represent zero or low vulnerability, as the PCVI value will be low ( Figure 8). In conclusion, the PCVI provides a great opportunity for a more informative CVI by including simultaneous probability information taken from different sets of vulnerability indices (BPVI and SEVI) and factors.

Vulnerability Maps
The coastal vulnerability maps obtained from three distinct aggregation frameworks are shown in Figure 9. Figure 9 shows the coastal hazards vulnerability maps developed using ACVI_0.5 (the numeric value 0.5 represents BPVI's weight) (Figure 9a), MCVI (Figure 9b), and PCVI (Figure 9c). The pattern of the three maps has a significant spatial distinction, which is especially prominent after comparing the deterministic CVI with stochastic CVI. The histogram in Figure 9d shows the percentage of census blocks classified under different vulnerability classes using three indices. The total census blocks classified under the moderate vulnerability class in ACVI_0.5 and MCVI maps are 61.95% and 67.26%, respectively, while the census blocks classified under the same vulnerability class in PCVI maps are 32.30%. In contrast, under the low and moderate low vulnerability classes in the case of the ACVI_0.5 and MCVI maps have fewer census blocks than the PCVI maps. The PCVI maps show 12.83% of census blocks are classified as low vulnerable, while the ACVI_0.5 and MCVI maps illustrate, respectively, 0.44% and 0.44% under the same vulnerability class. This result indicates that accounting for the importance of simultaneous vulnerability information in PCVI provides a lower degree of vulnerability than the other indices  (ACVI and MCVI) developed using a deterministic approach. Charleston County is the most vulnerable county in SC, regardless of the aggregation process, followed by Jasper County (Table S2 in Supporting Information S1). The statistical information of the census blocks under five vulnerability classes is described in Table S2 in Supporting Information S1. High vulnerability in the northeastern part of Charleston County is caused by the frequent occurrence of coastal hazard events and high vulnerability associated with the SEVI (Figures 4a and 4b).
The ACVI and MCVI show moderate and moderate low vulnerability in Horry County, while it has a notable shift to low and moderate low vulnerability based on the PCVI (Figure 9c).

CVI Evaluation
In this study, the CVI was evaluated with a composite index, which was compiled from bio-physical and socio-economic vulnerability components, by selecting three types of aggregation schemes. It is necessary to validate the outcomes of each aggregation approach because the explanatory power of CVI of different aggregation schemes is disputed by their spatial pattern. Therefore, a combined real-time flood map (Figure 10) of Hurricane Florence (2018) and Hurricane Matthew (2016) and NOAA cost of fatalities information (Figure 11) during these events were used to validate the CVI map. Approximately 158 census blocks were inundated during the two hurricanes. The percentage of the area flooded per unit area per event is determined in a range between 0.001% and 0.22% in those flooded census blocks. The NOAA reported total cost of fatalities during the hurricanes was about $35.37 million in the study area. The Natural Breaks method-based threshold of 0.21% is obtained for the binary classification of the reference flood data sets. The cost of fatalities was found in a range between $0.01 and $4.32 million throughout the census block ( Figure 11). The Natural Breaks method-based threshold of $0.52 million is obtained for binary classification of the cost of fatalities. Figure 12 shows the confusion matrices of the CVI indices developed in the current study. The results of the sensitivity analysis of the ACVI are shown in Figure 12a for the flood data sets and Figure 12b for the cost of fatalities. The result of the confusion matrices for PCVI provides 0.71 accuracy for flood data sets and 0.7 for the cost of fatalities. MCVI provides accuracy matrices 0.63 and 0.69 for flood data sets and cost of fatalities, respectively ( Figure 12). The numeric parts after the underscore in the ACVI label in Figure 12 indicate the weight of BPVI for ACVI development. For example, ACVI_0.1 refers to the CVI obtained using the BPVI's weight = 0.1 and SEVI's weight = 0.9 in Equation 5. Changing the weights of the biophysical vulnerability group in a range between 0.1 and 0.9 results accuracy matrices between 0.63 and 0.69 for the validation of flood data sets where the lowest accuracy is obtained for ACVI_0.1 and the highest accuracy is captured by the ACVI_0.4. From Figure 12b, when the weight of BPVI is decreased in ACVI_0.1, which also means that increasing the importance of socio-economic vulnerability factors, the accuracy of the ACVI_0.1 to capture information of coastal flood vulnerability data sets and cost of fatalities is enhanced. However, for the ACVI, the sensitivity analysis does not provide a distinct variation in the accuracy of presenting flood data sets. Still, comparing the explanatory power of disaster outcomes with ACVI regardless of changing weights, PCVI depicts a higher accuracy in both cases, which can lead to an informative and enhanced characteristic of the PCVI method. The sensitivity analysis suggests that instead of spending time picking the best set of weights for additive, either by using objective or subjective methods, one can use the probabilistic information for CVI development. The sensitivity analysis also suggests that characteristics of copula to generate joint probabilities can inform simultaneously about both sub-index: BPVI and SEVI feature, without compromising their importance. In other words, while evaluating the confusion matrix of the CVI indices, the PCVI method may be preferred because it balanced the multivariate nature of the socio-economic and biophysical vulnerability with reasonable accuracy in integrated coastal vulnerability analysis.
The reference data sets used in this study for CVI evaluation have event-specific disaster outcomes. When the vulnerability data set was compiled, it was gathered using historical long-term flood records to predict spatial vulnerability. As a future climate change impact, SLR-induced inundation may increase biophysical vulnerability. Hence, the false positive samples in the CVI evaluation may correspond to vulnerabilities in the near scenarios.
The observed data set was evaluated to determine the existence of the simultaneous vulnerability, in terms of flood hazard and the cost of coastal hazard. The data sets showed evidence of 12 census blocks that were jointly vulnerable in terms of flood and cost of coastal hazard vulnerability. The PCVI identified 10 census blocks correctly out of 12 observed census blocks, while the ACVI and the MCVI accurately identified 8 and 9 census blocks respectively. Therefore, the joint vulnerability was determined based on the observed data sets and was accurately depicted by the PCVI method, followed by the MCVI and ACVI methods.

Discussion
CVI as an indicator and index-based measure for integrated coastal vulnerability able to identify potential vulnerable places and communities of the coastal area. CVI map has become a widespread solution for coastal management and planning for identifying vulnerable places. CVI models are gaining significant importance for hazard mitigation and management planning, but their explanatory powers remain disputed in terms of interpreting the disaster outcomes for ranking and prioritizing the vulnerable places in coastal areas. The local governments, coastal planners, and floodplain managers will get benefit by identifying the vulnerable location of the SC coast and being able to re-evaluate existing hazard mitigation and recovery plans from the developed CVI maps in the current study. The approach adopted in the current study can also be followed for any other type of natural hazard vulnerability assessment.
CVI has become advantageous for summarizing and visualizing the vulnerability of the socio-environmental systems of the coastal communities; however, the robustness and accuracy of these indices remain unclear. In search of an informed and enhanced index, in the current study a convergent validation among the alternative CVI models is carried out where PCVI shows promising results balancing the importance of vulnerability information from both biophysical and socio-economic vulnerability groups. The ACVI method is commonly adopted in multi-criteria decision making and the CVI obtained in such an approach has an uncertainty in choosing weights. The MCVI is not sensitive to changing the weight assigned to the sub-index. However, in absence of sufficient vulnerability characterizing factors and information, expert judgment and opinions play roles in getting subjective weights for CVI development. The objective weighting approach is useful when a good amount of vulnerability characterizing factors are available.
According to the CVI development framework, the major stages are factors quantification and aggregation, and the CVI development and evaluation. In the previous study,  outlined factors of quantification and aggregation for fine-scale vulnerability assessment. In the current study, the CVI development and evaluation for integrated vulnerability analysis is extended at the census block level and enhanced using the earlier progress made by . The approach applied in factor quantification includes both ways of aggregation and disaggregation: fine and large scale can be used for US national-level CVI development. Overall, the CVI can benefit the integrated coastal vulnerability studies in other regions by using the development framework and evaluation approach.
There is a substantial spatial variation in the degree of vulnerability in the CVI maps, which arises from using different CVI development approaches. Spatial inconsistency in CVI also results in underestimating and overestimating the degree of vulnerability. The uncertainty of the factor's quantification stage may underestimate the degree of coastal vulnerability. Some of the sources of uncertainty may be associated with insufficient hydroclimate data sets to characterize the biophysical vulnerability, the spatial interpolation techniques for estimating the missing values during the factor quantification, and so on (Anandhi & Kannan, 2018).
The vulnerability function in the deterministic approach keeps the simultaneous information of the sub-index but underestimates the importance of the sub-index's vulnerability. The stochastic approach in copula has several advantages over a deterministic approach which are establishing the multivariate relationships among the characterizing factors, calculating the importance of simultaneous vulnerability of the sub-indices, and then using the joint probability to get the PCVI. If any unforeseen vulnerability characterizing events or instances goes beyond the factor quantification in the deterministic approach, it is not possible to reflect the unquantified degree of vulnerability through the deterministic CVI. In this regard, the PCVI is recommended here since the bivariate relationship is established considering a stochastic joint probability approach and uncertainty analysis. Therefore, PCVI can be adopted here for long-term coastal management while the ACVI and MCVI frameworks should be frequently updated with any unforeseen vulnerability characteristics.
Copula-based approach also reveals the bivariate relationship between the biophysical and socio-economic factors using the CDF plot and vulnerability surface plot. When the number of sub-indices increases, vine copula for multivariate copula analysis can be used for CVI development. Besides that, using two events: hurricanes Florence (2018) and Mathews (2016) for the validation of CVIs is not adequate. In future studies, proper consideration of a series of pre-storm CVI and post-disaster outcomes can be used. The spatial regression analysis also can be tested for CVI evaluation (Rufat et al., 2019) with more rational measures of disaster outcomes to see if the CVIs are spatially associated with disaster outcomes or not. In case the indices are not associated, an improvement of the factor selection process is also recommended.

Conclusions and Future Works
This study evaluated three methods of aggregation of coastal vulnerability assessment by developing the ACVI, MCVI, and PCVI, which account for spatially distinctive results and considers both biophysical and socio-economic factors using the SC coast as a case study. The major goal of vulnerability assessments is to give insight that will enhance the adaptation process and improve the adaptive capacity of the coast. As a result, a vulnerability index seeks to reduce several complex and interacting parameters, which are represented by joint socio-economic and biophysical vulnerability factors of a coast, to a form that is more meaningful, informative, enhanced, and useful as a management tool. A coastal management plan with vulnerability analysis can strengthen early disaster recovery from a natural hazard. In this study, the PCVI was developed using joint probability analysis, which searched for a wide range of copula functions and local optimization to reduce uncertainty in copula parameters. The Roch-Alegre copula was selected based on the goodness of fit measures. Other aggregation methods, such as ACVI and MCVI, did not have the same explanatory power of vulnerability as PCVI. The confusion matrix in the vulnerability evaluation showed that the PCVI achieved F1 scores of 0.42 and 0.54, with accuracy of 0.71 and 0.7, respectively, for the reference maps of flood and socio-economic vulnerability assessment. The importance of concurrent socio-economic and biophysical vulnerability information was well preserved in the copula-based PCVI. MCVI demonstrated F1 scores of 0.45 and 0.43 and an accuracy matrix of 0.63 and 0.69, respectively, for flood hazard and socio-economic vulnerability assessment data sets, respectively. The ACVI method showed good agreement with coastal flood vulnerability evaluation. However, ACVI showed a low F1 score and accuracy matrix for validation with the cost of fatalities. All CVI maps showed that the northeastern part of Charleston County is highly vulnerable because the biophysical vulnerability drivers of hazard and socio-economic conditions of this area are located at the upper margin of vulnerability. This coastal vulnerability study can provide decision support information for coastal infrastructure development, including ports and harbors management, treatment plants, infrastructural development to withstand the growing population, and economic activities of the coastal area.

Data Availability Statement
All data, codes, and developed shapefiles and analyzed geospatial data that are used for estimating PCVI in the study is publicly available, shared in a public database, and are licensed under a Creative Commons Attribution 4.0 International license .