Sensitivity of Storm-Induced Hazards in a Highly Curvilinear Coastline to Changing Storm Directions. The Tordera Delta Case (NW Mediterranean)

Abstract

Extreme coastal storms, especially when incident in areas with densely urbanized coastlines, are one of the most damaging forms of natural disasters. The main hazards originating from coastal storms are inundation and erosion, and their magnitude and extent needs to be accurately assessed for effective management of coastal risk. The use of state-of-art morphodynamic process-based models is becoming standard, with most being applied to straight coastlines with gentle slopes. In this study, the XBeach model is used to assess the coastal response of a curvilinear sensitive deltaic coast with coarse sediment and steep slopes (intermediate-reflective conditions). The tested hypothesis is that changes in wave direction may cause large variations in the magnitude of storm-induced hazards. The model is tested against field data available for the Sant Esteve Storm (December 2008), obtaining an overall BSS (Brier Skill Score) score on the emerged morphological response of 0.68. Later, the 2008 event is used as baseline scenario to create synthetic events covering the range from NE to S. The obtained results show that storm-induced hazards along a highly curvilinear coast are very sensitive to changes in wave direction. Therefore, even under climate scenarios of relatively steady storminess, a potential shift in wave direction may significantly change hazard conditions and thus, need to be accounted for in robust damage risk assessments.

1. Introduction

The impact of extreme storms on the coast is one of the costliest forms of natural disasters (Kron [1]; Bertin et al. [2]). In heavily urbanized coastal areas, such as the Mediterranean (in general) and the Catalan coast (in particular), where properties, infrastructures and businesses are located close to the shoreline, this kind of event usually results in the damage or destruction of exposed assets (Jiménez et al. [3]). These effects are the integrated consequences of two main storm-induced coastal hazards: inundation and erosion. In this context, an accurate assessment of the magnitude, location and extension of these hazards is becoming an essential part of the risk management process (e.g., Ciavola et al. [4,5]; Van Dongeren et al. [6], Jimenez et al. [7]; Plomaritis et al. [8], Harley et al. [9]) and, in this sense, the use of process-oriented models to forecast storm-induced morphodynamic changes under given scenarios is now standard (e.g., Roelvink et al., [10]; McCall, [11]; Van Dongeren et al. [12] and references therein, Dissanayake et al. [13]). Most of the studies on testing state-of-art morphodynamic process-based models have addressed cases characterized by straight coastlines and gentle slopes (i.e., conditions close to the comfort zone of the models) (e.g., McCall, [11]; Harter and Figlus [14]). However, applications to estimate costal hazards in hihgly curvilinear environments (e.g., deltaic cuspate coasts) have seldom been tested (e.g., Roelvink et al., [15]; Valchev et al., [16]; and Dissanayake et al. [13]). Furthermore, the effect of testing models based on surf-beat (i.e., the infragravity wave band) on steep slopes and coarse sediment has been recently undertaken mainly in 1D applications (e.g., Vousdoukas et al. [17]; Elsayed and Oumeraci [18]) but rarely so in fully 2DH (2-dimensional, depth-averaged) simulations.

Within this context, the magnitude of storm-induced hazards on a highly curvilinear coast by using XBeach is assessed in the present study. The relevance and main aim of this work is twofold: first, from a general standpoint, to test the use of Xbeach on a highly curvilinear coast characterized by coarse sediment reflective beaches, and second, from the local standpoint, to analyze the sensitivity of an already identified hotspot, the Tordera Delta (NW Mediterranean) (Jiménez et al. [7]), to assess storm impacts for different storm direction scenarios. Thus, the largest recorded storm in the area is used as base case scenario. It occurred in December 2008 and had the typical incoming direction of current climate conditions where eastern (E) incoming storms dominate (e.g., Mendoza et al. [19]). Existing storminess projections under climate change scenarios for the Western Mediterranean do not predict any increase in wave height (e.g., Lionello et al. [20]; Conte and Lionello [21]), but some projections identify potential changes in wave direction (Cases-Prat and Sierra [22,23]). Due to this and to the great sensitivity of cuspate coastlines to wave direction resulting from their curvature (e.g., Slott et al. [24]; Johnson et al. [25]), the study aims to assess the potential effects of changing wave direction on extreme storm-induced hazards for the Tordera Delta. The hypothesis to be tested is that changes in wave direction may cause large variations in the magnitude of storm-induced hazards. Other studies have included the sensitivity to incoming storm direction in their assessments, such as those by Mortlock et al. [26] in Australia, or de Winter and Ruessink [27] in Holland.

The article is arranged as follows: the second section introduces the study site and the data used, describes the Sant Esteve 2008 event, which is used as the base case storm-scenario, and presents the methodological part, i.e., the used morphodynamic model and the comparative assessment framework descriptions; the third section presents obtained results; and finally, the discussion and concluding remarks are presented in the fourth section.

2. Materials and Methods

2.1. Study Area

The Tordera Delta is located in the NW Mediterranean Sea approximately 60 km northwards of Barcelona (Figure 1). It is a coarse sand delta (i.e., sediment in the range between 0.8 and 1.6 mm) with an aerial surface of approximately 4.2 km² at the end of a small river basin of approximately 879 km² (Vila and Serra [28]). The Tordera river is dry during most part of the year, with long dry summers and episodic discharges after heavy rainfall (Martin-Vide and Llasat [29]). Coastal storms and heavy rainfall events are usually uncorrelated at the area and, in fact, during the simulated storm, significant river outflow was present after the storm peak, with a phase delay of 1 day (Sanchez-Vidal et al. [30]). It has a cuspate configuration with two well differentiated parts at each side of the river mouth (Figure 1). The northern part is fronted by the S’Abanell beach, which is oriented towards the E. It has a steep nearshore bathymetry without any relevant geomorphological features in shallow waters The northern hinterland presents a higher topography, except for the zone closest to the river mouth outlet. The southern part is fronted by the Malgrat de Mar beach and is oriented towards the S, which serves as natural protection of a lower hinterland. Shallow water bathymetry is characterized by the presence of a longshore bar running parallel to the coast from the river mouth to the SW, which encloses a plateau at 4 m water depth.

The delta coastline has been eroding during the last several decades as the net result of the littoral drift and the decrease of the Tordera river sediment output, with maximum measured retreats of approximately 120 m Jiménez et al. [31]. The combination of a progressively narrowing beach protecting a low-lying hinterland and this being mainly occupied by campsites makes this area a hotspot for storm-induced hazards (Jiménez et al. [7]) with different consequences depending on storm characteristics and beach morphology at the time of impact (see e.g., Sanuy et al. [32]). As Jiménez et al. [31] noted, under former accretive-stable deltaic conditions only extreme storms were able to exceed the capacity of protection provided by wide beaches, but under present medium/long-term erosion conditions, smaller storms have become able to exceed the dissipation capacity of the narrower beaches, increasing the frequency of storm-induced problems (e.g., Jiménez et al. [3]).

Wave climate at the NW Mediterranean is characterized in the wave rose in Figure 1 (Global Ocean Waves, GOW dataset hindcast, Reguero et al. [33]). Storm events are defined as events with significant wave height (Hs) exceeding 2 m during a minimum of 6 h (Mendoza et al. [19]). The main incoming wave direction at the site is NE and E, with some events arriving from the S, especially during spring. Within this two main groups, some residual events, can also be found. Thus, in propagated conditions (nearshore area) storm conditions can be characterized with waves in the range NE-S being the two extremes of the range the most frequent situations. Nonetheless, some studies rise concerns at the Catalan coast about future climate-induced changes in storm direction, particularly a frequency transfer from the current main cluster (NE-E) towards the secondary one (S) (Casas-Prat et al. [22,23]). The area is micro-tidal, with storm surges having a limited role in storm-induced processes due to its relative weight when compared to the wave component. Notably, surges are uncorrelated with waves, are most frequently under 25 cm, with some extreme events showing maximum recorded surges around 50 cm (Mendoza and Jiménez [34]).

2.2. Data

The wave data used in this study to characterize the Sant Esteve storm and validate the models were measured by a directional wave buoy located off the Tordera Delta at approximately 70 m water depth (Figure 1) belonging to the XIOM (Xarxa d’Instruments Oceanogràfics i Meteorològics) network, which is no longer operative (Bolaños et al. [35]). The Tordera Buoy was a Datawell waverider, sampling during 20’ every hour, and thus providing sea states and statistics with an hourly time-step. It covers the period from 1984 to 2013.

Wave spectra from deepwaters and wind fields used to force the model chain were provided by Puertos del Estado, from the WAM (WAve Model, version 4, European Centre for Medium-Range Weather Forecasts, Reading, UK) and HIRLAM-AEMET (HIgh Resolution Local Area Modelling version 7, Agencia Estatal de METeorología, Madrid, Spain) models respectively. Both datasets have a temporal time-step of 1 h. Nearshore water levels were obtained from the same source, provided by the HAMSOM model (HAMburg Shelf Ocean Model, barotropic version, Center for Marine and Atmospheric Sciences. Institute of Oceanography, Hamburg, Germany) with the same temporal resolution.

Storm-induced topographic changes have been quantified by using LIDAR-derived topographies obtained before (16 October 2008) and after the storm impact (17 January 2009) by the Institut Cartogràfic i Geológic de Catalunya. The data were provided as high-resolution digital elevation models (DEMs) with a 1-m grid step, a maximum vertical precision of 2–3 cm and overall RSME of 6 cm (Ruiz et al. [36]).

The bathymetry of the study area has been obtained by combining an offshore grid with a spatial resolution of 0.28’ derived from the GEBCO (General Bathymetric Chart of Oceans) database [37], and a finer inner topography while nearshore bathymetry was built by combining the LIDAR-derived topography and multi-beam bathymetric data provided by the Ministry of Agriculture, Fish, Food, and Environment, covering the whole area from the +2 to the −50 m with a 5 × 5 m resolution. Multiple bathymetries were available from different years, including 2006 and 2010, and these information was merged to properly fit the 2008 LIDAR shoreline and link the emerged topography with the submerged bed.

2.3. The Sant Esteve 2008 Storm

The storm of reference used in this study was a V-class (extreme) event, according to the classification of storms in the NW Mediterranean of Mendoza et al. [19]. This storm, known as the Sant Esteve storm, occurred on 26 December 2008 in the northern part of the Catalan Sea. It was created by the presence of a low-pressure center located over the Balearic Sea, with a minimum pressure of 1012 hPa, along with a high pressure center over northern Europe (1047 hPa). This is one of the typical mechanisms of cyclogenesis in the Mediterranean (e.g., Trigo et al. [38]), and it is the most common situation originating extreme storms along the Catalan coast (Mendoza et al. [19]). Under these conditions, the action of very strong NE winds in the Gulf of Lyon (wind velocities up to 20 m s⁻¹ were recorded at the coast) generated a wave field with a clear spatial pattern, with Hs values and power content decreasing from north to south along the Catalan coast (Jiménez [39]; Sánchez-Vidal et al. [30]). Thus, according to the data recorded by the Palamós buoy (see location in Figure 1), the storm lasted 73 h (Hs > 2 m) reaching a Hs of 7.5 m at the peak of the storm and a maximum wave height (Hmax) of 14.4 m. The associated return period of this storm was approximately 125 years according to the data on extreme climate obtained by Puertos del Estado [40] for this buoy and in light of data previous to the storm. The storm progressively lost strength as it moved south and, off the study site, the values recorded by the Tordera buoy showed a storm with a duration of 66 h, reaching a Hs at the peak of the storm of 4.65 m and Hmax of 8.0 m. Mean wave direction during the storm was E, which corresponds to the dominant direction during extreme storms in the area (Mendoza et al. [19]). Figure 2 shows the recorded values of wave parameters during the storm by the Tordera buoy (see location in Figure 1).

The impact of the storm produced a significant coastal morphodynamic response in the form of erosion and overwash for most beaches along the northern part of the Catalan coast (Plana-Casado, [41]; Jiménez et al. [42]; Durán et al. [43]). Moreover, the magnitude of the storm was so considerable that many benthic ecosystems were also significantly affected (Sánchez-Vidal et al. [30]; Teixidó et al. [44]; Pagès et al. [45]).

Storm-induced topographic changes in the surroundings of the Tordera river mouth are shown in Figure 3. The observed response was different on both sides of the river, with the largest erosion taking place in the northern part, the s’Abanell beach. This beach is oriented to the East, nearly perpendicular to storm waves, and thus, the beach presented a generalized erosive behavior along its total extension (2.4 km). The volume of sediment eroded from the subaerial part of the beach was approximately 66,000 m³, with a beach-averaged erosion rate of approximately 30 m³/m and a maximum value of approximately 80 m³/m at its northernmost part (section SB-1 in Figure 3). Storm-induced wave overtopping occurred along the entire beach, and in its southernmost part, close to the river mouth, overwash deposits up to 6 m³/m were detected (section SB-3 in Figure 3). These volume changes resulted in a beach-averaged shoreline retreat of 11 m, with a maximum recession of approximately 30 m (Plana-Casado, [41]; Jiménez et al. [42]).

From the river mouth to the south, the coast experienced a different morphodynamic response. This section, the Malgrat de Mar beach, is oriented to the S, resulting in a large obliquity to E incoming waves during the storm. Just south of the river mouth a small post-storm accretion spot of approximately 7000 m³ was detected. This seems to be related to the alongshore deposition of material eroded from the northern part. South of this area, the coastline shows a nearly generalized moderate erosion together with significant overwash deposits in the subaerial part of the beach (Figure 3). The spatially averaged erosion of the emerged beach was approximately 10 m³/m (one third of that observed for the northern beach) whereas the averaged overwash accumulation was approximately 7.5 m³/m (Plana-Casado, [41]; Jiménez et al. [42]).

2.4. Models

Simulations have been done by using a model setup composed of the SWAN model (Simulating Waves Nearshore, version Cycle III v41.01, Delft University of Technology, Deltares, The Netherlands) (Booij et al. [46]; Ris et al. [47], TU Delft [48]), which propagates waves to the coast, and of the XBeach model (eXtreme Beach behavior, Kingsday version, Deltares, Delft, The Netherlands) (Roelvink et al. [10]), which is used to assess two storm-induced coastal hazards: inundation and erosion.

SWAN (Simulating Waves Nearshore) is a third generation wave model which is based on the wave action balance equation. It simulates short-crested wind generated waves by incorporating wave–water interactions, wind growth and dissipation processes such as whitecapping, bottom friction, and depth-induced breaking. For more detailed insight into these mechanisms controlling energy and wave propagation processes the reader is referred to the SWAN manual and to SWAN scientific technical documents at http://swanmodel.sourceforge.net/. The implemented model version uses the Komen et al. [49] formulation to calculate whitecapping. This version permits counteracting the previously reported under-predictions of significant wave height and wave period in areas characterized by fetch-limited conditions under the influence of transient winds (see e.g., Bolaños [50]; Pallares et al. [51]), which are typical conditions for the NW Mediterranean coast.

The model has been implemented using a nested grid configuration (Figure 1). A coarse grid with a total extension of approximately 80 × 70 km is used to transfer offshore wave conditions to the study area. The bathymetry grid has a spatial resolution of 0.28′ whereas the wind field grid has a spatial resolution of 5′. This coarse setup is fed with wave spectra at 15 positions distributed along the offshore boundaries of the grid (Figure 1). The inner fine grid covers a domain of approximately 20 km × 26 km with a spatial resolution of 0.06’. This grid has been generated to properly reproduce the existing sharp changes in the bathymetry between intermediate-shallow waters due to the large steepness of the lower shoreface. This can permit a better simulation of wave propagation in the region close to the XBeach coastal grid (external boundary at 20 m water depth). A limitation of the SWAN model is its inability to simulate storm surges. To characterize storm surges in the study area we use water level predictions obtained with the HAMSOM model implemented by Puertos del Estado (Ratsimandresy, et al. [52]) at three locations close to the study site. The storm surge contribution to the total water level in the study area is of low magnitude and significantly smaller than wave-induced runup during storm conditions (e.g., Mendoza and Jiménez [34]). The validated SWAN model has been used to propagate waves from deep to shallow waters during the entire storm duration. Since a simultaneous time series of measured wave data was available within the modelled domain, it was possible to obtain transfer coefficients for wave conditions (wave height and direction) for any point in the grid with respect to buoy location. With this information, the recorded wave conditions at the buoy have been transferred to selected grid points to be used as input data for Xbeach modelling.

Once the forcing conditions during the storm have been propagated from deep water to nearshore, the remaining task is to propagate these conditions to the coast and to assess the magnitude of storm-induced hazards, i.e., erosion, overwash and overtopping, which is done with the XBeach model. XBeach is an open-source 2D depth averaged model which solves wave propagation, flow, sediment transport and bed level changes for varying wave and flow boundary conditions (Roelvink et al. [10]). It solves the time-dependent short wave action balance on the scale of wave groups, which allows for the reproduction of directionally spread infragravity motions (so called surf-beat) along with time-varying currents. The frequency domain is represented by a single representative peak frequency, assuming a narrow banded incident spectrum. Shallow water momentum and mass balance equations are solved to compute surface elevation and flow. Additionally, to solve the contribution of short waves to mass fluxes and return flows, XBeach uses the Generalized Lagrangian Mean formulation (Roelvink et al. [10]).

Sediment transport rates are calculated from the spatial variations in depth-averaged concentration, which are calculated from advection-diffusion equations with a source-sink term based on an equilibrium sediment concentration. The equilibrium concentration takes into account both the contribution of the suspended and bed loads by means of the Soulsby-Van Rijn formulation (Soulsby [53]) with a limitation of the maximum stirring velocity based on the Shields number at the start of the sheet flow (McCall et al. [11]). For deeper insight into the XBeach description, setup and equations, see Roelvink et al. [10].

The model has been implemented by using a curvilinear grid with variable cell size in both alongshore and cross-shore directions (see Figure 4). The extension of the mesh is approximately 1.5 km in the cross-shore direction, with cell size ranging from 5~6 m at the offshore boundary (20 m depth) to 0.7–0.8 m at the swash zone. In the alongshore direction the model has an extension of 4.5 km and the cell size ranges from 25 m at the lateral boundaries down to 2–3 m around the river mouth. The grid was obtained by means of the Delft3D-RGFGRID module, imposing consecutive maximum cell-size changes ~5–10% to ensure proper smoothing, while maintaining valid orthogonalisation. The final result has 669 alongshore by 568 cross-shore nodes. The model wave boundary conditions consist of wave characteristics specified at four different locations along the offshore boundary (see Figure 1), with a time-step of 1 h, which is the same resolution of the original data used to force SWAN. The use of four different locations aims to capture the difference in wave conditions on both sides of the river due to the different coastline orientation. Water level variations during storms to be simulated are directly introduced in the model from HAMSOM simulations, also with a time-step of 1 h. XBeach model computations (i.e., hydrodynamics and morphodynamcis) are performed with a temporal resolution of 1 s (dtbc = 1).

The morphology of the study area, characterized by coarse sand and steep reflective beaches, makes XBeach modelling a demanding exercise in terms of predicting beach morphodynamic response during storms (see e.g., Vousdoukas et al. [17,54]). Notably, coarse sediment environments are characterized by a lower frequency of the avalanching processes, a greater importance of the bed load over the suspended transport, and higher importance of mechanisms such wave asymmetry, or water infiltrations and groundwater effects. All these are by default configured in XBeach to work for fine sand environments, and must be revised and modified for its application at the study site, modelled with a D50 of 1.3 mm and a D90 of 1.9 mm.

To properly reproduce morphodynamic changes, both the surf-beat and the non-hydrostatic modes of the XBeach model were initially tested. The surf-beat model was observed to under-predict overwash when using typical XBeach-grid resolution for straight and mild-slope coasts. This under-prediction is also caused by the lower contribution of infragravity waves to the total run-up in steep profiles (Wright and Short [55]). However, the surf-beat model accurately reproduced the alongshore patterns of erosion in the entire domain. In contrast, the non-hydrostatic model was observed to have better performance in reproducing wave-by-wave run-up but lower accuracy in reproducing the alongshore morphodynamic patterns of erosion and deposition, with a quite higher computational cost (since it requires denser grids). An additional difficulty was using XBeach with a curvilinear grid to properly reproduce alongshore processes driven by the change of coast orientation at both sites of the river mouth. The final adopted approach was to use the surf-beat model with a higher resolution than the typically used in straight and gentle-slope coasts, which considerably improved the overwash prediction. The average cross-shore resolution of 5.2 m (typically 20~25 m) at the offshore boundary, going down to 0.7 at the swash zone, ensured a better reproduction of wave propagation. This setting aimed to properly capture the abrupt changes in the bathymetry from the 20 m to the 5 m depth by using a larger number of cells. In addition, the average alongshore resolution around the river mouth, where the largest alongshore gradients are present, was increased up to 2.3 m (typically 5~10 m) (Figure 4).

2.5. Scenario Testing

Since the recorded Sant Esteve 2008 event had the characteristics of a V-class extreme storm in the area, with the largest recorded Hs and with a typical E direction, it was used as base case scenario (C0). From this, six additional scenarios were defined to test the effect of different incoming wave directions. The procedure consisted of simulating the exact same wave time series as the Sant Esteve 2008 storm (keeping the same Hs, Tp, and directional spreading) and only changing the mean wave direction. This approach permits maintaining the same storm wave intensity as the reference case at the offshore boundary but under a different direction. The tested conditions were two scenarios where wave direction was shifted 20° and 40° counter-clockwise to the North (C20− and C40−) and four scenarios shifting 20°, 40°, 60°, and 80° clockwise to the South (C20+, C40+, C60+, and C80+). Thus, seven different scenarios, including the baseline condition, were simulated and compared to assess the differences in storm-induced hazards under incoming directions ranging from ~60° N (C40−) to ~180° N (C80+). The scenarios have been chosen to cover the tipical range of incoming conditions at the −20 m depth with directional spans of 20°.

The magnitude of storm–induced hazards was quantified in different control sectors along the study area to capture main factors potentially affecting the beach response to different incoming directions. Thus, the morphodynamic response was analyzed in five sectors: two 250 m long sectors northward of the river mouth (SN1, SN1), and three sectors southward of the river mouth (Figure 4). These S sectors are as follows: (i) SS1, a 200-m-long stretch between the river mouth and an existing rigid structure at the shoreline; (ii) SS2, a 200-m-long sector, southward of the mentioned existing structure; and (iii) SS3, a 500-m-long stretch at the southernmost end. At each sector, three variables are used to characterize morphodynamic changes: erosion volume (m³/m), overwash volume (m³/m) and profile retreat (m). All of the variables are calculated within sectors from XBeach gridded output, from which sector-averaged values and standard deviations are derived. Erosion volumes are computed in the inner part of the beach, from the subaerial part down to the −2 m level, which roughly determines the water depth where main inner profile changes occur. Overwash volumes are computed as deposited sediment volumes in those parts of the subaerial beach where vertical growth is detected. The profile retreat is measured at three elevations at the beachface (1, 1.5, and 1.75 m above mean water level).

To characterize inundation, just two sectors were selected, one for the region at the north of the river mouth and another to the south (Figure 4). The selected variable to characterize inundation was the inundation surface (Ha) over different thresholds of inundation depth: 0.05, 0.25, 0.5, and 1 m.

3. Results

3.1. Base Case Scenario (C0). The Sant Esteve Storm

This base case scenario corresponds to the model validation using the recorded Sant Esteve 2008 Storm event. The SWAN model was validated with wave conditions recorded by the Tordera wave buoy and by comparing the measured and modelled wave conditions (Hs, θ), shown in Figure 2. As can be seen, the model reproduces well wave parameters during the storm, especially during the peak of the event, when simulated variables almost coincide with recorded ones. The obtained root mean squared error (RMSE) of the model during the entire duration of the storm is 0.53 m in Hs and 12.5 degrees in θ, with the largest contributions to the error taking place during the relaxation phase of the storm, whereas error values at the storm peak are significantly lower (Figure 2). Since the largest storm-induced morphodynamic changes occurred during the storm peak, we can accept that the use of the SWAN model to simulate wave propagation in the study area is acceptable for the purposes of this research.

To calibrate XBeach, the model’s results were compared with LIDAR measurements of the emerged beach. The calibration of the model was performed by adopting a double approach: (i) by optimizing the Brier Skill Score (BSS), which quantifies model performance by comparing model output to the real post-storm LIDAR measurements of the emerged profile; (ii) by performing a qualitative assessment of the modelled features, such as alongshore and cross-shore patterns of bed level changes, magnitude, and location of the overwash deposits and validation of the inundation reach according to the available qualitative information on the event. The used BSS to characterize model predictive skill takes into account the measurement error (∆Z_e) as in Harley and Ciavola [56], and thus, the BSS score is given by:

$$\mathrm{BSS}=1-({\displaystyle \sum }{\left(\left|{\mathrm{z}}_{\mathrm{mf}}-{\mathrm{z}}_{\mathrm{mod}}\right|-\Delta {\mathrm{Z}}_{\mathrm{e}}\right)}^2/({\displaystyle \sum }{\left({\mathrm{z}}_{\mathrm{mf}}-{\mathrm{z}}_{\mathrm{mi}}\right)}^2$$

(1)

where z_mf is the final LIDAR measured bed level, z_mod the final modelled bed level, and z_mi the initial bed level. Here, ΔZ_e is considered as the LIDAR measurement error (i.e., RSME of the overall LIDAR product), which is 0.06 m. According to Sutherland et al. [57], the classification of models’ performance based on the BSS score can be considered to be very good for values over 0.4 and excellent for values over 0.5–0.6. Due to the morphology of the study area, which induces a differentiated morphodynamic response along the coast (Figure 5 and Figure 6), three BSS were calculated: (i) a global BSS, which is calculated for the entire area; (ii) a local BSS at the north of the river mouth; and (iii) a local BSS southward of the river mouth. Since the BSS assessment can only be performed for the emerged profile (where pre- and post-storm topographic data exist), a qualitative assessment of the modelled submerged profile was also performed. To this end, we analyzed the final shape of the modelled submerged profile taking into account the expected typical morphodynamic response under storm conditions at the site.

Although different combinations of model parameters resulted in BSS scores over 0.4 for the emerged profile, the qualitative assessment highlighted that some of them produced excessive deposition volumes in the submerged part. The final setup parameters, which resulted in an overall BSS of 0.68 (measured at the northern beach and first 600 m of the southern beach) and a meaningful qualitative simulation of the predicted submerged profile and alongshore bed level change patterns, are shown in

Table 1. Calibration parameters. Range of tested values and final parameter setup.

Parameter	Tested Values	Description	Final Set-up 1
gamma	0.55–0.7	Breaker parameter in Baldock or Roelvink formulation (default = 0.55)	0.7
delta	0–0.5	Fraction of wave height to add to water depth in wave breaking formulations (default = 0)	0.5
facAs	0.2–0.7	Calibration factor time averaged flows due to wave asymmetry (default = 0.1)	0.6
facSk	0.2–0.7	Calibration factor time averaged flows due to wave skewness (default = 0.1)	0.6
wetslp	0.3–0.8	Critical avalanching slope under water (dz/dx and dz/dy) (default = 0.3)	0.7
gwflow	0 and 1	Turn on groundwater flow (default = 0)	1
sedcal	0.1–1	Sediment transport calibration coefficient per grain type (default = 1)	0.1

¹ based on best BSS score and qualitative assessment.

along with parameter description and tested ranges.

The implications of the selected values are as follows. First, the increase in wave-attack on the coast improves the behavior of the model in terms of reproducing observed overwash deposits (gamma and delta). Second, the contribution of avalanching to bed level changes was limited, taking into account the steepness of the site (wetslp). Third, the non-linearity effect on sediment transport for steep profiles (facAs and facSk) may be taken into account as reported in Elsayed and Oumeraci [18]. Fourth, sediment particle mobilization may be limited (sedcal). With respect to this, other authors have already reported an excess of modelled sediment suspension because the shear stress values required to initiate particle motion are higher than those predicted by using the Shields curve, as noted in Elsayed and Oumeraci [18] or McCall [11]. Fifth and finally, the groundwater module was turned on (gwflow) because the role of infiltration is more significant in coarse sediment environments, such as the Tordera Delta, where grains are close to gravel-size. The value of the permeability factor has been estimated using the Kozeny–Carman formula as described in Carrier [58]. For the grain size in the study area, the value of the permeability factor was estimated to be 0.0058 m/s. The results of the storm-induced morphological changes simulated with the final adjusted set of model coefficients are shown in Figure 5 and Figure 6. As it can be seen, the measured changes in beach elevation at both sites of the river mouth are well reproduced by the model, with a BSS of 0.68 when the entire area is considered. The model also properly reproduces the differentiated response at both sites of the river mouth, mimicking the effects of change in coastal orientation with respect to the storm wave direction and the differences in coastal morphology.

Figure 7 shows the modelled depth-averaged wave-induced circulation during the peak of the storm, where a different circulation pattern is observed at both sides of the river mouth. Mean XBeach output (average conditions at 1 h output time-step) is used to average the current velocity and sediment transport during 4 h around the maximum peak of the event. At the northern part, where the beach is fully exposed to storm waves (i.e., a relatively small obliquity during the storm) and without any submerged morphological features, wave-induced circulation shows a typical quasi-uniform longshore current structure along the beach until the river mouth at the southern end. At the southern part, the coast is partially sheltered from storm waves (i.e., a large wave obliquity during the storm) and there is an alongshore bar running parallel to the shoreline with varying crest levels. This bar delimits a shallow shelf of approximately 4 m water depth to the shoreline. In this area, the induced circulation pattern during the peak of the storm shows a longshore current directed towards the south close to the shoreline and a local inversion of the current towards the N over the shallow shelf (Figure 7).

With respect to morphodynamic changes (Figure 5 and Figure 6), the model predicts, for the northern part, a nearly continuous erosion along the beach without significant alongshore gradients in sediment transport, and with significant sediment mobility down to about −5 m depth. The model properly reproduced the observed increase in erosion magnitude from B2 (P4) to B1 (P3), as well as most of the overwash deposits, which increase towards the south as the height of the berm decreases. In summary, the model reproduced well the observed variability in beach erosion and overwash along this northern site with a local BSS of 0.75. Southwards of the river mouth (Malgrat Beach), induced sediment transport and beach erosion present significant alongshore variations, with the largest erosion taking place just southwards of an existing structure located at the northern part of the sector, A2 (P2), which should act as a local boundary condition. This local effect can be seen in Figure 7 where a gradient in the longshore sediment transport, starting at the rigid structure north of A2, is detected. Longshore transport rates are larger than in the north, and mainly concentrated down to −3 m water depth. The model reproduced observed beach topographic changes with large erosion and overwash deposits due to a relatively low beach berm (A2, P2). At the southernmost part, erosion is only taking place at the upper level of the profile, whereas part of the material is deposited around zero level (A1, P1). The obtained local BSS for this sector was 0.60.

3.2. The Effects of Wave Direction on Storm-Induced Hazards

As was previously mentioned, once the morphodynamic model was calibrated, it was used to analyze the sensitivity of the area to changes in wave direction during storm impacts. Simulated morphodynamic changes and inundation for cases C20− (~80° N) to C40+ (~140° N) at northern and southern parts of the study area are shown in Figure 8, whereas the variation of integrated control variables for each sector can be seen in Figure 9 and Figure 10 for inundation and morphodynamic changes respectively.

With respect to inundation, the beach northward of the river mouth (Figure 9a) experiences an increase of the predicted inundation surface from C40− to C20+, when the surface reaches its maximum extension (5.47 Ha with 1 Ha over 0.5 m depth). Wave direction in this scenario corresponds to a nearly normal wave attack on the local coastline orientation. As incoming wave direction continues shifting towards the south, the predicted inundation surface progressively decreases, reaching values for the C80+ case similar to those observed for the C40− case. At the southern coast, the obtained pattern is significantly different (Figure 9b). Thus, for wave direction scenarios dominated by NE components (C40−, C20−) no significant inundation is observed, and this is consistent with a large sheltering from highly oblique wave incidence. For wave directions shifting from the base scenario to the south (C0 to C80+), the predicted inundated surface increases, reaching a maximum value for scenarios C60+ and C80+ of approximately 74 Ha (for an inundation depth above 0.05 m) and 33 Ha (for an inundation depth of 0.5 m), the C80+ case. Taking into account the total inundation of the study area, we can state that the magnitude of the inundation hazard significantly increases as wave direction shifts to the south. The southern coast is the most affected in terms of the magnitude of expected changes due to the local low-lying topography.

Regarding morphological changes, the northern beach (S’Abanell) shows a relatively low sensitivity to wave direction for scenarios C40− and C20+ when erosion volumes and profile retreats are considered (Figure 10). For this range of wave directions, the local wave-induced circulation is characterized by a southward directed longshore current along the beach (see Figure 7 and Figure 11) which turns north when the wave direction is C20+ (Figure 11). As the incoming wave direction turns southwards, a stronger north-directed alongshore current is induced with velocity increasing along the beach (Figure 11) which results in an increasing erosion and profile retreat. Figure 11 also shows hatching areas in the velocity field which originate in those areas where wave incidence is orthogonal to the nearshore bathymetry, characterized by a heterogeneous bar in front of the coast. This beach erosion increase is particularly observed at the southernmost end of the section (SN1), just northwards of the river mouth, where an existing revetment acts as a boundary condition (barrier) for northwards directed transport. This overall modelled behavior is consistent with local field observations, where the northern part of the beach (out of the domain) experiences a significant sediment deposition under the impact of southern storms, bringing sediment from the south.

Overwash deposits along the northern section present a variation pattern with wave direction consistent with modelled inundation. The largest overwash verifies in the southern end of this sector (SN1) due to its lower beach berm. Thus, maximum overwash deposition verifies at SN1 under C20+ and C40+, which were the scenarios producing the largest inundations. As wave obliquity increases (scenarios C60+ and C80+), overwash significantly decreases, which is also in agreement with the observed inundation decrease under these conditions.

Southwards of the river mouth (Malgrat de Mar), the analysed sectors show a differentiated response due to the local effect of the existing building in the coastline which acts as a boundary condition for sediment transport. The northernmost end, SS1, just southward of the river mouth, experiences accretion (i.e., positive profile displacement) under C40− and C20− storms, which would correspond to the deposition of sediments brought by the southward directed longshore transport from the northern beach. As wave direction shifts to the south (from C20+ to C80+) this area becomes slightly erosive. This morphodynamic behavior is consistent with the observed switch in alongshore current direction for scenarios C20+ and C40+ (Figure 11). The central sector SS2, located just southwards of the mentioned obstacle, shows a different morphological response with storm direction to that observed in SS1. Thus, the largest erosion rate verifies under C0 due to combined effect of the orientation of the coast, a gap in the submerged bar and the presence of the hard structure at the northern end of the sector. This can be seen in Figure 7, where a significantly large gradient in sediment transport rates under C0 is observed. As wave direction shifts towards south, erosion rates significantly decrease, reaching a nearly null integrated value under C40+ to C80+ scenarios. The largest erosion under scenarios C20− to C0 in this sector, in comparison with the rest of the southern coast, seems to reflect the effect of the revetment on enhancing local downcoast erosion, similar to the well-known flanking effects in seawalls (e.g., Kraus and McDougal [59]). The southernmost sector, SS3, is far from the local effect of the revetment, so much so that the aforementioned downcoast erosion enhancement is not detected and, under highly oblique storms, erosion rates are significantly lower. Total volumetric changes in this sector present a relative low sensitivity to changes in wave direction for southern wave scenarios (C20+ to C80+).

Regarding overwash, the southernmost section SS3 presents a wave-direction influence similar to the influence observed for inundation, with overwash rates growing as wave direction turns to the South. The sectors closer to the river mouth SS2 and SS1 present a similar response with a maximum for C20+ and C40+ scenarios, when the wave directions are close to local orthogonality (Figure 8), and they are significantly larger for SS2 (~12 m³/m) than SS1 (~2 m³/m).

4. Discussion and Conclusions

In this study, the potential effects of changing wave direction for the storm-induced hazards on a highly curvilinear coarse sandy coastline have been assessed. This sensitivity test has been selected because although storminess projections under climate change scenarios for the Western Mediterranean do not predict any increase in wave height (e.g., Lionello et al. [20]; Conte and Lionello [21]), some existing projections identify potential changes in wave direction (Casas-Prat and Sierra [22,23]). These changes in wave direction may have significant implications for coastal sediment transport and coastal stability, as has been confirmed for the interannual changes influenced by El Niño (e.g., Barnard et al. [60]). Moreover, regarding cuspate coastlines such as the study area, their greater sensitivity, due to their curvature, results in even more significant implications (e.g., Slott et al. [24]; Johnson et al. [25]).

The tested hypothesis is that changes in wave direction may cause large variations in the magnitude of storm-induced hazards. This effect has also been addressed in other studies such as those by Mortlock et al. [26] and de Winter and Ruessink [27], which specifically analyzed the effects of changes in wave direction on the storm-induced hazards in the SE Australian and Holland coasts respectively. To this end and to isolate the influence of wave direction, we used a recorded long-return period storm as a base case scenario and we built test scenarios just by changing wave direction while maintaining the other wave parameters as recorded during the base storm (wave height and period).

In any case, tested conditions have not been designed to be used as climate change induced projections, as that may require the proper forecasting of regional wave conditions under given climate scenarios (e.g., Casas Prat and Sierra [23]). These have to be considered from the perspective of coastal risk management, in which a set of possible conditions are analyzed to characterize coastal vulnerability and resilience to inform risk management under uncertainty (see e.g., Hinkel et al. [61] for an application of this perspective to analyze sea level rise). In the study area, the current storm wave conditions depend on direction, with largest wave height and power being associated with NE-E waves, whereas S storms are less frequent and present a smaller associated power (e.g., Sánchez-Arcilla et al. [62]; Mendoza et al. [19]). To assess the potential variability on storm-induced hazards, tested scenarios were built by just changing wave direction while the remaining recorded parameters (representative of a worst case scenario, according to recorded conditions) were maintained.

This analysis has been performed by using the SWAN and XBeach models to simulate storm-induced hazards. Both models were calibrated by using data recorded during the impact of an extreme storm recorded in December 2008, which is used as the base case scenario. Although it is desirable to use more than one event to properly calibrate/validate the models (e.g., Ranashinge [63]), data availability during storm conditions was restricted to this event. However, on the positive side, it has to be considered that this storm was the largest event recorded in the area and representative of extreme storms with a very long return period (Mendoza et al. [19]) under current climate conditions. The SWAN model was very successful in simulating wave conditions during the development phase of the storm up to the pass of the peak of the storm, with the larger differences between measured and simulated waves being detected during the relaxation phase of the storm, when most of the induced changes had already occurred. The default parametrization of the XBeach model had to be adapted for application at the side to represent the effects the coarse-sand environment. Sediment transport was limited by using the sedcal parameter, avalanching was limited by increasing the critical slope, wave asymmetry was increased as suggested in literature for steep slopes (Elsayed and Oumeraci [18]) and groundwater effects were included. Gamma and delta wave breaking parameters were also tuned (Table 1). Calibrated parameters setup for XBeach in the study area led to a BSS score of 0.68 in spite of the out-of-comfort tested conditions (i.e., highly curvilinear coast, steep beach, coarse sediment). Although the predictive skill was very good for the northern and southern beaches, the model performance was better in the northern domain (BSS = 0.75) than in the southern one (BSS = 0.60), since this last area presented a significantly larger obliquity to wave direction during the storm, and a more complex bathymetry.

The obtained results show a very high sensitivity of storm-induced processes, i.e., inundation and erosion, to changes in storm wave direction. With respect to inundation, expected changes in hazard magnitude are very significant, especially in the southern part of the study area, since its morphology is characterized by a lower berm, and its low-lying unprotected hinterland makes this area sensitive to storm flooding (Jiménez et al. [7]). Thus, as storm waves turn from the base case (C0) to the south, the inundated surface along this southern beach dramatically increases due to its direct exposure to the south. On the contrary, a potential shift of wave direction to the N will have a positive impact on inundation in this area, since it will be more sheltered from wave action. At the northern beach, the largest increase in inundation hazard verifies under C20+ and C40+ scenarios when waves face nearly orthogonally to the coastline, although due to local morphology, the affected surface is much lower than in the southern beach. The hinterland of the study area is mostly occupied by agriculture land and, in the outer fringe just behind the shoreline, by campsites. In this sense, to transfer the potential change in hazard magnitude to changes in damage risk, it should be important to consider not only the change in direction but also its seasonality. Thus, risk may vary dramatically between the summer season (when the campsite facilities are used by visitors) and the rest of the year when only installations will be affected (e.g., Merz et al. [64]). An analysis of the risk associated with storm-induced inundation for different storm conditions can be seen in Sanuy et al. [32].

Similarly, storm-induced morphodynamic changes are more sensitive to directional changes on the southern beach, where the magnitude of the changes is larger. The beaches at the south of the river mouth present a larger spatial variability than those in the north due to the presence of a local boundary condition in form of a revetment at the shoreline. This revetment, which modifies local longshore transport, significantly enhances downcoast erosion under storm conditions. This induces a southwards directed longshore sediment transport while simultaneously promoting the accumulation of upcoast sediment. This contrasting behavior is particularly observed in the base case scenario which seems to represent the optimum conditions for longshore sediment transport in the area, thus inducing the largest changes in the surroundings of the structure.

In a particular case, under C40− and C20− scenarios, when wave direction turns north, the beach sector just south of the river experiences an important sediment accumulation due to the apparently efficient transfer of sediment from the northern beach across the river mouth and the partial barrier effect of the existing revetment.

The magnitude of the erosional response along the two control sectors in the northern beach is similar, although a higher variability is detected in the area closest to the river mouth. In general, there is a slight increase in erosion rates as wave direction turns south. This variation should be indicative of the role of longshore sediment fluxes during storm conditions. Thus, as the controlled northern area is just besides the river mouth, where there is another structure acting as a boundary condition, the increase in longshore sediment transport as waves turn S (scenarios from C20+ to C80+) will increase sediment losses, which will be transported further to the north. This behavior is currently observed in the northernmost part of this beach (out of the control zone in Figure 10) which experiences sediment accumulation under the impact of southern storms.

As expected, changes in the magnitude of overwash deposits follow observed changes in inundation, i.e., they increase as wave direction turns to the south, with maximum values around C20+ and C40+. The exception to this is the predicted changes in the southernmost sector, which present the largest overwash for C80+ conditions. The spatial variability in the northern beach is significantly lower than in the south, with small variations in magnitude across the tested range. Moreover, and reflecting the observed differences in inundation, the magnitude of overwash deposits is much higher in the southern sector.

Finally, and as a concluding remark, this analysis has shown that storm-induced hazards along a highly curvilinear coast are extremely sensitive to changes in wave direction. This means that even under a climate scenario of relatively steady storminess (wave power and frequency), a potential shift in wave direction may significantly change hazard conditions and, in consequence, need to be accounted for in robust damage risk assessments. To this end, an analysis such as the one presented here also permits an assessment of how coastal geomorphology modulates induced changes. In the study area, the low-lying nature of the southern beach and its orientation with respect to the current dominant storm direction make this area much more sensitive to directional changes. This is especially relevant from the coastal management standpoint because this area has been already identified as a hotspot for storm impacts under current conditions. The use of detailed process-based models has permitted the identification and quantification of the drastic increase in sensitivity when anthropogenic perturbations are present along the coast. These perturbations act as boundary conditions modifying local hydrodynamics and associated transport. For the case study analyzed here, the obtained results clearly identify the hazardous potential of the existing revetment in the southern beach, which has also been identified under current conditions, suggesting that its removal will soften the estimated morphodynamic response.