Assessment of coastal flooding and associated hydrodynamic processes on the south-eastern coast of Mexico, during Central American cold surge events

Coastal flooding in the northern Yucatán Peninsula is mainly associated with storm surge events triggered by high-pressure cold front systems. This study evaluates the hydrodynamic processes of the Chelem lagoon, Mexico and the flooding threat from cold fronts for the neighbouring town of Progreso. A 30-year water-level hindcast (excluding wave set-up) was performed because of the lack of long-term tide gauge records. In order to assess the relative contribution from wave set-up and residual and astronomical tides to total flooding, the two worst storm scenarios in terms of maximum residual tide (Event A) and maximum water level (Event B) were simulated. Numerical results suggest that during Event A the wave set-up contribution reaches 0.35 at the coast and 0.17 m inside the lagoon, and these values are smaller for Event B (0.30 and 0.14 m, respectively). Results of the effect of the tidal phase on wave set-up and residual sea level show that (i) the wave set-up contribution increases during ebb tide and decreases during flood tide at the Chelem inlet, (ii) the residual tide is larger (smaller) near low (high) or receding (rising) tide, and (iii) maximum flooding occurs when the storm peak coincides with rising or high tide. The numerical results confirm the important role of wave set-up on the assessment of coastal flooding in micro-tidal coastal environments.

1 Introduction

The Yucatán Peninsula coast, located in south-eastern Mexico, is prone to coastal flooding due to both its geographical location and geological characteristics. On the one hand, extreme meteorological phenomena such as Central American cold surge (CACS) events (Appendini et al., 2018) and tropical cyclones are common in this area (Posada-Vanegas et al., 2011; Meza-Padilla et al., 2015; Rey et al., 2016). On the other hand, the region is characterized by a wide and shallow continental shelf, a low-lying coast and the presence of semi-enclosed back-barrier water bodies (lagoons, shelter ports, and wetlands). Thus, forcing agents and the coastal features enhance the physical vulnerability of this region.

Tropical cyclones occur during summer and early autumn months and are responsible for the worst coastal flooding events in the Yucatán Peninsula. However, these types of events have a low occurrence along the northern Yucatán coast, which has been exposed to only 25 cyclones in the past 150 years (until 2001), i.e. an average of only 0.16 events per year (Rosengaus-Moshinsky et al., 2002). CACS (locally knows as Nortes) events occur during late autumn and winter, with an annual mean ranging from 16 events (Reding, 1992) to 24.5 (Appendini et al., 2018), depending on how a CACS is defined. Following Reding (1992) and Schultz et al. (1998), the present study defines a CACS as an anticyclonic movement of a cold mass of air that originated in North America, which penetrates equatorward to at least 20^∘ N latitude. However, fewer efforts have been devoted to investigating the hazards associated with CACS events even though the annual frequency of CACS events is higher than hurricanes in the north-west of the Yucatán Peninsula. Thus, this study focuses on the effect of CACS events on coastal flooding.

The intensity of the waves and storm surges associated with CACS events depends on the magnitude of the shear stress on the sea surface, the fetch, and the duration of the events. The storm surge is enhanced by the wind shear stress on the sea surface and perturbations in the atmospheric pressure (Lin and Chavas, 2012). Since the inverse barometer effect on storm surge is small during low-pressure storm systems (Massey et al., 2007), the storm surge is mainly driven by wind stress, especially in shallow coastal waters (Flather, 2001). Moreover, considering that CACS are high-pressure systems, the storm surge is essentially driven by the direct wind effect. Depending on the shape of the coast (concave or convex), the coastal bathymetry, the extent of the continental shelf, as well as the direction and duration of the incident wind, a higher or smaller storm surge will occur. Additionally, seasonal variations in the mean sea level and pressure gradients induced by littoral currents and other factors (e.g. steric basin-scale anomaly, astronomical annual tidal component) may cause an increase in sea level, with maxima for this region in September–October (Zavala-Hidalgo et al., 2003).

Wave set-up can be important for the accurate estimation of the extreme water levels (Vousdoukas et al., 2016). During the wave-breaking process the kinetic energy is converted, to a great extent, to a quasi-steady potential energy, generating a water surface gradient to balance the onshore component of the momentum flux due to the presence of waves (Dorrestein, 1961; Longuet-Higgins and Stewart, 1963). Consequently, an increase in the water level along the shoreline is generated (Smith et al., 2000; Dodet et al., 2013), as well as wave-induced currents. Moreover, when inlets or port entrances are present, these processes play an important role in modifying the inlet and lagoon hydrodynamics and morphodynamics, which are regularly forced by the astronomical tide and fresh water input from springs. In fact, given that wave breaking in the inlet ebb shoal or further inside the main channel depends on the water depth, both the wave-breaking-induced acceleration and the wave set-up are tidally modulated (Olabarrieta et al., 2011). However, the wave set-up contribution is often neglected due to the computational time cost required for its modelling, especially in real-time forecasting of hurricane storm surges when prompt simulation results are required for preparing evacuation plans. Lin et al. (2012) simulated a set of over 210 extreme surge events with the ADCIRC model coupled with the SWAN wave model and found that the wave set-up accounted for less than 1.5 % of the surge for four locations around New York Harbor. They suggested that the limited wave set-up in this area may be related to the fact that large ocean waves break before entering the New York Harbor, and also because the near-shore wave breaking may not be captured by the large-scale computational mesh of models. These authors also found that the run time for the surge–wave coupled simulations was 1 order of magnitude larger than for the simulations, accounting only for the surge. For the Yucatán coast, flood hazard studies have not taken into account the wave set-up contribution and hence the present study aims to investigate its relative importance in coastal flooding during CACS events.

According to Merz et al. (2007), a flood hazard is defined as the probability of the induced potential damage caused by a flood over a determined area and period. The latter depends on several factors, such as the maximum water level reached, the flow velocity, the flood duration, the speed of sea-level rise, and the flood frequency. Flood hazard assessments are normally performed based on historic or synthetic flood data (Lin et al., 2010; Zachry et al., 2015). However, since wind reanalysis data sets have become available, such as the North American Regional Reanalysis (Mesinger et al., 2006) and the Climate Forecast System Reanalysis (CFSR; Saha et al., 2010), they have been used to force hydrodynamic models to generate sea-level reanalysis. Storm surge hazard due to hurricanes along the Yucatán coast has previously been studied (Posada-Vanegas et al., 2011; Meza-Padilla et al., 2015; Rey et al., 2016), as has the flood risk caused by hurricanes (Rey et al., 2016). Given that the wind fields during tropical cyclones are underestimated in these reanalyses (Swail and Cox, 2000), the sea level during these events is also underestimated by hydrodynamic models.

Long records of raw sea-level data are scarce for the northern coast of the Yucatán Peninsula (barely 5 years of high-quality historical raw data are available for the Progreso area: 1979–1984). Therefore, a sea-level hindcast was developed as part of this work in order to overcome this lack of data. The aim of this study is to assess flood hazards caused by CACS events in the northern Yucatán Peninsula with emphasis on evaluating the relative contribution of storm surge and wave set-up and the role of the astronomic tidal phase on these processes in Progreso, Yucatán.

2 Study area

The study area is located on the northern coast of the Yucatán Peninsula, Mexico (Fig. 1). The study focuses on (i) the town of Progreso, which is the most urbanized, populous, and economically important coastal city on the northern coast of Yucatán, and (ii) the back-barrier of Chelem lagoon located behind Progreso. The maximum ground elevation in Progreso is 2.1 m above Mexican Geoid GGM06, within a low-lying coast (average elevation in Merida, the state capital located 27 km inland, is only 8 m). The peninsula averages 10 m above sea level with only a small prominent sierra in the centre, where the maximum altitude reaches 150 m (Stringfield and LeGrand, 1974). The northern Yucatán coast is mostly sandy (85 % of its length), and 67 % is formed by coastal lagoons and barrier islands (Cinvestav, 2007).

Figure 1 Location map indicating the study zone and the town of Progreso.

The climate in the zone is characterized by three seasons: (i) warm and dry (March–May), (ii) rainy (June–October) and (iii) cold with occasional showers (November–February) associated with the CACS passage (Medina-Gómez and Herrera-Silveira, 2009; Schmitter-Soto et al., 2002). Predominant winds are associated with sea and land breezes from the NE and SE, which are more frequent and intense during spring (Figueroa-Espinoza et al., 2014). During winter months, strong northerly winds interact with the maritime tropical wind from the Caribbean, causing the CACS events (Schultz et al., 1998), which are distinct from the north-easterly sea breeze because of their duration (Reding, 1992), atmospheric pressure, and air temperature. The mean annual precipitation in this region varies between 444 and 1227 mm. Progreso is located in the driest area of the peninsula (INEGI, 2002).

In terms of bathymetry and energy dissipation, the karstic continental shelf is exceptionally wide (up to 245 km), shallow, and has a mild slope $-\mathrm{1}/\mathrm{1000}-$ (Enriquez et al., 2010). Therefore, the wind–wave energy is low (Lankford, 1976) near the coast, mainly due to bottom friction dissipation, while the storm surges are amplified. Long-term wave climate analysis in this zone was performed by Appendini et al. (2013) by means of a 30-year wave hindcast for the Gulf of Mexico and Caribbean Sea, showing that the continental shelf dissipates storm waves from distances far offshore, of the order of tens of kilometres. This, together with the mixed tide with a diurnal dominance and a small neap (spring) range of 0.1 m (0.8 m) (Cuevas-Jiménez and Euán-Ávila, 2009), provides the coastal zone with a low-energy regime during practically all circumstances (Salles et al., 2013), except at the eastern part of the peninsula, where the wave energy increases due to a narrower continental shelf (Appendini et al., 2012).

Appendini et al. (2012) reported a net westward longshore sediment transport along the entire northern Yucatán coast, ranging between 20 000 and 80 000 m³yr⁻¹, except west of Holbox (Fig. 2), where longshore transport direction is reversed at this location. The dominant westward longshore transport suggests an extremely sensitive shoreline to artificial littoral barriers (e.g. groins and jetties).

Figure 2 Computational domains and topo-bathymetry for the hydrodynamic and coupled models (polygon in red) as well as for the spectral wave model (entire Gulf of Mexico).

The northern coast of the Yucatán Peninsula is the first land mass to interact with the CACS events after they have crossed the Gulf of Mexico, where the wind speed can reach up to 30 m s⁻¹, according to measurements from the National Data Buoy Centre (NDBC) stations and Schultz et al. (1997). During the passage of CACS events, both erosion and flooding processes might occur, but given the low-lying coast, flooding represents a greater hazard (Mendoza et al., 2013), and given their relatively frequent occurrence, their socio-economic impact is high. For instance, during the passage of CACS events, both port and oil industry activities are affected in the southern Gulf of Mexico, which translates to economic loss. For example, Cold Front 4, which happened on 23 October 2007, caused major flooding and coastal structural damage in Villahermosa, in the nearby state of Tabasco, and caused an economic loss estimated at USD 2.45 billion (López-Méndez, 2009). Another case of flooding induced by a CACS occurred on 28 November 2013 in the northern coast of Yucatán, where the Dzilam de Bravo community was flooded, causing damage to around 80 houses and coastal structures (Cob-Chay et al., 2013).

3 Numerical model

3.1 Model description

The hydrodynamic model used for this study was MIKE 21, developed by DHI Water & Environment, which resolves the two-dimensional shallow water equations – the depth-integrated incompressible Reynolds-averaged Navier-Stokes equations invoking the assumptions of Boussinesq and Hydrostatic pressure (DHI, 2014a). Wetting and drying are included in the hydrodynamic model following the work of Zhao et al. (1994) and Sleigh et al. (1998). This model has been successfully used in recent scientific studies (Strauss et al., 2007; Appendini et al., 2014; Meza-Padilla et al., 2015).

The wave model used to compute the wave conditions and associated radiation stresses was the MIKE 21 third-generation spectral wave model. This model has been used for several spectral wind–wave modelling applications (Strauss et al., 2007; Appendini et al., 2013, 2015). For more detailed information about source terms, governing equation, time integration, and model parameters, readers are referred to Sørensen et al. (2004) and to the scientific manual documentation for the spectral wave model (DHI, 2014b).

3.2 Model set-up

3.2.1 Hydrodynamic model set-up

The hydrodynamic model was used in order to obtain 30-year currents and sea-level hindcast, not accounting for waves due to the high computational cost. The 30-year sea-level hindcast was developed as the basis for the extreme-level analysis, which is not possible from measurements due to the lack of long-term tidal gauge records. Unfortunately, the computational cost prohibits the modelling of coupled waves and hydrodynamics for such a long period (30 years). For instance, for a given period of 3 weeks, and the computational domain used for the hindcast shown in Fig. 2 (polygon in red), the computational time for the uncoupled (no waves) model was 12 h but up to 2 weeks for the coupled model. Therefore, coupled modelling was only considered for selected cases in order to include wave set-up and wave–current interaction (see Sect. 3.3).

The boundary conditions for the uncoupled model were treated following Enriquez et al. (2010):

• The Yucatán channel boundary was forced with a mean profile of the Yucatán current, which is constant in time and varies in space based on the results reported by Abascal et al. (2003). Part of this current has been attributed to mesoscale eddies, which are observed in the eastern Caribbean basin, the Cayman Sea, and western Caribbean passages (Athié et al., 2011).

• The Gulf of Mexico boundary was forced with the astronomical tide varying along the boundary; values were extracted from the global tide model (Andersen, 1995), which represents the major diurnal ($K_{\mathrm{1}},O_{\mathrm{1}},P_{\mathrm{1}}$, and Q₁) and semi-diurnal tidal constituents (M₂, S₂, N₂, and K₂) with a spatial resolution of 0.25 × 0.25^∘ (DHI, 2014c).

• The Campeche (western) boundary was considered open (sea level equal to zero).

• The southern boundary (land) was forced with a constant Yucatán aquifer discharge of 2.701 × 10⁻⁴m³s⁻¹m⁻¹ as reported by Weidie (1985).

On the surface the model was forced with wind and pressure fields from the CFSR database, which has a global atmospheric resolution of ∼ 38 km (T382) with 64 levels extending from the surface to 0.26 hPa. The global ocean resolution is 0.25^∘ at the equator, extending to 0.5^∘ beyond the tropics, with 40 levels from the surface to a depth of 4737 m. The National Centers for Environmental Prediction (NCEP) has created time series products at hourly temporal resolution by combining either (1) the analysis and 1 to 5 h forecasts, or (2) the 1 to 6 h forecasts, for each initialization time. When using these data products, it has to be kept in mind that only the 00:00, 06:00, 12:00, and 18:00 UTC fields are actual analyses, while the in-between hourly data are model forecast events (Saha et al., 2010, 2014).

Given that the spatial resolution of the CFRS grid is not regular, and the hydrodynamic model only takes wind and pressure data varying in space from a regular grid, CFSR wind and pressure fields were linearly interpolated from a T382 Gaussian grid resolution to a regular grid with spatial resolution of 0.3125^∘, which is coincident with the longitude of the T382 grid and close in latitude for the Gulf of Mexico. We assumed this resolution to be adequate to reproduce the CACS storm surge based on the work of Appendini et al. (2013), who showed that the resolution of NCEP–NCAR (National Centre for Atmospheric Research), ECMWF ERA-Interim (European Centre for Medium Range Weather Forecasts – European Reanalysis), and the North American Regional Reanalysis is sufficient for a wave model of CACS over the Gulf of Mexico. Indeed, given that CFSR data are superior to the above NCEP reanalyses due to (a) a finer resolution, (b) an advanced assimilation scheme, and (c) atmosphere–land–ocean–sea-ice coupling, it is expected to be a good compromise for this application. Moreover, the hourly resolution of CFSR allows this data set to capture extremes, such as the storm peak, which other reanalyses may miss, according to Sharp et al. (2015). These authors found a good correlation between the hourly CFSR data set and both onshore and offshore in situ measurements for the UK. For instance, NCEP FNL (Final), ECMWF ERA-Interim, and NCEP–NCAR provide data at 6 hourly intervals, which might not be able to capture storm peaks (Jørgensen et al., 2005). Therefore, when using these wind fields as forcing in hydrodynamic models, the maximum flooding areas may be underestimated.

The MIKE 21 hydrodynamic model uses a dynamic time step to optimize simulation speed while ensuring numerical stability. Hence, the time step may change during the simulation (large time step under calm conditions, smaller time step when flow becomes stronger). The user is allowed to set the minimum and maximum time steps in the model set-up. The actual dynamic time steps used are found to be in the range of 5 to 7.5 s. Since the time step for the CFSR is 1 h (3 orders of magnitude longer than the hydrodynamic model time step), the hydrodynamic model interpolates the CFSR data linearly to its own time step.

The bathymetry was extracted from the ETOPO1 database and was complemented with higher-resolution bathymetric data from 9 km long transects every 4 km along the coast. In addition, high-resolution topography (1 m spatial resolution) from a 2011 lidar survey of the entire town of Progreso was used. After a calibration process comparing model results with sea-level measurements during three CACS events (not shown) in Progreso, the bottom friction was defined using a constant Manning coefficient of 0.02, which, according to Arcement and Schneider (1989), corresponds to the average mean grain size (d₅₀) of the Yucatán sand beaches reported by Mendoza et al. (2013). For the horizontal eddy viscosity (Smagorinsky formulation) a constant coefficient of 0.28 was applied. The wind friction (Cd) was estimated based on the Garratt (1977) formulation modified by Lin and Chavas (2012) based on Powell et al. (2003) and further calibrated in this study, which consisted of varying (increasing and decreasing) the Lin and Chavas (2012) values and selecting the best combination that resulted in the smallest water sea-level error. The Cd used varies linearly with the wind speed, as suggested in other studies (Bryant and Akbar, 2016). The mesh selected for the hydrodynamic model is the result of a sensibility analysis of the domain size (not shown) that determines the size (Blain et al., 1994; Morey et al., 2006; Kerr et al., 2013) at which the model adequately reproduced the sea level recorded by a tide gauge at Progreso.

3.2.2 Wave and coupled model set-up

Once the sea-level hindcast was performed, two of the most extreme flooding events were identified from this data set. For the two identified events (A and B mentioned below), the wave set-up contribution was taken into account by coupling the hydrodynamic model and the spectral wave model, i.e. coupled model.

The same computational domain (Fig. 2, polygon in red) was used for both the hydrodynamic and the coupled models. The hydrodynamic model kept the same type of boundary conditions mentioned in Sect. 3.2.1. Given the unknown swell wave conditions for the ocean boundaries, a wave model was implemented for the entire Gulf of Mexico (Fig. 2) to reproduce distant wave climate accurately in order to use this as the forcing at the ocean boundaries in the coupled model. The wave model set-up was based on the study by Ruiz-Salcines (2013), who calibrated the model for mean and extreme conditions in the Gulf of Mexico and Caribbean Sea. This calibration was also used in Appendini et al. (2017).

3.2.3 Selection of CACS and simulated cases

In this study the term “storm surge” is used when referring only to the meteorological contribution to the total sea level; otherwise, we refer to “residual tide”, which may contain storm surge, tide–storm surge interaction, harmonic prediction errors, and timing errors (Horsburgh and Wilson, 2007). In other words, the residual tide is the total water level from the sea-level hindcast (without wave set-up) minus the astronomical tide.

From the 30-year sea-level hindcast, the largest yearly events were identified (excluding the sea level generated by tropical cyclones, given that they are underestimated), considering two criteria: (a) the largest residual tide that, being the less predictable tidal component, is relevant to flood hazard prediction, and (b) the sea surface elevation associated with each event from (a) astronomical tide plus residual tide, which is commonly considered the main parameter used to estimate the flood hazard. Clearly, the analysis of the residual sea level is crucial, since this is affected by greater uncertainties than the astronomical tide. While astronomical tide is characterized by periodic oscillations, the residual tide has high variability (Mel et al., 2014). Following the CACS identification method proposed by Reding (1992), all 30 events belonging to (a) were found to be CACS events. From these 30 CACS events, the two largest were selected for simulation and analysis: (a) the event with the largest residual tide (Event A, whose peak occurred during receding tide), which hit the peninsula from 12 March 1993 at 16:00 to 13 March at 23:00, and (b) the event with the largest sea surface elevation (Event B, whose peak occurred during rising tide), which occurred from 25 December 2004 at 15:00 to 26 December at 09:00.

Event A presented not only the largest residual tide but also the highest wind intensity and duration of the total 30-year hindcast. Bosart et al. (1996) and Schultz et al. (1997) called this event “the 1993 super storm cold surge, also known as the storm of the century”, which originated over Alaska and western Canada, and brought northerlies exceeding 20 m s⁻¹ and temperature decreases up to 15 ^∘C over 24 h into Mexico and Central America. Schultz et al. (1997) studied this CACS event in detail due to its exceptional intensity over the Gulf of Mexico, the important role that convection played in the incipient cyclogenesis, the planetary-scale antecedent conditions, and the merger of two short-wave troughs in the westerlies contributing to the extreme cyclogenesis.

Figure 3 NDBC station 42001 data: wind speed (10 m above sea level) and direction for events A (a) and B (b).

At the peak of the storm in Event A, the astronomical tide and residual tide were −0.35 and 1.14 m, resulting in a total sea level of 0.79 m. During the peak of the storm of Event B, the astronomical tide and residual tide were +0.44 and 0.72 m, resulting in a total sea-level elevation of 1.16 m, i.e. 0.37 m higher than the total sea level during Event A at 3 km offshore of Progreso. Figure 3 presents wind speed and direction from the 42001 NDBC buoy (see location in Fig. 2) during events A and B. Throughout Event A, the predominant wind direction (azimuth) was 315^∘, whereas for B it was 340^∘, i.e. closer than normal to the coast. Moreover, the maximum wind speed was similar for both events, but the duration of wind speeds higher than 20 m s⁻¹ was longer for Event A (11 h compared to 3 h for Event B). This suggests that the duration of the storm is a predominant factor in the generation of storm surges.

In order to investigate the relative contribution of the (i) storm surge and (ii) wave set-up during the flooding episodes, four different scenarios were implemented for events A and B:

• Case 1 (C1): this configuration considered only the hydrodynamic model (see details in Sect. 3.2.1).

• Case 2 (C2; storm surge): the hydrodynamic model was forced at the surface with pressure and wind fields from the CFSR database. Only the storm surge contribution was evaluated.

• Case 3 (C3; wave set-up): the hydrodynamic model was forced only with the radiation stresses obtained from the coupled model (wave–current interactions; see Sect. 3.2.2), obtaining only the wave set-up contribution and wave-induced currents.

In both cases, C2 and C3, the ocean boundaries surrounding the mesh formed by the polygon in red in Fig. 2 were open, and the coastline boundary was closed.

• Case 4 (C4; total sea surface elevation, TSSE, and total currents): the coupled model was used to investigate the contributions from the storm surge, wave set-up, astronomical tide, and the Yucatán current, for events A and B, to assess flood-prone areas in Progreso. The boundary conditions were set as mentioned in Sect. 3.2.2.

In addition, the role of the astronomic tidal phase in the flood hazard and the Chelem lagoon hydrodynamics was investigated. For this, three additional numerical experiments were performed as shown in Table 1.

These four “tide scenarios (TS; see Table 1)” were also used with Case 1 and Case 4 to study the variation in residual tide and maximum flood as a function of the astronomic tidal phase for Event A, respectively.

Table 1 Tide scenarios.

The dashed lines in Fig. 4 show the forcing tide for the scenarios mentioned above, which varied in phase but had the same amplitude. The first peak of TSSE for TS3 (time t1 in Fig. 4) was taken as a reference for varying the phase in the other scenarios: TS1 is 12 h ahead of TS3. The phase shift for TS2 and TS4 was set so that the water level is zero at time t1, during receding and rising tides, respectively. Both t1 (flood) and t2 (ebb) were the times used for assessing the wave set-up contribution inside the Chelem lagoon.

Figure 4 Tide forcing for the four scenarios at Progreso, for Event A (dashed lines); time series of TSSE at Progreso (continuous lines) associated with each scenario. For each tide scenario see Table 1.

3.3 Model validation

The hydrodynamic model results were validated with data from a tidal gauge in Progreso (21.3033^∘ N, −89.6667^∘ W), to perform a 30-year (1979–2008) sea-level hindcast for the northern Yucatán Peninsula. Figure 5 shows the measured sea level and hydrodynamic simulation results for the two storm events with the greatest residual tide in the 5-year tide gauge record. In general, there is good agreement for the sea surface elevation during such events. For the event shown in the top panel of Fig. 5, the Pearson correlation is 0.78 and the root mean squared error (RMSE) is 0.1 m, which correspond to 20.9 and 16.6 % of the measured and modelled sea-level range, respectively. For the other event shown in the lower panel of Fig. 5, the Pearson correlation is 0.87 and the RMSE is 0.17 m, which correspond to 16.6 and 18.3 % of the measured and modelled sea-level range, respectively. Based on the aforementioned model–data comparison the model validation is considered acceptable. For instance, for the SLOSH model, which is the official model used by the National Hurricane Centre to provide real-time hurricane storm surge data (Massey et al., 2007), the accuracy of the predicted surge heights is ±20 % when the tropical cyclone is adequately described (Jelesnianski et al., 1992).

Figure 5 Hydrodynamic model validation with tide gauge data from Progreso port in Yucatán. (a) Validation for a CACS event in 1979. (b) Validation for a CACS event in 1982.

In addition, both the CFSR winds and the wave model results (significant wave height and peak period) for events A and B were validated with measurements from the NDBC stations 42001, 42002, 42003, and 42055 (not shown). Figure 6 shows the validation for Event A with buoy 42001 data. Significant wave heights, H_s, and peak periods, T_p, from the wave model exhibited good Pearson's correlations (0.90 and 0.79) as well as a RMSE of 0.68 m and 1 s for H_s and T_p. Furthermore, a good Pearson's correlation (0.91) between the CFSR wind reanalysis and wind measurements from the same NDBC station was found (left bottom panel in Fig. 6).

Figure 6 Wave model and CFSR wind validation for Event A: (a) significant wave height, (b) wave peak period, (c) wind speed 10 m, and (d) wind direction.

4 Results

4.1 Sea-level hazard assessment

To assess the sea-level hazard, an extreme analysis was performed with the 30 largest yearly CACS events selected in Sect. 3.3: (a) the 30 largest annual events in terms of residual tide and (b) the sea surface elevation associated with each event from (a). In both cases the extreme analysis was performed using model data from a point situated 2 km offshore (5 m depth), where the Progreso tide gauge is located. Both data sets were fitted to the generalized extreme values (GEVs) distribution probability function H (Ho et al., 1976; Jenkinson, 1969), using the following equation:

$$\begin{array}{}\text{(1)}& H\left(x,\mathit{\mu },\mathit{\psi },\mathit{\xi }\right)=\exp {\left\{-\left(\mathrm{1}+\mathit{\xi }{\displaystyle \frac{x-\mathit{\mu }}{\mathit{\psi }}}\right)\right\}}^{-\mathrm{1}/\mathit{\xi }},\end{array}$$

where x is the data to fit, μ is the location parameter, ξ is the shape parameter and ψ is the scale parameter.

By means of the maximum likelihood estimation method the following function parameters were found: (a) μ=0.469, ξ=0.189, and ψ=0.105 when considering the events with the largest residual tide (where Event A is the largest), and (b) μ=0.526, $\mathit{\xi }=-\mathrm{0.295}$, and ψ=0.261, when considering the sea-level elevation associated with events from (a) (where Event B is the largest).

On the one hand, when considering only the residual tide in the extreme analysis, the resulting return periods for events A and B were 67 and 7 years (Fig. 7a), respectively, given that the residual tide for Event A was larger than for Event B (1.14 m vs. 0.72). On the other hand, when considering the total sea surface elevation, i.e. including the residual and astronomical tide, the return periods for events A and B become 3 and 78 years (Fig. 7b), respectively. This is due to the fact that Event A, for which the residual sea level was the largest of the 30-year hindcast, happened near low tide, while Event B occurred during spring high tide. Therefore, including the astronomical tide in the extreme analysis is crucial for estimating the return periods. Moreover, in Sect. 4.2 and 4.3, and for locations near semi-enclosed coastal bodies – as is the case for Progreso – the importance of not only considering the tidal phase but also the local wind direction when determining the flood hazard inside the lagoon is discussed. Wind set-up in this semi-enclosed body can be of major significance (Carniello et al., 2005), and is crucial for the aim of this study, since the town of Progreso has its southern limits bordering the Chelem lagoon.

Figure 7 Yearly maximum residual tide heights adjusted to the GEV distribution probability function for Progreso (2 km offshore). The solid curve is the fitted curve. The dots correspond to the 30 events with (a) the yearly maximum residual tide, and (b) the yearly maximum sea-level elevation. The dashed curves are the 95 % confidence limits.

4.2 Contribution of storm surge and wave set-up to flooding for events A and B

Inside the lagoon, the wind set-up (and hence the wind direction and storm duration) and the hydrodynamics at the inlet play an important role. Figure 8 presents maps of the highest wind stress and wave height values for events A and B. During Event A (storm shown in Fig. 6) wave heights of 4–5 m occurred 2.5 km offshore. In contrast, the same wave height values occurred 5 km offshore for Event B. Moreover, during Event A, the wind direction and wind stress (top left panel), propagated the waves (top right panel) from west to east inside the lagoon, leading to a large wave height (around 0.9 m) within the lagoon in the south-eastern part of Chelem. During Event B, the wave height (bottom right panel) inside the lagoon was weaker than during Event A, in part due to the weaker wind stress (lower left panel) over the lagoon. These variations caused significant differences in the induced wave set-up and wind set-up inside the lagoon for each event as shown in Fig. 9.

Figure 8 Highest wind stress values during Event A (a) and Event B (b), and highest wave height values for Event A (c) and Event B (d).

Figure 9 Main contributions for maximum flooding for Progreso during events A and B: (a) maximum storm surge for Event A, (b) maximum storm surge for Event B, (c) maximum wave set-up for Event A, and (d) maximum wave set-up for Event B. The black line in panel (c) is a transect discussed in Sect. 4.4.

The maximum values of the storm surge and wave set-up for events A and B are shown in Fig. 9. The top panels (left and right) show the model results of the maximum storm surge. During normal weather conditions there is a predominant littoral current along the northern coast of the Yucatán Peninsula from east to west (Enriquez et al., 2010), which is attributed to sea breezes (Torres-Freyermuth et al., 2017), wind events from the south-east, and in part to the Yucatán current which floods the Yucatán Shelf from the east (Abascal et al., 2003). However, during CACS events, both (i) winds (from the north-west and north) produce a large and northerly shear stress on the sea surface, and (ii) pressure gradients due to atmospheric pressure perturbations drive water towards the peninsula. As a consequence, the predominant longshore current switches direction from west to east and leads to an increase in the sea level along the northern Yucatán coast and hence inside the coastal lagoons due to the orientation of the coast (see Fig. 1). It is evident from Fig. 9a and b that the storm surge during Event A (the event with the longest duration) was larger than Event B both outside and inside the lagoon. This suggests that the duration of the storm is an important factor in the generation of storm surges. It can also be seen from these two panels that the surge is greater in the eastern portion of the lagoon due to the direction of the wind stress – which accounts for a significant amount of the surge – and the corresponding wind set-up, in addition to the water volume inflow through the inlet.

Regarding the wave set-up (Fig. 9c and d), Event A presented larger values (0.35 m along the Progreso coast and roughly 0.17 m inside the Chelem lagoon) compared to Event B (0.3 and 0.14 m, respectively). These differences seem to be related to the tidal phase and offshore wave energy (Dodet et al., 2013) as well as to the wave direction (Guza and Feddersen, 2012) when the storm peak took place in each event. The wave set-up for Event A and Event B at the inlet contributed up to 19 and 14.5 % of the TSSE, respectively. This shows the importance of taking into consideration the wave set-up contribution in the flood hazard assessment of coastal lagoons.

4.3 Flood-prone areas for events A and B

From events A and B (Fig. 10), the most affected area on the sea side seems to be the stretch of coast from the Chelem inlet to Progreso Pier, partly due to its concave shape. Flooding in this area was larger during Event B (bottom panel) because this CACS event hit the peninsula during high tide and the TSSE was larger. Inside the Chelem lagoon, flooding is in general larger than along the open coast for both events. This is due to its limited capacity to regulate increased volumes of water flooding through the inlet as well as the wind set-up, which can be larger for certain wind directions than for the open coast. In fact, even if the sea surface elevation on the sea side is lower for Event A than for Event B (as shown in Sect. 3.3), a higher water level inside the lagoon, particularly in its eastern sector, was found during Event A. This is related to a limited tidal range and shallow water depth in the lagoon, resulting in a strong correlation between the storm surge levels and the wind stress, as well as a stronger wave set-up, resulting in larger flooded areas of Progreso. In terms of the specific flooding areas in the town of Progreso, the model wetting-and-drying algorithm performs skillfully inland, in particular in the eastern part of the back-barrier lagoon. The total number of flooded city blocks in Progreso during Event A was 157 (25 % of the town surface area), distributed as follows: 8 along the Progreso beach and 149 along the eastern lagoon shores. For Event B, the total number of flooded blocks were 110 (18 % of the town surface area), with 18 corresponding to the sea side and 92 on the lagoon side. Since assessing their impact on infrastructure is one the main objectives of this research, the block was chosen as the unit with which to show the inhabited flood-prone areas. In this regards, quantifying areas without inhabitants is beyond the aim of this study. These results show that the most flood-prone areas in Progreso are those located in the south-eastern area, bordering the Chelem lagoon, and not along the coast, mainly because of the semi-enclosed nature of the lagoon and the wind and wave set-up associated with storms.

Figure 10 Maximum flooding over Progreso during (a) Event A and (b) Event B.

Figure 11 Modulation of the wave set-up for the ebb/flood current at the Chelem inlet during Event A. (a) TSSE, storm surge and wave set-up for times t1 and t2. (b) V Total component velocities and maximum wave height for t1 and t2.

4.4 Hydrodynamics and wave set-up at the inlet of the Chelem lagoon

In order to further investigate the evolution of the wave set-up through the inlet during the two TSSE peaks during Event A (times t1 and t2 on the green solid line in Fig. 4), Fig. 11 shows the profiles of TSSE, storm surge, wave set-up, cross-shore flow velocity (V), and the wave height in the transect perpendicular to the coast. The 9 km long transect passes through the inlet and starts 1000 m inside the lagoon (i.e. the coastline is at x=1000 m), the end is offshore as shown in Fig. 9c (black line).

For time t1 the TSSE (top panel, black continuous line) was higher on the sea side than in the inlet channel. As a result, the sea-level slope at this time induced flood currents toward the inlet, which is in agreement with the negative (landward) cross-shore velocity (bottom panel, black solid line). The TSSE was lower than the storm surge height because the astronomical tide level was negative for both t1 and t2 (Fig. 4, green dotted line), where t1 was the time at which the TSSE reached its maximum value for Event A.

The maximum wave height for t1 (bottom panel, blue solid line) before the breaking point (where the waves breaks) was higher than 2.8 m. The peak wave period and the mean wave direction for time t1 were 6 s and 308^∘ (49^∘ with respect to the coast, from the north-west; not shown). The wave set-up reached 0.08 m inside the main channel and 0.07 m at the breaking point.

Table 2 Blocks of Progreso affected as a function of tidal phase during Event A.

For time t2, the TSSE inside the inlet channel was higher than offshore (12 cm difference; top panel, black dotted line), inducing an ebb flow (positive flow velocity V; bottom panel, black dotted line). Furthermore, the storm surge dropped significantly (more than 40 cm offshore; top panel, red dotted line) due mainly to weaker winds at the end of the storm. In terms of the wave climate, the maximum wave height for t2 decreased with respect to t1 (from 2.8 to 2.4 m), the peak wave period increased from 6 to 11 s – suggesting that longer swell waves reached the coast by the end of the storm – and the mean wave direction became more cross-shore (76^∘ compared to 49^∘ at t1). Under these conditions, the wave set-up reached 0.17 in the inlet channel and 0.14 m at the breaking point, which are significantly larger values than during the peak of the storm (t1). In fact, the wave set-up at the inlet channel represented 7 % of the wave height at the breaking point at t2 and only 2.8 % at t1 (Fig. 9).

The above analysis shows that the maximum wave set-up for Event A occurred at t2, when the mean wave direction reached its maximum northerly value and when the ebb reached its maximum. The wave set-up is controlled by the cross-shore radiation stress component and reaches its maximum when the incident wave direction is normal to the coast (Guza and Feddersen, 2012).

5 Discussion

5.1 Flooding at Progreso modulated by the tidal phase and the wind set-up inside the lagoon

The role of tidal modulation on the CACS events was studied by means of the numerical model. Since (i) Event A flooded larger areas than Event B, even if it occurred during low tide, and (ii) Event A could have occurred during other periods of the tide, a numerical experiment was carried out to assess the effect of a tidal phase on the flooding, wave set-up, and the residual sea level induced by CACS over Progreso (using the four tide scenarios shown in Fig. 4).

Figure 12 shows the maps of maximum flood (TSSE) corresponding to the four tide scenarios carried out for Event A (TS3 corresponds to the actual tidal phase that occurred during Event A). It is important to note that the time at which the maximum flood occurred is not the same on the sea side as inside the lagoon, nor is it the same for each simulation, similar to the criteria used with the inundation threat analysis in other studies (Zachry et al., 2015). The worst-case flood scenarios were for TS1 (high tide during the peak of the storm) and TS4 (rising tide near mean sea level during the peak of the storm), and TSSE were significantly lower with scenarios TS2 and TS3 (see summary in Table 2). This is partly due to the fact that during the storm – between 00:00 and 15:00 of 13 March 1993 (see Fig. 4, lower panel), i.e. during the period in which the local winds were stronger and able to produce large wind set-up inside the lagoon – the tide was high for TS1 and rising for TS4, while for TS2 and TS3 the tide was receding and near-low, respectively. Therefore, the wind and astronomical tide effects added up for TS1 and TS4. In turn, for TS2 and TS3, the tide was low and the residual sea level and wave set-up were the dominant factors (Fig. 13, top panel). Comparing TS1 and TS4, it shows that the highest TSSE on the sea side occurs with TS4, due to the eastward tidal currents during rising tide that contribute to the storm surge (not shown), while for TS1 the tidal currents were westward during receding tide and did not contribute to the storm surge. However, inside the lagoon, the maximum flood for TS1 occurred during stronger local winds (07:00), which generated a larger wind set-up, while for TS4 it occurred 4 h later, when the local wind intensity was receding, thus producing a smaller wind set-up. Therefore, when the wind direction is parallel to the main lagoon axis, a large wind set-up is generated at the eastern part of the lagoon. As for the results for TS2 compared to TS3, the flood levels occurred at the same time (07:00) but were higher for TS2, particularly inside the lagoon. This is mainly due to the fact that for TS2 the tide was receding and for TS3 the tide was near slack-low, which in turn produced a larger wave set-up inside the lagoon for TS2 (Fig. 13, top panel, right axis). This is in agreement with the findings of Dodet et al. (2013), who stated that during the ebb (receding tide), waves break over the ebb shoal, leading to stronger values of wave radiation stress than during the flood (see wave set-up for TS1 in Fig. 13), resulting in a larger wave set-up which propagates inside the lagoon.

Figure 12 Flooding maps with total sea surface elevation at Progreso under four tide scenarios for Event A: (a) maximum flooding for TS1, (b) maximum flooding for TS2, (c) maximum flooding for TS3, and (d) maximum flooding for TS4.

Figure 13 Residual sea level and wave set-up modulated by the tide. (a) Residual sea level and wave set-up for the four tide scenarios. (b) Wind speed at 10 m a.m.s.l. and H_s used as forcing agents for the numerical experiment with the the HW model.

In summary, the astronomical tidal phase during the passage of the storm is a very important factor for flooding, not only because of the tidal level itself, but also because of the interactions with the other contributors to the TSSE. For instance, Fig. 12 shows the city areas (blocks) affected for the different tide scenarios, showing that a storm with the characteristics of Event A would have been much more destructive if it had occurred during high or rising tide, as in TS1 and TS4. The astronomical tide also plays an important role in the total flooding, for instance in 2012 Hurricane Sandy inundated New York city at high tide, raising the water level to 3.5 m above mean sea level in the Battery Park area (southern tip of Manhattan), which exceeded the maximum water level during a hurricane in 1821 when water rose approximately 3.2 m at near low tide. If the 1821 event were to occur at high tide, a higher water level may be expected than the one observed during Hurricane Sandy (Woodruff et al., 2013). The flooded area increases when storm events occur at high tide (Rey et al., 2016; Zachry et al., 2015). The passage of CACS events, besides affecting water exchange with the sea and renewal dynamics (Viero and Defina, 2016) inside the Chelem lagoon, produce significant wind and wave set-up, characterized by non-linear interactions between meteorological forcings and the astronomical tide.

5.2 Residual sea level and wave set-up modulated by the tidal phase

Figure 13 shows the residual sea level (top panel, left axis) and wave set-up (top panel, right axis) for Event A under different tidal phases at the Chelem inlet as well as the wind speed offshore (low panel, left axis) and the significant wave height (H_s, low panel, right axis) used to force the numerical experiment. The residual sea level was obtained by means of a harmonic analysis using the T_Tide programme (Pawlowicz et al., 2002), from which the residual tide was determined by using the sea surface elevation for each scenario at the Chelem inlet from the hydrodynamic model using Case 1 boundary conditions described in Sect. 3.3. From these numerical experiments, it was found that the residual tide at the Chelem inlet is larger during low (TS3) or receding tide (TS2) but is smaller during high (TS1) or rising tide (TS4). The variation is non-linear and consistent with prior studies (Lin et al., 2012; Rego and Li, 2010), which attributed this behaviour to non-linear effects of the bottom friction and momentum advection on the surge due to the presence of the tide. This is also confirmed by Horsburgh and Wilson (2007), who, in a study on the coast of Great Britain, stated that surge peaks never occur during high water.

While the local wind contributes to the wind set-up inside the Chelem lagoon, the offshore wind (e.g. at location B located 161 km from Progreso; see Fig. 2) has a good correlation with the residual sea level (Fig. 13). For instance, the second peak in the offshore wind speed (location B) is also present as a second peak in all the residual sea-level scenarios as well as in H_s.

Model results from this study (Fig. 13) show that the wave set-up contribution inside the Chelem lagoon is tidally modulated as found in other studies for other sites (Smith et al., 2000; Smith and Smith, 2001; Olabarrieta et al., 2011; Dodet et al., 2013). Dodet et al. (2013) observed similar inlet hydrodynamic behaviour from data analysis and numerical modelling of the Albufeira lagoon (in Portugal). They stated that, during the ebb tide, currents cancel the intrinsic group velocity at the inlet, and waves are refracted by the ebb-jet current at the entrance of the inlet. Furthermore, this ebb-jet current caused the wave height at the inlet to increase, leading to more energetic wave breaking, and thus a greater wave set-up contribution inside coastal lagoons. Similarly, Gonzalez et al. (1985) observed from a case study at the Columbia River entrance that wave height increases during the ebb and decreases during flood tide. The opposing current retards the advance of a wave and can even block wave energy transport when the upstream component of the wave group velocity is equal to the current velocity, and the flood-induced current enhances the advance of the wave (Olabarrieta et al., 2011) but does not contribute significantly to the wave set-up.

Olabarrieta et al. (2011), who identified the effects of wave–current interaction on the circulation at Willapa Bay (Washington state), showed that the wave set-up inside estuaries increases with offshore energy, and Malhadas et al. (2009) suggested that wave set-up height inside the lagoon depends not only upon offshore significant wave height but also on tidal inlet morphology (mainly depth and length). These authors demonstrated by means of numerical solutions of simple idealized models that the deeper and shorter the morphology, the more the wave set-up is reduced. This is in fact the case in the present study. The Chelem inlet is 130 m wide, and the inlet channel is roughly 1.2 km long and 3 m deep (decreasing further along the channel), resulting in significant wave set-up inside the lagoon during strong storm events.

The above shows that the maximum residual tide and the maximum wave set-up did not occur at the same time at the inlet for any of the tidal phase scenarios used, but it does not seem to be determinant in the high flooding levels. Instead, what seems to play a more important role is the duration of these maximum values of residual tide and wave set-up. In fact, the highest values of TSSE at the inlet and those with longer duration are found in the TS4 tidal phase scenario (Fig. 4). This is mainly due to (a) a longer duration of significant residual tide and wave set-up, and (b) a rising astronomical tide, which translates to higher astronomical tidal levels and is associated with easterly tidal currents that contribute to the sea-level anomaly.

From this study, it was found that the most important contributions during CACS events, which pass over the peninsula, are the residual tide, the astronomical tidal phase, and the wind and wave set-up (inside Chelem).

5.3 Flood return periods

Regarding the non-linearity in the processes contributing to flooding and return periods, the extreme sea-level analysis showed that the astronomical tide has an important effect on determining the return periods for possible sea levels. In fact, the CACS events can occur at any tidal phase and thus the CACS residual tide is independent of any simultaneously occurring tidal phase. However, the interaction between residual and astronomical tides is an interesting subtle point to study: if this interaction is linear, the storm tide (TSSE without waves) probability would be roughly equal for each tidal phase and would be simply equal to the sum of the residual sea level and the astronomical tide, which means independence between the astronomical tide and the residual sea level. In that case, the use of joint probability methods could be used. These methods provide the chance to source variables by taking values at the same time and to create a scenario in which a flooding event may occur. This method is usually used for independent events (Chini and Stansby, 2012). For instance, Zhong et al. (2013) assumed independence between the astronomical tide and residual tide and estimated the storm tide probability of a joint probability method. However, hydrodynamic numerical experiments for Progreso in the present study suggested that the astronomical tide and residual sea level have a non-linear relationship, and thus applying the joint probability method for this case would not be adequate. A pragmatic and simplistic approach can consist of using a large data set of sea-level reanalysis (for instance 30-years in this case), assuming that a sufficiently large number of combinations of storm surges and astronomical tidal levels are present in the data set, and it can perform an extreme analysis of the sea level as a single variable. However, if quantifying the non-linear effect is the objective of a future study for Progreso, there are other options for estimating the storm tide probability given the height of the astronomical tide and phase when the storm surge arrives (assuming that the surge can happen at any time during a tidal cycle with equal likelihood), such as the one proposed by Lin et al. (2012). These authors developed an empirical function based on over 200 extreme tropical cyclone events (with both the storm tide and the storm surge simulated for the full range of tidal phases) for New York City. The implementation of this method for this zone is beyond the scope of this study, where the main objective is to study the hydrodynamic and flood-prone surface areas for Progreso during CACS events. However, this study shows the need to develop an empirical probability function based on the data for this area to estimate the probability of any storm surge for a given astronomical tidal phase as well as the wind direction and intensity with respect to the main lagoon axis.

6 Conclusions

This study has developed a 30-year sea-level hindcast using a hydrodynamic model forced by tides, mesoscale currents, as well as wind and pressure fields from hindcast data. Modelling results allowed extreme water levels to be identified and their probability of occurrence to be characterized using the GEV distribution function at Progreso port, Yucatán. Furthermore, the role of wave set-up was also investigated for two selected storm events, which correspond to the largest residual tide (Event A) and to the largest storm tide (Event B) identified during the simulated period. The analysis of the results shows the following:

• (a)

  The wave set-up and storm surge are large for events with strong wind intensity and long duration. The wave set-up at the inlet of the Chelem lagoon represented between 14.5 and 19 % of the TSSE, depending on the tidal phase (flood or ebb) in which the storm peak occurs. This contribution is higher during ebb than during flood tide, mainly because in the former the wave radiation stress terms are stronger due to current-induced wave breaking.

• (b)

  The local winds play an important role inside the Chelem lagoon, especially when the wind direction is parallel to the main axis of the basin, producing a large wind set-up. In terms of flooding, the most affected areas for both storm events (A and B) were along the eastern shores of the lagoon, due to significant wind and wave set-up, in particular for Event A. In this sense, the increasing water level inside Chelem lagoon during CACS events is not only due to the exchange of water through the inlet, but also because of wind and wave set-up over the lagoon, as well as non-linear interactions between these forcing agents and the astronomical tide.

Besides, the role of the tidal phase on the residual sea level, the wave set-up, and the total flooded area in Progreso were investigated based on numerical experiments varying the tidal phase for Event A. The results suggest the following:

• (c)

  The wave set-up and residual sea level are tidally modulated where the maximum wave set-up inside the Chelem lagoon occurs during receding tide (ebb) when the ocean water level is near the mean sea level and the incident wave direction is almost perpendicular to the coast. The residual sea level is larger during low or receding tide and smaller during high or rising tide. However, as expected, maximum flooding occurs when the CACS peak coincides both with rising tide (sea level near zero) and high tide (TS4 and TS1 scenarios). Nevertheless, the maximum values for the residual sea level and the wave set-up did not occur at the same time in any of the numerical experiments. The tidal phase difference with respect to the storm arrival determined the flood duration and the maximum water depth reached for each scenario.

• (d)

  If the largest CACS residual sea level (Event A) had occurred during high spring tide, the percentage of blocks flooded in the city would have increased from 25 to 60 %. This implies the need to accurately estimate the probabilities of residual and tidal levels in conjunction with local winds and wave set-up for a reliable estimation of coastal flood hazard caused by CACS events. This requires the definition of empirical probability functions specific for the area based on the astronomical tidal amplitude and phase, storm surge, and set-up due to both wind and waves.