Hurricane Flood Hazard Assessment for the Archipelago of San Andres, Providencia and Santa Catalina, Colombia

Despite the low occurrence of tropical cyclones at the archipelago of San Andres, Providencia, and Santa Catalina (Colombia), Hurricane Iota in 2020 made evident the area vulnerability to tropical cyclones as major hazards by obliterating 56.4 % of housing, partially destroying the remaining houses in Providencia. We investigated the hurricane storm surge inundation in the archipelago by forcing hydrodynamic models with synthetic tropical cyclones and hypothetical hurricanes. The storm surge from synthetic events allowed identifying the strongest surges using the probability distribution, enabling the generation of hurricane storm surge flood maps for 100 and 500 year return periods. This analysis suggested that the east of San Andres and Providencia are the more likely areas to be flooded from hurricanes storm surges. The hypothetical events were used to force the hydrodynamic model to create worst-case flood scenario maps, useful for contingency and development planning. Additionally, Hurricane Iota flood levels were calculated using 2D and 1D models. The 2D model included storm surge (SS), SS with astronomical tides (AT), and SS with AT and wave setup (WS), resulting in a total flooded area (percentage related to Providencia’s total area) of 67.05 ha (3.25 %), 65.23 ha (3.16 %), and 76.68 ha (3.68%), respectively. While Hurricane Iota occurred during low tide, the WS contributed 14.93 % (11.45 ha) of the total flooded area in Providencia. The 1D approximation showed that during the storm peak in the eastern of the island, the contribution of AT, SS, and wave runup to the maximum sea water level was −3.01%, 46.36%, and 56.55 %, respectively. This finding provides evidence of the water level underestimation in insular environments when modeling SS without wave contributions. The maximum SS derived from Iota was 1.25 m at the east of Providencia, which according to this study has an associated return period of 3,234 years. The methodology proposed in this study can be applied to other coastal zones and may include the effect of climate change on hurricane storm surges and sea-level rise. Results from this study are useful for emergency managers, government, coastal communities, and policymakers as civil protection measures.

Despite the low occurrence of tropical cyclones at the archipelago of San Andres, Providencia, and Santa Catalina (Colombia), Hurricane Iota in 2020 made evident the area vulnerability to tropical cyclones as major hazards by obliterating 56.4 % of housing, partially destroying the remaining houses in Providencia. We investigated the hurricane storm surge inundation in the archipelago by forcing hydrodynamic models with synthetic tropical cyclones and hypothetical hurricanes. The storm surge from synthetic events allowed identifying the strongest surges using the probability distribution, enabling the generation of hurricane storm surge flood maps for 100 and 500 year return periods. This analysis suggested that the east of San Andres and Providencia are the more likely areas to be flooded from hurricanes storm surges. The hypothetical events were used to force the hydrodynamic model to create worst-case flood scenario maps, useful for contingency and development planning. Additionally, Hurricane Iota flood levels were calculated using 2D and 1D models. The 2D model included storm surge (SS), SS with astronomical tides (AT), and SS with AT and wave setup (WS), resulting in a total flooded area (percentage related to Providencia's total area) of 67.05 ha (3.25 %), 65.23 ha (3.16 %), and 76.68 ha (3.68%), respectively. While Hurricane Iota occurred during low tide, the WS contributed 14.93 % (11.45 ha) of the total flooded area in Providencia. The 1D approximation showed that during the storm peak in the eastern of the island, the contribution of AT, SS, and wave runup to the maximum sea water level was −3.01%, 46.36%, and 56.55 %, respectively. This finding provides evidence of the water level underestimation in insular environments when modeling SS without wave contributions. The maximum SS derived from Iota was INTRODUCTION Tropical Cyclones (TC) present a major hazard for countries in tropical areas (Ortiz-Royero, 2012;Martell-Dubois et al., Considering the exposure of the SPSC to TC, the archipelago is exposed to hurricane SS flooding, and assessing the inundation hazard is relevant for decision making and planning. Despite the relevance for assessing flood levels, there are only a few historical events to perform a robust statistical description of the TC climate in this area. This limitation has been solved regionally utilizing synthetic TC (Emanuel et al., 2006(Emanuel et al., , 2008 enabling robust statistical analysis to characterize the TC climatology, or using hypothetical TC (Zachry et al., 2015;Rey et al., 2019a) to assess flood worst-case scenarios. In this study, we followed the method from Ruiz-Salcines et al. (2021), using synthetic events (Emanuel et al., 2006(Emanuel et al., , 2008 to develop a model-based hurricane SS flood hazard using a hydrodynamic model. As in previous studies (Lin et al., 2010;Meza-Padilla et al., 2015;Appendini et al., 2019;Marsooli and Lin, 2020), this flood hazard assessment methodology does not rely on historical data on TC or surges, but rather generates synthetic storms that are derived based on the physics of TC (Emanuel et al., 2006(Emanuel et al., , 2008, enabling the assessment of SS flood and their return periods in areas with low TC occurrence, as in the archipelago of SPSC. As such, a large number of synthetic TC were generated over the Caribbean Sea, to use those passing 250 km from the middle point between the San Andres and Providencia islands, and used to force a hydrodynamic model to determine hurricane SS and their associated flood in the archipelago. Additionally, given the scarcity of major hurricanes in the synthetic hurricane dataset, we created hypothetical category V hurricane events for each island to assess the flood worst-case scenario (Zachry et al., 2015;Rey et al., 2019a;Ruiz-Salcines et al., 2021). Considering the above, this study aims to assess the inundation threat from TC for the archipelago of SPSC by means of a hydrodynamic model, which was forced with synthetic tropical cyclones and hypothetical hurricanes, the former to determine the probability for different hurricane SS levels and the latter to provide hurricane SS flood worst-case scenarios. Although the wave contribution to surges in insular environments is significant (Chen et al., 2017), in this study has been neglected due to computational constraints to simulate hundreds of TC events. However, the contribution of the astronomical tide (AT), SS, the wave setup (WS), and wave runup to the maximum seawater level during the pass of hurricane Iota in 2020 for Providencia and Santa Catalina is investigated to know the underestimation of seawater level when waves and tides are not considered. The description of materials and methods as well as description of the study area is introduced in section "Materials and Methods." The results and discussions are conducted in section FIGURE 1 | Study area. (A) Topographic elevation and bathymetric depth, computational domain (unstructured mesh), Hurricane Iota track (gray color). Ground elevation data for the islands of (B) San Andres, and (C) Providencia, and Santa Catalina. Santa Catalina is the smallest island located in the north part of Providencia. Modified from Rey et al. (2019b).
Frontiers in Marine Science | www.frontiersin.org "Results and discussions, " and section "Conclusion" presents the main conclusions.

Description of the Study Area
The Colombian Caribbean Sea extends from 8 • to 13 • N latitude and from the country's border with Panama in the southwest (SW) at longitude 79 • W to the Guajira Peninsula in the northeast (NE) at longitude 71 • W. The insular Colombian Caribbean is formed by the archipelago of SPSC, the Roncador, Quitasueño, Serrania, and other cays, adjacent islets, and deedwater coral reefs (Otero et al., 2016). The study area is the archipelago of SPSC (Figure 1), composed of the San Andres Island, Providencia, and Santa Catalina and a group of lesser islands, atolls, and coral reefs that are home to a variety of marine flora and fauna (Geister, 1973), making this area the most important tourist destination in Colombia, with roughly 500,000 visitors a year (Ortiz Royero et al., 2015). Figure 1A shows the ground level, which is referred to as the Mean Low Water Spring (MLWS), and the computational domain for the hydrodynamic model. San Andres is located 190 km east of Nicaragua, and roughly 480 km northwest of the Colombian coast. Providencia is 90 km north of San Andres, and Santa Catalina is separated from Providencia by a 150 m wide channel. The total surface areas for San Andres, Providencia, and Santa Catalina are approximately 2,678 ha, 2,065 ha, and 119 ha, respectively (Rey et al., 2019b). Based on topographic LIDAR data with a spatial resolution of 5 m for the archipelago of SPSC (surveyed from 2010 by the Oceanographic and Hydrographic Research -CIOH, of the General Maritime Direction -DIMAR), the San Andres Island ( Figure 1B) is characterized by a low-lying coastal plain in the eastern side with altitudes lower than 5 m, whereas in the western side there is a narrow low-lying area. In the inner sector, there are topographic elevations up to 91 m. The black polygon in Figure 1A represents the bathymetric data collected by CIOH between 2015 and 2017, surveyed with single beam echo sounders in shallow waters and multi-beam echo sounders in deep waters. Outside the polygon, the bathymetric data was complemented with the General Bathymetric Chart of the Oceans (GEBCO) database (International Hydrographic Organization, Intergovernmental Oceanographic Commission) (IHO-IOC, 2018). Figure 1B shows that the insular shelf extends outside the coral reef northeast of San Andres, while the southeast area has a narrow shelf, with corals very close to the coastline or absent. In Providencia Island (Figure 1C), the main low-lying coastal plain area is in the northeast, north, and northwest side. The channel between Santa Catalina and Providencia has average depths of 3 m. The inner sector of Providencia and Santa Catalina have altitudes up to 358 m and 135.6 m, respectively. At the northeast, east, and southeast of Providencia, the insular shelf goes beyond the coral reef.
The Colombian climate is characterized by dry and wet seasons. The former occurs from December to May, and the latter in the rest of the year, interrupted by a relative minimum in June and July known as the Indian summer (Andrade and Barton, 2000). Hurricanes occur from June to November with a peak in activity in September. During most of the wet and dry seasons (September to April), cold fronts generate swells that can create considerable damage to coastal structures (Ortiz et al., 2014). On average there are six cold fronts per year in the Colombian Caribbean zone, but years such as 2010 have presented 20 events, and are relevant for coastal planning as they generate the largest waves in the Colombian Caribbean Sea (Ortiz-Royero et al., 2013), except for the east of the Guajira Peninsula and the archipelago of SPSC where extreme waves are related to TC (Ortiz-Royero, 2012;Otero et al., 2016). During average conditions, the most energetic waves in the continental Colombian Caribbean coast are in the regions of Cartagena, Barranquilla, and Santa Marta while in insular part are in San Andres and Providencia, where stronger trade winds from the northeast generate local sea waves .

Tropical Cyclone Datasets Used
To characterize the TC climatology we used synthetic events generated using the method proposed by Emanuel et al. (2006Emanuel et al. ( , 2008, consisting of a downscaling technique to allow a better representation of TC in locations with a scarcity of historical events (Emanuel and Jagger, 2010), including gray swan events, which are not be predicted based on history but can have a high-impact and can be foreseeable using physical knowledge together with historical data (Lin and Emanuel, 2016). The TC model involves storm tracks from genesis to lysis. The genesis starts with random seeding of warm-core vortices with peak wind speeds of 12 m/s. These vortices may decay or become TC depending on the ocean and atmospheric environmental conditions. The surviving vortices are then moved by a beta-advection model (Marks, 1992) and intensified based on the Coupled Hurricane Intensity Prediction System (CHIPS) (Emanuel et al., 2004), according to environmental factors based on the National Center for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis (Kalnay et al., 1996), referred hereafter as NCEP. The resulting synthetic contained storm parameters such as date, time, position, maximum winds, central pressure, and storm size, which is measured by the radius of maximum wind (R mw ). For this study, we employed a large number of synthetic tracks under NCEP conditions for the late twentieth century (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) over the Caribbean Sea. The events included 1000 tracks passing within 250 km of the middle point between the San Andres and Providencia islands (81.55 • W, 12.95 • N). The synthetic hurricane dataset was split into TC categories, from tropical storms to category V hurricanes. These events were used to force a hydrodynamic model to generate SS and determine inundation levels.
The use of a synthetic hurricane dataset allows us to perform robust statistics for TC climate in the area. However, the dataset has only few major hurricanes, with most events being minor hurricanes and tropical storms. As such, the synthetic database limits the determination of the worst-case flood scenario in case of major hurricanes striking the study area. Therefore, we additionally used hypothetical events (to differentiate from synthetic events), which have straight-line trajectories, considering different possible angles of approach to a given site, and different combinations of storm sizes, and forward speeds (Ruiz-Salcines et al., 2021). The hypothetical events we modeled considered category V hurricanes approaching the archipelago of SPSC. The final events had a constant maximum wind speed intensity of 95.17 m/s, a constant forward speed of 5.87 m/s, and a constant R mw of 56.3 km during their lifetime, for tracks with 11 different approach angles to the area of interest (Figure 2). At the regional level, Hurricane Mitch in 1998 reached category V (79.73 m/s) before making landfall in Honduras, while at a local level, Hurricane Iota in 2020 made landfall in Providencia as category V (72 m/s). Nevertheless, we selected the storm parameters for the hypothetical events based on a conservative perspective of the climatology for this site. Due to the distance between San Andres and Providencia (90 km), we used different hypothetical hurricane datasets for each island. Each of the eleven directions considered three parallel tracks separated by a distance equal to the R mw , with the central track passing through the center of the island, giving a total of 33 hypothetical events for each island. Figure 2 shows the hypothetical tracks for San Andres and Providencia, the area of movement for hurricane tracks, the computational domain for calculating the hydrodynamic, and the domain for the wind and pressure fields (also used for the synthetic events). The hydrodynamic model was forced with the hypothetical dataset in each island to generate SS to calculate the inundation areas.
To assess the flooding from Hurricane Iota in 2020 we used storm parameters from the best track dataset from the National Hurricane Center and Central Pacific Hurricane Center (NHC, 2020) to create the input wind fields for the hydrodynamic model.

Hydrodynamic Model Setup
For the hydrodynamic simulations, we used the two-dimensional model MIKE 21 HD FM (DHI, 2014a), which solves the momentum, continuity, temperature, salinity, and density equations with turbulent closure scheme equations. It is based on the incompressible Reynolds averaged Navier-Stokes equations (RANS) under the Boussinesq and hydrostatic pressure approximation (DHI, 2014a). The spatial discretization of the equations was based on the finite volume scheme. The model uses a dynamic time step to optimize simulation speed while ensuring numerical stability and considers a barotropic density with a varying Coriolis force according to the domain. The wetting and drying algorithm is included following the work of Zhao et al. (1994) and Sleigh et al. (1998). The HD FM model used an unstructured mesh (Figure 1) with a resolution of 40 km in the offshore areas far from the islands, increasing its resolution to 80 m in the vicinity of the archipelago and with up to 30 m resolution in the nearshore and flood-prone areas.
Regarding the boundary conditions on the HD FM model, when using the synthetic hurricane dataset, the hydrodynamic model was only forced with wind and pressure fields to account only the TC SS. The computational domain and ocean boundaries were set to zero water level and the boundary along the Nicaraguan coast was set to land. The wind fields were generated at the gradient level by the parametric wind model embedded in the SLOSH hydrodynamic model (Jelesnianski et al., 1992) using the storm characteristics from synthetic and hypothetical TC datasets. This parametric wind model is used by the National Hurricane Center (US) to create wind fields and force the hydrodynamic SLOSH model to forecast hurricane SS inundation in the Gulf of Mexico and the Atlantic ocean (Zachry et al., 2015). As shown by Chen et al. (2019), hybrid wind fields (parametric winds embedded in reanalysis data) provides a better estimation of the tropical cyclone far field winds. While this is not possible for synthetic events, it could be done for historical events such as hurricane Iota (used in this study). Nevertheless, the far field winds are not expected to influence the maximum flood levels in the study area, which is located within the R mw , thus, the use of parametric wind models is considered more accurate.
The formulation for the parametric wind model is given by where V (r) is the gradient wind velocity, at a radius r in m/s, R mw is the radius of maximum wind in km and V max is the maximum sustained 1 min speed in m/s. For TC SS modeling purposes, the wind gradient level requires to be translated to surface level to account for surface friction. The wind correction was done based on the method proposed by Lin and Chavas (2012), and used in other studies (Ruiz-Salcines et al., 2019;Rey et al., 2020). As the model provides symmetric wind fields at gradient level, the translation to surface level provides asymmetric wind fields used to force hydrodynamic models. For the MIKE 21 storm surge simulations, we adjusted the 1-min. wind to a 10-min. average by a reduction factor of 0.893 (Powell et al., 1996). For a given height, wind speed averaged over 10-30 min speeds will typically be about 80% of the maximum 1-min wind speeds (Resio and Westerink, 2008).
The pressure fields were calculated from the pressure distribution model of Holland (1980) as follows, using the storm characteristics from synthetic and hypothetical TC datasets.
where P c is the central pressure, P n the ambient pressure, and B is the Holland shape parameter obtained from Eq. 2 where ρ is the air density and e is the base of the natural logarithm.
Regarding Hurricane Iota, we used the approximation given by Silva et al. (2002) (Eq. 4), as the R mw is not provided by the NHC (2020).
where p c represents the central pressure of the TC. This equation was found by Chang et al. (2015) as the most accurate when compared to other formulations. When using the hypothetical hurricane dataset, each simulation was run with an initial water level of 0 m and 0.354 m. The latter is the equivalent to the Mean High Water Level (MHWL) as in Rey et al. (2019a) and was calculated from tidal records in San Andres Island covering from January 1997 to February 2015. The tide series have a relative reference and several gaps, which were filled using the Fourier series approximation. Therefore, the gaps were filled before calculating the MHWL. The tidal form factor, or Courtier index (F), obtained was 1.45 and calculated according to the following equation: where K 1 , O 1 and M 2 and S 2 are the diurnals (1) and semidiurnals (2) tide constituents, respectively. In the study area, the tides have a mixed regime, with a predominantly semi-diurnal component (0.25 < F < 1.50), having two high tides and two low tides every day, as defined by Jigena-Antelo et al. (2015). Because of the tide series gaps, an approximation, based on Fourier series was applied to predict the values of the lost data (Alvera-Azcárate et al., 2005;Bauer et al., 2017). The prediction of missed values of the surface water level (including the residual tide) was reconstructed with the following equation: where: x0 = 2.38; A = 0.085; a = 30678 This is a procedure that can be applied repeatedly in a cycle until the series converges. Unfortunately, there are no tide gauge records neither in San Andres nor in Providencia during the passing of TC to validate the model results. Therefore, we employed the same model setup as in Rey et al. (2019aRey et al. ( , 2020, consisting of a constant eddy viscosity of 0.28 under the Smagorinksy formulation and varying bottom friction based on the Manning coefficients for various categories of land as in previous studies (Mattocks and Forbes, 2008;Zhang et al., 2012). The wind friction varied with the wind speed, as in Rey et al. (2018), with a constant value of 0.001255 for wind speed below 7 m/s and a constant value of 0.002425 for wind speed above 25 m/s, and a linear variation between those values.
For modeling the flood-prone areas derived from Hurricane Iota in Providencia and Santa Catalina, we used four configurations for the model setup: (a) the first to determine SS levels, consisting on forcing the HD FM model with the wind fields from Hurricane Iota over a constant water level set to zero in the computational domain and on the ocean boundaries, (b) the second to account for AT and SS contributions to water levels, where the previous configuration was modified to include tides in the ocean boundaries, as obtained from the tide global model (Andersen, 1995), which represents the major diurnal (K1; O1; P1, and Q1) and semi-diurnal tidal constituents (M2, S2, N2, and K2) with a spatial resolution of 0.25 • × 0.25 • . Besides the contribution of AT and SS to the water levels, wave runup including both setup and swash (Dorrestein, 1961;Longuet-Higgins and Stewart, 1963;Sallenger, 2000;Stockdon et al., 2007) may be important, particularly for insular environments where the bathymetry is usually quite steep near the coast. Usually, the 2% exceedance level for vertical wave runup (R 2 ) is used as an indicator of this value (Ruggiero et al., 2001). Therefore, as a first approximation to include the wave contribution to water levels at the coast, we included (c), a third configuration to account for WS contribution by adding wave radiation stresses to configuration (b). To do so, we coupled the HD FM model with the MIKE 21 third-generation Spectral Wave (SW) model (DHI, 2014b). For information regarding source terms of the SW model, governing equations, time integration, and model parameters, readers are referred to Sørensen et al. (2004) and the manual documentation (DHI, 2014b). However, this coupled HD FM -SW model does not take into account the swash contribution to seawater levels. Consequently, a fourth configuration (d) was implemented to solve in detail the coastal processes in the surf zone, such as the wave propagation in shallow water, and the wave runup. As such, we implemented the X Beach 1D model (Roelvink et al., 2010). For validation purposes, we used a bathymetry profile in the eastern part of Santa Catalina (north of Providencia) whose initial point (x = 0.0 m) corresponds to the location of the AW600 sensor (Figure 3), whose recorded data were used as a boundary condition. The last point of the profile (x = 7196 m) ends in the coast, passing through two sensors in x = 3024 m (AW1000) and x = 4008 (RBR), which were used to validate the model. The model was forced with wave parameters (significant wave height Hs, peak wave period Tp, mean wave direction and spectral conditions) and a sea water level recorded by the AW600 sensor.
For model calibration, we focused on bed friction parameters as this parameter does not change from mean to storm conditions (bed friction due to seafloor configurations). Specifically, we calibrated the values of the bed friction coefficient C f , which has been used in the same model for similar sea floor configurations (Roelvink et al., 2021), getting values of 0.01 for the points before x = 3600 m and after x = 7090 m, and values of 0.9 in between, where coral rubble can be found (Díaz et al., 2000). According to the seafloor configuration for other profiles, we FIGURE 3 | Computational domain for the X Beach model. extrapolated these results to find the zones of the domain where each C f value applies.
Using the non-hydrostatic mode, from the X Beach profile we extracted the free sea surface elevation every second in the points where the AW1000 and RBR sensors are located. We validated the model by comparing simulated Hs to observations. After validation, we used the same configuration to run the model in profiles located in the eastern and western parts of Providencia. To do so, we forced the X Beach model with wave parameters and water levels obtained from the SW and HD FM models, respectively. The HD FM model results were processed following the National Hurricane Center (NHC) methodology to obtain the Maximum Envelope of Water (MEOWs) and Maximum of MEOWs (MOMs) (Zachry et al., 2015). The MEOW represents the maximum SS resulting from a set of hypothetical storms of a given category (from tropical storm to category V) and forward speed with varying sizes (related to the R mw ), initial water level, and unique track direction. The MOMs are composed of the maximum SS at each grid/cell element based on the category of storms from all the storm directions, storm sizes, forward speeds, and AT levels (zero or high). It means that MOMs are separated by category and initial water level, retaining the maximum SS value in each grid cell. MEOWs for the synthetic tropical cyclone dataset were not computed because events for a given storm category are allowed to move randomly. There is not a group of parallel synthetic tracks moving in the same direction to reach the coast as in the hypothetical dataset. Therefore, only MOMs for all the categories (tropical storm-category V) for zero tide level were calculated for the synthetic dataset. For the hypothetical dataset, MEOWs and MOMs for hurricane category V were computed for zero and high tide levels.

Extreme Storm Surge Statistical Analysis
The synthetic tropical cyclone dataset is composed of 1000 events. However, we used only events with equal or higher intensity than tropical storms inside a polygon centered in the middle point between San Andres and Providencia and a radius of 250 km. In this sense, only 738 events were considered. The total number of events per category were 366, 256, 49, 36, 22, and 9 for tropical storms, and categories I, II, III, IV, and V, respectively. The surge heights return periods were calculated by combining the SS cumulative distribution function (CDF) and the annual storm frequency (Lin et al., 2010). The SS CDF is estimated based on exceedance statistical methods. Since extreme events usually exhibited a heavy tail, we modeled the tail of the SS CDF using the peaks-over-threshold method (POT) with a generalized Pareto distribution (GPD) and maximum likelihood estimation (Coles, 2001). The threshold values were selected with the criterion of finding the best GPD curve fitting the empirical data.

RESULTS AND DISCUSSION
This section shows the 2D MIKE 21 and 1D Xbeach inundation results for Hurricane Iota at Providencia and Santa Catalina, synthetic and hypothetical tropical cyclones inundation maps, as well as the statistical analysis of the synthetic TC SS for the archipelago of San Andres, Providencia and Santa Catalina.
Tide gauge records for San Andres have a mixed tide regime with semidiurnal dominance with strong neap-spring variability, and a tidal range of 1.01 m, varying from 1.90 m during neaps, to 2.91 during springs. The MHWL was calculated based on the Fourier series approximation. Because of the lack of tide gauge records during the pass of TC near San Andres, it was not possible to validate the HD FM model results. However, we used the same calibration parameters used with the hydrodynamic model in previous studies (Rey et al., 2019a(Rey et al., , 2020, and we consider this model results reasonable as a first approximation for assessing the inundation threat from tropical cyclones at the study site. Figure 4 shows hydrodynamic modeling results for the maximum envelope of water depth (i.e., water level above the terrain) induced by Hurricane Iota for Providencia and Santa Catalina, considering the SS contribution (Figure 4A), SS plus AT (Yu et al., 2019), which compose the storm tide (Figure 4B), and storm tide plus WS ( Figure 4C). The total flooded area considering only SS was larger than storm tide since Iota passed over Providencia during low tide. The simulation considering storm tide plus WS showed a 14.93 % increase of the flooded area compared to the case when considering only storm tide (an increase of 11.45 ha), similar as in Wu et al. (2018), which used an idealized, alongshore uniform domain, finding that waves have a major influence on the maximum inundation distance inland, with contribution up to 16% of the inundation distance, depending on the different storm characteristics. This significant contribution of the WS to the flood-prone area is because of the bathymetry, which is quite steep beyond the coral reef. Waves make the dominant contribution to the total surge in coastal areas with steep slopes such as near reefed islands or levees. For instance, during typhoon Soudelor in 2015 the wave setup reached up to 0.5-1 m along the northeastern coast of Taiwan (Chen et al., 2017). But even on shallow slopes, the wave contribution could be significant, as an example waves added about 0.6-1.2 m to the surges during Hurricane Katrina in the coast of Louisiana, United States (Resio and Westerink, 2008). Due to computational constraints and the increased computational time when including WS, synthetic and hypothetical events simulations did not include WS, despite the increased flooded area observed under Iota when including WS.
Hurricane Iota passed 11 km north of Providencia. Thus, the maximum winds were from offshore north of the island in the forward direction of movement. However, over the island, the winds were in the backward direction of Iota's track, resulting in the highest flooded water depths northwest of the island. If Iota had passed on the south of Providencia with a distance similar to the R mw , the hurricane SS would have been larger in the east of the island. Table 1 shows the land area (in hectares, ha) for Providencia and Santa Catalina, as well as the total flooded area induced by hurricane Iota, and flooded area for different water depth intervals and associated percentage related to the total area of each island. The modeled SS induced by Iota on San Andres was negligible (not shown), as Iota passed roughly 100 km to the north.
The previous analysis did not consider the swash contribution to flooding. Therefore, to solve the hydrodynamic details in the surf zone, we modeled the wave propagation in shallow water and the wave runup employing the X Beach 1 D model. Figure 5 presents the Hs and tide levels used as forcing (AW600, Figure 5E) and to validate the X Beach model (AW100 and RBR, Figures 5F,G). Free sea surface levels were recorded at the three sensors between 6/03/2021 at 11h and 15/03/2021 at 15h. The profile and position of the sensors are presented in Figure 5D. The correlation coefficient between the modeled and measured Hs by the AW100 and RBR sensors are 0.93 and 0.80, respectively; and the RMSE are 0.012 m and 0.0025 m, respectively. These results show the capacity of the model to adequately reproduce the hydrodynamics process in the surf zone.
Once the X Beach 1D model was validated, we used the model setup for modeling the coastal hydrodynamic in other profiles at Providencia during the pass of Hurricane Iota, specifically at the eastern and western part of the island, where the flood analysis in 2D shows the largest flood-prone areas ( Figure 4C). We assessed the contribution of AT, SS, and R 2 to the maximum water level reached, R high (Sallenger, 2000) (Figure 6), which has vertical and horizontal displacements, varying on each time step. R high is modulated by SS ( Figure 6B). Based on the R high values at hourly intervals, we estimated the 2% exceedance water level for vertical R high (R high2 ) during the simulation period (42 h). R high is simulated with the X Beach model, and SS and the storm tide with the HD FM model. R 2 was obtained by subtracting storm tide from R high2 . R low is composed of storm tide and WS (Sallenger, 2000), simulated with the coupled HD FM -SW model. The WS was estimated by subtracting the modeled storm tide from R low . Percentage contributions to R high2 during the storm peak by AT, SS, WS and R 2 were −3.01 %, 46.3 %, 17.17 %, and 56.55 %, respectively ( Figure 6C). The negative value is because the AT levels for the storm peak period were low. The spatial evolution of these contributions as well as Hs for the storm peak, using as a reference of the maximum R high2 (2.22 m) show how the attenuation of Hs is compensated with the increase of SS and WS ( Figure 6D). However, the increase of R high2 is remarkable, with values above 2 m, and this elevation is evident in the swash excursion of around 500 m (Figures 6E,F) inshore. These results confirm oral testimonies of inhabitants and show  the importance of wave setup and runup contributions in the flooding process. Finally, Figure 6G shows the percentiles for the simulation period, the storm peak, and the temporal evolution of R high where the exceedance at 90% is 1.925 m during the storm peak, also confirming the flooding and infrastructure loss in the island during the pass of Hurricane Iota. A similar analysis as the previous one is shown for a profile in the western side of Providencia (Figure 7). The maximum of R high2 occurred 1 h later than in the eastern part of the island. Due to the location of the profile and the predominant eastern wave direction, the percentage contributions of the AT, SS, WS, and R 2 to R high2 for the storm peak were −3.56 %, 80.42%, 12.07%, and 23.14%, respectively, which suggests that the SS was the predominant contribution to the maximum seawater level reached on the western side of the island. Therefore, while in the eastern part of Providencia, the predominant contribution to R high2 was by R 2 , in the western was by the SS. The contributions to flood from R 2 and SS are mainly conditioned to the coastal environment such as the bathymetry and configuration of the coast as well as the storm characteristics. Figure 8 shows the classification of the synthetic tropical cyclones by category, where each panel shows the entire track of the color-coded events with the category at each track location.  Please note that the classification by category is done by using the highest category attained when the event pass inside the red circle, which has a center in the middle point between San Andres and Providencia and a radius of 250 km. As expected, when increasing the hurricane category, the number of events decreases.  Histograms of the simulated SS generated by the synthetic tropical cyclones at San Andres and Providencia islands are shown in Figure 9. The SS values used to create the histograms were extracted for specific points located at sea and inside the computational domain. The coordinates are 81.707 • W and 12.5721 • N (tide gauge location), and 81.3537 • W, 13.3608 • N for San Andres and Providencia, respectively. The histogram peaks are approximately at SS level of 0.1 and 0.12 m for San Andres and Providencia, respectively, and rapidly decrease as the values of SS increase, showing a large tail that extends to the maximum SS of 1.68 m and 1.46 m, respectively. Figures 10, 11 show synthetic TC MOMs from tropical storm to category V hurricane for zero tide level at San Andres, Providencia, and Santa Catalina, respectively. Each category MOM produces higher SS inundation compared to the corresponding lower category MOM. The city blocks (black lines on the maps) were generated by the National Administrative Department of Statistics of Colombia (DANE, 2018). The results show that the highest exposure area to SS flooding is in the northeast part of San Andres Island, particularly the areas surrounding the harbor and the north part of the island, which is also the most populated area. Other areas do not seem to be exposed to floods, mainly because of their physical characteristics. For instance, while the northeast part of the island has a flat insular shelf and coral reef, which dissipate the waves but amplifies SS (Flather, 2001), the southeast of the island has a narrow shelf with a scarce presence of coral reefs, resulting in negligible SS despite suffering significant wave impact (Ortiz Royero et al., 2015;Bernal et al., 2016). However, the coral barrier reef seems to work as a natural wall that decreases the impact of hurricanes in those islands. The role of the coral reef in inundation, though somehow studied (Gallop et al., 2020), will be a topic of future research for those islands. The total flooded area and associated percentage related to the total area of San Andres are 37.08 ha (1.38%), 52.97 ha (1.97%), 59.55 ha (2.22%), 78.70 ha (2.93%), 95.86 ha (3.57%), and 148.93 ha (5.56 %) for tropical storm, and categories I, II, III, IV, and V, respectively.
The northeast of Providencia Island (Figure 11) is the most prone to be flooded because of the flat insular shelf, coral reef, and low topography. Whereas there are no settlements in that area, there is an airport exposed to flooding. For northwest of Providencia, as well as southeast of Santa Catalina, settlements exist in flood-prone areas. The total flooded area and associated percentage related to the total area of Providencia are 3.29 ha (0.15%), 27.32 ha (1.32%), 52.48 ha (2.54%), 58.29 ha (2.82%), 66.49 ha (3.21 %), and 69.30 ha (3.35 %) for tropical storm and categories I, II, III, IV, and V, respectively. These values for Santa Catalina are 0.83 ha (0.7 %), 0.92 (0.78%), 1.64 ha (1.38%), 2.22 ha (1.87%), 2.92 ha (2.46%), and 3.62 ha (3.05%) for tropical storm and categories I, II, III, IV, and V, respectively. The percentage of the flooded area for Providencia derived from the Iota SS is similar to the synthetic hurricane category IV MOM (Figure 11e), and the one generated for Iota SS in Santa Catalina is larger than the synthetic hurricane category V MOM by 8.81% (Figure 11f). Therefore, MOMs based on hypothetical events are important to characterize and cover all susceptible flood-prone areas from TC, which is complicated with synthetic events. More events would be needed.
As an alternative to the lack of enough major hurricanes in the synthetic dataset to assess all the potential exposed areas to hurricane floods, we used hypothetical tracks. This allows assessing the inundation threat for a given area to provide valuable information for evacuation plans. Figure 12 shows the MOMs for hypothetical category V hurricane during (a) low tide and (b) high tide, as well as the flood-prone areas (maximum water depth envelope) induced by the maximum SS resulting from the 738 synthetic tropical cyclones (c) for San Andres. Figures 12d-f shows the same information as described before for Providencia and Santa Catalina. Whereas the hypothetical hurricane category V MOM maps for (a), (b), (d) and (e) in Figure 12 are useful for showing the potential hurricane SS floods, (c) and (f) may not be called MOM because is composed of several synthetic events with different categories, the MOM is composed for events with the same category. However, the information provided in (c) and (f) is useful to determine the probability for different hurricane SS levels. For San Andres, the MOMs show that the maximum water depth is generated around the harbor. This maximum is a result of the diminished flushing capacity through the inlet (Rey et al., 2018). As result, the surrounded area such as the south of the harbor is flooded, this area is characterized as having low altitudes (less than 1 m) and being covered mainly by mangroves. Fortunately, there are  just a few houses in that area, at the head of the harbor. The total flooded area and associated percentage related to the total area of San Andres are 354.83 ha (13.24 %), 392.23 ha (14.64%), and 164.63 ha (6.14 %) for the hypothetical hurricane category V MOM for zero tide level, hypothetical hurricane category V MOMs for high tide, and the flood-prone area induced by the synthetic TC SS, respectively. When comparing the flooded area for the hypothetical hurricane category V MOM for zero tide level tide (Figure 12a) with the one for high tide, the latter ( Figure   12b) is larger up to 9.53%. It is recommended for future studies to create MEOWs and MOMs based on hypothetical tracks for all the categories (tropical storm to Category V) since synthetic TC dataset is not enough to reproduce the hurricane SS inundationprone areas. MEOWs and MOMs are used for evacuation plans from a conservative perspective rather than to forecast a single event (Zhang et al., 2013). The flooded area induced by the nearworts-case hypothetical MOM category V for low tide (Figure  12a) is larger up to 53.60 % than the one induced by the synthetic TC events (Figure 12c). As major hurricanes are scarce in the synthetic TC dataset, the use of hypothetical hurricanes provides the opportunity to assess the flooding from higher category events. However, it is important to note that the maximum hurricane SS does not depend only on the hurricane intensity, but also on the storm size, forward speed, direction of approach, and the configuration of the basin and the coastal bathymetry (Resio and Westerink, 2008;Lin et al., 2010).
Regarding the total flooded area (Figures 12d-f) and associated percentage related to the total area of Providencia, the values are 91.12 ha (4.41%), 98.57 ha (4.77%), 69.84 ha (3.38 %) for the hypothetical hurricane category V MOM for low tide, and high tide, and for the flood-prone area induced by the synthetic events, respectively. For Santa Catalina, these values are 6.79 ha (5.72 %), 8.89 ha (7.49 %) and 3.62 ha (3.05 %) for the hypothetical hurricane category V MOM for low tide, and high tide, and the flood-prone area induced by the synthetic events, respectively. The flooded area induced by hypothetical hurricane category V MOM for high tide is 7.55%, and 23.62% larger than for low tide for Providencia, and Santa Catalina, respectively.
The main contributors to coastal flooding during the pass of tropical cyclones are the AT, the SS (wind setup plus pressure setup), the wave runup (swash plus wave setup) and the freshwater (river, rainfall). MEOWs and MOMS do not account for the inclusion of wave runup neither freshwater to the total water depth in each hydrodynamic simulation and thus in the MEOWs and MOMs products (Zachry et al., 2015). Additionally, the dynamically forced tidal signal is not included, only a constant water level is considered. Once known these limitations, products from this study may be used adequately.
Another geospatial tool often used for assessing coastal floods is the Geographic Information System (GIS) models. However, the use of hydrodynamic models is preferable to the GIS models for modeling coastal floods (Seenath et al., 2016;Ruiz-Salcines et al., 2021). The latter is primarily based on topography, and it does not take into account the physics of the fluid dynamics (e.g., bottom friction, fluid flow direction, structural barriers). Therefore, in certain cases, it is possible that the use of GIS over-estimates the flood levels and consequently causes overengineered schemes for defenses and inflated estimates of socioeconomic impact. For instance, if the GIS model had been used in this study, all the east parts of San Andres Island from north to south would have been resulted exposed to inundation because the topographic levels are less than 5 m over that zone. However, in the southeast part of the Island, the shelf is narrow with steep slope. Therefore, under these conditions, the SS induced coastal flood is negligible, which is unknown by a GIS model. Figure 13 shows the hurricane SS return periods for San Andres and Providencia at the same locations used for creating the SS histograms (Figure 9). The thresholds of 0.18 m and 0.14 m were used to define the tail for statistical fitting with the GPD, for San Andres and Providencia, respectively. In order to know the Hurricane Iota SS return period at the location of the extreme analysis at Providencia, we extracted the maximum SS of Iota from the computational domain for that point, which is 1.25 m. Thus, looking into the SS return period curve for Providencia, this value has an associate return period of 3,234years ( Figure 13B). Given the low probability of the Hurricane Iota SS under present climate, future studies should consider the effect of climate change on tropical cyclones Marsooli and Lin, 2020) to know if storm surge values as the one for Hurricane Iota will be more frequent in Providencia and Santa Catalina at the end of the century.
Based on the water depth induced by the synthetic TC events, we estimated the 100-year (Figures 14a,b) and 500-year (Figures 14c,d) return periods for San Andres, Providencia, and Santa Catalina. From each element of the computational domain on the flooded area, we extracted the water depth associated with 100-year and 500-year return periods. This analysis was done based on the empirical estimated water depth return levels. It means that the water depth at each element was not fitted to a canonical probability distribution. The water depth heights return periods were calculated by combining the water depth heights CDF and the annual storm frequency. Whereas water depth levels associated with 100-years return periods cover partially the east part of these islands, the 500-years return periods increase significantly the flooded area around the harbor and the airport for San Andres and Providencia, respectively. Flooded areas outside of the 500-years flood return period map, such as the ones shown in previous figure (Figures 12c,f), have larger return periods. It means that the water depth levels on the north of San Andres have return periods larger than 500 years.
The use of synthetic or hypothetical datasets to do flood hazard assessments depends on the objective of the study. As such, synthetic events provide a robust dataset for statistical analysis to characterize present and future climate trends, whereas hypothetical events allow assessing the inundation threat, which is very useful for evacuation plans from a conservative perspective (Rey et al., 2019a), and also for estimating the near-worst-case flood scenario, including climate change (Ruiz-Salcines et al., 2021). For coastal environments with steep slopes, the contribution of waves may be larger or equal than SS to flooding (Stockdon et al., 2007;Chen et al., 2017) while in wide and shallow continental shelf with a mild slope (e.g., continental shelf for Yucatan in the Gulf of Mexico) this behavior is inverted (Rey et al., 2020).

CONCLUSION
This study presents a characterization of the hurricane SS inundation threat for the archipelago of San Andres, Providencia, and Santa Catalina, which is a site with scarce historical events. SS was investigated using a hydrodynamic model forced with synthetic tropical cyclones and hypothetical hurricanes. The former are mathematical events derived using hurricane physics useful to determine the probability for different hurricane SS levels and the latter are events with category V hurricanes during all their lifetime and with different track direction toward the archipelago to establish the near-worst-case hurricane SS flood scenarios. The use of synthetic events allowed us to perform statistical analysis of hurricane SS heights at San Andres, Providencia and Santa Catalina. This analysis suggested that the eastern part of San Andres is more likely to be flooded from hurricanes SS than the north, and for Providencia the east is more likely to be flooded than the northwest. The hypothetical category V hurricane MOMs for high tide show the worst-case flood scenario from hurricanes SS and a constant water level associated with the MHWL. In this sense, the flooded area and associated percentage related to the total area of San Andres, Providencia, and Santa Catalina are 392.23 ha (14.64%), 98.57 ha (4.77%), 8.89 ha (7.49 %), respectively. This is because of the role that plays the flat insular shelf and coral reef, which would be the topic of a future study. Additionally, we modeled the Hurricane Iota (2020) inundation for Providencia with 2 D and 1 D models. For the former we considered contributions from SS, AT, and WS. The flooded area derived from these three contributions was 75.93 ha, which represents 3.68% of the island. The contribution of the wave setup was 14.93 % (11.45 ha) of the total flooded area modeled. For the latter, modeling results with X-beach (1 D) in the eastern side of the island showed that the AT, SS, and R 2 contribution to the maximum water level reached -3.01%, 46.36%, and 56.55 %, respectively during the storm peak. This information provided an approximation of the underestimation of the seawater level when modeling only SS without including the wave contribution through the WS and swash in insular environments. The maximum SS derived from Hurricane Iota for the east of Providencia was 1.25 m, which according to this study has an associated return period of 3,234-years. Therefore, future studies should consider the effect of climate change on tropical cyclones to know if the probability of storm surge values as the one for Hurricane Iota may be increased in Providencia Island with the current global warming trends. This methodology may be applied to other coastal zones and extended to consider the effect of climate change on hurricane SS and the sea level rise, which together with the increase of settlements along coastal zones will increase the flood risk by the end of the century. This study provides tools for decision-makers to implement initiatives to make those islands more resilient.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.  . This work was funded by COLCIENCIAS, DIMAR-CIOH, and partially by the Centro Mexicano de Innovación en Energia del Oceano, CEMIE-Oceano (OLE-1). AFO and JPR acknowledge the data provided by the project Convenio 003 de 2020 CORALINA-UNAL, "Estudio para definir el riesgo por inundación en la isla de San Andres, ante huracanes en las categorías más probables".