Spatiotemporal Storm Impact on the Northern Yucatan Coast during Hurricanes and Central American Cold Surge Events

: We investigate the storm impact associated with historical events in the northern Yucatan Peninsula. The study area is prone to coastal ﬂooding due to both its geographical location and low-lying areas. Extreme events associated with tropical cyclones and Central American cold surge (CACS; locally known as Nortes ) are ubiquitous in this region, and coastal development in the study area has exacerbated the erosion of the sand beach-dune system. This study aims to assess the impact on the northern coast of Yucatan associated with di ﬀ erent types of storms and to investigate the role of the dune in its spatial variability. Nearshore hydrodynamics, associated with hurricanes (Gilbert: 14 September 1988; Isidore: 22 September 2002) and energetic Nortes ( Norte A: 12 March 1993; Norte B: 25 December 2004), were computed using a numerical model. The beach and dune characteristics were extracted from a LIDAR ﬂight with a spatial resolution of 1 m conducted in 2011. Furthermore, the extreme water levels and the spatiotemporal variability of the storm-impact regime (swash, collision, overwash, or inundation), along a 41.5 km stretch of coast, were derived using both runup parametrizations and the modeling results. On the one hand, the predominant storm impact regimes for Hurricanes Gilbert and Isidore were inundation and overwash, respectively. The ﬂood that propagated from east to west in the northern Yucatan was due to westerly-directed hurricane tracks. On the other hand, for the Norte events, the predominant impact regimes were collision and overwash for Nortes A and B, respectively. This di ﬀ erence in the impact regime between Norte events can be ascribed to tidal di ﬀ erences. Moreover, during the passages of Nortes A and B, the ﬂood was propagated from west to east in the northern Yucatan, consistent with cold-front paths. The results suggest that the western part of the study area presented a stronger impact regime due to the dune degradation caused by coastal infrastructure and settlements established in those areas. This work highlights the important role of sand dunes in providing natural coastal protection during Norte events.


Introduction
Extreme storm events, such as Central America cold surge and tropical cyclones, play an important role in coastal flooding [1,2] and beach erosion [3,4] on the northern Yucatan coast. The low-lying coast and wide continental shelf of the northern Yucatan Peninsula, as well as the presence of semi-enclosed back-barrier lagoons and wetlands, enhance the physical vulnerability to coastal flooding during extreme storm events [1]. The net potential alongshore sediment transport is to the west, ranging between 20,000 and 80,000 m 3 yr −1 along the northern Yucatan Peninsula [5]. Beach erosion in the Yucatan coast is critical at some locations [6] and is mainly associated with longshore transport gradients induced by coastal infrastructure (e.g., jetties, groins, breakwaters, offshore port in Progreso) [3,7,8]. Moreover, coastal development has contributed to the degradation of dunes, making the coasts less resilient to storm events.
The storm impacts on barrier islands present a significant spatial variability mainly due to the alongshore-variability of both beach morphology and the intensity of the ocean's forcing such as storm surge, waves, and wave runup [9][10][11], which make the berms and dunes more vulnerable to erosion and overtopping because the waves attack the higher part of the dune profile when the water levels increase along the coast [10]. In fact, minor hurricanes or tropical storms can lead to overwash and breaching if no protection is afforded by dunes, even when the storm passes within some distance of a given beach. For instance, Santa Rosa Island, located in northwest Florida, suffered beach erosion and widespread overwash and breaching during Hurricane Katrina (2005) even though the storm made landfall 145 miles to the west [12]. Sallenger developed a simple model [9], based on the relative elevation between the extreme water levels induced by storms and the height of the foredune that categorizes storm-induced impacts associated with different levels of erosion on barrier islands. For each category, net sand-transport directions and processes are argued to be unique, and the limits between them represent thresholds where processes and magnitudes of impact change significantly [9]. This approach provides a framework to forecast the impact of a storm and has been validated [10] and implemented in several studies showing satisfactory results at regional [13,14] and local levels [15,16]. Assessing the coastal response to tropical and extratropical storms [17,18] may be more useful for decision-makers than studying the intensity of these natural phenomena.
This research aims to assess the coastal response of the northern Yucatan coast to Central American cold surge (CACS) and hurricane events. We focus on the spatiotemporal variability of the storm impact owing to differences in storm tracks and beach morphology features.

Study Area and Historical Events
The study site is located on the northern coast of the Yucatan State, comprising a 41.5 km stretch of coast between Yucalpetén and Telchac ports ( Figure 1). The coastal geomorphology corresponds to a barrier island with both natural inlets and artificial channels, which allow sea and fresh-water exchange with lagoons and wetlands where the predominant freshwater input is from springs formed by fractures of the coastal confining layer in the large coastal Yucatan aquifer [19]. This dynamic exchange leads to water and soil salinity variability, decrease of soil pH, and excessive siltation that increases the wetlands' physiologic stress level and induces changes in vegetation distribution [20]. The land cover classification for the study zone is variable with forest, shrublands, grasslands, and water bodies, according to the 30 m resolution dataset from the National Geomatics Center of China dataset [21], as shown in Figure 1c. A detailed description of the study area characteristics is found in [1] and [16].
Of particular interest is the port of Progreso, which is the largest and most economically important coastal city of the State of Yucatan, located at the east of the Yucalpetén harbor. At the south and west of Yucalpetén is the Chelem Lagoon, which extends along the coast bordering the inland side of the sand barrier over which Progreso and smaller towns (Chelem, Chuburná) are settled. Considering the location of Progreso over the sand barrier, the town is also threatened by inundation from the lagoon as storm events may raise water levels in the Chelem Lagoon [1].
interaction with strong northerly winds from a CACS event, and subsequently left the Peninsula as a tropical storm. Norte A corresponds to a CACS event that reached the study area during spring low tide, occurred at the same time than "the 1993 superstorm cold surge", also known as "the storm of the century" [22,23]. Norte A brought northerly winds exceeding 20 ms -1 and a temperature drop of up to 15 °C over 24 h into the Gulf of Mexico. Norte B occurred during spring high tide. The maximum wind speeds for Norte B were similar to Norte A, but the period with wind speeds larger than 20 ms -1 was longer for Event A (11 h compared to 3 h for Event B). A detailed characterization of the wind and wave climate associated with these events can be found in [1].

Methods
To investigate the spatiotemporal evolution of the storm impact on the northern coast of the Yucatan State associated with different storm events (i.e., hurricanes and Nortes), we employed the storm impact model proposed by [9]. This model compares the extreme water level with the dune morphology to forecast the impact of a storm as it approaches the coast. The storm impact was evaluated along the stretch of coast between Yucalpetén and Telchac at different cross-shore transects. The dune morphology was characterized every 500 m using LIDAR information (1 m resolution) from 2011. Numerical models and empirical parameterizations were employed to The selected storm events for this study are Hurricanes Gilbert (September 1988) and Isidore (September 2002), and two energetic CACS events named Norte A (March 1993) and Norte B (December 2004). The tracks for Hurricanes Gilbert and Isidore are shown in Figure 1. Gilbert reached the east of the Yucatan Peninsula as Category 5 and left the Peninsula passing over Progreso port as Category 3, according to the Saffir-Simpson hurricane wind scale (SSHWS). Hurricane Isidore made landfall at Telchac port as Category 3, remained roughly for 36 h over the Peninsula due to the interaction with strong northerly winds from a CACS event, and subsequently left the Peninsula as a tropical storm. Norte A corresponds to a CACS event that reached the study area during spring low tide, occurred at the same time than "the 1993 superstorm cold surge", also known as "the storm of the century" [22,23]. Norte A brought northerly winds exceeding 20 ms −1 and a temperature drop of up to 15 • C over 24 h into the Gulf of Mexico. Norte B occurred during spring high tide. The maximum wind speeds for Norte B were similar to Norte A, but the period with wind speeds larger than 20 ms −1 was longer for Event A (11 h compared to 3 h for Event B). A detailed characterization of the wind and wave climate associated with these events can be found in [1].

Methods
To investigate the spatiotemporal evolution of the storm impact on the northern coast of the Yucatan State associated with different storm events (i.e., hurricanes and Nortes), we employed the storm impact model proposed by [9]. This model compares the extreme water level with the dune morphology to forecast the impact of a storm as it approaches the coast. The storm impact was evaluated along the stretch of coast between Yucalpetén and Telchac at different cross-shore transects. The dune morphology was characterized every 500 m using LIDAR information (1 m resolution) from 2011. Numerical models and empirical parameterizations were employed to estimate extreme water levels at each transect. Detailed information on the methods is provided below, where Section 2.2.1 describes Sallenger's model and Section 2.2.2 shows the characterization of the dune morphology for the northern Yucatan coast. The procedure to obtain the extreme water levels are presented in Sections 2.2.3 and 2.2.4, where the numerical (hydrodynamic and wave) models and the parametric wave runup model are described, respectively.

Storm-Impact Scaling Model
The storm impact scale proposed by [9] predicts the storm impact on barrier islands by coupling the fluid forcing (R high , R low ) during the pass of storm events and the dune morphology (D high , D low ). R high is defined as the water level obtained as the sum of the astronomical, the storm surge, and the 2% of exceedance runup (R 2% ). The wave runup includes both setup and swash contributions [24]. R low represents the elevation below which the beach is continuously subaqueous and is defined as the sum of storm surge, astronomical tide, and the wave-induced setup (η). Thus, the difference between R high and R low is the 2% given by the exceedance swash amplitude included only in R high . D high is defined as the dune crest height or, when the foredune is not present, the beach berm. D low is the elevation of the base of the dune (toe). Therefore, in the absence of the foredune ridge D low = D high .
Considering the variation of R high and R low with respect to D low and D high , four impact regimes can be identified [9]: (1) The swash regime is the condition during storms in which the wave runup is confined to the foreshore region (R high < D low ). The beach foreshore erodes and sand is moved offshore and returned to the beach during more calm conditions [10]. (2) The collision regime occurs when the maximum water level exceeds the base of the dune (D high > R high > D low ). The runup collides with the base of the dune causing erosion that may be more long-lasting than foreshore erosion as dunes are typically rebuilt through slower aeolian processes [10]. Eroded sand can be transported offshore or longshore at the same time and usually does not return to re-establish the dune. (3) The overwash regime occurs when R high > D high , implying that the runup overtops the dune (or berm crest), leading to erosion of the dune. The sand is transported landward and not readily returned to the seaward side of the island under post-storm conditions. However, other processes, such as aeolian transport may return overwashed sand to the beach foreshore [25]. The overwash regime contributes to a barrier island's net landward migration. (4) The inundation regime occurs when R low > D high and the sea level increase is enough to submerge a barrier island completely. The beach and dunes are completely and continually subaqueous. Limited observations suggest that the eroded sand undergoes a net landward transport on the order of 1 km. Therefore, to refine and expand the characterization of the inundation regime, more observations are needed [9].
Given that dune morphology and extreme water levels associated with extreme events vary in the alongshore direction, spatiotemporal variability of the impact regime is expected during the pass of storm events. However, the dune geometry can change during a storm, which in turn can change the impact regime. Therefore, the storm-impact scaling model is relevant to forecast the initial impact of a storm when the geometry change is not significant (D high does not change to the extent that a new threshold is exceeded). If that change is significant, the result would be an increased impact level [9].

Beach-Dune Morphology
The digital elevation terrain (DET) used in this study is based on the LIDAR topographic data acquired in 2011 for a coastal stretch of 1.5 km wide along the 365 km of coast in the Yucatan state.
The spatial resolution of the DET is 1 m 2 . A total of 83 cross-shore profiles were extracted from the DET, spaced every 500 m, along the 41.5 km stretch of the Yucalpetén-Telchac coast ( Figure 1). All the profiles have the origin at the shoreline (X sl , Z sl ) and were extended up to 60 m inland as shown in Figure 2, where (X sl , Z sl ), (X dt , Z dt ), (X dc , Z dc ) are the cross-shore position at the shoreline, dune toe and dune crest, respectively. The shoreline, in this case, was defined as a line which separated the water from the land during the pass of the LIDAR flight. This line is roughly at Z sl =0.25 m, referring to the GGM06 Mexican geoid. Following [10], the width of the beach was defined as the horizontal distance between the shoreline (X sl , Z sl ) and the base of the dune (X dt , Z dt ) or, in the absence of a dune, the crest of the berm (X sl -X dt ). Dune characteristics and beach slope were estimated for the 83 transects. to the GGM06 Mexican geoid. Following [10], the width of the beach was defined as the horizontal distance between the shoreline (Xsl, Zsl) and the base of the dune (Xdt, Zdt) or, in the absence of a dune, the crest of the berm (Xsl-Xdt). Dune characteristics and beach slope were estimated for the 83 transects.

Waves and Nearshore Hydrodynamics
The offshore wave conditions were modeled using the MIKE 21 third-generation spectral wave (SW) model [26]. For information regarding source terms of the SW model, governing equations, time integration, and model parameters, readers are referred to [27] and the manual documentation [26].
The bathymetry was extracted from the ETOPO1 database [28] and was complemented with higher-resolution bathymetric data from 9-km long transects every 4 km along the coast ( Figure 3). In addition, the DET along a stretch of the Yucalpetén-Telchac was used. The model setup was based on previous studies [29,30], which calibrated the model for mean and extreme conditions in the Gulf of Mexico and the Caribbean Sea. Given the direction of motion for Hurricanes Gilbert and Isidore, each profile experiences a different wave height at each time step. To calculate the runup, for each profile a wave height time series was extracted from the computational domain at 10 m depth, as shown in Figure 3b.

Waves and Nearshore Hydrodynamics
The offshore wave conditions were modeled using the MIKE 21 third-generation spectral wave (SW) model [26]. For information regarding source terms of the SW model, governing equations, time integration, and model parameters, readers are referred to [27] and the manual documentation [26].
The bathymetry was extracted from the ETOPO1 database [28] and was complemented with higher-resolution bathymetric data from 9-km long transects every 4 km along the coast ( Figure 3). In addition, the DET along a stretch of the Yucalpetén-Telchac was used. The model setup was based on previous studies [29,30], which calibrated the model for mean and extreme conditions in the Gulf of Mexico and the Caribbean Sea. Given the direction of motion for Hurricanes Gilbert and Isidore, each profile experiences a different wave height at each time step. To calculate the runup, for each profile a wave height time series was extracted from the computational domain at 10 m depth, as shown in Figure 3b.
The hydrodynamic module of MIKE 21 is a two-dimensional model developed by the Danish Hydraulic Institute (DHI) Water & Environment [31], which resolves the depth-averaged (2D) incompressible Reynolds-average Navier-Stokes equations invoking the assumptions of Boussinesq and hydrostatic pressure. This hydrodynamic (HD) model solves the continuity and full momentum equations. The wetting and drying algorithm is included following the work of [32] and [33]. For more information about the HD model, readers are referred to the scientific manual documentation [31]. The hydrodynamic module of MIKE 21 is a two-dimensional model developed by the Danish Hydraulic Institute (DHI) Water & Environment [31], which resolves the depth-averaged (2D) incompressible Reynolds-average Navier-Stokes equations invoking the assumptions of Boussinesq and hydrostatic pressure. This hydrodynamic (HD) model solves the continuity and full momentum equations. The wetting and drying algorithm is included following the work of [32] and [33]. For more information about the HD model, readers are referred to the scientific manual documentation [31].
In the computational domain (Figure 3a), bathymetric data and wind fields used as forcing in the HD model are the same as for the SW model. The values used for bottom friction and horizontal eddy viscosity (Smagorinsky formulation) are the same used in [1]. The HD model was forced with astronomical tides extracted from the global tide model [34], varying in time and along the boundary between Cancun and Cuba and between Cuba and Florida (Figure 3a) as in [1].
For Norte events, the HD model was forced with pressure and wind fields from the climate forecast system reanalysis (CFSR) dataset [35]. For Hurricanes Gilbert and Isidore the CFSR was not used because it underestimates the wind fields [36] so they were instead computed as the sum of two components: i) the axisymmetric wind field, which was provided by the Sobey and Young (SB) parametric wind model at gradient level [37] and ii) the environment background field calculated from the storm's translational velocity. These two components were computed using as input storm characteristics (i.e., track, intensity, and size) from the North Atlantic Historical Hurricane Database (HURDAT) [38]. The formulation for the SB parametric wind model is given by In the computational domain (Figure 3a), bathymetric data and wind fields used as forcing in the HD model are the same as for the SW model. The values used for bottom friction and horizontal eddy viscosity (Smagorinsky formulation) are the same used in [1]. The HD model was forced with astronomical tides extracted from the global tide model [34], varying in time and along the boundary between Cancun and Cuba and between Cuba and Florida (Figure 3a) as in [1].
For Norte events, the HD model was forced with pressure and wind fields from the climate forecast system reanalysis (CFSR) dataset [35]. For Hurricanes Gilbert and Isidore the CFSR was not used because it underestimates the wind fields [36] so they were instead computed as the sum of two components: (i) the axisymmetric wind field, which was provided by the Sobey and Young (SB) parametric wind model at gradient level [37] and (ii) the environment background field calculated from the storm's translational velocity. These two components were computed using as input storm characteristics (i.e., track, intensity, and size) from the North Atlantic Historical Hurricane Database (HURDAT) [38]. The formulation for the SB 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. The use of wind fields for storm surge and wave modeling requires to translate the axisymmetric component of the wind field and the environment background field component from gradient level to surface level to take into account the surface friction experienced by the storm. For that, an empirical surface wind reduction factor (SWRF) of 0.85 [39] and an inflow angle based on the formulation of [40] was applied to the axisymmetric component. For the background wind field, a wind reduction factor of 0.55 and a counter-clockwise rotation angle of 22 • over the direction was employed as suggested in [41] and used by [39]. Corrections to wind speeds are needed to provide a more accurate asymmetric shape to the storm [42]. Similar to previous studies [41,43], the pressure fields were calculated for Hurricanes Gilbert and Isidore from the pressure distribution model of Holland [44] and the storm characteristics from the HURDAT database. The wind friction (Cd) for the four storm events is the same employed in [1].
One of the main issues for modeling the Yucatan coastal hydrodynamics is the lack of observations to validate the HD model. No sea surface elevation records were available in the study area for hurricane events and only a 3 m watermark registered at Telchac port during the pass of Hurricane Isidore was used to validate the HD model. It was assumed that the HD model generates adequately water levels for Gilbert using the same model setup as for Isidore.
Extreme water levels (R high, R low ) associated with hurricane and CACS events require knowledge of the tide, storm surge, wave setup, and wave runup associated with each event during its landfall. Tide gauge data was not available during the four storm events here investigated. To overcome this issue, the R low was estimated using the HD model coupled with the SW (hereafter called HD-SW) model following the model setup presented in [1]. R high was estimated by simulating the interaction between the astronomical tide and the storm surge using the HD module [31]. The wave runup and setup were calculated from simulated offshore wave conditions [26] using the empirical parameterization of [16] as explained below.

Parametric Wave Runup Predictions
The wave runup is the sum of the wave setup and swash oscillations [10]. The maximum elevation of the runup has been found to be a function of the deep-water wave height H 0 and wave period T 0 and the foreshore beach slope β 0 [24]. However, previous studies [16,45] have pointed out the suitability of a parametric model independent of the beach slope for the estimation of the vertical runup elevation (e.g., [46]). In this study, the wave runup parameterization of [16] was used for each of the four events.
This model was developed for the northern coast of Yucatan based on numerical simulations of nonlinear wave transformation [47] using a 30 years wave hindcast [48] in which Nortes and tropical cyclones were presented. The parameterization is a generalized expression for the prediction of the 2% exceedance value of runup (R 2% ) as a function of H and Z and a simple linear relationship was fitted to the hyperbolic tangent fit parameters, a and b. The generalized expression of R 2% is given by where H is the water wave height at 10 m deep, Z is the astronomic tidal level, a = 1.615Z + 1. Given that the water level from the HD model for the four events took into account not only the astronomic tidal level (Z) but also the storm surge, which was significant for the four events, the values of a and b for HWL were chosen. This parameterization was developed using 600 simulations from a non-hydrostatic nonlinear shallow water equations numerical model, Simulating Waves till Shore (SWASH) [47], implemented in this region. The correlation coefficient between numerical results and the runup parameterization was 0.80 during Norte events, which was considered suitable for this study.

Results and Discussion
This section shows the results for the spatiotemporal evolution of H s , R 2%, R low , and R high for the four simulated storm events.

Hydrodynamic Forcing
Water level time series, estimated by the HD and HD-SW models, during the crossing of Hurricane Isidore at Telchac port are shown in Figure 4. The spatiotemporal evolution of the extreme water level along the stretch of coast for the four cases are shown in the following figures. The HD and the HD-SW models reproduced a peak water level of 2.28 m and 2.39 m (referred to the GGM06 Mexican Geoid), respectively. This 11 cm difference is related to the wave setup contributions to the seawater level. Given the lack of water level measurements, the hydrodynamic model was validated with a watermark of 3 m (Figure 4), surveyed at Telchac port during the passage of Hurricane Isidore, and referred to the GGM06. The difference between the maximum water level reproduced by the HD-SW model and the watermark is 0.61 m, which may be related to the swash amplitude contributions, not considered in the numerical model. Given the lack of observed data for validating the HD-SW model for Hurricane Gilbert, the same calibration parameters as for Hurricane Isidore were employed.

Results and Discussion
This section shows the results for the spatiotemporal evolution of Hs, R2%, Rlow, and Rhigh for the four simulated storm events.

Hydrodynamic Forcing
Water level time series, estimated by the HD and HD-SW models, during the crossing of Hurricane Isidore at Telchac port are shown in Figure 4. The spatiotemporal evolution of the extreme water level along the stretch of coast for the four cases are shown in the following figures. The HD and the HD-SW models reproduced a peak water level of 2.28 m and 2.39 m (referred to the GGM06 Mexican Geoid), respectively. This 11 cm difference is related to the wave setup contributions to the seawater level. Given the lack of water level measurements, the hydrodynamic model was validated with a watermark of 3 m (Figure 4), surveyed at Telchac port during the passage of Hurricane Isidore, and referred to the GGM06. The difference between the maximum water level reproduced by the HD-SW model and the watermark is 0.61 m, which may be related to the swash amplitude contributions, not considered in the numerical model. Given the lack of observed data for validating the HD-SW model for Hurricane Gilbert, the same calibration parameters as for Hurricane Isidore were employed.    Given the storm tracks, the associated waves were propagated from east to west during the hurricanes' passage and from west to east during the CACS events. Therefore, spatiotemporal variability is different for each case. Hurricane Gilbert and Norte B were the events with the strongest and weakest wave energy, respectively.

Maximum Storm Impact along the Coast
The storm impact at the northern Yucatan coast is assessed by comparing morphology (D high , D low ) measurements with water level estimates. Figures 6 and 7 show the spatial distribution of D high , D low , the maximum values of the simulated R high (astronomic tide + storm surge + wave runup) (top panel) at each profile, and the corresponding storm impact regime for the simulated hurricanes and Norte events. The mean values of D high and D low along the 41.5 km stretch of coast are 2.53 m and 1.77 m, respectively, with a spatial variability (standard deviation) of 0.42 m and 0.41 m, respectively. The maximum values of R high for Gilbert and Isidore were 6.99 m and 3.98 m, respectively. The predominant storm impact regime was inundation for Gilbert and overwash for Isidore. Considering that hurricanes reaching the Yucatan peninsula are expected to increase their intensity by the end of this century [49], the conservation of the dune as a system defence play an important role to reduce the storm impact regime.

Maximum Storm Impact along the Coast
The storm impact at the northern Yucatan coast is assessed by comparing morphology (Dhigh, Dlow) measurements with water level estimates. Figure 6 and Figure 7 show the spatial distribution of Dhigh, Dlow, the maximum values of the simulated Rhigh (astronomic tide + storm surge + wave runup) (top panel) at each profile, and the corresponding storm impact regime for the simulated hurricanes and Norte events. The mean values of Dhigh and Dlow along the 41.5 km stretch of coast are 2.53 m and 1.77 m, respectively, with a spatial variability (standard deviation) of 0.42 m and 0.41 m, respectively. The maximum values of Rhigh for Gilbert and Isidore were 6.99 m and 3.98 m, respectively. The predominant storm impact regime was inundation for Gilbert and overwash for Isidore. Considering that hurricanes reaching the Yucatan peninsula are expected to increase their intensity by the end of this century [49], the conservation of the dune as a system defence play an important role to reduce the storm impact regime.
These results show the first approximation of the response of the Yucatan coast to tropical cyclones. For the present study, the wave runup parametrization of [16] was used; it was developed based on the numerical modeling of selected cases in a 30 years wave hindcast [48] in which Nortes and tropical cyclones were considered. However, the wave hindcast information [48] employed by [16] underestimated the hurricane-induced waves because of the use of wind reanalysis [50][51][52] as input. Furthermore, [16] did not take into account the storm surge in the modeling of the wave runup, suggesting that future efforts should focus on modeling hurricane-induced wave runup for the Yucatan coast. On one hand, the use of parametric wind models [42,53,54] may improve the modeling of hurricane's winds, which may be used as forcing for modeling hurricane-induced waves [26,55]. On the other hand, modeling of hurricane storm surge is needed for modeling levels during the pass of hurricanes. These two ocean's forcings may be used as input to a phase-resolving nonlinear nonhydrostatic model such as SWASH [47], for a better estimate of the wave runup. The modeling of the wave runup using numerical models is beyond the aim of this research. However, given that the continental shelf for Yucatan is wide (up to 245 km) and shallow, with a mild slope of -1/1000 [56], the wind-wave energy is expected to be low, and flooding can be dominated by storm surge during extreme events. Therefore, the Yucatan coast is more prone to be flooded from hurricane storm surge [2] than from hurricane-induced wave runup. The maximum values of Rhigh for Norte A and Norte B were 2.34 m and 2.32 m, respectively. Consequently, the predominant spatial regime was similar for both, mainly collision and overwash. Although the Rhigh for the hurricanes was found to be significantly larger than for the Nortes, the latter presents a more frequent overwash regime (up to 24.5 per year) [57].  (middle panel) and Isidore (lower panel) was calculated as swash = R high < D low ; collision = D high > R high > Dl ow ; overwash = R high > D high ; and Inundation = R low > D high .
These results show the first approximation of the response of the Yucatan coast to tropical cyclones. For the present study, the wave runup parametrization of [16] was used; it was developed based on the numerical modeling of selected cases in a 30 years wave hindcast [48] in which Nortes and tropical cyclones were considered. However, the wave hindcast information [48] employed by [16] underestimated the hurricane-induced waves because of the use of wind reanalysis [50][51][52] as input. Furthermore, [16] did not take into account the storm surge in the modeling of the wave runup, suggesting that future efforts should focus on modeling hurricane-induced wave runup for the Yucatan coast. On one hand, the use of parametric wind models [42,53,54] may improve the modeling of hurricane's winds, which may be used as forcing for modeling hurricane-induced waves [26,55]. On the other hand, modeling of hurricane storm surge is needed for modeling levels during the pass of hurricanes. These two ocean's forcings may be used as input to a phase-resolving nonlinear non-hydrostatic model such as SWASH [47], for a better estimate of the wave runup. The modeling of the wave runup using numerical models is beyond the aim of this research. However, given that the continental shelf for Yucatan is wide (up to 245 km) and shallow, with a mild slope of -1/1000 [56], the wind-wave energy is expected to be low, and flooding can be dominated by storm surge during extreme events. Therefore, the Yucatan coast is more prone to be flooded from hurricane storm surge [2] than from hurricane-induced wave runup. The maximum values of Rhigh for Norte A and Norte B were 2.34 m and 2.32 m, respectively. Consequently, the predominant spatial regime was similar for both, mainly collision and overwash. Although the Rhigh for the hurricanes was found to be significantly larger than for the Nortes, the latter presents a more frequent overwash regime (up to 24.5 per year) [57]. The highest spatial impact regime at each profile during the pass of Norte A (top panel) and Norte B (lower panel) was calculated as swash = Rhigh < Dlow; collision =Dhigh > Rhigh > Dlow; overwash =Rhigh > Dhigh. Note: after km = 35 some profiles did not get wet for Rlow nor outputs of the HD model. Then, the main contribution to the water level for those profiles was the wave runup. Because of the limitation on the interpolation of the topographic data to the mesh, the origin (shoreline) for those profiles are higher than 0.25m. Figure 7. Elevations of D hihg and D low , maximum R high, and associated R low reached at each profile for Norte A and Norte B (a). The highest spatial impact regime at each profile during the pass of Norte A (top panel) and Norte B (lower panel) was calculated as swash = R high < D low ; collision = D high > R high > D low ; overwash = R high > D high . Note: after km = 35 some profiles did not get wet for R low nor outputs of the HD model. Then, the main contribution to the water level for those profiles was the wave runup. Because of the limitation on the interpolation of the topographic data to the mesh, the origin (shoreline) for those profiles are higher than 0.25 m.

Spatiotemporal Variability of Storm Impact
The maximum values of R high for Norte A and Norte B were 2.34 m and 2.32 m, respectively. Consequently, the predominant spatial regime was similar for both, mainly collision and overwash. Although the R high for the hurricanes was found to be significantly larger than for the Nortes, the latter presents a more frequent overwash regime (up to 24.5 per year) [57]. Figures 8 and 9 show the temporal variation at each profile location of the wave runup (R 2%, (t, y), top row), R low (astronomic tide + storm surge + wave setup, second row), R high (third row) , and the impact regime (fourth row) for the hurricanes and CACS events, respectively. It is clear that since the highest water levels occurred during Gilbert, the strongest impact (inundation) on the coast was for that event.

Spatiotemporal Variability of Storm Impact
of the simulated period) along the entire barrier island. For Isidore the predominant regime was overwash (12 h, from 12:00 to 24:00) along the coast during the storm peak, except for some locations where the regime impact was inundation given the low dune topography. The predominant temporal impact regime for Norte A and Norte B were similar (overwash) during the storm peak but the duration for Norte A (4 h, from 07:00 to 11:00) was shorter than for Norte B (9 h, from 04:00 to 13:00), especially in the western region of the study zone (Yucalpetén and Progreso).  Regarding the predominant impact regime of both hurricanes, it was observed that while Gilbert was associated with inundation, Isidore produced overwash. This is mainly because the wind intensity for Gilbert was higher than for Isidore and consequently its associated storm surge, as well as the runup, were also higher. As a result, Gilbert opened several breaches in the western part of the study zone. For Norte A, the predominant impact regime during the storm peak was mainly collision (with some brief overwash periods in some locations) and for Norte B overwash periods were frequent. This difference is ascribed to tidal differences between such events: spring low tide during Norte A and spring high tide during Norte B. Besides, as expected, the propagation of waves, storm surge and as a consequence, the flood over coastal areas depend on the direction of motion of the storms: during the pass of Hurricane Gilbert and Isidore the first areas to flood or overwash were those around Telchac and the flood propagated to the west until reaching Yucalpetén. The opposite occurred during the pass of the Norte events.
The differences among the wave climate, storm surge, total water level and impact regime induced by hurricanes and Norte events are clear. However, the Norte events have a high occurrence Figure 9. Temporal variation at each profile location of the wave runup (R 2%, (t, y), top row), R low (astronomic tide + storm surge + wave setup, second row), R high (third row), and the impact regime (four row) for Norte A (left column) and Norte B (right column).
The maximum wave runup for Gilbert was 1.58 m, which is the maximum value allowed due to runup saturation in the runup parameterization of [16], and 1.56 m for Isidore. For Norte A and Norte B the runup was 1.35 m and 1.25 m, respectively.
The flood propagated from east (Km = 40) to west (km = 0) for both hurricanes and from west to east for the Norte events. The duration of flood for Gilbert was around ten hours (from 12:00 to 22:00 of the simulated period) along the entire barrier island. For Isidore the predominant regime was overwash (12 h, from 12:00 to 24:00) along the coast during the storm peak, except for some locations where the regime impact was inundation given the low dune topography. The predominant temporal impact regime for Norte A and Norte B were similar (overwash) during the storm peak but the duration for Norte A (4 h, from 07:00 to 11:00) was shorter than for Norte B (9 h, from 04:00 to 13:00), especially in the western region of the study zone (Yucalpetén and Progreso).
Regarding the predominant impact regime of both hurricanes, it was observed that while Gilbert was associated with inundation, Isidore produced overwash. This is mainly because the wind intensity for Gilbert was higher than for Isidore and consequently its associated storm surge, as well as the runup, were also higher. As a result, Gilbert opened several breaches in the western part of the study zone. For Norte A, the predominant impact regime during the storm peak was mainly collision (with some brief overwash periods in some locations) and for Norte B overwash periods were frequent. This difference is ascribed to tidal differences between such events: spring low tide during Norte A and spring high tide during Norte B. Besides, as expected, the propagation of waves, storm surge and as a consequence, the flood over coastal areas depend on the direction of motion of the storms: during the pass of Hurricane Gilbert and Isidore the first areas to flood or overwash were those around Telchac and the flood propagated to the west until reaching Yucalpetén. The opposite occurred during the pass of the Norte events.
The differences among the wave climate, storm surge, total water level and impact regime induced by hurricanes and Norte events are clear. However, the Norte events have a high occurrence in the southern coast of the Gulf of Mexico, including the study zone, with an annual mean ranging from 16 events [58] to 24.5 [57], depending on how CACS are defined. Therefore, it is important to take into account the coastal response to CACS events under scenarios of dune degradation. For instance, if Norte A had occurred during spring high tide the predominant impact regime may have switched from collision to overwash or even inundation.
The impact regime of the coast to storm events depends on both the peak sea water level and the dune height. The higher the dune, the lesser is the impact the coast suffers during the pass of storm events. Most extreme conditions occurred during Hurricane Gilbert when the inundation regime occurred in the entire study area for around 10 h. High storm impacts during hurricane conditions are mainly associated with more intense winds than dune morphology. Thus, the dune height elevation did not play a role in coastal protection with such extreme conditions. The most affected areas during both hurricanes and Nortes occurred near Progreso (km = 0 to km = 15), where population density is larger. Less impact occurred in those areas with well-preserved dunes due to less anthropization (settlements and coastal structures). However, during the pass of Gilbert and Isidore, Telchac experienced the highest values of R high because this area experienced the maximum winds during the crossing of those hurricanes.

Conclusions
We investigated the storm impact associated with historical events in the northern Yucatan Peninsula using LIDAR information, a numerical hydrodynamic model, and an empirical runup parameterization. The analysis of the results suggests that the predominant impact regimen associated with hurricanes and Norte events ranged from overwash to inundation and from collision to overwash, respectively; i.e., less frequent hurricane events have a greater storm impact than recurrent Norte events. Spatiotemporal evolution of waves and flood along the stretch of coast follows the storm path and is from east to west for the two hurricane events; it is the opposite during the Norte events. The western part of the study area presents the higher impact regime during Norte events owing to the less preserved coastal dunes in those areas. This is highly correlated with urban expansion in this part of the coast. Therefore, the implementation of coastal dune protection initiatives is needed to make the Yucatan coast more resilient to Nortes and moderate hurricane events.  Acknowledgments: W.R. thanks DHI Water & Environment for facilitating a student license of MIKE 21 SW and HD, Gonzalo U. Martín-Ruiz, Juan Alberto Gómez-Liera and José López-González for their help with computer resources and technical support, as well as Cecilia Enriquez-Ortiz for providing the bathymetric data of the Chelem lagoon. The authors gratefully acknowledge the support of NVIDIA Corporation for the donation of the Tesla K40 GPU used for this research. Gabriela Medellín is acknowledged for discussion on the use of wave run up parameterization. Daniel Behm is acknowledged for reviewing this paper. Finally, authors acknowledges to the three anonymous reviewers that greatly contributed to the improvement of this paper.