Human Impact on the Spatiotemporal Evolution of Beach Resilience on the Northwestern Yucatan Coast

The northern Yucatan peninsula is prone to coastal flooding and erosion owing to its low-land elevation and high exposure to storms. Therefore, it is important to assess the capability of the beach-dune system to resist, recover, and adapt from storms in the context of coastal development and climate change. This work aims to investigate the role of human impacts on the spatiotemporal evolution of the Coastal Resilience Index (CRI) in the area. The study is conducted on a prograding micro-tidal beach located on the vicinity of coastal structures. Beach and dune morphometrics, characteristic beach parameters, and maximum shoreline recession were determined from the analysis of beach profiles undertaken along a 2-km straight of coastline during the 2015–2020 period. Moreover, the maximum extreme water levels were estimated using in situ data and numerical models. This information is employed to assess the alongshore and temporal variability of the beach resilience. The results suggest that the beach and dune morphology present alongshore and temporal variability due to the human impacts associated to the dune degradation and the presence of coastal structures. The analysis shows that coastal resilience has been increasing over the past years but presents significant alongshore variations. High CRI values were found at beach transects presenting low anthropogenic impact, whereas low resilience was observed at transects with a degraded dune or located in the vicinity of coastal structures regardless of presenting high progradation rates. The observed beach response during the passage of recent tropical storms is consistent with the CRI values.


INTRODUCTION
Low-laying coastal areas are prone to natural and anthropogenic perturbations associated to population growth, extreme events, and climate change (Wong et al., 2014). Sea level rise due to climate change is expected to reach between 0.5 and 1.0 m by the end of the century (Nicholls et al., 2014). The population growth in coastal areas has been accelerated over the past decades owing to increasing economic activities and is expected to continue accelerating during the current century (Neumann et al., 2015).
Coastal ecosystems (beaches, dunes, mangrove, seagrass, reefs, and wetlands) provide numerous ecosystem services including habitat to species and carbon dioxide sequestration (Pendleton et al., 2012). Furthermore, they can provide natural coastal protection by inducing wave energy dissipation during extreme events (Narayan et al., 2016), and the foredune acts as a natural barrier which provides ultimate coastal protection against coastal flooding (Sallenger, 2000). However, coastal development has threatened such ecosystems, often triggering a greater exposure to extreme events and sea level rise (Neumann et al., 2015), reducing their resilience.
Beach morphology plays a key role in the coastal protectionacting as a positive feedback-by adjusting to the prevailing hydrodynamic conditions to dissipate more effectively the wave energy (Hoefel and Elgar, 2003). Stationary changes in the pattern of cross-shore transport, associated with wave conditions, are responsible for beach accretion-erosion cycles (Aubrey, 1979). Nevertheless, the beach capacity to resist storms, without losing its functionality, can be heavily undermined by the presence of coastal structures (groins and breakwaters) which interrupt the littoral transport. This can become a major problem at locations with a large net alongshore sediment transport. Furthermore, dune degradation also has important implications on resisting flooding and erosion during storm events (Durán and Moore, 2013;Feagin et al., 2019).
The coastal resilience has been employed to provide a framework to propose sustainable solutions (Masselink and Lazarus, 2019). The National Academy of Sciences defines the resilience as "the ability of a system to prepare, absorb, recover, and more successfully adapt to an adverse event" (National Research Council, 2012). The resilience concept has been addressed from disciplinary and inter-disciplinary perspectives. Differences between the stable equilibrium (i.e., bounce back) and the dynamic equilibrium (i.e., change of state) approaches for the study of resilience have been pointed out by Holling (1996). For the study of beach morphology changes due to storms, the dynamic equilibrium approach is more suitable and hence is adopted in the present work.
Previous studies have proposed methodologies to assess the beach-dune resistance during extreme events. For instance, the storms impact scale proposed by Sallenger (2000) correlates the dune parameters with the extreme water levels to identify the coastal resistance. This approach has been successfully employed to conduct studies at different costal locations around the world (Stockdon et al., 2007;Medellín et al., 2016;Franklin et al., 2018; among many others). Sallenger's storm impact assumes a morphological response based on the storm impact scale (i.e., swash, collision, overwash, or inundation). However, previous works have pointed out some inconsistencies with field observations. For instance, Cohn et al. (2018) found that the dune can increase during the collision regime, whereas Sallenger (2000) predicts erosion. More recently, Dong et al. (2018) introduced the Coastal Resilience Index (CRI), based on the relationship between beach/dune morphometrics and storm-induced extreme water levels and beach erosion, which incorporates the role of the beach morphology (beach width and beach elevation) and beach erosion (maximum shoreline recession) on coastal protection. The CRI is considered a more integrated approach, and thus was employed in the present work given the availability of beach-dune morphology data with a high spatial and temporal resolution.
The aim of this work is to investigate the role of human impact on costal resilience on a microtidal prograding beach located in Sisal, on the northwestern Yucatan peninsula's coast. This approach allows us to improve our understanding on the elements that determine a resilient beach in the study area to further identify strategies to improve its resilience. The outline of this paper is the following. First, material and methods are described, including the description of the study area, data sources, and the framework employed in this work. Section "Results" estimates the key elements required to determine the trajectory of the CRI in the study area. Finally, discussion and concluding remarks are given.

Study Site
The northern Yucatan Peninsula coast is located in the Gulf of Mexico (Figure 1A), and is characterized by a wide and shallow continental shelf, karstic geology, low-land elevation, and the presence of a barrier island ( Figure 1B). Micro-tidal diurnal regime (0.4 m spring amplitude) and low-energy wave conditions (Hs = 1 m) are ubiquitous in this region (Torres-Freyermuth et al., 2017). In the Yucatan Peninsula, tourism activities have triggered a threefold population growth between 1990 and 2010, more considering the floating population. The northern Yucatan peninsula has been experiencing chronic beach erosion at some locations, ascribed to both the degradation of the dune that protects the coast during storm events and the impact of coastal structures on the littoral transport (Appendini et al., 2012). This problem has been exacerbated by sea level rise, which was of 2.5 mm/year during the second half of the past century (Zavala-Hidalgo et al., 2010).
In addition, winter storms, known as Central American Cold-Surge (CACS) or Norte events, associated to the cold-front passages occurring in the Gulf of Mexico between October and April (Medina-Gómez and Herrera-Silveira, 2009), are frequently present in the northern Yucatan Peninsula Coast, with a yearly average of 22 events . Wave conditions in the study area during Norte events are characterized by NNW relatively energetic (Hs > 2 m) swell waves (T > 6 s). Storm surges associated to such events are typically 0.5 m (Rey et al., 2018). On the other hand, less frequent local storms causing wind squalls (Turbonadas) are more likely to occur during the springsummer months. Turbonadas drive strong winds and energetic sea waves (Hs = 2 m and Tp = 3 s), and the associated storm surge during these events can reach up to 1.0 m. Hurricane landfalls in the northern Yucatan peninsula are not very frequent (only two major hurricanes in 32 years), owing to the east-west shoreline orientation facing the Gulf of Mexico to the north ( Figure 1A). Besides, mean wave conditions, associated to sea-breeze events (Torres-Freyermuth et al., 2017), are responsible of a westward net sediment transport on the order of 30,000-60,000 m 3 /year (Appendini et al., 2012). These conditions enhance the erosion downdrift of coastal structures along that stretch coast (Medellín et al., 2018).
This study is conducted in Sisal beach, located on a barrier island in the vicinity of the Sisal's Port and pier ( Figure 1C). The urbanized area has grown mainly toward the wetlands south of the barrier, where land reclamation by population is a common practice in the region. Furthermore, the beach and foredune have been affected due to the presence of the structures and touristic activities that have conducted to the degradation of the dunes ( Figure 1C).

Field Data
Beach profiles were measured along 20 equally spaced cross-shore transects, encompassing a 2-km straight of coast ( Figure 1C). Beach surveys were undertaken on average every 15 days, with the highest and lowest temporal resolution of 7 and 30 days, respectively. The beach transects cover from the subaerial beach to a water depth of approximately 1.5 m. A Differential Global Positioning System (DGPS) was employed in Real Time Kinematics (RTK) mode, carrying the rover receiver on a backpack, for undertaking the 116 beach surveys analyzed for this work. The first and last surveys correspond to 08/05/2015 and 05/03/2020, respectively. Control points are measured at the beginning and at the end of each survey to measure the antenna elevation. Furthermore, a bathymetric profile along transect P05 ( Figure 1C) was measured 10 km offshore till a 10 m water depth on August 2019 with the aim of propagating the measured wave conditions to the shore. The measured elevation data are referenced to the MEX97 geoid 1 . Further information on the monitoring program and the interannual and seasonal beach variability at this site can be found in Medellín and Torres-Freyermuth (2019).
Water levels and offshore wave conditions were measured in the study area; the former by means of an ultrasonic tide gauge installed inside the port of Sisal (Figure 1C 2 ) and the latter by means of a bottom-mounted Acoustic Doppler Current Profiler deployed 10 km off the coast at 10 m water depth ( Figure 1B). Data gaps in the records were filled with reanalysis information from Wave Watch III 3 . Hindcast and measured data at the ADCP location have an R 2 value of 0.79 (Medellín and Torres-Freyermuth, 2019).

Model Simulations
Extreme water levels were estimated by means of numerical models. A maximum dissimilitude algorithm (e.g., Camus et al., 2011a,b;Guanche et al., 2013;Medellín et al., 2016) was employed to select 300 representative cases based on the combination of the measured mean sea level and offshore wave conditions from the 4-year dataset. The bathymetry information and the beach profile, representative of natural conditions (P05), were employed for the model simulations. The selected wave conditions were propagated from 10 m water depth to the coast by means of the coupling of wave transformation (SWAN; Booij et al., 1999) and hydrodynamic (SWASH; Zijlema et al., 2011) models. The measured or reanalysis wave data provide SWAN with offshore (h = 10 m) wave conditions that were propagated to 4 m water depth, followed by SWASH model simulations of the nearshore wave transformation including wave breaking and wave run-up. The run-up time series were further employed to estimate the high water level (R HIGH = WR+MS+AT), which includes the elevation due to the runup of 2% of exceedance (WR) for given values of storm surge (MS) and astronomical tide (AT) elevations for the simulated cases. This information was employed to reconstruct the time series of R HIGH for the 4-year period.

Assessment of Coastal Resilience
For the study of resilience, it is important to specify both the system and disturbances of interest (i.e., resilience of what to what) (Carpenter et al., 2001). This allows to identify the perturbations (response), the system resistance, threshold values, scales, and trajectories (Piégay et al., 2020). We adopt the CRI, proposed by Dong et al. (2018), focused in the beach and dune system and its response to storms. A description of each of the required components for its assessment is provided below.

Morphology Indicators
Relevant morphologic parameters are used to identify the characteristic parameters and assess the spatiotemporal evolution of the beach-dune system. The protective width PW, protective elevation PE, and the dune crest elevation DE (Figure 2A) were identified following Dong et al. (2018). Both, natural and human perturbations are implicitly incorporated in the evolution of the morphometric parameters. Each parameter was calculated for each transect (e.g., Figure 1C) and for all surveys.
The CRI requires knowledge on characteristic beach morphology parameters PW 0 , PE 0 , and CF 0 , where the dune crest freeboard CF 0 = DE 0 − (MS + AT) depends on the MS and AT. Here, we estimated the characteristic parameters from the spatial (20 beach transects) and temporal (116 surveys) average of the individual morphology parameters (Supplementary  Figures 1, 2).

Disturbances and Perturbations
Coastal disturbances that induce beach/dune morphology changes can be either natural or anthropogenic, and instantaneous (shocks) or long-lasting (stresses). The response of the system to such disturbances (e.g., shoreline and dune changes) is known as perturbations. In the study site, human disturbances include the presence of coastal structures (jetty and pier) and the degradation of the dune. For instance, the port's jetty ( Figure 1D) is responsible of shoreline perturbations associated to large progradation rates (7 m year −1 ) updrift of the jetty. Also, a pier (Figure 1F), located 1.8 km east from the port, induces beach morphology response on adjacent transects (P01, P02, and P03 in Figure 1C). Dune degradation, associated to touristic activities, is also ubiquitous at some transects (i.e., P01, P02, and P09 in Figure 1C). All these perturbations are implicitly incorporated in the morphometric parameters, with increasing influence near the structures.
Natural disturbances in the study site are often related with Norte events, local storms, and less frequent tropical storms. To identify the maximum storm-induced shoreline recession (MSR), beach width and beach volume differences between subsequent surveys were computed for the study period (Supplementary Figure 3A). The MSR corresponds to the maximum negative difference between campaigns from all beach transects. On the other hand, tide gauge information and numerical model results were employed to determine the extreme water levels. The maximum mean sea level (MS+AT) max is obtained from the tidal record. Similarly, the maximum high water level (R HIGH ) max can be extracted from the time series reconstructed from the numerical model simulations (Supplementary Figure 3B).

Coastal Resilience Index
Larger beach width, subaerial volume, and dune height contribute to increase the resilience of beaches to extreme events Frontiers in Marine Science | www.frontiersin.org FIGURE 2 | (A) Beach parameters and water levels required to estimate the Coastal Resilience Index according to Dong et al. (2018), spatiotemporal evolution of (B) protective width, (C) dune elevation, (D) subaerial beach volume, and (E) subaerial beach elevation between May 2015 and March 2020 along a 2-km straight of coast. Elevation in C and E is referenced to the geoid MEX97. (Dong et al., 2018). The coastal resilience factors proposed by Dong et al. (2018) are adopted in this work and are defined as, where a is the protective elevation factor, b is the volume density, c is the protective width factor, d is the crest freeboard factor, and e is the wave runup factor. The latter (e) has been modified to make it more consistent with Sallenger's model and to allow us to investigate the CRI associated to a given return period. The percentage of fine sediment s is set to 0.2. These factors are computed using the morphology parameters and perturbations obtained as described in previous sections. Based on the CRI values, Dong et al. (2018) defined the resilience as low (CRI < 1.5), medium (1.5 < CRI < 2.0), and high (CRI > 2).

Morphology Indicators
High-resolution beach surveys allow us to capture the spatiotemporal beach evolution during the 5-year study period along the 2-km stretch of coast. The beach protective width presents the highest increasing trend at transects located near the jetty (P15-P20) ( Figure 2B). A net increase in the protective width PW is observed for transects P6-P20, whereas more stable conditions are shown at P1 and P3-P5. The beach width evolution at P2, located east from the pier, is similar to the one observed at P20, but with a reduced magnitude. The dune elevation DE shows the larger values away from the structures ( Figure 2C). The maximum DE growth was observed at P05 (Figure 2C and Supplementary Figure 1C). Significant DE growth is also observed for P03, P07, P08, P10, and P14.
The spatiotemporal evolution of the beach volume ( Figure 2D) presents a similar behavior than the beach width ( Figure 2B). These two parameters determine the protective elevation PE (Figure 2E) which is strongly dependent on the DE as noticed by the similar spatiotemporal evolution.

Characteristic Beach-Dune Parameters
The temporal average (denoted by . ) and standard deviation of each parameter are computed for each transect. The PW and its variability increase toward the port's jetty where largest high progradation rates occur (Supplementary Figure 2A). On the other hand, central transects, located between the port and the pier, present less variability implying larger stability (e.g., P07) to increase again near the pier. The characteristic protective width PW 0 can be estimated by the spatial and temporal average of all transects, finding a value of PW 0 = 38 m which is close to the one observed at central transects.
The PE presents the lower values near the jetty (Supplementary Figure 2B) as a result of a very wide beach with low elevation (no dune). More uniform conditions with similar variability occurred between P01 and P17, except for P09 located on a heavily impacted beach area that has inhibited the dune recovery (Supplementary Figure 1C). The characteristic protective elevation PE 0 is estimated to be 0.8 m, which is similar to the mean values of most of the transects (Supplementary Figure 2B). The DE also shows the lower values near the jetty and the highest values associated to low progradation rates (Supplementary Figure 2C). A characteristic dune elevation of DE 0 = 2.2 m is employed to determine the characteristic dune crest freeboard CF 0 =0.9 m.

Storm Perturbations
The maximum storm-induced shoreline recession during the study period corresponds to a Norte storm event that arrived to the study area on 5 February 2016 (Figure 1E and Supplementary Figure 3A). The significant wave height at 10 m water depth reached 2.5 m (Figure 3A). The event was followed by two milder Norte events with maximum H s = 1.5 m. Beach profiles conducted before (03/02/2016) and after (11/02/2016) the storm sequence show that the maximum shoreline recession MSR = 17 m ( Figure 3B) occurred at the beach transect located near the Port's jetty. The associated subaerial beach volume loss at this transect, after the storm sequence, was 15 m 3 /m.
Numerical model results and tide gauge data were used to identify the maximum mean water (MS+AT) and high water levels (R HIGH ). Both, tide gage measurements and modeling results, show maximum values for the 4th May 2016 storm event ( Figure 1G and Supplementary Figure 3B). A local storm, known as Turbonada, was generated by strong winds driving energetic waves ( Figure 3C) and a large storm surge ( Figure 3D). Calm conditions (Hs = 0.5 m) preceded the storm arrival with H s = 1.8 m. The inundation was strongly related with the mean sea level increase (MS+AT = 1.3 m) during the peak of the storm. Numerical model simulations predicted a (R HIGH ) max =2.3 m during this event which was responsible of dune overwash at some locations ( Figure 1G). The maximum values of the measured/simulated water levels and beach erosion (MS+AT, R HIGH , and MSR) during the period of the study are employed to assess coastal resilience in the study site.

Coastal Resilience Trajectories
Storm perturbations, characteristic beach-dune parameters, and morphometrics evolution allow us to estimate trajectories of the CRI. Supplementary Figure 4 shows the spatiotemporal evolution of the different factors that contribute to coastal resilience. The protective elevation factor a, crest freeboard factor d, and wave runup factor e show lower values in the vicinity of coastal structures and at the degraded dune transects (e.g., P09). On the other hand, the volume density and protective width factors present higher values east (updrift) from the jetty and groin. The crest freeboard and wave run-up factors are strongly dependent on the dune elevation and control the CRI.
The CRI, obtained by the summation of the five factors, shows significant spatiotemporal differences (Figure 4). In general, high CRI values are located at the mid-transects away from the structures (Figure 4A). The differences in the temporal evolution can be seen more clearly at selected alongshore locations ( Figure 4B). The time series of the CRI from all selected transects started with negative values of CRI due to the absence of a dune at the beginning of the monitoring program. The resilience increases as the foredune and beach grow at each transect partially due to beach progradation. Transect P05, presenting a high dune growth (Supplementary Figure 1C), shows a sustained CRI increase despite the low progradation rate (Supplementary Figure 1A). Near the jetty (P20) and the pier (P01-02), the CRI presents low values despite the net increase in the PW (Supplementary Figure 1A). At the transect with the most degraded dune (P09), surrounded by transects with high resilience, the CRI remains low during the study period ( Figure 4B).

DISCUSSION AND CONCLUSION
Field observations and numerical model simulations were employed to investigate coastal resilience on a prograding beach located on the northwestern Yucatan coast. The results show an overall increase of the CRI over the past 5 years. The analysis suggests significant spatial gradients of beachdune morphometrics along the 2-km stretch of coast, and hence different CRI trajectories were observed. The low DE growth observed at some locations can be attributed to human impacts such as the influence of the Port's jetty (P20) and dune degradation (P09). Thus, low resilience occurs at alongshore locations impacted by human activities, whereas high resilience is observed at the most preserved cross-shore transects.
To validate our results, we employed beach surveys collected after the study period, corresponding to transects P05 and P20 that were undertaken before and after tropical storms Gamma and Delta (2-12 October 2020) ( Figure 3E). The low resilient transect P20, located near the jetty, receded more than 20 m due to the sequence of storms ( Figure 3F) and to the lowest beach elevation and presence of the jetty (Figures 2C,E). These observations are consistent with the low CRI (Figure 4) at this location (CRI < 1). On the other hand, transect P05, with a high CRI (CRI > 3), shows small changes during the same period. Hence, the CRI seems to provide a good proxy to predict the beach-dune response during storms. Coastal flooding on the barrier island during these events occurred on the wetlands border due to the heavy rain and saturated water table, emphasizing the need to investigate the resilience across the barrier island.
The present study highlights the important role of beachdune conservation for beach resilience. The beach evolution toward a more resilient state was observed at beach transects with low human impact. The alongshore location with the largest protective width PW, associated to sediment impoundment, was significantly less resilient than those transects with high dune elevation DE, low human impact, and low progradation rates. This demonstrates the important role of dune conservation for building resilient coasts on this region.
Future work should be devoted to improving the CRI. For instance, additional aspects such as the presence of submerged coastal ecosystems (seagrass and coral reefs) or submerged morphology features (sand bars) play an important role on wave energy dissipation and must be incorporated. Another drawback in the current CRI is that, like in other coastal vulnerability models (e.g., Gornitz, 1990;Sallenger, 2000; among many others), it does not incorporate the socio-economic aspects and the strong interaction between the large confined aquifer and the coastal ocean, which are fundamental for a more integral study of coastal resilience. Thus, future efforts should also be devoted to integrate ecological and socio-economic components to quantify the coastal resilience. Finally, the present work shows the importance of high spatial and temporal resolution data to assess the CRI trajectories at a local scale. However, extending this study to a regional scale requires the use of remote sensing information, for the derivation of digital elevation models (Almeida et al., 2019) and the assessment of water levels and wave conditions, in combination with numerical models.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
AT-F conceived the study, assisted in some of the field data collection, and wrote the first draft. GM conducted all the topographic surveys and data analysis, conducted the numerical model simulations, and assisted in discussion and writing. PS assisted with the discussion, analysis of the results, and writing. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
Technical field support was provided by José López González and Juan Alberto Gómez Liera (ADCP data collection and bathymetric profile), and IT support by Gonzalo Uriel Martín Ruiz. Tidal data were provided by the National Tidal Service (Servicio Mareográfico Nacional). We acknowledge NOAA for making WWIII wave hindcast data available. Two reviewers provided valuable comments that helped to substantially improve the manuscript. Finally, we thank Delft University of Technology for making the development of SWAN and SWASH models possible.