Modeling of Coastal Erosion in Exposed and Groin-Protected Steep Beaches

Praia de Faro is a sandy beach located in the barrier island system of Ria Formosa, on the southern Portuguese coast [Figs. 1(a–c)]. The barrier splits the Atlantic Ocean on the front side from a coastal lagoon on the backside [Fig. 1(d)]. The study site is an elongated peninsula orientated 130° approximately (measure from the north). This beach is subject to significant urban development pressure, and the oceanfront is partially stabilized with rocks/walls or naturally protected by a dune. The dune morphology varies alongshore with higher robustness and elevation [7–8 m above mean sea level (MSL)] at the western portion of the study area [F6 in Fig. 1(d)], while at the central and eastern part, the dune is lower (6–7 m above MSL) and weaker. The central part, F4 and F5 in Fig. 1(d), is periodically overtopped during spring tides or storm conditions (Ferreira et al. 2019; Matias et al. 2010). The site presents a subtidal terrace and steep beach face, with an average slope of around 0.10 (Almeida et al. 2012; Vousdoukas et al. 2012a). In previous studies, such as Ferreira et al. (2016), this site has been classified as reflective during calm conditions and intermediate during energetic conditions with the formation of a longshore bar. A beach berm can be normally found, except after very energetic storms, with variable widths ranging from less than 15 m to more than 40 m (Ferreira et al. 2016). The study area responds rapidly to storm events and variations of the wave forces, and the beach can regain a large part of the sediment after storm events in days/weeks. Moreover, the site is characterized by the presence of multiple highly dynamic beach cusps at the lower-beach face that interact with the more persistent upper-beach face cusps (Vousdoukas et al. 2012a). Sediments are medium to very coarse, moderately well sorted sands with a median (d₅₀) of around 500 µm and d₉₀ around of 2,000 µm (Vousdoukas et al. 2012b).

The town of Quarteira is located 10 km northwest from Praia de Faro [Fig. 1(c)]. The analyzed sector consists of a set of three sandy beaches with a total longshore length of 900 m. The main orientation of the coastline is 120°N, and the average beach slope is 0.10. Sediment grain size is slightly finer than in Praia de Faro, and d₅₀ and d₉₀ are 485 and 900 µm, respectively. The three beaches are laterally limited by 150 m-long rocky groins [Fig. 1(e)]. These groins make the three beaches behave like “manmade pocket beaches” subject to beach rotation as a function of wave direction. During the dominant conditions, the groins maintain the sediment in the system; however, during very energetic conditions, the sediment can fall outside of the system. Beach nourishments have been performed in the area to guarantee a reasonable beach width for bathing conditions, with a total of 360,000 m³ placed in 1998 (Pinto et al. 2018). At the backside, the beach is limited by a long promenade with an elevation ranging from 6 to 8 m. While Quarteira represents a relevant touristic destination in Portugal, the beach morphodynamic has been poorly investigated.

Two storms impacting the area have been considered for analysis. Storm Emma, in March 2018, was a severe storm that traveled toward the northeast (Ferreira et al. 2019). The maximum significant wave height, Hs, during this storm was 6.55 m, with a peak period and direction of 13 s and 240°N, respectively [Fig. 2(a)]. Moreover, the timing of the storm matched with a spring tide, exacerbating the impact of the storm. The wave characteristics were measured at the Faro buoy and the water levels were extracted from the tidal gage of Huelva (Spain). Large damages were reported after the storm, especially in Praia de Faro (Ferreira et al. 2019). In December 2019, two consecutive storms hit the area, Storm Daniel (December 16) and Storm Elsa (December 19–20). These storms traveled across the North Atlantic from east to west, and their effects were widely felt in several Western European countries. Storm Elsa created Hs up to 5.15 m and a 11-s peak period at the Faro buoy [Fig. 2(b)]. This storm coincided with neap tides, probably reducing the negative erosive effects of the storm. Maximum Hs and peak period during Storm Daniel were 3.85 m and 9 s, respectively.

Based on the analysis of Pires (1998), Storm Emma corresponds to a 16-year storm (Ferreira et al. 2019), Elsa is a 5-year storm, and Daniel is a 1-year storm. Moreover, at Praia de Faro, during Storm Emma, the collision regime was observed in sections where the dune system is higher, while in sections where the dune is weaker and lower, overwash occurred. Storm Elsa caused mainly the swash regime; the collision regime was only observed in sections with a limited berm, and only the swash regime occurred during Storm Daniel. Therefore, for the southern Portuguese coast, Storms Emma, Elsa, and Daniel can be defined as high-energy, mid-energy, and low-energy events, respectively.

To assess the ability of XBeach to model morphological changes during those storms, several topographic data sets were collected (Table 1). At Praia de Faro, after the peak of Storm Emma, on March 2, 2018, five profiles, F2–F6 [Fig. 1(d)], were surveyed with a DGPS (differential global positioning system). Within the COSMO program (PROGRAMA COSMO n.d.), an unmanned aerial vehicle (UAV) survey and a bathymetric survey were conducted at Praia de Faro in October 2018. The surveys covered the entire study site, including the nearshore area (13.5 m below MSL), beach face, dune system, and the urbanized area. The mean vertical error of the October 2018 survey along five profiles ranged between 0.14 and 0.20 m, considering the DGPS survey as the benchmark. Finally, another survey was conducted on December 20, 2019, to measure the elevation of four profiles (F1, F2, F3, and F5) after Storm Elsa.

On the Quarteira site, three surveys were conducted. First, on May 28, 2019, a Mavic 2 Pro UAV, was used to obtain a digital elevation model (DEM) and orthophoto map of the study area (Table 1). A total of 45 ground control points were marked at the fixed structures and at the beach to build the DEM. The elevation of the DEM was compared against four cross-shore profiles measured, and the estimated by root mean square error (RMSE) for the entire area was 0.09 m. The built DEM presented a resolution of 0.014 m, while the orthophoto had less than 1 cm of resolution. Second, two consecutive surveys were performed by using DGPS along nine profiles [Fig. 1(e)], on December 17, 2019, before Storm Elsa, and December 20, 2019, after Storm Elsa.

Two complementary sources were utilized to extract the elevation of the sea bottom in offshore areas of the model domains. The Agência Portuguesa do Ambiente (APA) surveyed the nearshore areas of the Quarteira site in 2018. The processed information covers from −1 to −8 m MSL. Moreover, a regional bathymetry of the entire southern Portuguese continental shelf was extracted from MIRONE (a Windows MATLAB-based framework tool) (Luis 2007) to cover the deepest region of the grid models. Both data sets have a resolution of 10 m.

The comparison of the different profiles measured along Praia de Faro [Fig. 1(d)], together with those reported in Ferreira et al. (2019) and Garzon et al. (2020), allowed defining a fully developed-berm and a weak-berm profile, which are representative of this site. The first one, with a beach face slope of 0.13, had a well developed berm around the 4.4 m elevation, and the dune toe was found at 5 m [Fig. 3(a)]. The second profile displayed a beach face slope of 0.12, and the dune toe was located at 4.6 m. The berm of this profile was less marked than the previous profile and was therefore named the weak berm [Fig. 3(b)]. The vertical datum of the aforementioned elevation data were MSL as well.

The Praia de Faro 2D grid had an extension of 3,000 m longshore and 3,900 m cross-shore [Fig. 1(c)]. The numerical grid, with a variable longshore and cross-shore grid cell spacing, was built using the OpenEarth tool (Van Koningsveld et al. 2010). The grid optimization and interpolation were made based on the bathymetry data from the COSMO program survey (PROGRAMA COSMO n.d.), and the continental shelf referred to MSL. The minimum cross-shore and alongshore resolution in the sub-aerial beach were 2 and 10 m, respectively. Using satellite imagery available via Google Earth, parking lots, infrastructures, and building locations were identified and superimposed on the grid. At those locations, grid cells were set as a nonerodible layer. Therefore, these grid cells cannot be eroded or destroyed. This has important implications for Praia de Faro, where the seawall impedes erosion to reach the urbanized areas during severe storms (Ferreira et al. 2019). The 1D simulations were performed in cross-shore profiles extracted from the 2D grid at specific locations (F2, F3, F4, F5, F6), and thus, they maintained the same cross-shore resolution of the 2D grid model.

In Quarteira, two 2D grids were built using the OpenEarth tool. These grids shared geometry, domain, and resolution but differed in the intertidal and the dry region elevation to account for two different initial topographies. Regarding the topography, two data sets were used: (1) the 2D unmanned aerial survey of May 28, 2019 (hereafter the Quarteira May grid), and (2) interpolation of the nine profiles measured on December 17, 2019 (hereafter the “Quarteira December grid”), into a 2D elevation model of the dry beach combined with the elevation of the groins measured in the May 2019 survey. The dimension of the grids was 1,100 m longshore by 5,800 m cross-shore [Fig. 1(c)]. The grids presented a variable resolution with a minimum resolution of 2 and 5 m in the cross-shore and longshore directions, respectively. The longshore resolution was reduced to better capture the groin geometry. The nearshore bathymetry measured by APA and the regional bathymetry of the continental shelf were merged to cover the offshore and nearshore areas. Rocky groins and infrastructures determined as nonerodible in the model were identified from the orthophoto. The grid model elevation was referred to MSL.

This study used the XBeach X 1.23.5526 version with lateral Neumann conditions in the nonlinear shallow water equations and cyclic wave boundary conditions. The “single_dir” option was selected to simulate the propagation of the directionally spread short-waves group in the 1D and 2D models. In the 2D approach, the mean wave direction was intermittently computed using the stationary solver within XBeach, and then it propagated the wave energy along these directions. Thus, it preserved the groupness of the wave, leading to higher forcing on the infragravity waves (Roelvink et al. 2018). In the 1D approach, “single_dir” used a single-direction bin, considering waves reaching normally to the coast and ignoring refraction (Deltares 2018). To support transferability, many of the XBeach free parameters that can impact the ability of the model to predict morphological changes were either maintained in their default value or substituted by values reported in the literature (Table 2). For instance, in the absence of wave measurements in the surf zone, the breaker parameter, “gamma,” was set to 0.56, as suggested by Callaghan et al. (2013). The parameters “gamma” (maximum ratio wave height to water depth) and “beta” (the breaker slope of the roller) were set to 2.364 and 0.138, respectively, following the values reported by Do et al. (2018), a study that investigated the modeling of dune erosion. To reduce the number of parameters to calibrate, the Manning’s coefficient associated with bed roughness was based on land-cover classification for sandy sediment (Garzon and Ferreira 2016; van der Lugt et al. 2019). Also, a morphological acceleration factor of 10 was applied following previous studies, such as McCall et al. (2010).

A SWAN model (Booij et al. 1999) was created for the southern Portuguese coast to propagate the wave conditions from the Faro buoy to the study sites [Fig. 1(b)]. The grid was curvilinear to match the variations of the bathymetry, with a varying resolution (from a few km to hundreds of m). SWAN outputs from the nonstationary computations were extracted at the offshore boundary of the XBeach computational grids, between water depths of 25 and 30 m [Fig. 1(c)]. The wave conditions considered in the XBeach simulations for each storm are displayed in Fig. 2. In the SWAN computations, water levels did not experience variations, and model parameters were maintained by default.

As previously expressed, three model parameters (facua, bermslope, and wetslope) that may have a strong influence on the morphodynamics were investigated. The parameter “facua” is highly relevant in the morphological module, as it mainly governs the net cross-shore sediment transport (Elsayed and Oumeraci 2017). As XBeach only simulates short-wave energy averaged over a wavelength (wave phase is not considered), the sediment advection velocity u_a responsible for stirring the sediment and transporting it to the shore must be approximated. Van Thiel de Vries (2009) proposed an expression to indirectly calculate that velocity as a function of a wave skewness parameter (S_k), a wave asymmetry parameter (A_s), the root mean square velocity (u_rms), and facua following the form:In shallow areas, with highly nonlinear waves, higher values for u_a are expected since the difference between S_k and A_s increases, and consequently there is larger onshore sediment transport due to the wave nonlinearity occurs.

Another parameter investigated affecting the cross-shore sediment transport was bermslope. Under the assumption that the intertidal beach slope remains relatively stable in intermediate and reflective beaches, Roelvink and Costas (2017) proposed a simple mechanism to address the challenge of simulating intertidal beach slopes in long-term simulations. Later, Roelvink et al. (2019) extended the method to two dimensions. In XBeach the sediment transport (S) is computed aswhere u_sed = depth-averaged sediment advection speed; h = water depth; and c = depth-average concentration:

where S_swash = corrected transport; f_swash = transport factor (15); z_b = bottom elevation; x = cross-shore distance; and the last term is the equilibrium slope near the waterline. This corrected transport term is only applied in the swash region, defined by a limit where wave height is larger than water depth (Roelvink et al. 2019).

During the collision regime, dune face erosion or slumping is predominantly triggered by a combination of infragravity swash runup on the previously dry dune face and the critical wet slope. In the model, this process is considered by an avalanching mechanism triggered when the infragravity incursions reach the dune front and the dune front becomes wet. The transition of the wet and dry grid cells is controlled by a user-specified water depth (“h_switch”). If the critical slope between two adjacent grid wet cells is exceeded, sediment is exchanged between these cells in the amount needed in order to bring the slope back to the critical slope. As a result, sediment is brought from the dry dune into the wet profile, where it is transported further seaward by the undertow and infragravity backwash (Roelvink et al. 2009):

The maximum dune erosion rate can be limited (“dz_max”). In this study, dz_max and critical dry slope used the default values, while the “h_switch” used the lowest limit within the recommended values (0.01), since the coarse material of the sites enhances water infiltration and soil saturation and thus reduces soil cohesion and, ultimately, dune resistance.

By varying the parameter values according to a set of combinations, a sensitivity analysis of the horizontal retreat and XBeach calibration and validation were performed (Fig. 4). Note that the default value of facua was not considered in this evaluation, since previous authors recommended higher values of facua for steep profiles (Vousdoukas et al. 2012b). The values proposed for bermslope were chosen based on values found in the literature (Roelvink et al. 2019) and the average beach face slope of Praia de Faro. For wetslope, the evaluated values covered the range of ±0.1 around the default value (Table 3). Fig. 4 illustrates the running test sequence.

First, Storm Emma was used to assess the model sensitivity, in terms of horizontal retreat of two 1D profiles: a fully developed berm and a weak berm (Figs. 3 and 4), by applying combinations of the parameters listed in Table 3. In total, 36 simulations were performed. For each run, the erosion indicator selected was the relative horizontal displacement computed at the three elevations depicted in Fig. 3 (3 m, berm crest, and dune toe) as the relative displacement with respect to the displacement simulated by using the parameter combination with the lowest facua and wetslope, and bermslope off. This combination was chosen as a benchmark because it was expected to simulate the largest retreat. Second, two 1D cross-shore models, Profiles F4 and F6 [Fig. 1(d)], were selected to reproduce the morphological changes induced by Storm Emma in Praia de Faro (Fig. 4). The profiles were selected due to the different morphologies that they exhibited for the same exposed beach. The model was calibrated using the parameters previously studied (Table 3), and the rest of the main settings are indicated in Table 2. Third, Storm Emma was used to validate the 2D model of Praia de Faro and to compare model discrepancies between the 1D and 2D approaches (Fig. 4). Moreover, further investigations are presented to better understand the role of the gravity and infragravity wave modeling approach (2D “single_dir,” 1D “single_dir,” and 1D “multi_dir”) in those discrepancies. Using Storm Emma as a calibration event, the model was validated for Storm Elsa (Fig. 4). For Storm Emma, the model was initiated on February 28, 2018, 00:00, and finished on March 2, 2018, 18:00 (Fig. 2). For Storm Elsa, the model started on December 18, 2019, 00:00, and finished on December 20, 2019, 18:00 (Fig. 2). Fourth, at the site of Quarteira, a 2D model calibration, using Storm Elsa, was carried out (Fig. 4); the 1D approach was not applied for this site. The parameters used for calibration are listed in Table 3 and the rest of the main settings in Table 2. After model calibration, the role played by the initial or prestorm topography was assessed by simulating Storms Elsa and Emma for both Quarteira May and Quarteira December topographies (Fig. 4). These two storms were included in the analysis to study the influence of storm energy on the model sensitivity to the input topography. In addition, the impact of Storm Emma on this site was qualitatively analyzed (Fig. 4).

Three typical and widely applied model skills were used: bias, RMSE, and Brier Skill Score (BSS). The bias is the difference, in meters, in central tendencies of the predicted or modeled elevation values, Z_modeled, and the measured elevation values, Z_measured, at each considered grid cell:

a positive bias means that the bed level is higher in the computed data set than in the measurements (erosion underprediction). The angled brackets indicate the average of the readings. The RMSE is the quadratic mean, in meters, of the differences between predicted values and measured values. The RMSE is determined asThe BSS provides an objective method for assessing the performance of morphological models. Conversely to the RMSE and bias skills, BSS is dependent on the profile morphology before the storm. The classification used (van Rijn et al. 2003) considers values between 0.8–1.0 as excellent, 0.6–0.8 as good, 0.3–0.6 as reasonable, 0–0.3 as poor, and below 0 as bad. The BSS can be computed as

where Z_initial = initial elevation.

Moreover, a new classification for the RMSE is proposed in this study, such that excellent represents a score lower than 0.25 m, good between 0.25 and 0.4 m, reasonable between 0.4 and 0.6 m, poor between 0.6 and 0.8 m, and bad larger than 0.8 m. This classification aims at complementing the van Rijn et al. (2003) classification. Furthermore, a coastal erosion indicator relevant for coastal risks such as the dune toe retreat (Ferreira et al. 2017) was computed to further analyze differences in model performance between the 1D and 2D schemes.

The 1D sensitivity analysis performed at the full-berm profile revealed that (1) when the wetslope was set to 0.2 [Fig. 5(a)], the horizontal retreat at 3 m was sensitive to bermslope and facua. The retreat decreased linearly with increasing values of facua for all bermslope conditions. Moreover, when the bermslope was deactivated, higher erosion was computed (for the same facua) than when it was activated. Furthermore, an increase of 0.02 in bermslope resulted in more than 20% less horizontal erosion. At the berm height, facua was similarly important in controlling erosion. On the other hand, differences between bermslope off and 0.10 were minimum, but bermslope set to 0.12 reduced the retreat between 15% and 20%. At the dune toe, variations in facua still had an almost linear relationship with horizontal erosion (except for bermslope 0.12), and the erosion can be reduced by 90% when setting facua to 0.3 when compared with 0.15. The parameter bermslope was partially important; a value of 0.10 provided almost a similar effect to deactivated, but 0.12 reduced the horizontal erosion, between 70% and 100% (null horizontal displacement). (2) When bermslope was deactivated [Fig. 5(b)], a higher wetslope reduced retreat; an increase of 0.05 produced 5%–10% approximately less erosion, regardless of the value of facua. At the berm height, the influence of wetslope was higher, and increases of 0.05 decreased erosion by 12%–18%. At the dune toe, wetslope played an essential role and with facua set to 0.15, wetslope of 0.25 and 0.3 can reduce the erosion 70% and 100%, respectively, in respect to wetslope equal to 0.2. (3) When bermslope was set to 0.10 [Fig. 5(c)], the model was not sensitive to changes on the wetslope at 3 m in height. Moreover, the sensitivity to wetslope at the berm height was lower than in the case when bermslope was deactivated. At the dune toe, the model was still sensitive to wetslope. Thus, with facua set to 0.15, wetslope of 0.25 and 0.3 can reduce the erosion by 40% and 70%, respectively, in comparison with wetslope equal to 0.2.

The weak-berm profile responded to changes at the tested parameters like the full-berm profile (Fig. 6), except at the dune toe, where the model still maintained a certain sensitivity to the three parameters but was less sensitive than the full-berm profile. While all runs performed on the weak-berm profile resulted in dune toe retreat, some runs performed on the full-berm profile simulated no dune retreat (Figs. 5 and 6).

A total of 36 simulations were performed for each profile at Praia de Faro (Table 3). For easier visualization of the calibration exercise, both the van Rijn et al. (2003) and the RMSE classifications were transformed into a color code, with excellent being represented by black and bad by white (Fig. 7). The skills BSS and RMSE of the two analyzed profiles, F4 and F6, were averaged and used to compare the performance of each simulation. The characteristics of the tested profiles were reasonably different (Fig. 8). F4 was shorter, with a beach face slope of 0.12, and it was considered nonerodible after 3,732 m (cross-shore distance). F6 had a similar slope (0.11), but its backshore was wider and connected to a 6.5-m MSL height dune (Fig. 8). The calibration results revealed that when bermslope was deactivated, facua required the highest value assessed (0.3) to classify the model results as excellent (Fig. 7), regardless of the value of wetslope. As facua was reduced, the performance of the model was reduced as well. When bermslope was set to 0.10, all simulations, in general, resulted in high scores, especially simulations with facua equal to 0.20 and 0.25. Erosion was underestimated for simulations with facua equal to 0.3 and overestimated when facua was set to 0.15 (not shown here). Finally, when bermslope was set to 0.12, it was observed that decreasing facua improved the skills of the model; the best score was found with “acua set to 0.15. Therefore, several combinations of facua, bermslope, and wetslope yielded excellent data for both skills, BSS and RMSE (Fig. 7; black cells represent the best performance). Among these 36 simulations, the skills of the simulation with the highest scores for bermslope off, 0.10, and 0.12 are presented in Table 4 as Simulations 1, 2, and 3, respectively. The simulations with bermslope activated obtained slightly better skills than the simulations with bermslope off. The simulation with the lowest scores (facua = 0.15; wetslope = 0.4; bermslope = off) is represented for comparison in Fig. 8(d), highlighting the sensitivity of the model for the cases assessed in this study.

As has been previously highlighted (Bugajny et al. 2013; Elsayed and Oumeraci 2017; Simmons et al. 2019; Splinter and Palmsten 2012; Vousdoukas et al. 2012b), the modeled horizontal erosion was sensitive to facua. In general, a linear relationship was found between increasing this parameter and decreasing the percentage of horizontal erosion in this study (Figs. 5 and 6). Similarly, van der Lugt et al. (2019) found that changes in erosion volume scaled linearly with facua. The model seemed more sensitive to changes in facua in the dune toe for both considered profiles at Praia de Faro. One of the reasons was that the actual magnitude of the horizontal retreat at the dune toe was lower than at 3 m or at the berm crest, and therefore small changes in the model results caused large differences in terms of percentage. This also explained the bigger sensitivity found in the full-berm profile, since in this profile, the magnitude of the horizontal retreat at the dune toe was smaller than at the weak-berm profile. For instance, variations in facua from 0.15 to 0.25 (wetslope = 0.2 and bermslope = off) can lead to changes at the dune toe retreat up to 2.2 m in both profiles. Thus, simulations with wetslope equal to 0.25 caused 50% and 30% less retreat than the benchmark case for the full-berm and weak-berm profile, respectively [Figs. 5(a) and 6(a)]. The model was also sensitive to wetslope, and increasing this parameter reduced the erosion, particularly in the dune toe [Figs. 5(b and c) and 6(b and c)], as was also observed by previous studies (Armaroli et al. 2013; Cho et al. 2019; Vousdoukas et al. 2012b). This was especially observed when bermslope was off. When bermslope was set to 0.10, the model was no longer sensitive to this parameter at 3 m height and slightly sensitive at the berm height, but it was still sensitive at the dune toe [Figs. 5(c) and 6(c)]. This would partially agree with Roelvink et al. (2019), since they stated that bermslope could replace the effect of wetslope. Also, an increase in bermslope produced an enhanced onshore sediment transport and the erosion decreased with bermslope set to 0.12 [Figs. 5(a) and 6(a)] as Roelvink et al. (2019) found in their study.

The 1D calibration proved that several model combinations of facua, wetslope, and bermslope produced excellent results for the evaluation of coastal erosion at steep beaches (Fig. 7). These results also confirmed that steep profiles required larger values of facua than did the default value to compensate for the onshore transport induced by the incident-band swash processes, which are not included in XBeach (Elsayed and Oumeraci 2017; Roelvink and Costas 2017; Vousdoukas et al. 2012b). Moreover, it is important to highlight that within the limits proposed on the RMSE-based and van Rijn et al. (2003) classifications, the sensitivity to wetslope was low (Fig. 7). In the modeled profiles, the performance of XBeach was only evaluated up to 3–4 m MSL, and the dune erosion was not considered (lack of measured data), where the avalanching and slumping processes controlled by wetslope are more significant. The addition of the bermslope parameter and the reduction of facua resulted in an excellent model prediction in the tested steep profiles with different berm morphologies. Thus, when bermslope was set to 0.10 and 0.12, facua values of 0.2 and 0.15 produced excellent results. Roelvink et al. (2019) also noticed that the combination of a moderate bermslope and low facua provided good results. Similary, Lashley et al. (2019) found that only when bermslope was activated, the model was able to reproduce the steep poststorm dune profile.

The same parameters were used to calibrate the 2D model in Praia de Faro, yielding several combinations of excellent results. As it was found in the 1D calibration, when bermslope was off and facua was set to 0.3, the model tended to provide accurate results (Table 4). However, when bermslope was included, excellent results were obtained for several combinations. This confirmed that the inclusion and use of bermslope were also adequate for a 2D XBeach model for steep beaches. Moreover, the use of bermslope reduced the variability between profiles, providing more robust results alongshore (Fig. 9). The poststorm profiles presented a relatively uniform behavior despite the existence of a prestorm alongshore variability, due to the presence of beach cusps (Vousdoukas et al. 2012a). It is important to note that the model parametrization providing the most accurate results used values that were close to the values recommended in the XBeach manual for facua and wetslope. This can avoid the model from behaving abnormally (e.g., excessive profile flattening around the waterline) as a consequence of unusual values on those model parameters, as observed in Fig. 9 and previously reported by other authors, such as Simmons et al. (2019). Furthermore, the value of bermslope chosen here was similar to the beach face slope found in the study area (Figs. 9 and 12). As the present work only focused on destructive processes induced by storms, these settings might not be appropriate to simulate constructive morphological processes in the long term, as found by Kombiadou et al. (2021).

The comparison between the 1D and 2D hydrodynamic and morphological outputs for the high-energy Storm Emma demonstrated that the higher erosion simulated by the 1D model was related to the higher infragravity wave height simulated in the 1D model (Fig. 10). Different model domains (1D or 2D) and different wave spreading approaches (“single_dir” and “multi_dir”) resulted in discrepancies in the simulated infragravity energy, Fig. 10(b). Roelvink and Reniers (2012) declared that the swash processes in the infragravity band play an essential role in the avalanching mechanism, one of the main factors in dune erosion. Thus, the authors hypothesized that the larger energy of the infragravity band of the 1D model might lead to an enhanced erosion of the dune with respect to the 2D model, as observed in Figs. 8 and 10(d). Roelvink et al. (2018) also stated that the 1D model with “single_dir” simulated higher runup in steep beaches than a 2D model using “single_dir” as well, confirming the results depicted in Fig. 10. In fact, they compared model results against field observations and noted that the 2D model was able to predict the runup while the 1D overestimated the measurements. On the other hand, the beach face seemed similarly eroded at the end of the storm with both the 1D and 2D approaches [Figs. 8 and 10(d)].

The good performance of the model in Praia de Faro to simulate the erosion caused by the mid-energy Storm Elsa demonstrated that the results obtained in a calibration process, for one storm, can be applied to successfully validate a second storm, even with different severity. This statement agrees with Simmons et al. (2019), who found that a second storm modestly improved the calibration results, and suggested that coastal practitioners should focus more on collecting data in more locations rather than collecting data for several storms. Nevertheless, the findings of this site cannot be transferred to another site located just 10 km away (Quarteira). XBeach applied to Quarteira was also successfully validated, quantitatively and qualitatively, for Storms Emma and Elsa, but the parametrization providing good results in Praia de Faro failed in predicting the erosion in Quarteira. Conversely to the Praia de Faro model, bermslope did not provide satisfactory results in Quarteira. The model faced large problems to simulate morphological changes around the groins, and highly overestimated the erosion. When bermslope was off, it can be observed that the value of facua in this site was considerably lower than in the Praia de Faro model. It is not totally understood why the model required lower onshore transport, since the profiles also presented a steep beach face slope. Some of the reasons can be associated with the presence of the groins. Bugajny et al. (2013) carried out several 1D XBeach model simulations on unprotected, protected, and heavily protected coasts and noticed a decrease in facua with an increase in the engineering protection, reporting values of facua of 0.5, 0.3, and 0.1, respectively. Further investigations of the role of facua in protected coasts are suggested to shed some light on this process. Moreover, the good skills of the model simulating erosion under the impact of two very different storms (high-energy and mid-energy) demonstrated these calibration parameters produce adequate results for both hydrodynamic conditions.

Elevation measurements immediately before the storm in the surf zone area and the dry beach are not always available and model experiments might have to use data only collected weeks or months prior to the storm. This gap in time might lead to some errors or uncertainties when evaluating morphodynamic models. The comparison between the outputs from the Quarteira May and December grids proved that the model was largely more sensitive to input parameters than initial topography, mainly if the initial different morphologies still represented similar beach volumes, which is in line with other studies, such as Armaroli et al. (2013). This was especially evident for the high-energy storm (Emma) when compared with the mid-energy storm (Elsa).

When implementing the XBeach model on an unexplored steep beach, the first efforts should be focused on collecting field data to validate and calibrate the morphodynamic model. This has special relevance in reflective beaches, since all wave processes occurring in the swash zone of these types of environments are not totally reproduced by models, namely XBeach (surfbeat). Thus, a default model parametrization could provide completely erroneous results. In the case of the lack of data for model validation and calibration, modelers working on exposed steep profiles should pay attention to facua and bermslope. In these beaches, values of facua close to the default value only perform properly when bermslope is activated, and acceptable results can be obtained by using several combinations of these two parameters. A tentative value for bermslope can be the beach face slope of the profiles. Moreover, if a dune system is present, the wetslope value must be carefully chosen but it might be close to the recommended value in the XBeach manual. On the other hand, if the reflective site presents engineering protection structures, such as groins, the impact of the wave asymmetry sediment transport is lower and, consequently, values of facua must be reduced. The use of the bermslope in these protected sites must be carefully assessed, as they might largely overestimate the erosion. Furthermore, the model would not be required to be calibrated against several storms, for both landscapes (exposed and groin-protected sites), which is an advantage for designing and predictive purposes.

While the 1D and 2D models simulate similar erosion in the beach face (using the same model parametrization), the erosion predicted at the upper beach and in the dune by a 1D cross-shore model was higher than the erosion computed on a 2D grid. Therefore, if the dune erosion is a major concern, different model parametrization is required in 1D and 2D models. This is especially important for the design of coastal interventions, namely if dealing with risk. Thus, at an initial phase, both approaches should be tested against observations to find the optimum parametrization for the 1D and 2D models. If the computational power is a major limitation, then multiple 1D simulations are preferable to a 2D model. Moreover, “single_dir” is more efficient than the “multi_dir” approach (Roelvink et al. 2018).

A systematic collection of elevation data can be very challenging. Thus, numerical models usually rely on a static initial topo-bathymetry, that is the elevation of the grid is not updated. In highly dynamic environments such as sandy beaches, this might result in some limitations, adding some uncertainties to the predicted morphological changes. However, even for different prestorm profile morphologies, if they maintain approximately the same amount of sand volume, the predicted erosion (poststorm shoreline and berm position) on those different morphologies is similar, mostly for highly energetic storms. This would imply that the topography of the model does not need to be periodically updated, particularly when taking into consideration the impact of storms with high energy or return period. Furthermore, steep beaches exhibit a rapid poststorm recovery response, and they are able to gain a large part of the eroded material in the order of weeks, supporting the idea that the static initial elevation approach would be sufficient to obtain reliable model results.