Natural Hazards and Earth System Sciences Wave-induced extreme water levels in the Puerto Morelos fringing reef lagoon

Wave-induced extreme water levels in the Puerto Morelos fringing reef lagoon are investigated by means of a phase-resolving non-hydrostatic wave model (SWASH). This model solves the nonlinear shallow water equations including non-hydrostatic pressure. The one-dimensional version of the model is implemented in order to investigate wave transformation in fringing reefs. Firstly, the numerical model is validated with (i) laboratory experiments conducted on a physical model ( Demirbilek et al. , 2007) and (ii) field observations ( Coronado et al. , 2007). Numerical results show good agreement with both experimental and field data. The comparison against the physical model results, for energetic wave conditions, indicates that highand low-frequency wave transformation is well reproduced. Moreover, extreme water-level conditions measured during the passage of Hurricane Ivan in Puerto Morelos are also estimated by the numerical tool. Subsequently, the model is implemented at different along-reef locations in Puerto Morelos. Extreme water levels, wave-induced setup, and infragravity wave energy are estimated inside the reef lagoon for different storm wave conditions (Hs > 2 m). The numerical results revealed a strong correlation between the offshore sea-swell wave energy and the setup. In contrast, infragravity waves are shown to be the result of a more complex pattern which heavily relies on the reef geometry. Indeed, the southern end of the reef lagoon provides evidence of resonance excitation, suggesting that the reef barrier may act as either a natural flood protection morphological feature, or as an inundation hazard enhancer depending on the incident wave conditions.


Introduction
Fringing reefs are ubiquitous on tropical regions and are thought to provide natural shore protection during extreme wave events.The wave breaking at the shallow reef barrier induces radiation stress gradients, which are balanced by the water level increase inside the reef lagoon (i.e.setup).On the other hand, long-or infragravity (IG) waves, travelling with the short-wave group, are able to propagate into the reef lagoon with less dissipation.Resonant conditions develop if the incident group frequency is similar to the natural resonant frequency within the basin (Mei, 1983).Previous studies have shown that this physical mechanism can increase the coastal flooding hazard (e.g.Péquignet et al., 2009), and hence its contribution to the wave-induced water level variability needs to be considered.
The understanding of wave transformation in this environment has motivated several mathematical (Symonds et al., 1995) and experimental investigations (Seelig, 1983;Gourlay, 1994;Demirbilek et al., 2007).Prior studies (e.g.Massel and Gourlay, 2000) employed linear wave theory to estimate the wave setup inside the lagoon.Furthermore, spectral wave models, like SWAN (Booij et al., 1999), have demonstrated to be a valuable tool for understanding the role of waves on the reef-lagoon circulation (Lowe et al., 2009;Marino-Tapia et al., 2011).However, due to the limitations on the parameterizations of wave transformation processes (e.g.wave breaking, reflection and triad interactions), these models cannot fully resolve the wave transformation in this environment (Marino-Tapia et al., 2011).Therefore, the accurate modelling of IG waves transformation on fringing reefs requires the use of phase-resolving nonlinear wave models which solve triad interactions and wave breaking.
With the recent advances in computational technology, the use of numerical models based on the nonlinear shallow water equations (Kobayashi et al., 1989;Watson and Peregrine, 1992) and the Boussinesq approximation (Wei et al., 1999;Chen et al., 2003;Lynett, 2006;Cienfuegos et al., 2006) has become available for modelling nearshore wave transformation including wave breaking.Only recently, such numerical models have been applied to study wave transformation on fringing reefs (e.g.Nwogu and Demirbilek, 2010;Sheremet et al., 2011;Yao et al., 2012;Roeber and Cheung, 2012).Nwogu and Demirbilek (2010) validated a Boussinesq-type model using the data set of Demirbilek et al. (2007) and demonstrated the capability of the model for describing nonlinear processes.Sheremet et al. (2011) investigated the capability of nonlinear spectral wave models, both deterministic and stochastic, for the study of wave transformation on reef environments.However, the wave breaking model (i.e.Janssen and Battjes, 2007) requires calibration of three free parameters.Another approach for simulating breaking waves is the application of numerical models based on the Reynolds-averaged Navier-Stokes (RANS) equations (e.g.Lin and Liu, 1998), which have been validated for the study of breaking waves (Torres-Freyermuth et al., 2007;Pedrozo-Acuna et al., 2010) and low-frequency wave transformation (Torres-Freyermuth et al., 2010;Torres-Freyermuth and Hsu, 2010;Lara et al., 2011) in different environments.In a recent study by Lanza et al. (2012), the RANS model application over fringing reefs demonstrates its potential for describing reef hydrodynamics.However, the model application in large-scale domains (> 3 km) is still impractical due to computational constraints.
Limitations in classical nonlinear shallow water models have been overcome by incorporating wave dispersion in multi-layered nonlinear shallow water models.Zijlema et al. (2011) introduced SWASH (Simulating WAves till SHore), a numerical model based on the nonlinear shallow water equations including non-hydrostatic pressure.This model has been demonstrated to be relatively efficient in the study of rapidly varying flows including wave breaking, waveinduced setup, and nonlinear wave interactions (Zijlema et al., 2011).Therefore, the SWASH model should be suitable for simulating the wave transformation on fringing reefs.
In the present study, the one-dimensional (1-D) version of the SWASH model is employed to estimate wave-induced extreme water levels (i.e.setup and IG waves) in Puerto Morelos, Mexico.The numerical model allows us to evaluate (i) the role of the reef morphology in coastal protection; and (ii) a conservative estimate of extreme water levels inside the reef lagoon.It is well known that directional spreading decreases the energy transfer from high-frequency (hf-) to low-frequency (lf-) wave components owing to the detuning from resonance (e.g.Herbers and Burton, 1997).Thus, in the present study both wave-induced setup and lf-wave contribution would represent the worst case scenario.Although twodimensional processes (e.g.wave diffraction and refraction) and wind effects are relevant for modelling wave-induced circulation in these environments, they are out of the scope of this study and hence were not considered.This paper is organized as follows.In Sect.2, the reef morphology and storm wave conditions in Puerto Morelos are briefly described.Then, the model description and its validation with laboratory (Demirbilek et al., 2007) and field observations (Coronado et al., 2007) are presented in Sect.3. Section 4 presents the model application to study the dependence of wave setup, IG wave, and extreme water-levels, inside the reef lagoon, on the wave forcing/reef geometry.Finally, concluding remarks are given in Sect. 5.

Study area
The Puerto Morelos fringing reef lagoon is located at the north east of the Yucatan Peninsula, in Mexico, 30 km south of Cancun (Fig. 1).The coast of Puerto Morelos is protected by a 4 km barrier reef with variable water depth h r (5.3< h r <0.0 m), creating a reef lagoon of variable width W (550< W <1500 m).The lagoon is relatively shallow with an average depth of 3-4 m and maximum depth h l of 8 m at its southern inlet (Coronado et al., 2007).The lagoon is connected to the ocean by two gaps/inlets at the north and south, 6 and 8 m deep, respectively.Table 1 presents the reef-lagoon characteristics at the different along-reef locations depicted in Fig. 1.
The wave climate in Puerto Morelos is dominated by easterly swell coming from the Caribbean Sea with occasional influence from the northerly swell generated by winter storms in the Gulf of Mexico.In this area, waves greater than 2 m are considered as high energy storm conditions (Coronado et al., 2007;Marino-Tapia et al., 2011).The system is micro-tidal, with average semi-diurnal oscillations of 0.40 m.

Mathematical formulation
We use the open source model SWASH (http://swash.sourceforge.net),developed at Delft University of Technology (Zijlema et al., 2011).The model solves the twodimensional depth-averaged shallow water equations in non-conservative form.However, for simplicity the onedimensional governing equations, used in this study, are presented here: where t is time, ζ is the free surface measured from the still water level, d is the still water depth, h = ζ + d is the total depth, u is the depth-averaged flow velocity in x direction, q b is the non-hydrostatic pressure at the bottom, g is the gravitational acceleration, c f is the friction coefficient, ν t is the horizontal eddy viscosity due to wave breaking and subgrid turbulence, and w s and w b are the vertical velocity at the surface and the bottom.The wave breaking is accounted for in the model by considering the similarity between breaking waves and bores.To initiate the wave breaking process, steep bore-like wave fronts need to be tracked.When the steepness exceeds a fraction of the speed of the wave front such as the non-hydrostatic pressure is then neglected and remains so at the front of the breaker.The parameter α > 0 determines the onset of the breaking process.A value of α = 0.6 seems to work well for different wave conditions (The SWASH team, 2012).This is in close agreement to Longuet-Higgins and Fox (1977), who showed that α = 0.585 is the limiting slope for waves.The finite difference numerical model discretizes the vertical direction by means of a fixed number of layers, and the Courant number is the criterion to control and adjust dynamically the time step.For more details, interested readers should refer to the paper by Zijlema et al. (2011) where model equations, coefficient values, and boundary and initial conditions were first introduced.

Model validation
The SWASH model has been validated for different ranges of model applications (Zijlema et al., 2011), but has not for simulating wave transformation on a fringing reef.Therefore, model validation is here presented using results from laboratory experiments (Demirbilek et al., 2007) and field observations (Coronado et al., 2007).

Laboratory experiments (Demirbilek et al., 2007)
Laboratory experiments were carried out at the University of Michigan wind-wave flume (35 m long, 0.7 m wide, and 1.6 m high), equipped with a non-absorbing plunger-type wavemaker installed at one end and capable of generating irregular waves (H s < 10 cm and 0.4< f < 10 Hz).An idealized 1:64 model of a fringing reef (Fig. 2) was built inside, consisting of a 1:12 beach followed by a 4.8-m-long reef flat and a composite slope reef face similar to the one used by Seelig (1983).Free-surface elevation time series were measured using nine capacitance wave gauges located at different cross-shore locations (see Fig. 2).Different wave conditions, including irregular and bichromatic waves, for various water levels, were run by Demirbilek et al. (2007).
In this study, irregular waves, assuming a JONSWAP spectrum with γ = 0.3, for different wave conditions (H s and T p ) and the same water level h r were selected for the model validation (see Table 2).For validation, the numerical model was forced with the free-surface elevation time series measured at the offshore sensor WG1.The computational domain was set  to be 14-m-long with a regular grid of x = 0.01 m.The bottom friction coefficient c f was set equal to 0.002 and the wave breaking parameter α equal to 0.6 for all simulated cases.The peak frequency f p of the wave spectrum at the offshore location was employed in order to define the frequency threshold between high-frequency (f > 0.5f p ) and low-frequency (f < 0.5f p ) wave components (see Fig. 3).The high-(low-) significant wave height is calculated as four times the standard deviation of the high-(low-) passed water level fluctuation.The cross-shore variations of the hfand lf-significant wave heights (H s hf and H s lf , respectively) for Test 29 are shown in Fig. 4. The hf-wave components (f > 0.5f p Hz) experience strong energy dissipation due to wave breaking at the reef crest (Fig. 4a).On the other hand, the low-frequency wave (f < 0.5f p Hz) is partially transmitted into the reef lagoon and increases its energy with respect to the offshore conditions as it propagates shoreward.It is concluded that the numerical model predicts satisfactorily the measured value at the onshore location WG9 in Test 29 (Fig. 4a and Table 3).
With respect to wave-induced setup, the numerical model accurately predicts the relatively constant mean water level above the flat reef-top (e.g.Fig. 4b).Some discrepancies in the measured and predicted mean water level are observed at the breaking point (i.e.WG6).These differences have been also reported in other modelling studies (i.e.Nwogu and Demirbilek, 2010) employing the same data set.
Table 3 presents the summary of results for the simulated cases of measured and predicted wave setup η, H s hf , and H s lf at the onshore sensor (Table 2).In general, a good agreement between measured and predicted values of water level statistics is found, with a mean relative error of 7 % for all cases.Larger discrepancies (relative errors of 29 %, 13 %, and 13 %) between measured and predicted water level statistics (H s hf , η, and H s lf ) are found in the less energetic case (i.e.Test 26), whereas relative errors of only 6 % are found for other cases (e.g.Test 29).By conducting a sensitivity analysis, with different values of c f , we found no significant differences when the c f value is halved (i.e.c f = 0.001), whereas increasing the c f value by one order of magnitude (c f = 0.01) decreases both the hf-and lf-wave height by 18 % and 32 %, respectively, whereas the wave setup remained unchanged.Extreme water levels were measured during the nearby passage of Hurricane Ivan in September 2004.The significant wave height H s and peak wave period T p at LPM3 were 3.8 m and 14 s. Figure 5 shows the free-surface elevation time series recorded at LPM3 and the corresponding lowfrequency (f < 0.5f p ) wave component during one burst interval.Significant sea-swell dissipation occurred at LPM1 due to wave breaking at the barrier reef.The maximum wave setup (η = 0.36 m) was recorded at LPM1 when the lf-significant wave height inside the reef lagoon was important (H s lf > 0.6 m).Field observations suggest that, under extreme wave conditions, the tidal level modulates the wave setup inside the reef lagoon by determining the breaking type over the reef barrier (not shown), whereas the H s lf inside the lagoon is mainly controlled by the offshore wave energy at LPM3.The wave setup and IG wave energy were considered for estimating the long period water level fluctuations inside the lagoon and were obtained by averaging the free-surface elevation time series over a time interval equal to 2T p (e.g.Seelig, 1983) where T p (i.e.T p = 1 f p ) is the peak wave period at LPM3. Figure 6 (lower panel) shows the instantaneous water level and its long-term water level fluctuation at LPM1.Moreover, the water level probability distribution curve during Hurricane Ivan is also calculated (Fig. 6, upper panel), where the extreme water level η 2 % is here defined as the value above the still water level that is exceeded only two percent of the time.An across-reef transect extracted from the bathymetry, passing through the measuring locations (LPM3 and LPM1) and extending further to the shore, is employed in the numerical model.The computational domain is 3000 m long with a regular grid of x =1 m.The numerical model is forced with the free-surface time series measured at LPM3, and the still water level is referenced to the mean water level at the most offshore location LPM3 (i.e. 9 cm).The bottom friction in the numerical model is employed as the calibration parameter (c f =0.001), which is an order of magnitude smaller than the one reported at this site by Coronado et al. (2007).Thus, the calibration parameter seems to compensate other processes not accounted for in the modelling such as the wind setup, which can have an important contribution to water level increase during hurricane conditions.The significant wave height H s , lowfrequency significant wave height, wave setup, and extreme water level η 2 % are computed at the LPM1 location and compared with observations (see Table 4).The numerical model slightly overpredicts (underpredicts) the wave-induced setup (low-frequency wave energy) at this location, whereas the extreme water level is underpredicted by only 6 %.However, an overall good agreement between measured data at LPM1 and numerical results is observed.

Numerical implementation
Following the results of the previous section, the model is applied to different cross-reef transect extracted from the Puerto Morelos fringing reef lagoon bathymetry (Fig. 1).Despite the fact that two-dimensional effects are important, the one-dimensional mode is employed in order to simulate a broad range of wave conditions with affordable computational time.This allows us to investigate the IG growth mechanism and the coastal flooding hazard induced by the wave alone.For this task, cross-reef bathymetry profiles were extracted at different locations (see Fig. 7) and employed for the 1-D simulations.The profile locations were chosen in order to capture the strong spatial variability of the reef-lagoon geometry.The transects were oriented orthogonal to the shoreline for given locations (Fig. 1) that represent the southern area (P0), the main lagoon (P1), the northern reef-barrier gap (P2), and the northern end (P3).The reef profile characteristics, including water depth at the reef crest (h r ), maximum depth inside the lagoon (h l ), and reef-lagoon Table 5. Range of the reef geometry and the wave conditions parameters employed for the simulated cases.

Basin geometry
Wave conditions width (W ), at each location are summarized in Table 1.The fore-reef slope presents values close to 1 : 50.As shown in Fig. 7, the reef geometry presents a strong spatial variability; therefore differences in the wave transformation processes would be expected between locations.The numerical model is forced at 20 m water depth with (synthetic) irregular waves time series based on different wave parameters (H s and T p ) and assuming a JONSWAP spectrum with a peak enhancement factor of 3.3.The ranges of wave conditions employed for the simulated cases are shown in Table 5.The computational domain is 3.5 km long with a uniform grid system and cell size of x = 1 m.All simulations are 3600 s long and assume a constant tidal level during the computation period.The wave statistics are calculated from the last 1024 s of the computed signal of each simulation.The time step, t, is automatically adjusted during the computation in order to satisfy the stability constraints.The computational time for each simulation is of the order of few minutes on an Intel Xeon (quad-core) 2.53 GHz computer.The friction coefficient was set as a constant c f = 0.001 which is the value estimated during the model calibration in the study area (Sect.3.2.2).

Wave-induced setup
The wave setup is calculated at the intersection between the still water level (z = 0) and the cross-reef profile.Figure 8a shows the wave-induced setup at four along-reef locations (Fig. 7) for different wave conditions (Table 3).In general, the wave setup increases with increasing wave energy density (H 2 s T p ) at all locations.For low energy values (H 2 s T p < 1 × 10 2 m 2 s) the wave setup is relatively uniform along the reef-lagoon, with slightly lower values at P0.However, a strong dependence of wave-induced setup on the water depth at the reef crest h r is observed for the more energetic cases (H 2 s T p > 1 × 10 2 m 2 s), with maximum values occurring at locations with the higher reef crest P0 and P3), whereas smaller setup estimates are presented at the northern gap (P2, see Fig. 8a).This suggests that barrier reef only provides natural protection against wave-induced setup under less energetic wave conditions.On the other hand, for very large waves, the reef barrier induces a larger setup due to the increase in the stress gradients associated with the wave breaking process.Consistent with previous studies (Seelig, 1983), results show that the reef-lagoon width W does not play an important role on the wave-induced setup.

Infragravity wave energy
The other important contribution to the water level is the IG waves defined here as f < 0.5f p (Roelvink and Stive, 1989).The standard deviation σ of the low-passed signal is taken as an indicator of the amount of energy in this frequency band.Figure 8b presents the IG energy as a function of incident wave conditions.For relatively low wave energy conditions (H 2 s T p < 1 × 10 2 m 2 s), the IG energy is significantly higher at the northern gap (i.e.P2) with respect to the other locations, in particular P3.Under such conditions, the IG energy increases within the fringing reef lagoon due to nonlinear energy transfer from higher frequencies.Moreover, reef-induced wave reflection seems to be significantly smaller than at the other locations where the reef barrier is present.For the more energetic cases (H 2 s T p > 1×10 2 m 2 s), the IG energy is highly correlated with the lagoon width, with larger values occurring at the locations with the smaller W values (i.e.P0 and P1).At some instances, the IG energy at P0 is twice as large as at any other along-reef location.Notably, the reported oscillations in low-frequency energy at P0 (Fig. 8b) would suggest that resonance occurs at this location.The resonant modes can be estimated following Kowalik and Murty (1993) as  where n = 0, 1, 2, .. is the n-th mode of oscillation.The first three periods of oscillation were estimated for the different cross-reef profiles.The oscillation periods at the five locations are presented in Table 6.
According to Table 6, the T n values in P0 are within the range of the incoming (IG) wave periods and hence it is very likely that resonant conditions will occur at this location.It would appear that in Puerto Morelos the width of the reeflagoon has a greater influence on the resonant periods than the water depth (W = 460 m at P0; see

Extreme water levels
Finally, the extreme water level η 2 % at each location was calculated, following the methodology described in Sect.3.2.2, for the different wave conditions (Fig. 9).In general, the extreme water level is highly correlated with the IG energy contribution at most locations (Fig. 8a) with larger values at P2 and P0 for the less (H 2 s T p < 1×10 2 m 2 s) and more energetic wave conditions (H 2 s T p > 1 × 10 2 m 2 s), respectively.Therefore, the numerical results suggest that IG energy is very important for determining the hazardous areas behind this reeflagoon system.The cross-reef transect located at the northern end (i.e.P3) presents the lower extreme water levels for all range of wave conditions.

Conclusions
Wave transformation on fringing reefs is simulated using a one-dimensional nonlinear shallow water equations model.The numerical model is validated with both laboratory experiments and field observations.The model is in agreement with laboratory measurements of wave-induced setup and low-frequency significant wave height over the reef flat, whereas significant differences are observed for H s hf during mild wave conditions.The model is further applied to simulate wave transformation and wave-induced water levels inside the Puerto Morelos reef lagoon.Extreme water levels during Hurricane Ivan were satisfactorily predicted inside the reef lagoon.The numerical simulation of storm wave conditions suggests that wave-induced setup increases with decreasing depth at the reef-barrier for the more energetic wave conditions, while the IG energy increases with decreasing reef-lagoon width under such conditions.The northern end presents the lower values of IG energy and extreme water level for all range of wave conditions.It is very likely that IG wave resonance happens at the southern end of the reeflagoon for highly energetic waves.Thus, the reef geometry may increase or decrease the flood hazard in coastal areas at Puerto Morelos depending on incident wave conditions and reef geometry.Non-uniformity of extreme water levels suggests that two-dimensional effects must be incorporated in the modelling of this system and deserve further investigation.

Fig. 1 .
Fig. 1.Location of the study area and reef-lagoon bathymetry.The numbered circles represent the locations of the cross-reef profiles extracted for the simulations.The crosses represent wave measurement locations by Coronado et al. (2007).

Fig. 6 .
Fig.6.Water level probability distribution at LPM1 for incident H s = 3.82 m and T p = 14 s.The probability distribution curve is calculated from the moving average (dark, solid line) of the instantaneous free-surface elevation (light, thin line).
Fig. 8. (a) Wave setup and (b) infragravity-wave energy as a function of wave conditions for all simulated cases at the intersection between the still water level and the cross-reef profile for the locations shown in Fig. 1.

Fig. 9 .
Fig. 9. Extreme water level as a function of incoming wave energy at the intersection between the still water level and the cross-reef profile for the locations shown in Fig. 1.

Table 1 .
Cross-reef profile characteristics at different locations.

Table 3 .
(a) Measured and (b) predicted wave-induced setup (η), hf-(H s hf ) and lf-(H s lf ) significant wave height at the onshore sensor WG9.Units in cm.

Table 4 .
Model-data comparison of significant wave height (H s ), low-frequency significant wave height (H s lf ), mean water level (η), and extreme water level (η 2 % ) at LPM1.Wave conditions at LPM3 are H s = 3.8 m and T p = 14 s.

Table 6 .
The first three oscillation periods for each transect.

Table 1 )
. This is in agreement with results presented by Pomeroy (2011). www.