Run-up parameterization and beach vulnerability assessment on a barrier island : a downscaling approach

We present a downscaling approach for the study of wave-induced extreme water levels at a location on a barrier island in Yucatán (Mexico). Wave information from a 30-year wave hindcast is validated with in situ measurements at 8 m water depth. The maximum dissimilarity algorithm is employed for the selection of 600 representative cases, encompassing different combinations of wave characteristics and tidal level. The selected cases are propagated from 8 m water depth to the shore using the coupling of a third-generation wave model and a phase-resolving nonhydrostatic nonlinear shallow-water equation model. Extreme wave run-up, R2 %, is estimated for the simulated cases and can be further employed to reconstruct the 30year time series using an interpolation algorithm. Downscaling results show run-up saturation during more energetic wave conditions and modulation owing to tides. The latter suggests that the R2 % can be parameterized using a hyperbolic-like formulation with dependency on both wave height and tidal level. The new parametric formulation is in agreement with the downscaling results (r= 0.78), allowing a fast calculation of wave-induced extreme water levels at this location. Finally, an assessment of beach vulnerability to wave-induced extreme water levels is conducted at the study area by employing the two approaches (reconstruction/parameterization) and a storm impact scale. The 30-year extreme water level hindcast allows the calculation of beach vulnerability as a function of return periods. It is shown that the downscaling-derived parameterization provides reasonable results as compared with the numerical approach. This methodology can be extended to other locations and can be further improved by incorporating the storm surge contributions to the extreme water level.


Introduction
The assessment of beach vulnerability in low-lying areas is important for coastal managers and decision makers.Furthermore, these coastal systems are particularly sensitive to climate change effects, such as mean sea level increase and storm intensification (Wong et al., 2014).Thus, it is anticipated that low-lying areas will experience more severe coastal flooding and beach erosion during the following decades.
The beach vulnerability can be estimated by comparing extreme water level elevations to those of the beach morphology features (Sallenger, 2000).For instance, the storm impact scale proposed by Sallenger (2000) couples the fluid forcing and the beach morphology by examining the relationship between the dune height and the water level due to the storm surge, wave setup, and extreme run-up.This approach was validated in Stockdon et al. (2007) for a stretch of coast in North Carolina.Stockdon et al. (2007) employed lidar-derived measures of pre-storm dune and berm elevation and hurricane-induced water levels to hindcast the potential storm impact regime to the landfalls of Hurricane Bonnie (1998) andHurricane Floyd (1999), which were further compared to the observed response.More recently, long-term observations were employed together with a run-up parameterization in order to determine the return periods correlated to the storm impact scale on the coast of Oregon (Serafin and Ruggiero, 2014).
A great effort has been devoted to the development of methodologies for storm surge estimation (e.g., Lin et al., 2010;Irish et al., 2011).However, less attention has been given to the development of reliable approaches for the estimation of wave-induced run-up.Wave run-up is often calculated using parameterizations based on field observations (e.g., Ruessink et al., 1998;Ruggiero et al., 2001Ruggiero et al., , 2004;;Stockdon et al., 2006;Senechal et al., 2011).However, their performance is questionable during extreme wave conditions (Stockdon et al., 2014).Furthermore, run-up parameterizations are strongly dependent on the beach morphology features, tidal level elevation, and wave forcing conditions.Therefore, a universal run-up parameterization is not available and site-specific parameterizations should be developed.
The advent of nonlinear phase-resolving wave transformation numerical models allows the simulation of wave run-up in a wave-by-wave basis.Different approaches have been developed with different degree of sophistication, including nonlinear shallow-water equation (NLSWE) models (e.g., Kobayashi and Wurjanto, 1992;Raubenheimer and Guza, 1996;Zijlema et al., 2011), Boussinesq-type models (e.g., Wei et al., 1999;Chen et al., 2003), Reynolds-averaged Navier-Stokes models (e.g., Lin and Liu, 1998;Losada et al., 2008), and large eddy simulation (Christensen, 2006;Zhou et al., 2014) models.The capabilities of more sophisticated approaches for addressing the study of small-scale processes demand higher computational cost.Non-hydrostatic NL-SWE models (e.g., SWASH (Simulating WAves till SHore)) allow us to overcome some of the limitations in classic NL-SWE models by incorporating wave dispersion in the simulations.This numerical approach has been employed for the study of extreme water levels on a fringing reef lagoon (e.g., Torres-Freyermuth et al., 2012), wave run-up on beaches (e.g., Ruju et al., 2014;Guimarao et al., 2015), and infragravity shoreline dissipation (e.g., de Bakker et al., 2014).Furthermore, the potential for the run-up parameterization has been shown in Brinkkemper et al. (2013).
The characterization of small-scale coastal processes such as run-up from long-term data sets requires the implementation of a statistical method.Downscaling of wave conditions for the study of nearshore processes is possible through data reduction using a rigorous statistic approach (Camus et al., 2011a;Guanche et al., 2013).Camus et al. (2011a) employed a hybrid downscaling methodology to transfer wave climate to coastal areas.They show that interpolating simulated results using a radial basis function provides good estimates to characterize a complete year of hourly sea states and that the decrease of the error is negligible considering a subset of cases.This methodology was further extended to reconstruct time series of stability parameters on vertical breakwaters by Guanche et al. (2013).Therefore, combining numerical mod-els with statistical methods provides a means to characterize coastal dynamics and reduce the computational cost.
The aim of this work is to present a methodology for the assessment of wave-induced vulnerability at a location on a barrier island located in Yucatán (Mexico).The methodology combines numerical, statistical, and probabilistic methods for the estimation of wave-induced water levels associated with return periods.Moreover, this approach allows the derivation of a site-specific run-up parameterization for the study area.The outline of this paper is the following.Firstly, the study area location and characteristics are described in Sect. 2. The methodology for downscaling wave conditions in order to obtain a 30-year run-up hindcast is presented in Sect.3. Section 4 presents the reconstruction of the extreme water level time series, the derivation of a new run-up parameterization, and an assessment of beach vulnerability in the study area.An uncertainty analysis regarding the use of wave hindcast information is presented in Sect. 5. Finally, concluding remarks and future work are presented (Sect.6).

Study area
Dzilam de Bravo is located on a barrier island in the northern Yucatán Peninsula (Fig. 1).The coastline is fronted by a 200 km wide continental shelf with a very mild (1 : 1000) beach slope (Enriquez et al., 2010).The tidal regime is micro-tidal and wave conditions in the study area are dominated by local sea breezes and mesoscale meteorological events (cold fronts) known as Nortes (Appendini et al., 2013).Furthermore, less frequent hurricane events can also affect the study area.According to Mendoza et al. (2013), the Yucatán coast is more vulnerable to flood than to erosion during the impact of storms.Dzilam de Bravo is characterized by submarine dune fields (Cuevas et al., 2013) that induce a complex nearshore wave transformation.Cuevas et al. (2013) characterized the submarine dunes by means of sub-bottom seismological profiles, finding at Dzilam de Bravo a mean dune height ranging from 0.8 to 1.0 m and a mean dune wave length of 98-120 m, predominantly moving northwestward.These sedimentary deposits (see Fig. 2) might play an important role in wave energy dissipation, providing a natural protection from storms to this site.
The beach profile at Dzilam de Bravo was measured using the Differential Global Positioning System (DGPS) and extended landward with terrestrial lidar information acquired in 2011.Moreover, the beach profile was further extended offshore to a water depth of 10 m assuming an equilibrium profile according to Dean (1991).Wave information at 8 m water depth is available from a 30-year hindcast  for the Gulf of Mexico and the Western Caribbean Sea (Appendini et al., 2014).These data were estimated by means of a third-generation spectral model forced with wind data from the North American Regional Reanalysis (NARR) (Mesinger et al., 2006)   calibrated/validated in deep waters with wave buoys (Appendini et al., 2013) and altimeter information (Appendini et al., 2014).Data measured with an acoustic Doppler current profiler (ADCP) located near the study area, between Chuburná and Yucalpetén (see Fig. 1) at approximately 8.5 m water depth, were available for a 2.5-year period (June 2010-December 2012).Comparison of the in situ data and wave hindcast information (NODE12972, Fig. 1) presents good correlation between the model and observations for H s 1 m (see Fig. 3).The model underestimates H s for values between 1.2 and 1.7, whereas hindcast data overestimate observations for H s 1.7 m.The differences can be ascribed to the relative coarse resolution of the wind data employed to generate the 30-year wave hindcast by Appendini et al. (2013Appendini et al. ( , 2014) ) for resolving wave generation by local winds and the lack of high-resolution bathymetry available for this area in the ETOPO1 (Amante and Eakins, 2009).However, an overall good agreement is observed between model and data.Thus, in this study a wave hindcast node located at ap-proximately 10 m water depth in front of Dzilam de Bravo (NODE11583, Fig. 1) was selected as the offshore boundary condition (H s , T p , θ ).

Methods
We extended a methodology to downscale wave information to the nearshore as proposed by Camus et al. (2011a) and Guanche et al. (2013) for the assessment of storm impact on barred beaches.The methodology is as follows.Firstly, a subset of wave conditions is selected from the 3-hourly 30-year wave hindcast.Then, the selected sea states, with the corresponding tidal level, are propagated from deep waters to the shore by employing numerical models.Extreme run-up is computed and further employed for calculating a 30-year run-up hindcast by means of interpolation.Subsequently, the 30-year run-up information was employed in order to derive a run-up parameterization for the study area.Finally, the storm impact for different return periods can be obtained using both numerical results and the new parameterization.This methodology assumes that the dissimilarity in offshore wave conditions leads to comparable dissimilarity in the run-up data.This assumption is supported by the strong correlation between run-up and offshore wave conditions reported in previous research (e.g., Stockdon et al., 2007).It would be expected that only small differences in dissimilarity could arise when run-up height saturates under extreme conditions (e.g., Senechal et al., 2011).

Selection of wave conditions
The available 30-year wave hindcast (Appendini et al., 2014) consists of a total of 87 664 sea states (H s , T p , and θ ), one every 3 h.Due to the computational effort involved to downscale the complete data set for all sea states, it is desirable to obtain a representative subset.A comparison of selection algorithms applied for the analysis of wave climate is presented in Camus et al. (2011b).They found that the subset of wave conditions obtained by implementing the maximum dissimilarity algorithm (MDA) was representative of the variety of sea states and therefore appropriate for downscaling wave climate.
The aim of the MDA, described in detail in Camus et al. (2011a, b) and Guanche et al. (2013) for coastal engineering applications, is to identify the most dissimilar subset of multivariate vectors (i.e., wave parameters) in a database.Therefore, the extracted subset of M vectors represents the diversity of the data set consisting of N n-dimensional vectors.In this study, the multivariate data include significant wave height, H s , peak period, T p , mean wave direction, θ m , and mean sea level, Z m .Wave parameters were obtained from the wave hindcast, while the time series of sea level corresponds to the astronomical tide prediction for this area (www.predmar.cicese.mx)during the same time period.
Following the procedure described in Camus et al. (2011a) and Guanche et al. (2013), the multivariate data at deepwater are defined as where N corresponds to the 87 664 sea states from the 30year wave hindcast.The first step in the methodology described in Camus et al. (2011a) is to normalize the vector components so they can be evenly weighted in the similarity criterion, defined by the Euclidean distance.Special consideration should be made for the circular variable (direction) when it is adapted to a linear scale, since it is recorded in a continuous scale where 0 and 360 • are identical.Therefore, the circular distance should be implemented in that case where the maximum distance in the circle is equal to π .The sample data consisting of N dimensionless vectors are defined as from which a set of M vectors D 1 . . .D M is selected by means of the MDA.
The selection starts by transferring one vector from the data sample to the subset D.Then, the rest of the M − 1 vectors are selected calculating the dissimilarity between each of the remaining elements in the database and the elements in the subset, transferring the most dissimilar one to the subset, considering the MaxMin version of the algorithm as proposed by Camus et al. (2011a).This procedure is repeated iteratively until the M elements are selected.
For instance, having a subset of R(R ≤ M), the dissimilarity among vector i of the data sample N − R and the j vectors of the R subset is determined as Then the dissimilarity between vector i and subset R, d i,subset , is obtained as (4) Now, having calculated the N − R dissimilarities, the following data to be selected have the maximum d i,subset .In this work, the Euclidean distance was computed using the Dis-tanceMatrix algorithm developed by Fasshauer (2007), modified for the case of the directional parameters considering the circular distance as described in Camus et al. (2011a) and given by the following expression: (5) The final step is to denormalize the selected subset using Here, a total of M = 600 sea states were selected that adequately represent the whole sample, and their distribution uniformly covers the area of the input data as well as its borders (Fig. 4).It is also worth considering that the selected sea states are well distributed along the time series of wave parameters and sea level (Fig. 7).

Propagation of selected wave conditions
Wave propagation from 8 m water depth to the shoreline was performed employing the coupling of a spectral wave model (Simulating WAves Nearshore, SWAN; Booij et al., 1999) and a phase-resolving nonlinear non-hydrostatic model (SWASH; Zijlema et al., 2011).The SWAN model is a third-generation wave model for coastal regions, based on a Eulerian formulation of the discrete spectral balance of action density, which accounts for wind generation, whitecapping, triad and quadruplet wave-wave interactions, bottom friction, and wave-induced wave breaking (Booij et al., 1999).The SWASH model employs the nonlinear shallow water equations, including terms for non-hydrostatic pressure, which makes the model suitable for simulating wave transformation due to nonlinear wave-wave interactions in both surf and swash zones, wave-current interaction, wave breaking, and wave run-up.Wave breaking is included in the model based upon the bore formation concept.Flooding and drying of grid cells is important for a correct run-up simulation.In this model, no special features are needed to model dry cells accurately when the time step is chosen correctly, as flooding never happens faster than one grid size per time step (Zijlema et al., 2011).
The SWAN model is run in stationary one-dimensional mode (mesh size 1 m) along the section from 8 to 4 m water depth and forced with a JONSWAP distance was computed using the spectrum at the offshore boundary.The wave energy spectrum at 4 m depth calculated by the SWAN model is employed as the seaward boundary forcing for the SWASH model (Fig. 5).The SWASH domain extends from 4 m water depth to the shoreline with a mesh size of 0.1 m.The initial time step is 0.025 s with a maximum Courant number of 0.5.Simulations were sampled for 2 170 s, after 530 s of spin-up time.

Model data analysis: extreme water level calculation
The instantaneous water level elevation, η(t), relative to mean sea level was extracted from the SWASH simulations for each sea state propagated as the height of the bottom profile at the location of the wet-dry interface with respect to time (Fig. 6a).This location was tracked as the first grid point in which the water depth was less than 0.005 m in order to obtain a continuous time series.Subsequently, the extreme run-up was calculated from the run-up maxima (R) following the work by Stockdon et al. (2006) as the 2 % exceedance value (Fig. 6b).Additionally, the mean value of the wave run-up time series (< η >), which corresponds to a super-elevation of the mean water level due to the presence of waves known as the wave setup (Longuet-Higgins and Stewart, 1964), was obtained for each case.Following Sallenger (2000) and Stockdon et al. (2007) we define the extreme water levels R high = R 2 % + Z and R low = < η > + Z for each simulated case.

Reconstruction of the extreme water level (R high ) time series: radial basis function (RBF) interpolation
The extreme water levels, R high and R low , associated with each of the 600 selected sea states were employed to reconstruct the 30-year-long time series.Notice that the storm surge is not included in Z for this work but can be incorporated.The extreme water level time series reconstruction is performed by means of an interpolation technique based on the RBFs.This is an exact interpolation technique, given that the interpolated surface always passes exactly through the data points, and is suitable for multivari- ate scattered data interpolation.Franke (1982) tested the performance of about 30 methods for scattered data interpolation, finding that the best and second best were methods based on RBFs.This method has been previously implemented in diverse applications such as the reconstruction of topographic surfaces based on coordinate data (Hardy, 1971) and, more recently, for the downscaling of wave parameters (Camus et al., 2011a;Guanche et al., 2013).Following the method presented in Camus et al. (2011a) and considering that X i = {H si , T pi , θ mi , Z i }; i = 1, . . ., N represents each sea state in the 30-year-long time series and . ., M represents each one of the M = 600 selected cases associated with a value of R high , and the interpolation function is given by where is the radial basis function, ||.|| denotes the Euclidean norm, p( and is a Gaussian function defined as where the points D j , j = 1, . . ., M are the centers of the RBF approximation and c is a shape parameter that must be carefully selected since it has a strong influence on the accuracy of the solution (Rippa, 1999).The interpolation based on RBF was performed by means of an algorithm developed by Fasshauer (2007) which incorporates the algorithm proposed by Rippa (1999) for the selection of an optimal value for the shape parameter c.The selection is performed by minimizing the root mean square error (RMSE) of a data fit based on a radial interpolant for which one of the centers was left out ("leave-one-out cross validation" approach).
The coefficient b of the monomials and the coefficient a of the RBF are obtained by the interpolation conditions (Camus et al., 2011a): where D p,j are the real functions defined by the calculated extreme water level values R high which correspond to the M selected sea states (D j ).Subsequently, the R high time series can be reconstructed for the 30-year period by means of the RBF as follows (see bottom panel of Fig. 7): R high,i = RBF R high D j , R high,j (j = 1, . . . M) , X i , (10) where i = 1, . . ., N.
Similarly, the wave-induced mean water level time series, R low = < η > + Z, is reconstructed following the same methodology with M = 600.Camus et al. (2011a) compared the reconstructed time series for a range of M values (i.e., M = 25, 100, and 1000) against simulated time series of N = 8784, finding that the error obtained in the estimation of wave parameters is almost negligible considering only M = 100 cases.Therefore, they recommend that for the specific application of transformation of wave climate from deep to shallow waters, 100 ≤ M ≤ 200 is an adequate number of cases.Guanche et al. (2013) validated their interpolated values with those calculated analytically, finding that the reconstructed series with more than 100-200 cases out of 500 000 cases reached values of less than 1 % error, and with 500 cases the error made is almost negligible.
In this work, a sensitive analysis on the dependency of R high to the number of cases employed for the reconstruction was conducted for M = 50, 100, 200, 300, 400, 500, 600.The mean and standard deviation of each time series were computed, finding that for M > 300 cases the variability of these statistic parameters is insignificant (not shown).

Run-up parameterization
The reconstructed 30-year extreme water level time series provides a mean to correlate R 2 % to offshore wave conditions.The run-up is obtained by subtracting the astronomical tide Z from R high for each case.Different run-up parameterizations from the literature (Ruggiero et al., 2001;Senechal et al., 2011) were employed and calibrated using the downscaled data.Furthermore, due to the observed modulation of run-up by tides (e.g., Guedes et al., 2011), a new parameterization of run-up and setup was derived for the study area as a function of the tidal level and wave conditions employing the 30-year R 2 % data.

Extreme value analysis of R high and R low
In order to incorporate the probability to a given extreme water level, the annual maxima of R high and R low were fitted to the generalized extreme value (GEV) distribution of Jenkinson (1955), which has been widely employed for modeling extremes of natural phenomena.The GEV distribution is given by where µ and σ are the location and scale parameters, respectively, and the shape parameter k determines which extremevalue distribution is represented: Gumbel (k = 0), Fréchet (k > 0), and Weibull (k < 0).WAFO-group (2000) toolbox was used for the GEV model, in which the parameter estimation methods used are the maximum likelihood method (Prescott and Walden, 1980) and the probability-weighted moments method (Hosking et al., 1985).The latter was selected for parameter estimation since it is more suitable for  (Sallenger, 2000;Stockdon et al., 2007).
small samples (N = 15, 25) (Hosking et al., 1985) and resulted in a better goodness of fit to the annual maxima data.The estimated parameters for R high yearly maxima are the shape parameter k = 0.3057 with a 95 % confidence interval, and location and scale parameters estimated to µ = 1.5739 and σ = 0.1238.Regarding the R low annual maxima data, the estimated parameters are k = 0.3184, µ = 0.7247, and σ = 0.0896.

Storm impact scale
A storm impact scale for barrier islands that considers the magnitude of fluid forcing (storm-induced water levels) relative to beach morphology (dune/berm elevation) was first proposed by Sallenger (2000).The model defines four storm impact regimes (Table 1) depending on the relative relationship between the sand dune/beach berm elevation and the storm-induced water levels.The elevation measures are defined as R low (the sum of storm surge, astronomical tide, and wave setup), R high (the sum of storm surge, astronomical tide, and R 2 % ), D high (dune crest), and D low (dune toe).
The results of the extreme value analysis of R high and R low for different return periods were correlated with the beach morphology features (i.e., D high and D low ) derived from the lidar data.

Results
The 30-year time series of R high (Fig. 7) and R low (not shown) reconstructed by means of the RBF interpolation of the simulated cases are further employed for the parameterization of run-up and setup and an assessment of the vulnerability of the beach by means of the storm impact scale proposed by Sallenger (2000).

Run-up parameterization
Beach vulnerability is often evaluated using run-up parameterizations (e.g., Stockdon et al., 2007;Serafin and Ruggiero, 2014).Therefore, the development of suitable parameterizations is important for vulnerability studies.The extreme run-up results, obtained from the reconstructed extreme water levels, are further analyzed in terms of their relationship with offshore wave parameters, beach conditions, and astronomical tide.For that purpose, the 5 % exceedance value of water level according to the astronomical tide Z was found for high (Z ≥ Z 5 % = 0.32 m) and low water level (Z ≤ Z 5 % = −0.32m), while for the mean water level the values considered were 0.05 ≥ Z ≥ −0.05 m.Analyzing the behavior of the R 2 % and < η > values it was clear that they are modulated/saturated by the tides/wave energy conditions.Firstly, a linear relationship for R 2 % is employed, as proposed by Ruggiero et al. (2001), which depends on deepwater wave parameters (H s , L 0 = g T 2 /2π ) and beach slope, S. Deepwater wave parameters correspond to the wave hindcast data, while the value of S = 0.09 was considered according to Brinkkemper et al. (2013).This parameterization describes remarkably well the behavior of the values corresponding to high water level (Z ≥ Z 5 % = 0.32 m, darker gray dots in Fig. 8a) and is very similar to the best fit (zero intercept) to the data associated with high water (solid line).Furthermore, a relationship obtained in a previous study (Brinkkemper et al., 2013) performed on the same area but employing the results of only five simulations, corresponding to energetic wave conditions associated with high water level (Fig. 8a, dash-dot line), is almost identical to the one obtained on this study for all values corresponding to high water level (Fig. 8a, solid line).Even though the r 2 value of the linear fit (zero intercept) is satisfactory (r 2 = 0.87) for the case of high water level, these linear relationships are only valid up to a value (1.1 R 2 % 1.3 m) where saturation in maxi- mum run-up values is observed for more energetic conditions ((SH s L 0 ) 1/2 4 m).The linear fit performed to the whole data set (dotted line) showed a lower r 2 value and a higher RMSE (see Table 2).Regarding the setup values, the linear parameterization proposed by Stockdon et al. (2006) is employed, which includes the deepwater wavelength, L 0 (T 0 ), and the foreshore slope, β f .This formulation lies on the upper limit of the data (Fig. 8b, dashed line), better describing the values associated with high water even though it differs from the linear fit (zero intercept) to those values.The r 2 corresponding to the linear fit associated with high water values is greater than that obtained by considering the whole data set.Similarly to run-up values, saturation and a dependency on water level are observed for setup values.The saturation value for setup associated with high water level is around 0.35 m.
Due to the observed saturation of run-up and setup values associated with the water level (astronomical tide), a more suitable hyperbolic-like relationship is employed.This saturation was previously examined by Brinkkemper et al. (2013) following the relationship proposed by Senechal et al. (2011).It was found that for H s 1.5 m the saturation value of R 2 % varies accordingly with water level (Fig. 9).Therefore, a hyperbolic tangent relationship was fitted through the method of least squares separating the data in high (Z ≥ Z 5 % = 0.32 m), mean (0.05 ≥ Z ≥ −0.05 m), and low (Z ≤ Z 5 % = −0.32m) water levels.For each set of data a hyperbolic tangent fit was performed (Fig. 9), obtaining acceptable values of r 2 and RMSE (Table 3).The relationship proposed by Senechal et al. (2011) falls on the upper limit of the data corresponding to high water level.The tanh argument of the relationship presented in Senechal et al. (2011) and that obtained in this study for high water level are almost the same (0.39H 0 and 0.4H 0 ) and very similar to that obtained by Brinkkemper et al. (2013) of 0.5H 0 .However, the saturation value in the case of Senechal et al. ( 2011) is much higher (2.14 m) than that obtained by Brinkkemper et al. (2013) of 1.62 m, which is exactly the same as the one obtained in the present study.
Thus, in order to obtain a generalized expression for the prediction of the 2 % exceedance value of run-up (R 2 % ) and setup (< η >) as a function of H 0 and Z, a simple linear relationship was fitted to the hyperbolic tangent fit parameters, a and b (Table 3), with r 2 values of 0.99 and 0.95, respectively, in the case of the run-up fit and 0.97 and 0.63 for setup fit parameters.The generalized expression obtained for R 2 % is given by where H 0 is the deepwater wave height, a = 1.615Z + 1.098, and b = −0.297Z+ 0.476 m.A similar expression was obtained for the case of the setup, with a = 0.23Z + 0.27 and b = 0.15Z + 0.46.The run-up and setup were calculated with the linear and the generalized expressions, the latter depending only on the deepwater wave height and the astronomical tide.The values obtained through the linear parameterization showed a greater dispersion for more energetic conditions with respect to the values obtained using the hyperbolic parameterization, due to the saturation of wave run-up not accounted for on the linear relationship, showing a correlation of 0.73 with respect to the reconstructed run-up values (Fig. 10a).Regarding the hyperbolic parameterization, the parameterized values were compared to the reconstructed values obtaining, for the case of the R 2 % , an r 2 value of 0.78 considering the whole set and an r 2 = 0.86 considering only the waves approaching from the north (NNW, N, NNE): 22.5 > θ > 337.5 (Fig. 10b).Regarding the setup, < η >, the correlation obtained for the whole data set is 0.75 and is 0.86 for the data associated with waves arriving from the north (not shown).For both cases, R 2 % and < η >, the dispersion is greater for smaller values, R 2 % 0.5 m and < η > 0.1 m, and for waves arriving with an angle 22.5 < θ < 337.5.
The overall modulation of run-up and setup (not shown) due to wave height and astronomical tide is captured by the generalized hyperbolic parameterization and better illustrated in the reconstructed vs. parameterized run-up time series (Fig. 10c).Some deviations are observed but the general behavior is reproduced satisfactorily as compared to the linear parameterization where overestimation of run-up maxima is observed and the modulation of the run-up due to tides is not well captured (Fig. 10c, dashed line).

Extreme water levels
The reconstruction of the 30-year extreme water levels allows us to determine the corresponding return periods.In order to do that, a GEV distribution is fitted to the 30-year (reconstructed) R high time series.Figure 11 shows the return level extrapolation of R high for the 100-year return period with the 95 % confidence bounds (100 Similarly, the extreme analysis for R low can be performed.These parameters can be employed in order with associate the storm impact scale to a given return period (e.g., Serafin and Ruggiero, 2014).Additionally, the same procedure was followed for the parameterized values of R 2 % obtained by the generalized expression described in the previous section, incorporating the contribution of Z to obtain the R high and R low 30-year time series.The values of R high associated with a 5-, 10-, 50-, and 100year return periods obtained from the reconstructed time series are smaller than the ones predicted from the parameterized time series for all return periods (Table 4).However, the R low values predicted from the reconstructed time series are greater than the parameterized R low values (Table 5) for all return periods (5, 10, 50, and 100 years).

Beach vulnerability assessment on a barrier island
Using the 30-year-long time series of R high (missing surge) and R low (missing surge), obtained by means of both the RBF interpolation and the parameterizations, together with the topographic elevations (D high and D low ) obtained from the lidar, the storm impact regimes can be estimated.
Based on the return values of R high and R low (Tables 4  and 5) corresponding to a 5-, 10-, 50-, and 100-year return period, the associated storm impact regime "collision" (D high > R high > D low ) was found (Table 6).This regime would occur even if the storm surge were not considered to cause long-lasting erosion and the possibility of sediment not returning from offshore.However, considering a typical storm surge elevation of ≈ 0.5 m associated with the frequent cold fronts in the study area in addition to the mean return values of R high would result in water elevations that would exceed the dune crest for a return period of 10 years or more (Table 6,Fig. 11), leading to the "overwash" R high > D high storm impact regime, where run-up overtops the dune and the sand transported landward does not return seaward to the beach under post-storm conditions.The storm impact scale is a valuable tool for predicting coastal response to storms (with an accuracy that depends on the regime) as well as for analyzing the longshore variability of coastal change in a stretch of coast (Stockdon et al., 2007).However, in this work, the longshore variability is not evaluated since a single beach profile is considered.

Discussion
In order to evaluate the impact associated with driving the model with wave hindcast information, we employed the same methodology followed in this work but using the ADCP wave data.We selected 60 conditions from the available 3year measurement period in order to conduct beach run-up simulations.The same analysis was conducted using wave hindcast information during the same period.
Figure 12 shows a comparison between reconstructed time series of R 2 % derived from measured and hindcasted wave data for the same time period used in Fig. 3. Wave run-up obtained from hindcasted wave data is poorly reproduced for less energetic conditions but satisfactorily describes the upper envelope of wave run-up values with respect to measured wave data (Fig. 12a).Run-up estimates seem reliable for storm waves approaching from the north (NNW, N, NNE), despite differences in offshore wave height (see Fig. 3).These differences can be ascribed to the fact that the run-up calculations based on hindcast information are compensated by the slight overprediction of the peak wave pe- riod.Due to this compensation by the wave period, extreme run-up heights retrieved from model runs with the hindcasted data as input are in good agreement with those retrieved from the model runs with measurements as input (Fig. 12b).For mean wave conditions, however, the run-up is consistently overpredicted.A summary of the extreme run-up statistics is shown in Table 7 for all wave conditions and storm conditions only.The correlation coefficient for the whole period is very poor (r 2 = 0.43 and RMSE = 0.23) due to a limited wave hindcast resolution for resolving local processes (i.e., sea/landbreezes).However, the correlation increases significantly (r 2 = 0.80 and RMSE = 0.16) when constraining the analysis to storm waves associated with Nortes.For storm conditions, relative errors of the run-up statistics between hindcast and measured data are smaller than 20 % with a relative error of only 4 % for the maximum R 2 % .The latter suggests that the methodology employed in the present work is valid, since we are focused on the study of extreme events.However, the use of high-resolution wind fields for driving wave generation models is necessary for the study of wave run-up under mean conditions.

Conclusions
Extreme water levels on a barrier island, located on the northern Yucatán Peninsula, are investigated using a downscaling approach based on wave hindcast information.Wave run-up on the study area presents a dependency on offshore wave height and tidal elevation.A new parameterization which incorporates saturation and tidal modulation is derived from the downscaling information.Both downscaling results and the run-up parameterization provided similar results in terms of return periods and the storm impact at this location.The uncertainty analysis of the impact of employing wave hind-cast information suggests that it does not significantly affect the extreme water level analysis.Future work will be devoted to conducting the model calibration using run-up measurements and to the inclusion of the storm surge contributions in the extreme water levels.These two aspects need to be addressed in order to achieve a more reliable analysis of beach vulnerability in this area.

Figure 1 .Figure 2 .
Figure 1.Location map indicating the study site (Dzilam de Bravo) at the barrier island backed by the wetlands and the mainland of Yucatán Peninsula, the position of the hindcast nodes used for the simulations, and the ADCP location in between Chuburna and Yucalpeten ports.
Figure 3. (a) Time series section of observed and hindcast significant wave height and (b) Q-Qplot showing the comparison between ADCP wave data and modeled hindcast data at a location close to the study area.Solid line in (b) indicates perfect correlation.

Figure 4 .
Figure 4. Significant wave height (H s ) against wave period (T p ) and mean direction corresponding to the 30-year wave hindcast data, and astronomical tide (Z), showing the distribution of the selected cases using the MDA algorithm for M = 600.

Figure 5 .
Figure 5. Beach profile at Dzilam showing the section corresponding to simulations performed by the SWAN model (10-4 m depth) and the section for SWASH simulations (4 m depth to shore).Dotted line represents the mean sea level.

Figure 6 .
Figure6.Section of water level elevation time series relative to mean sea level (η(t)) extracted from the wet-dry boundary of SWASH simulations for every sea state propagated to shore, indicating tide Z, run-up maxima, R, and setup at the shoreline < η > (left panel).The 2 % exceedance value was extracted from the cumulative distribution function (cdf) of the R values (right panel).

Figure 7 .
Figure 7. Time series of the wave hindcast data (H s , T p , θ ), sea-level (Z), and interpolated run-up value (R high = R 2 % + Z), indicating the selected/simulated cases.

Figure 8 .
Figure 8. Linear fit to (a) R 2 % values associated with high water level (Z ≥ Z 5 % = 0.32 m, darker gray dots) and comparison to Ruggiero et al. (2001) and Brinkkemper et al. (2013) parameterizations; (b) setup values (all and high water levels) compared to Stockdon et al. (2006) parameterization.The dotted line corresponds to the linear fit performed to the entire set of data.

Figure 10 .
Figure 10.Reconstructed R 2 % values obtained from the RBF interpolation correlated against R 2 % values obtained from (a) linear vs. hyperbolic tangent fits; (b) hyperbolic tangent fit to all R 2 % values vs. those corresponding to waves approaching from the north (337.5 ≤ θ ≤ 22.5); (c) comparison of a time series section of reconstructed vs. parameterized R 2 % values.

Figure 11 .
Figure 11.Return value extrapolation of (a) R high and (b) R low in the GEV model with the 95 % confidence bounds.

Table 1 .
Storm impact scale regimes

Table 2 .
Zero intercept linear regression to R 2 % values with respect to (SH s L 0 ) 1/2 and < η > values with respect to β f (H 0 L 0 ) 1/2 , considering the whole set of data and those associated with high water level (HWL) (Z ≥ Z 5 % = 0.32 m), correlation (r 2 ), and root mean square error (RMSE) in meters.

Table 6 .
Storm impact regimes corresponding to the 5-, 10-, 50-, and 100-year return values of R high and R low , considering a value of D high and D low of 2.27 and 0.8 m, respectively, and a typical value of storm surge (≈ 0.5 m).

Table 7 .
Extreme run-up statistics for measured and hindcasted wave conditions during a 3-year time period(2011)(2012)(2013).The error analysis and statistics correspond to all wave conditions and only storm conditions.