Modeling of coastal erosion in exposed and groin-protected steep

10 Process-based models are suitable tools for reproducing storm-driven erosion. However, their performance has been mainly 11 examined on mild-slope sandy beaches and their use on steep beaches represents still a challenge. Here, XBeach experiments were 12 combined with topographical measurements collected for two storms (16-yr and 5-yr return period) to obtain a reliable model. The model 13 parameters “facua” (parameterized wave asymmetry and skewness sediment transport component), “bermslope” (upslope transport term 14 for semi-reflective beaches), and “wetslope” (critical avalanching submerged slope) were utilized for calibration and validation. The 16-yr 15 storm simulations on an exposed beach revealed that whether “bermslope” increased, “facua” must be reduced, and vice versa, to properly 16 simulate erosion. Adding “bermslope” provided excellent results for these storms when using “facua” and “wetslope” values close to the 17 recommended values. In a groin-protected site, XBeach was successfully calibrated and validated for the tested storms using these 18 parameters, although with different values. These experiments demonstrated that the appropriate use of these parameters can satisfactorily 19 simulate morphological changes on steep beaches for different hydrodynamic conditions and coastal settings (exposed and groin-20 protected


22
Sandy coasts are among the most populated areas worldwide and they host a large number of socio-economic activities.

23
However, these environments are susceptible to the impact of coastal storms, with storm surges and waves generated during energetic 24 events causing severe erosion and shoreline retreat. Moreover, this problem might be exacerbated by rising sea levels and changes in 25 storminess. Under these threats, predicting erosion due to extreme oceanic events is essential to improve coastal management and 26 2 implement mitigation measures (e.g. early warning systems) that contribute to avoiding or minimizing their effects in socioeconomic 27 activities and ecosystems services. Among the different alternatives to support engineers and decision-makers, morphodynamic numerical 28 models are increasing their popularity. For instance, the open-source process-based model XBeach (Roelvink et al. 2009) has been applied 29 and validated in many sandy beaches worldwide impacted by severe tropical and extratropical storms. This model solves wave breaking, 30 surf and swash zone processes, dune erosion, and overwashing in a one-dimensional or a horizontally two-dimensional computational grid 31 under the assumption of a saturated surf zone, mainly occurring in dissipative beaches.

32
Dissipative beaches present a mild slope in the intertidal region and wave processes are dominated by skewed and asymmetric 33 short-crested waves, undertow, and infragravity waves. Several modeling studies have focused on these types of environments (e.g. slope is forced to stay close to a given slope, "bermslope". Therefore, an onshore transport term that is proportional to the difference 53 between the actual slope and a prescribed "bermslope" is added in the swash zone. While this transport term was initially thought for 54 3 improving the reliability of the XBeach model in the long-term simulations, Roelvink et al. (2019) demonstrated that this new addition has 55 a beneficial effect on the profile evolution during storm events on steep beaches. Moreover, they suggested that "bermslope" could 56 minimize the importance of other slope parameters implemented in the model such as "wetslope", which establishes the critical bed-slope 57 of the wet profile before the initiation of avalanching. Cho et al. (2019) declared that XBeach is more sensitive to changes in "wetslope" 58 values in a steep profile than in a mild profile, and hence, the selection of this parameter should be carefully considered in these beaches.

59
In terms of computational effort, using a 1D approach to calibrate a model instead of a 2D model would be more efficient since 60 it would allow for a more rapid assessment of all free model parameters, especially for highly complex, process-based models. The 61 downside of this approach is that the findings of a 1D calibration model might not be directly applied to a detailed 2D model. Cross-shore 62 profile models have the inherent limitation of longshore uniformity, and they are not capable of capturing the effects of longshore transport 63 gradients. Conversely, the 2D models incorporate longshore variations, and they are not limited to straight-line coastal systems, being able designing an EWS; yet, the downside of such an approach was not fully investigated.

68
While process-based models seek to explicitly represent the crucial physical dynamics, in practice, they include semi-empirical 69 parameterizations to improve their efficiency. This results in a large number of free model parameters to calibrate, especially for coastal 70 morphologies where research efforts have been less intense (e.g. steep beaches). The constant evolution of the models implies the 71 necessity of continuous calibration and validation. Also, calibration parameters vary according to the specific characteristics of each site.

72
The present study has a twofold objective: firstly, to perform a sensitivity analysis of the main morphological parameters used for XBeach 73 calibration ("facua", "bermslope", and "wetslope"); and secondly, to use the above results to calibrate and validate XBeach for two storm 74 events with different severity at two steep beaches with different levels of human intervention. The results will contribute to obtaining 75 better model performances in such environments. Moreover, in cases where the available data for model calibration and validation are 76 limited, recommendations for the model implementation are indicated.

79
Praia de Faro is an sandy beach located in the barrier island system of Ria Formosa, on the southern Portuguese coast (Fig. 1I -80 III). The barrier splits the Atlantic Ocean on the front side from a coastal lagoon on the backside (Fig. 1 IV). The study site is an elongated 81 4 peninsula orientated 130º approximately (measure from the north). This beach is subject to significant urban development pressure and the 82 oceanfront is partially stabilized with rocks/walls or naturally protected by a dune. The dune morphology varies alongshore with higher 83 robustness and elevation (7-8 m above mean sea level, MSL) at the western portion of the study area (F6 in Fig. 1 IV), while at the central 84 and eastern part, the dune is lower (6-7 m above MSL) and weaker. The central part, F4 and F5 in Fig. 1  The town of Quarteira is located ten kilometers from Praia de Faro towards the northwestern direction ( Fig. 1 III). The analyzed 95 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 96 the average beach slope is 0.10. Sediment grain size is slightly finer than in Praia de Faro and d50 and d90 are 485µm and 900 µm 97 respectively. The three beaches are laterally limited by 150 m length rocky groins (Fig. 1 V). These groins make the three beaches behave 98 like "manmade pocket beaches" subject to beach rotation as a function of wave direction. During the dominant conditions, the groins 99 maintain the sediment in the system; however, during very energetic conditions the sediment can fall outside of the system. Beach 100 nourishments have been performed at the area to guarantee a reasonable beach width for bathing conditions, with a total of 360000 m 3 101 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.

102
While Quarteira represents a relevant touristic destination in Portugal, the beach morphodynamic has been poorly investigated.

105
Two storms impacting the area have been considered for analysis. Storm Emma, in March 2018, was a severe storm that traveled 106 towards the northeastern direction (Ferreira et al. 2019). The maximum significant wave height, Hs, measured at the Faro buoy during this 107 storm was 6.55 m, with a peak period and direction of 13 s and 240º N respectively (Fig. 2). Moreover, the timing of the storm matched 108 5 with a spring tide, exacerbating the impact of the storm. Large damages were reported after the storm, especially in Praia de Faro (Ferreira 109 et al. 2019). In December 2019, two consecutive storms hit the area, storm Daniel (December, 16) and storm Elsa (December, 19-20).

110
These storms traveled across the North Atlantic from east to west and their effects were widely felt in several western European countries.

111
Storm Elsa created Hs up to 5.15 m and 11 s peak period, at the Faro buoy (Fig. 2). This storm coincided with neap tides, probably 112 reducing the negative erosive effects of the storm. Maximum Hs and peak period during Daniel were 3.85 m and 9 s, respectively.

139
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 140 (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 141 marked than the previous profile being therefore named weak-berm (Fig. 3 b). The vertical datum of the aforementioned elevation data 142 was MSL as well. 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.

153
In Quarteira, two 2D grids were built using the OpenEarth tool. These grids shared geometry, domain, and resolution but 154 differed in the intertidal and the dry region elevation to account for two different initial topographies. Regarding the topography, two was 1100 m longshore by 5800 m cross-shore ( Fig. 1 III). The grids presented a variable resolution with a minimum resolution of 2 m and 159 5 m in the cross-shore and longshore directions. The longshore resolution was reduced to better capture the groin geometry. The nearshore 160 bathymetry measured by APA and the regional bathymetry of the continental shelf were merged to cover offshore and nearshore areas.
7 Rocky groins and infrastructures determined as non-erodible in the model were identified from the orthophoto. The grid model elevation 162 was referred to MSL.

163
This study used the XBeachX 1.23.5526 version with lateral Neumann conditions in the non-linear shallow water equations and 164 cyclic wave boundary conditions. The "single_dir" option was selected to simulate the propagation of the directionally-spread short waves 165 group in the 1D and 2D models. In the 2D approach, the mean wave direction was intermittently computed using the stationary solver 166 within XBeach and then, it propagated the wave energy along these directions. Thus, it preserved the groupness of the wave, leading to 167 higher forcing on the infragravity waves (Roelvink et al. 2018). In the 1D, "single_dir" used a single directional bin, considering waves As previously expressed, three model parameters ("facua", "bermslope", and "wetslope") that may have a strong influence on 183 the morphodynamics were investigated. The parameter "facua" is very relevant in the morphological module as it mainly governs the net 184 cross-shore sediment transport (Elsayed and Oumeraci 2017). As XBeach only simulates short wave energy averaged over a wavelength 185 (wave phase is not considered), the sediment advection velocity ua responsible for stirring the sediment and transporting it to the shore 186 must be approximated. Van Thiel de Vries (2009) proposed an expression to indirectly calculate that velocity as a function of a wave 187 skewness parameter (Sk), a wave asymmetry parameter (As), the root mean square velocity (urms) and "facua" following the form: In shallow areas, with highly nonlinear waves, higher values for ua are expected since the difference between Sk and As increases and 190 consequently, larger onshore sediment transport due to the wave nonlinearity occurs.

191
Another parameter investigated affecting the cross-shore sediment transport was "bermslope". Under the assumption that the where used is the depth-averaged sediment advection speed, h the water depth and c the depth-average concentration.
where Sswash is the corrected transport, fswash a transport factor (15), zb represents the bottom elevation and x the cross-shore distance, and 199 the last term is the equilibrium slope near the waterline. This corrected transport term is only applied in the swash region, defined by a 208 209 The maximum dune erosion rate can be limited ("dzmax"). In this study, "dzmax" and critical dry slope used the default values while the 210 "hswitch" used the lowest limit within the recommended values (0.01) since the coarse material of the sites enhances water infiltration and 211 soil saturation, and thus reduces soil cohesion, and ultimately dune resistance.

212
By varying the parameter values according to a set of combinations, a sensitivity analysis of the horizontal retreat and XBeach 213 calibration and validation were performed (Fig. 4). Note that the default value of "facua" was not considered at this evaluation since 214 previous authors recommended higher values of "facua" for steep profiles (Vousdoukas et al., 2012b). The values proposed to 9 For "wetslope", the evaluated values covered the range of ± 0.1 around the default value (Table 3). The running test sequence is illustrated 217 in Fig. 4.

218
Firstly, storm Emma was used to assess model sensitivity, in terms of horizontal retreat of two 1D profiles, a fully-developed 219 berm and a weak-berm ( Fig. 3 and Fig. 4), by applying combinations of the parameters displayed in Table 3. In total, 36 simulations were 220 performed. For each run, the erosion indicator selected was the relative horizontal displacement computed at the three elevations depicted 221 in Fig. 3 (3 m, berm crest, and dune toe) as the relative displacement with respect to the displacement simulated by using the parameter 222 combination with the lowest "facua" and "wetslope" and "bermslope" off. This combination was chosen as a benchmark because it was 223 Fig. 1 IV) were selected to reproduce the morphological changes induced by storm Emma in Praia de Faro (Fig. 4). The profiles were 224 selected due to the different morphologies that they exhibited for the same exposed beach. The model was calibrated using the parameters 225 previously studied (Table 3) and the rest of the main settings are indicated in Table 2. Thirdly, storm Emma was used to validate the 2D 226 model of Praia de Faro and to compare model discrepancies between 1D and 2D approaches (Fig. 4). Moreover, further investigations are 227 presented to better understand the role of the gravity and infragravity wave modeling approach (2D "single_dir", 1D "single_dir", and 1D 228 "multi_dir") in those discrepancies. Using storm Emma as a calibration event, the model was validated for storm Elsa (Fig. 4). Fourthly, 229 in the site of Quarteira, a 2D model calibration, using storm Elsa, was carried out (Figure 4); the 1D approach was not applied for this site.

230
The parameters used for calibration are displayed in Table 3 and the rest of the main settings in Table 2. After model calibration, the role 231 played by the initial or pre-storm topography was assessed by simulating storm Elsa and Emma for both Quarteira May and Quarteira 232 December topographies (Fig. 4). These two storms were included in the analysis to study the influence of storm energy on the model 233 sensitivity to the input topography. Also, the impact of storm Emma on this site was qualitatively analysed (Fig. 4). 234

235
Three typical and widely applied model skills were used: bias, RMSE and BSS. The bias is the difference, in meters, in central 236 tendencies of the predicted or modeled elevation values, Zmodeled and the measured elevation values, Zmeasured, at each considered grid cell. 237 242 The Brier Skill Score, BSS, provides an objective method for assessing the performance of morphological models. Conversely to 243 the RMSE and bias skills, BSS is dependent on the profile morphology before the storm. The classification used, from van Rijn et al. 247 where Zinitial represents the initial elevation. The 1D sensitivity analysis performed at the full-berm profile revealed that: 1) when the "wetslope" was set to 0.2 (upper panels 256 in Fig. 5), the horizontal retreat at 3 m was sensitive to "bermslope" and "facua". The retreat decreased linearly with increasing values of 257 "facua" for all "bermslope" conditions. Moreover, when the "bermslope" was deactivated, higher erosion was computed (for the same 258 "facua") than when it was activated. Furthermore, an increase of 0.02 in "bermslope" resulted in more than 20% less horizontal erosion. At 259 the berm height, "facua" was similarly important in controlling erosion. On the other hand, differences between "bermslope" off and 0.10 260 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 261 almost linear relationship with horizontal erosion (except for "bermslope" 0.12), and the erosion can be reduced by 90% when setting 262 "facua" to 0.3 when compared to 0.15. The parameter "bermslope" was partially important; a value of 0.10 provided almost a similar 263 effect to deactivated, but 0.12 reduced the horizontal erosion, between 70 and 100% (null horizontal displacement). 2) When "bermslope" 264 was deactivated (middle panels Fig. 5), higher "wetslope" reduced retreat; an increase of 0.05 produced 5-10% approximately less erosion, 265 regardless of the value of "facua". At the berm height, the influence of "wetslope" was higher, and increases of 0.05 decreased erosion by 266 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 267 erosion 70 and 100%, respectively, in respect to "wetslope" equal to 0.2. 3) when "bermslope" was set to 0.10 (lower panels in Fig. 5), the 268 11 model was not sensitive to changes on the "wetslope" at 3 m height. Moreover, the sensitivity to "wetslope" at the berm height was lower 269 than in the case when "bermslope" was deactivated. At the dune toe, the model was still sensitive to "wetslope". Thus, with "facua" set to 270 0.15, "wetslope" of 0.25 and 0.3 can reduce the erosion by 40 and 70%, respectively, in comparison to "wetslope" equal to 0.2.

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

276
A total of 36 simulations were performed for each profile at Praia de Faro (see Table 3). The characteristics of the tested profiles 277 ( Fig. 7) were reasonably different: F4 was shorter with a beach face slope of 0.12, and it was considered as non-erodible after the 3732 m 278 (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.  (Table 4), regardless of the value of "wetslope". As 283 "facua" was reduced, the performance of the model was reduced as well. When "bermslope" was set to 0.10, all simulations, in general, 284 resulted in high scores, especially simulations with "facua" equal to 0.20 and 0.25. Erosion was underestimated for simulations with 285 "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 286 observed that decreasing "facua" improved the skills of the model; the best score was found with "facua" set to 0.15. Therefore, several 287 combinations of "facua", "bermslope", and "wetslope" yielded excellent for both skills, BSS and RMSE (see Table 4 where a cell entirely 288 black represents the best performance). Among these 36 simulations, the skills of the simulation with the highest scores for "bermslope" 289 off, 0.10, and 0.12 are presented in Table 5 -i, -ii, and -iii respectively. The simulations with "bermslope" activated obtained slightly 290 better skills than the simulations with "bermslope" off. The simulation with the lowest scores ("facua" = 0.15; "wetslope" = 0.4; 291 "bermslope" = off) is represented for comparison in Fig. 7 IV), highlighting the sensitivity of the model for the cases assessed in this 292 study.  Table 5 -i, -ii, and -iii were plotted against the measurements taken after storm 296 Emma in Fig. 7 I -III (note that post-storm Emma measurements did not cover the dune in F6). The morphological changes simulated by 297 the 2D model using the same values for "facua", "wetslope" and "bermslope" as the 1D model were represented as well (Fig. 7 I -III). 298 The comparison between the 1D and 2D model approaches showed that the modeled erosion at the beach face was, in general, slightly 299 higher 300 in the 1D approach, especially in F4 (Fig. 7 I -III). However, this comparison indicated that the erosion of the dune in F6 was 301 considerably larger in the 1D model than in the 2D model, with dune toe retreat differences up to 7 m in some runs (e.g. Table 5 -ii). In 302 general, the average BSS and RMSE (F4 and F6) for the 1D and 2D models were similar (Table 5). Among all setting combinations, and 303 after plotting several 2D results (not shown here), the best morphological representation of the storm-induced erosion was provided by the 304 simulation with "facua" = 0.15, "wetslope" = 0.2, and "bermslope" = 0.12, see Fig. 8. It does not exhibit entirely the best skills (Table 5 -305 v), but these settings replicated more accurately the observed dune retreat (Garzon et al. 2021) while simulating also correctly the erosion 306 in the beach face. Furthermore, it is important to highlight two aspects: 1) when including "bermslope", the erosion in all profiles was 307 reasonably well predicted (Fig. 8), including F3, whose initial profile showed a lower sand volume (beach cusp trough) compared for 308 instance with F2 (beach cusp crest). Conversely, when "bermslope" was deactivated the modeled erosion was overestimated (~ 1 m of 309 vertical erosion) in F3 (Fig. 8). Thus, including "bermslope" produced, in general, more consistent results along the five profiles; and 2) 310 the simulation with "bermslope" deactivated displayed a milder slope below 0 m MSL than the one with "bermslope" (Fig. 8). The value 311 of "bermslope" was similar to the actual beach face slopes found in the profiles F2-F6 that ranged between 0.09 and 0.12 (Fig. 8). 312 To further investigate discrepancies on the modeled erosion between the 1D and 2D approaches for storm Emma, the 313 hydrodynamic (wave height, sea level and infragravity wave height) and morphological outputs were plotted (Fig. 9). The two models, 1D 314 and 2D (using single direction) simulated similar wave height (Fig. 9 a) and sea level (Fig. 9 c), but the infragravity wave height computed 315 on the 1D model was higher than on the 2D model (Fig. 9 b). When comparing the final profile ( Fig. 9 d), it can be concluded that dune 316 erosion was enhanced on the 1D model with regards to the 2D. Also, the outputs of the 1D model using both the single direction and 317 multiple direction options were contrasted in Fig. 9. The multiple direction approach computed slightly lower infragravity wave heights 318 than the single direction approach (but still higher than the 2D approach) as seen in Fig. 9 b, while the rest of the hydrodynamic variables 319 were similar (Fig. 9 a, c). When comparing the post-storm profiles in Fig. 9 d, it was observed that the three approaches simulated similar 320 13 erosion at the beach face (below 3 m above MSL) but large differences were found in the dune face. For instance, the dune retreat at 5 m 321 above MSL obtained on the 2D model, 1D approach -multiple directions and 1D approach -single direction was 0, 3 and 6 m 322 respectively, Fig. 9

d. 323
The spatial prediction of the seabed change for Praia de Faro is displayed in Fig. 10, where the red color represents erosion, blue 324 deposition and yellow minimal changes. Maximum vertical erosion of 2.5 m in the beach face was simulated but this erosion was not 325 uniform alongshore as the result of the alongshore variability induced by the presence of the beach cusps (Fig. 10). Also, Fig. 10 326 demonstrated the effect of the non-erodible layer to hinder the erosion in the urbanized area. Moreover, the transition between the erodible 327 and non-erodible regions was correctly simulated (Fig. 8 and Fig. 10).

337
The calibration was carried out on the Quarteira December grid and involved the same parameters presented before (Table 3).

338
The modeled profiles displayed a well-developed berm, especially on the eastern and central profiles of each pocket beach (Q1, Q2, Q4, 339 and Q5) and a beach face slope between 0.08 and 0.12 (Fig. 12). Model experiments (not shown here), revealed that conversely to Praia de 340 Faro, "bermslope" caused unrealistic overprediction of the erosion in the groin heads (updrift), and hence, this parameter was deactivated 341 for the rest of the simulations in this site. Thus, the best calibration established "facua" and "wetslope" set to 0.15 and 0.20 respectively, 342 and "bermslope" off. The BSS and RMSE scores for the eastern and central profiles of each pocket beach were classified as excellent, 343 based on the previously discussed classifications, as displayed in Table 6. It is important to point out that the model reduced its ability to 344 predict the morphology of the profiles immediately downdrift of the groins (Q3, Q6, Q9; see Fig. 12 and Table 6). 345

346
To evaluate the sensitivity of the model to the initial topography, storm Emma was simulated in the Quarteira site under the two 347 topography conditions, May and December of 2019. There were differences between the initial profiles in May and December (Fig. 13)  348 and for instance, the berm crest was notorious in the latter survey, while in May the transition from the berm to the beach face was 349 smoother, particularly in Q1, Q2, Q4 and Q5. Moreover, in the May survey, mainly in the eastern profiles of each pocket beach (Q1, Q4, 350 and Q7), the slope was milder and the volume of sediment above MSL was larger than in the December survey. The numerical settings 351 were equal to those used to calibrate storm Elsa. XBeach results for storm Emma showed that the outputs from the two considered grids,

352
May and December, were close (Fig. 13) and, for instance, the onshore limit of the eroded profile was similar for both initial conditions. 353 Even in profiles with clearly different morphologies such as Q1, Q4, and Q7, the final computed impact of the storm was almost 354 equivalent.

355
When simulating storm Elsa on the Quarteira May grid, XBeach was still able to fairly reproduce the erosion as well, even if the 356 profiles were different from the actual pre-storm morphology. However, the RMSE increased in almost all profiles suggesting a lower 357 performance. The largest differences in model performances were found in Q1, Q6, Q8, and Q9 (Table 6). In all profiles, the positive bias 358 increased for the May grid simulation with regards to the December grid, indicating that erosion underestimation was higher on the May 359 grid. The model ability to simulate erosion at profiles immediately downdrift of the groins (Q3, Q6, Q9) on the May grid was also low.

361
Visual inspections after storm Emma revealed that erosion barely reached the urbanized area at Quarteira and that morphological 362 changes were more significant in the eastern side of each pocket beach (downdrift). It was hypothesized that the typically eastward 363 longshore transport before the storm might have accumulated more material updrift of the groins. Thus, as observed in Fig. 13, if there 364 was more available sediment, the erosion was higher in these regions tending to reach a similar equilibrium profile in the area. The 365 modeled results of storm Emma for Quarteira showed that the groins had a clear influence on the morphodynamic of this site, and the 366 erosion (volume and berm retreat) was more significant immediately updrift of the groins (Q1, Q4, and Q7), with a maximum vertical 367 erosion of 2.5 m (Fig. 14). The definition of these structures as non-erodible at the model settings allowed successfully replicating the 368 morphodynamic of this site. Thus, as it was also observed in the Praia de Faro simulations, the non-erodible layer implemented in the 369 model behaved properly (Fig. 14). The upper section of the eastern and central beaches was not eroded, while the erosion at the western 370 beach extended almost to the promenade (Fig. 14). It, in general, matched the field (visual) observations, confirming the positive 371 performance of the model. found between increasing this parameter and decreasing the percentage of horizontal erosion in this study (Fig. 5 and Fig. 6). Similarly, 377 van der Lugt et al. (2019) found that changes in erosion volume scaled linearly with "facua". The model seemed more sensitive to changes 378 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 379 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 380 in terms of percentage. This also explained the bigger sensitivity found in the full-berm profile, since in this profile, the magnitude of the 381 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 382 ("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 383 "wetslope" equal to 0.25 caused 50% and 30% less retreat than the benchmark case for the full-berm and weak-berm profile respectively 384 (upper panels in Fig. 5 and Fig. 6). The model was also sensitive to "wetslope", and increasing this parameter reduced the erosion, 385 particularly in the dune toe (mid and lower panels in Fig. 5 and Fig. 6), as was also observed by previous studies (Armaroli et al. 2013;386 Cho et al. 2019; Vousdoukas et al. 2012b). This was especially observed when "bermslope" was off. When "bermslope" was set to 0.10, 387 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 388 dune toe (lower panels in Fig. 5 and Fig. 6). This would partially agree with Roelvink et al. (2019) since they stated that "bermslope" 389 could replace the effect of "wetslope". Also, an increase in "bermslope" produced an enhanced onshore sediment transport and the erosion 390 decreased with "bermslope" set to 0.12 (upper panels in Fig. 5 and Fig. 6) as Roelvink et al. (2019) found in their study. 391 The 1D calibration proved that several model combinations of "facua", "wetslope", and "bermslope" produced excellent results 392 for the evaluation of coastal erosion at steep beaches (Table 4). These results also confirmed that steep profiles required larger values of 393 "facua" than the default value to compensate for the onshore transport induced by the incident-band swash processes, which are not  (Table 4). In the modeled profiles, the performance of XBeach was only evaluated up to 3 -4 m MSL, and the dune erosion was not 397 considered (lack of measured data), where the avalanching and slumping processes controlled by "wetslope" are more significant. The 398 addition of the "bermslope" parameter and the reduction of "facua" resulted in an excellent model prediction in the tested steep profiles 399 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 400 16 results. Roelvink et al. (2019) also noticed that the combination of a moderate "bermslope" and low "facua" provided good results.

401
Similary, Lashley et al. (2019) found that only when "bermslope" was activated, the model was able to reproduce the steep post-storm 402 dune profile.

404
The same parameters were used to calibrate the 2D model in Praia de Faro yielding several combinations of excellent results. As 405 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 406 (Table 5). However, when "bermslope" was included, excellent results were obtained for several combinations. This confirmed that the 407 inclusion and use of "bermslope" were also adequate for a two-dimensional XBeach model for steep beaches. Moreover, the use of 408 "bermslope" reduced the variability between profiles, providing more robust results alongshore (Fig. 8). The post-storm profiles presented 409 a relatively uniform behavior despite the existence of a pre-storm alongshore variability due to the presence of beach cusps (Vousdoukas et 410 al. 2012a). It is important to note that the model parametrization providing the most accurate results used values that were close to the 411 values recommended in the XBeach manual for "facua" and "wetslope". This can avoid the model behaving abnormally (e.g. excessive 412 profile flattening around the waterline) as a consequence of unusual values on those model parameters, as observed in Fig. 8 and  413 previously reported by other authors such as Simmons et al. (2019). Furthermore, the value of "bermslope" chosen here was similar to the 414 beach face slope found in the study area ( Fig. 8 and Fig. 11). As the present work only focused on destructive processes induced by 415 storms, these settings might not be appropriate to simulate constructive morphological processes in the long term as found by Kombiadou

417
The comparison between the 1D and 2D hydrodynamic and morphological outputs for the high-energy storm Emma 418 demonstrated that the higher erosion simulated by the 1D model was related to the higher infragravity wave height simulated on the one-419 dimensional model (Fig. 9). Different model domains (1D or 2D) and different wave spreading approaches ("single_dir" and "multi_dir") 420 resulted in discrepancies in the simulated infragravity energy, Fig. 9 b. Roelvink and Reniers (2012) declared that the swash processes in 421 the infragravity band play an essential role in the avalanching mechanism, one of the main factors in dune erosion. Thus, the authors 422 hypothesized that the larger energy of the infragravity band of the 1D model might lead to an enhanced erosion of the dune in respect to 423 the 2D model, as observed in Fig. 7 and Fig. 9 d. Roelvink et al. (2018) also stated that the 1D model with "single_dir" simulated higher 424 runup in steep beaches than a 2D model using "single_dir" as well, confirming the results depicted in Fig. 9. In fact, they compared model 425 results against field observations and noted that the 2D model was able to predict the runup while the 1D overestimated the measurements.

426
On the other hand, the beach face seemed similarly eroded at the end of the storm with both the 1D and 2D approaches (Fig. 7 and Fig. 9  427 d.)

428
The good performance of the model in Praia de Faro to simulate the erosion caused by the mid-energy storm Elsa demonstrated 429 that the results obtained in a calibration process, for one storm, can be applied to successfully validate a second storm, even with different 430 severity. This statement agrees with Simmons et al. (2019) that found that a second storm modestly improved the calibration results and 431 suggested that coastal practitioners should focus more on collecting data in more locations rather than collecting data for several storms.

432
Nevertheless, the findings of this site cannot be transferred to another site located just 10 km away (Quarteira). XBeach applied to 433 Quarteira was also successfully validated, quantitatively and qualitatively, for storms Emma and Elsa, but the parametrization providing

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

463
While the 1D and 2D models simulate similar erosion in the beach face (using the same model parametrization), the erosion 464 predicted at the upper beach and in the dune by a 1D cross-shore model was higher than the erosion computed on a two-dimensional grid.

465
Therefore, if the dune erosion is a major concern, different model parametrization is required in a one and two-dimensional model. This is 466 especially important for the design of coastal interventions, namely if dealing with risk. Thus, at an initial phase, both approaches should 467 be tested against observations to find the optimum parametrization for the 1D and 2D models. If the computational power is a major 468 limitation, then multiple 1D simulations are preferable to a 2D model. Moreover, "single_dir" is more efficient than the "multi_dir" 469 approach (Roelvink et al. 2018).

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

480
Complex morphodynamic models are suitable to assess and investigate the storm impact in coastal areas. However, these models 481 can be sensitive to a large number of free parameters and require calibration and validation. The sensitivity and performance of XBeach 482 19 have been mainly investigated in mild slope beaches where the saturated surf zone condition is matched, while model behavior in steep 483 beaches has received less attention. Here, numerical experiments were combined with topographical measurements collected for two large 484 storms (16-yr and 5-yr return period) to obtain reliable settings for better model performance in steep beaches. These experiments 485 demonstrated that: 1) the model was sensitive to "facua" (parameterized wave asymmetry sediment transport component), "bermslope" 486 (upslope swash zone transport term), and "wetslope" (critical avalanching slope) when simulating high energy storms. However, if 487 "bermslope" was activated, the effect of "wetslope" was reduced in the beach face but still relevant in the dune face; 2) the model 488 calibration in an exposed beach (Praia de Faro) for a high-energy storm showed that when "bermslope" increased, "facua" must be 489 reduced, and vice versa, to properly simulate the erosion. Moreover, "bermslope" reduced the model results variability alongshore, 490 minimizing the effect of upper beach face cusps on the final model erosion; 3) with similar settings, 1D and 2D models simulated similar 491 erosion in the beach face, but the erosion in the dune increased in the 1D simulations; 4) after calibrating with one storm, the 2D model 492 ability to simulate erosion during two storms was classified as excellent. The value of "bermslope" can be related to the beach face slope 493 and it contributed to the utilization of values of "facua" and "wetslope" close to the default values; 5) in a groin-protected site, Quarteira, 494 the 2D model was also successfully validated, although required different settings when compared to the exposed beach. Also, the 495 predicted erosion in this site was not especially sensitive to initial beach topography. These findings demonstrate that these parameters 496 produce adequate results for both hydrodynamic conditions and coastal settings. This work provides new insights on how to improve the 497 modeling of coastal erosion processes in steep beaches and supports the implementation of morphodynamic models at exposed and groin-498 protected beaches.

500
All data, models, or codes that support the findings of this study are available from the corresponding author upon reasonable request.    is equal to 0.20. In the middle plots, "bermslope" was deactivated and in the lower plots, the "bermslope" was 0.10. 626 Fig. 6. Weak-berm profile sensitivity results displayed as erosion percentage against the value of "facua". In the upper plots, the 627 "wetslope" is equal to 0.20. In the middle plots, "bermslope" was deactivated and in the lower plots, the "bermslope" was 0.10.   Fig. 8. Five cross-shore profiles extracted from the two-dimensional model. The simulations correspond to Table 5 -i and Table 5 -v for 631 storm Emma.