Hydro-morphodynamic modelling in Mediterranean storms – errors and uncertainties under sharp gradients

This paper deals with the limits in hydrodynamic and morphodynamic predictions for semi-enclosed coastal domains subject to sharp gradients (in bathymetry, topography, sediment transport and coastal damages). It starts with an overview of wave prediction limits (based on satellite images and buoy records) in a restricted domain, namely the Mediterranean basin, followed by an in-depth analysis of the Catalan coast, one of its land boundaries. The morphodynamic modelling for such regions is next discussed, based on the impact of a characteristic storm. The driving wave and surge conditions produce a morphodynamic response that is validated against the preand post-storm emerged beach state, recovered from two lidar images. The quality of the fit is discussed in terms of the physical processes and the suitability of the employed modelling suite. From here an assessment of errors and uncertainties is presented, with the aim of establishing the prediction limits for flooding and erosion analyses, key elements for coastal engineering decisions.

The general circulation pattern in this area follows the Northern or Ligurian current flowing from Southern France towards the central Spanish Mediterranean coast. It is characterized by average speeds around 10 cm seg −1 , although presenting winter intensifications up to 50 cm seg −1 . Some pulsating events and mesoscale eddies may lead to higher current velocities that in all cases are seldom strong enough near the 5 coast to transport sediment particles that are, from the circulation perspective, stable at the coastal zone. Mean water levels in the area are dominated by storm surges which may reach up to 50 cm (Conte and Lionello, 2013). The astronomical tidal range, between 10 and 30 cm is quantitatively less important, making the area a micro-tidal environment.

10
This average Southwest flowing current is super imposed on near inertial motions, linked to small sized meteo events which dominate the fluctuations over shelf regions (Rippeth et al., 2002). The flow is mostly topographically steered, including the effect of the shoreline boundary (Grifoll et al., 2012). Near the coast the mean flow is Southwest wards also, but with reversals throughout the year, particularly under high energetic 15 events mainly from the East and Northeast sectors.
In this type of coastal domains, the sharp gradients (due to the topographic and bathymetric features) in the spatial patterns of wind, wave and circulation fields impose a tough challenge for numerical simulations. It requires using a mesh discretization fine enough to be able to solve these spatial variations, without introducing numerical errors 20 or diffusivity that degrades the quality of the calculated results.
To prove this point, the paper will start by presenting wave simulations, the nesting strategy and how the computational results have been calibrated with available observations. This allows establishing some performance limits for the hydrodynamic variables (Sect. 2). The paper continues with a description of a significant storm (in 25 hydrodynamic and morphodynamic terms) that has recently affected the Catalan coast and for which there are some observations available (Sect. 3). We shall present hydrodynamic results from deep water down to the surf zone and assess the resulting errors, linking them to the controlling physical processes. It follows an analysis of the morphodynamic simulations performed, presenting the morphodynamic modelling suite and how it can provide erosion and flooding results based on the driving hydrodynamics for the selected storm (Sect. 4). This is the basis to calculate the errors between observations and simulations for hydrodynamic variables and the uncertainties we find when obtaining morphodynamic predictions (Sect. 5).

5
After this there is a discussion session on the limits of hydro-morphodynamic predictions and the implications for engineering design and management decisions (Sect. 6), followed by some conclusions (Sect. 7).

Hydrodynamic modelling. Waves, currents and mean sea level
The wave code selected is the third-generation model SWAN (Simulating Waves 10 Nearshore), which is based on a finite differences discretization of the wave action balance equation (Booij et al., 1999). It incorporates wind wave growth, quadruplet and triad wave-wave interactions plus depth-induced breaking. In this paper, the selected physics consist of exponential wind wave growth computed by means of the Janssen (1991) formulation. Non-linear wave interactions are parametrized through 15 the DIA scheme (Hasselmann et al., 1985) and the LTA method (Eldeberky, 1996) at deep and shallow waters, respectively. Regarding sink terms, whitecapping dissipation corresponds to the Komen et al. (1984) formula and bottom friction dissipation is derived from the JONSWAP results (Hasselmann et al., 1973). Finally, depth-induced breaking is triggered through the Battjes and Janssen (1978) bore-based model. 20 Calculations are performed through a set of three nested domains that cover the Mediterranean Sea. The first one is the Mediterranean sea grid which has a resolution of 10 km × 10 km (0.09 The main aim of this nesting strategy is to solve the spatial variability of the region and to improve wave estimations near the coast, since the morphodynamic results will 5 depend directly on the quality of these nearshore values, which happen to show a larger hydrodynamic error than that found for the more offshore (open) domain. Because of that we have tested higher resolution wind models going down to a mesh size of 4 km and a wind input every one hour (Bolaños-Sanchez et al., 2007). In all cases the differences between the input and dissipation terms were not relevant for such limited coastal domains. Moreover, the fact that SWAN uses a semi implicit scheme has also proven to be advantageous in terms of selecting time and space steps.
Model boundary conditions are provided by the wind, mean sea level and current fields. Input wind fields are obtained from IFREMER blended wind fields (Bentamy et al., 2007)  frequency. Data comprises operational ECMWF wind fields blended with remote sensing observations (QuikSCAT scatterometer and SSM/I radiometers) through optimal interpolation. These wind estimates were compared with moored buoy data, both at deep and shallow waters and the results exhibit better agreement in open sea than in nearshore areas. The conclusion is that near the coast the observations were up to 20 2 m s −1 stronger and having a weaker onshore component. The input included daily averaged sea level and velocity fields with a horizontal grid resolution of 1/16 • (c.a. 6-7 km), as provided by MyOcean (Tonani et al., 2009). Simulations span from 20 December to 31 December 2008 (i.e. they are initiated six days prior to the storm peak on 26 December). This start point guarantees enough 25 time margin to avoid warm-up phenomena, excluding it as a possible uncertainty factor during the extreme event simulation. The model outputs are compared with different data sources ranging from remote (satellite) images (deep water) to buoy network data (shallow water). Usually, a storm event is characterized by the storm duration and wave Introduction height, peak period and direction. Fortunately, wave height records are available both at regional and local scale, becoming the "connecting" physical magnitude along the different stages of this study. The models' suitability for the study domain and some of the associated errors, when compared with observations, have been also analysed in Grifoll et al. (2013).

5
The first step in the error estimation is based on the relation between satellite preprocessed wave height data (Queffeulou, 2004) and the numerical output at deep waters. Up to 14 000 data are within the considered regional area, derived from the combined tracks of four satellites: ERS-2 (21.43 % of the total), ENVISAT (21.43 %),  and . The second step deals 10 with shallow water comparisons, based on the XIOM meteo-oceanographic network observations. In this case, hourly wave height, period and direction are recorded.  ). In the studied period, going from 22 December to 31 December 2008, there were intervals of average and extreme wave energy levels. The error analysis begins two days after the simulation initial time, be-20 cause the model needs less than 48 h to reach a steady state from a stationary cold start.

Storm simulation and analysis
In the last two decades, up to 200 damaging wave storms have been reported at the Catalan coast. Although this zone could be considered as well monitored in terms of 25 hydrodynamics (see Sect. 2), the lack of pre and post-storm beach profile data represent an important constraint in terms of storm impacts. However, in the selected storm Introduction (December 2008), pre-and post storm subaerial beach profiles from a LIDAR campaign were available, offering an opportunity for improving morphodynamic response assessment under extreme conditions. The studied event consisted in a Mediterranean tropical-like storm that affected the entire Catalan coast during almost 3 days (from 26 to 29 December). A high pressure 5 centre in northern Europe helped to locate a weak low pressure centre in front of the Catalan coast which kept on growing (due to the difference in air-water temperatures) generating strong eastern winds with speeds of about 50 km h −1 and gusts of up to 85 km h −1 . As a consequence, a big storm travelling from North to South affected the Spanish NE coastal zone, producing severe damages on many coastal infrastructures.  If only the wave height is considered, the storm had an associated return period of about 5 yr. However, its duration and the energy integrated contents arriving to the 20 coast makes the event extraordinary and can be classified with the maximum category that a storm can have in the Catalan coast (Sánchez-Arcilla et al., 2011b).
The impact of this storm was significant shoreline retreat in most of the beaches, with overtopping and flooding in some cases, especially at the northern sector, together with damages in many coastal and harbor structures. There were also three casualties and   Inverse uncertainty quantification through parameter calibration has been performed with the same dataset already presented in Sect. 2. The method consists in fitting a free parameter with the aim of minimising a quadratic cost function (Eq. 1), establishing a ratio of model error to the observed variance. If the cost, χ , is equal to zero, both recorded and computed time-series agree without "discrepancies". Increasing "costs" 5 indicate lower degrees of similarity. Moreover, the quadratic term ensures that large differences weigh heavily in the total χ value. Such a cost function can be written as: Where X and Y refer to the computed and observed signals respectively, N is the 10 number of data points and σ 2 Y is the variance of the observations. The whitecapping dissipation coefficient, normally used to balance wind input (Leckler et al., 2013), has been chosen as the free fit parameter. This choice is based on the assumption that whitecapping is one of the worse known processes in the wave action balance, despite recent efforts to improve its parameterization (Ardhuin et al., 2010). 15 The minimum χ cost obtained with altimeter wave height is 0.37 with a rate of whitecapping dissipation of 2.25 × 10 −5 . Costs between 0.4 and 1.0 could be interpreted as the range in which the variables are well-modeled and there is "enough" predictive skill (Holt et al., 2005).
This dissipation parameter has been obtained with a sample that includes both mod-20 erate and extreme conditions (33.77 % of the data records are above H s = 2 m) and, thus, could be considered to span the different storm stages (i.e. calm, growth and decay). Higher data densities are found in the interval of H s = 0.5 to 1.5 m, despite of statistics showing a wave height mean of 1.81 m and a standard deviation of 2.39 m.
In the same interval, when plotting model values vs. observations, an almost linear fit 25 is revealed (see red line in Fig. 3 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | a forecasting error from 0.5 m to 1.0 m would be within a confidence interval from 70 % (magenta lines) to 95 % (blue lines), respectively. The spatial distribution of the wave height differences has been analyzed with a Normalized Root Mean Square Error approach (NRMSE, Eq. 2). This statistical indicator is obtained clustering the data subset within a 0.5 • × 0.5 • grid and a time window of 5 20 min. In case of having more than one time point on the same grid point, a mean NRMSE is obtained instead (Fig. 4).
Where X and Y refer again to the computed and observed signals respectively. This 10 procedure highlighted that at deep and "open" waters (i.e. with no obstacles such as islands, mountains and peninsulas that modify the meteo-oceanographic flow) the proposed modelling sequence is able to reproduce wave height accurately with values ranging from 0.2 to 0.6. However, sheltered nearshore areas such as the Southern part of France, Northern part of Adriatic Sea, Aegean Sea and the Eastern part of 15 Cyprus confirm the more limited model performance (NRMSE near 1.0) in restricted domains when using the coarse mesh. This phenomena is partially alleviated by a suitable nesting stategy and by implementing the required nearshore physics.
For the sake of completeness, the XIOM buoy network has been compared with the finer mesh (1 × 1 km) results, but also refining the time step discretization and looking 20 for time trends. As expected, correlated wave height costs are found in shallow waters, exhibiting larger errors than for deep waters, with cost values of χ = 0.30 at Llobregat and similar at Blanes (χ = 0.38) and Cap Tortosa (χ = 0.40) buoy positions. At Blanes a 1.5 m difference between predicted and measured wave height was obtained during the first storm peak, although the second one presented a much better fit. In the south- Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | measured vs. 3.5 m modelled) while the second storm maximum is underpredicted. At the same place, it is also observed that a 3 h lag exists for the predicted storm peak. It has to be considered that Llobregat and Cap Tortosa show a smaller error than the Blanes location (the Llobregat buoy for the first storm peak presented "only" a discrepancy of 0.56 m). In these central and southern parts of the coast the storm impact 5 was lower. In all cases, the storm (wave height) growth and decay slopes were well reproduced by the model, which also captured well the storm onset and peak pattern. The peak period and average direction exhibit a similar behaviour, but with slightly higher cost levels (with a maximum value of χ = 0.55), underlining that the model reproduces better the wave pattern in the central coast (Llobregat buoy) than in the other 10 two coastal sectors (northern and southern). Regarding the peak period fit quality, it must be stressed that such variable depends on the directional wave spectral shape, which is more complicated to reproduce than the integral (bulk) properties of the sea state. This makes the fit assessment more difficult and uncertain, particularly when considering the directional veering and spread typical of semi-enclosed domains with 15 variable winds. This illustrates also a "physical" limitation of wave models, that need a minimum time to "react" from sharp wind direction shifts (Bolaños-Sanchez et al., 2007).

Hydro-morphodynamic impact: erosion and flooding
Once the hydrodynamic main driver (wind generated waves) has been quantified, 20 coastal impacts (erosion and flooding) can be assessed using a "response" model. Mataró beaches, both located at the central part of the Catalan coast ( Fig. 1), where the error levels for wave drivers are smallest. The Badalona beach is a typical Mediterranean urban beach with a seafront promenade in its backshore between +5 m and +6 m above MWL. The beach is characterized by coarse sediment, between 350 µm and 600 µm. The main storm impacts reported 5 were severe shoreline retreat and flooding due to run-up episodes.
The Mataró beach is located upstream of the Mataró harbor, in an accretive sector due to the barrier effect of the harbour. It features an alongshore revetment in its backside that protects a railway line from the impact of incoming waves. The beach width is variable, ranging from 10 m at the northern-most part up to 100 m at the south (closer to the harbor barrier). The sediment is also coarse, with median diameter of about 500 µm. The main reported storm impacts were overtopping beyond the revetment and a reshaping of the upper part of the beach, which showed an important growth at the southern part due to the impoundment produced by the harbor breakwaters under storm waves. 15 To start the analysis, prior to modelling tasks, the differences between pre and post storm LIDAR bathymetries are calculated, helping to determine which processes must be reproduced by the morphodynamic model. Hydrodynamic results from Sect. 3 constitute the first step, where the central buoy position (Llobregat buoy) provided the best fit. This buoy is close to the simulated Badalona domain and can be considered to be 20 reasonably near the "upstream" of the Mataró beach. The emerged beach evolution under the storm presented some general features (overall retreat, loss of sedimentary volume and flooding) but was also conditioned by local constraints such as the railway revetment at Mataró or the groyne at Badalona.
In previous sections it has been shown that the December 2008 storm was a two 25 peak event with a predominant Eastern direction. Waves coming from the East induce cross-shore directed water mass fluxes and a longshore current from Northwest to Southwest due to wave obliquity with respect to the shoreline. Cross-shore currents together with high run-ups (enhanced by the associated storm surge) are able to capture subaerial beach sediment. This material could be transported landward through overtopping or offshorewards driven by undertow currents. In addition, long-shore currents advect sand from NW to SW which is deposited at the emerged southern part of the domain, impounded by the barrier presence and, thus, resulting in accretion. This barrier effect is produced by the Badalona port, located at the southern end of the Badalona 5 sector; for the Mataró beach there exists a non-permeable groyne also at the southern part.
During the studied storm there was a significant precipitation volume, reaching daily mean values of 23.1 mm on 26 December and 16.6 mm on 27 December. Storm runoff contained fine material making it difficult to distinguish between the inland or offshore 10 origin of the sediments.
When looking at the detailed hydro-morphodynamics, some differences between the two cases emerge. The Badalona sector (Fig. 5a)  . There is a 2 km long, wide and open beach just outside the analyzed domain (northern part) constituting a sediment source that is transported towards the Badalona sector by longshore currents. This material is retained by a groyne (see green line) located at about 100 m from the North border of the coastal stretch. In addition, crossshore currents drive material towards the offshore part of the profile. From +5.0 m to 20 +2.5 m erosion dominates, relocating the sediment at the shoreface (i.e. from +2.0 m to 0 m). On the contrary, deposition occurs at the central B-C zone where the sewer outfall acts as a sediment trap. In contrast to the northern sector, in this zone crossshore currents drive sediment from the offshore to the inshore. The beach slope at B-C is more reflective than at the B-N zone. Thus, wave overtopping and run-up are higher 25 at B-C, helping to store sediment at the most inshore part of the profile. Finally, the southern part acts as a sink zone, receiving sediment from longshore currents and the sewer outfall but remaining sheltered by the Badalona port. eroding fluxes that take material from the subaerial beach towards the submerged beach, mainly by undertow currents. Such an offshore flux can be advected by longshore currents, that take it from the narrower northern sector (M-N) to the southern side (M-S), storing the sand next to the groyne. The stream mouth (see magenta line) shows a locally important accretion that should come from the river solid discharge, 10 trapped by the low elevation zone present there.
Numerical modelling has provided further insight into how hydrodynamic drivers reshape the beach after a storm event. For this purpose, a locally adapted XBEACH version (Roelvink et al., 2009), termed PREMOS, has been selected to reproduce nearshore hydrodynamics and sediment transport for a range of conditions. This model 15 solves nearshore waves and currents by means of a wave action balance equation (similar to SWAN), coupled with a 2DH Navier Stokes shallow-water system (Stelling and Duinmeijer, 2003). This coupling provides circulation and bed shear stress fields which are used to feed an advection-diffusion equation computing the 2DH sediment concentration patterns. This allows estimating sediment fluxes based on the continuity 20 of sediment mass equation with suitable source/sink terms. From here we can calculate erosion/accretion at every grid node, including the possibility to trigger an avalanching mechanism (van Thiel de Vries, 2009) on the wet-dry (shoreline) beach boundary wherever the slope requires it.
The two studied beaches (Badalona and Mataró) have been discretized with an ir- Where z b,m is the post-storm bathymetry from LIDAR, z b,c is the model output grid and z b,0 is the pre-storm bathymetry from LIDAR. A BSS value of 1 means that the simulated and measured observations fully agree. As these two datasets diverge, the 15 BSS is reduced even reaching negatives values. The threshold value that indicates a prediction as "acceptable" is 0.4. For both studied beaches the LIDAR campaign only covered the emerged part. Hence, the bathymetry has been digitized from the most updated available nautical charts. Due to this restriction, metrics have been calculated only where active cell 20 points are present. Non-erodable points are not taken into account, avoiding BSS values "masked" with non-active points (since zero changes would artificially increase the BSS index). Previous experiences at the Ebro delta (where Cap Tortosa buoy from the XIOM network is located) have resulted in a BSS of 0.44, which can be considered as acceptable (Gràcia et al., 2013b).
Offshore boundary conditions for the hydrodynamic drivers required wave and sea level fields, since the general circulation has been verified to be too small to transport 1707 Introduction 3) every 20 min. Level boundary conditions are generated through blending Barcelona port sea level measurements (provided by Puertos del Estado) and assimilated local atmospheric pressure gauges, with a time sampling of 300 s. Storm surges did not influence significantly the sea level height, except on 5 the 25 December when there was a mean decrease of 8 cm and the second half of 26 December when a mean increase of 12 cm occured. As noted in Figs. 5b and 6b, the simulations have reproduced the observed erosion and accretion patterns. From a qualitative point of view, the model has even been able to capture the influence of local constraints such as the Badalona groyne or the existing 10 revetments. However, the sedimentary input from land, although observable, has not been taken into account due to the lack of reliable data.
With these settings a BSS of 0.36 has been obtained for Badalona beach which can be considered as "acceptable". Better results have been achieved at Mataró beach, despite being a larger area (BSS = 0.38). The main reason of these apparently low index 15 values stems from the inaccuracies in the computed erosion/accretion magnitudes. Figures 5c and 6c show that the simulations are "smoother" than measured values, with a clear underprediction, reaching maximum errors of −0.8 m in the Badalona case, and even worse, 1.5 m in the Mataró case. The reason of these discrepancies will be addressed in the next sesion. 20

Errors and uncertainties
As shown above, the observed hydromorphodynamic patterns have been successfully reproduced, with a relative error level of around 20 %; however, the resulting error in morphodynamics is of order 1.0 m (difference between observations and simulations) which leads to a relative error level of 50 %. The error is propagated in a non-linear  (Cavaleri et al., 2010;Alomar et al., 2010) is the main source of error for wind waves (the key hydrodynamic driver) and, thus, the first step to bound this sequential error would be to deal with the wind field uncertainty. For this purpose a linear correlation between wind velocity measurements and observations, analogous to the one presented in Sect. 3 (i.e. altimeter wind modulus vs. 5 wind field modulus) has been obtained (not shown here). Although the wind field analysis against altimeter observations is out of the scope of this paper, it can be shown that when increasing the wind velocity by 10 % the fit to observations gets "better" (in average) for the area that spans the whole Mediterranean (i.e. D-1 domain, see Fig. 1). This is in accordance with previous results (Bertotti and Cavaleri, 2011) and shows that 10 a 15 % error in wind forecasting is not unusual, especially when sharp gradients occur. An alternative to improve the hydrodynamic results is to perturb the wind field modulus (Alomar et al., 2009). To analyse the sensitivity to the various options that are computationally feasible, the following subset of cases have been selected: no perturbation (A); 20 % increase (B), 10 % increase (D) and 5 % (F) increase; and, correspondingly 20 % 15 (C), 10 % (E) and 5 % (G) decrements. The same computational approach exposed in Sects. 3 and 4 has been here employed.
For the sake of simplicity, hydrodynamic analyses will be presented at shallow waters only, comparing the XIOM buoy network wave height data with the SWAN output from the D-3 domain (i.e. 1 × 1 km mesh). A sample comparison of the observations from 20 two buoys (Llobregat and Blanes locations) vs. the model output using Taylor diagrams (Taylor, 2001) appears in Fig. 7. This diagram summarises observations and model agreement, establishing a graphical relationship among standard deviation (scaled by the number of samples), centered root mean square difference (CRMSD, scaled by the number of samples) and correlation coefficient. The model performance can be 25 assessed from the distance between the observations reference point (LLOB for Llobregat and BLAN for Blanes buoy, respectively) and the model points (see red points for Blanes and blue ones for Llobregat). Note that the reference point position in the x axis is related to the standard deviation intrinsic to each buoy. In other words, since higher waves were recorded and more variability exists at the Blanes location, the standard deviation is slightly larger. Interestingly, each buoy presents distinctive trends but with common traits that support the existence of a link between deep and shallow water wave conditions (through the propagation physics). If the wind modulus field is decreased, better performance 5 (i.e. minimum distance) is found at Blanes, but the opposite occurs at Llobregat. The model skill became worse for extreme cases: B and C did not show any improvement. On the same line, E (−10 %) case at Blanes and F (+5 %) case at Llobregat show shorter distances than the "unperturbed" model. These results reinforce the concept that a perturbation ranging from 5 to 10 % could represent more accurate wind fields in 10 terms of the resulting waves.
In both cases, model and observations could be considered as "well correlated" with mean coefficients near 0.95 at Llobregat and 0.97 at Blanes. Likewise, CRMSD is higher at Blanes (mean value of 0.45) than at Llobregat (0.41). Moreover, distances among the different cases are greater at Blanes than at Llobregat, suggesting that 15 Blanes is more "sensitive" to perturbations. This could be attributed to the fact that the Blanes location is more exposed to Eastern wave storms (more energetic and thus, more capable of producing morphodynamic impact) than the Llobregat location.
Cap Tortosa results (not shown here) are similar to the ones from Blanes: minimum distances are found by decreasing the wind modulus in 10 % (E case). Correlations are 20 lower than the ones described above (0.94), but also in a positive way, with CRMSD (0.4) and standard deviation (1.01) showing the lowest values. These metrics could be explained as the result of the "milder" wave conditions that characterize the Tarragona sector (see Fig. 2), suggesting that simulations for "moderate waves" show better skill than for more extremes (high or low) ones. 25 In previous Catalan Coast analysis (for wave storm periods), it was shown that hydrodynamic drivers play a fundamental role and their accurate forecasting becomes crucial (Gràcia et al., 2013b). Consequently, it was expected to find a wide range of BSS values when perturbing sea level and wave conditions. Directional spectra coming from NHESSD Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the SWAN model simulations altered the metrics but, surprisingly, sea level alterations did not trigger any significant change in the response. This can be due to the reserve beach height (berm) level in some profiles and to the crest level of the promenade or revetment in others (when there was hardly any beach available). In Sect. 4 it was noted that PREMOS results were smoother than the observations (see Figs. 5c and 6c). This could be partially solved by increasing the energetic level of the hydrodynamic drivers, enhancing sediment fluxes that would lead to a sharper morphodynamic evolution. As it can be seen from Fig. 8, in Mataró (see purple line) a moderate increase of drivers (+5 %, F scenario) improves BSS from 0.38 to 0.44. This behaviour reinforces the critical role played by hydrodynamic forcing, whose best model 10 skill is obtained by increasing the wind strength in a moderate way. However, there are still exceptions such as the G case (−5 %), where a decrease of hydrodynamic action improves the subaerial beach accretion in Mataró South (see Fig. 6b), going from +6 m to +4 m and resulting in more balanced metrics. Therefore, when the wind modulus (and, thus, the wave height) is augmented, both 15 the "correct" and "incorrect" patterns would tend to be amplified. The proposed balance concept for the fit quality could be assessed in terms of the bias (see green line). As highlighted by the green line, the bias is proportional to the drivers magnitude; higher BSS are, thus, not necessarily correlated with lower global mean error. For very energetic levels, e.g. B scenario (+20 %), some simulations may even become 20 morphodynamically unbounded and, thus, unstable. Similar conclusions can be derived from the Badalona beach case (see red line), but in this case morphodynamic patterns are more complex due to local constraints (i.e. groyne, sewer outfall and Badalona port). If the wind modulus is decreased by 5 % (G case), the balance between "good" and "poor" patterns leads to higher BSS (from 0.36 25 to 0.40). Note that it is another case where bias increases (see blue line) and BSS increases, indicating that the overall pattern reproduction is improved while the model error is not. It could be explained in terms of a fragile balance between local variables and the overall energetic level. For instance, when the wind input is reduced (−10 % or −20 %, i.e. E and C cases), the observed morphodynamic patterns (see Sect. 4) are not properly developed, resulting in a BSS around 0.2. Moreover, when the wind strength is decreased unrealistically, the sediment fluxes coming from the Northern part of the domain start to bypass the groyne (see red line), reversing the behaviour in the central part [B-C], Fig. 5a and c, from erosion to accretion.

6 Discussion
Uncertainty reduction, as seen in previous sections, depends critically on the accuracy of the forcing fields but is also conditioned by all terms appearing in the governing equations (the balance expression from previous sections). In this balance, the physics, the corresponding parameterizations and the numerical discretizations play all important 10 roles.
As shown in Sect. 3, strong bursts of wind momentum transfer were captured by the wave model, albeit with large errors in areas near the land-ocean border (e.g. Southern France, see Fig. 4). The implemented wave growth parameterizations are based on "moderate" wave conditions, because extreme event data for semi-enclosed domains are scarce and, thus, not well included in the calibrated formulations. Moreover, as sharp gradients appear, the stationarity assumptions are no longer valid and nonlinearity gains in importance.
On the same line, as spatial and temporal scales become shorter (semi enclosed domains), dominant processes have a significantly shorter "life", making the simulated 20 results too dependent on numerical inertia and boundary conditions errors. Since, the time scale of hydro and morphodynamics may differ by one order of magnitude or more, this introduces additional metric-related uncertainty (see Sects. 3 to 5). Higher resolution implies also getting closer to the prediction limits and a higher sensitivity to "numerics" and physics. The Mataró domain (1700 m alongshore) is four times big- Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the same wind input changes "affect" Badalona more than Mataró according to the employed metrics (see Fig. 8).
Sharp gradient modelling (again typical of semi-enclosed domains), requires to retain the maximum amount of information among the different nested or coupled models, establishing "links" that control information exchanges. For this reason, directional wave 5 spectra have been selected as the link between the hydrodynamic and morphodynamic models, instead of using spectral integral parameters (e.g. Hs, Tp etc) that may mask the true driving term characteristics (e.g. two crossing wave trains transporting sediment in opposite directions that are "condensed" into a single average direction that is not physically correct). In all cases due attention should be paid to the consistency 10 of "internal" boundary conditions (wind, wave and current values mutually compatible) and the amount of information (energy) imposed at such boundaries, much more critical for restricted domains than for open sea areas.
Regarding the morphodynamic response under combined drivers, it is a qualitatively well-known field but that is still in its infancy in modeling or quantitative terms. Adding storm surges to incident wave action (see Sect. 4) should have deep and significant morphodynamic consequences, although in some runs they may be numerically limited. A perturbation of +0.1 m in sea level due to surges or uncertainty may provoke a markedly different beach response and/or flooding consequences. This depends on berm or structural crest levels (e.g. the promenade at the Badalona case and the revet-20 ment at the northern part of Mataró case). The coupling between mean sea level and wind-generated waves is particularly important for determining gradients and it is these gradients that control the sedimentary response.
The simulated bottom and coastline evolution will depend on the overall sediment budget, which means considering the land discharge (in this case the sewer outfall at 25 Badalona and the river stream at Mataró). They have not been, however, included in the analysis, due to the lack of quantitative observations but they are expected to have a significant (even if local) effect on the morphodynamic evolution under storm events, when torrential rain discharges in a few hours most of the land-derived sediments throughout the year (see Sect. 4,Figs. 5a and 6a). In order to quantify this discharge, a run-off model coupled with the morphodynamic model could be an interesting solution but also a rather complex one.
The morphological evolution will also depend on the initial bathymetry. In the presented cases, as it usually happens in coastal studies, the emerged beach could be 5 characterized (e.g. using LIDAR) in a more quantitative manner than the submerged part of the profiles, that had to be at least in part interpolated. This introduces another uncertainty in the computations, although in the presented cases the interpolated submerged beach behaved (qualitatively) according to expectations from previous knowledge. The interpolation provides, in addition, the smooth bathymetry that is required to avoid instabilities in the numerical runs, which requires often to artificially "smooth" surveyed isobaths to ensure numerical stability, particularly when sharp gradients are present. Because of that it cannot be stated, in general, that high resolution bathymetry (with abrupt elevation changes) would significantly improve the results with respect to a "smoothed" one. 15 The connection between the emerged and submerged parts of the beach through the swash zone is another important hindrance of present numerical approaches. It depends on non-linear processes such as wave assymetry, skewness and non-linear moments of bottom orbital velocity, including the turbulent component. As it has been noted in Sect. 4, van Rijn (2007) is one of the formulations that takes these processes 20 partly into account. Some of these processes has been included in PREMOS, but in a simplified form and without considering the new knowledge recently derived from large scale laboratory data (Alsina and Cáceres, 2011;Sánchez-Arcilla et al., 2011a).
In addition to this land boundary condition, the beach evolution will depend on the lateral sediment fluxes, which constitute another uncertainty source. These lateral con-25 ditions determine the real sand volume available in a littoral cell and, therefore, the amount that can be "mobilised". The uncertainty can be, however, partially alleviated with nested simulations that move far "enough" these poorly known sediment boundary fluxes. In all the above remarks there is the question of multiple sediment sizes, which is the normal situation in real beaches but that is seldom considered in the simulations, due to the poor performance of transport formulations for a wide granulometric distribution. Effects such as heterogeneous mixing and bottom armoring are still not considered in state of the art models, due to the limited knowledge about the underlying physics.

5
In summary, present models combine multiple limitations that, in addition to natural variability, result in uncertainties at multiple scales which are still hard to quantify, both individually and as part of a joint assessment.

Conclusions
Hydro-morphodynamic predictions in restricted domains show larger errors than simu-10 lations for open sea cases, where temporal and spatial scales show smaller gradients, easier to capture by numerical models. Wave action is the main driver for morphodynamic evolution and the errors in spectral wave parameters near the coast are about twice what has been found in the middle of the Mediterranean basin. The spectral shape has also significant effects on sediment transport since, for instance, the exis- 15 tence of more than one spectral peak results in long waves or crossing wave trains, both of which may modify critically the associated morphological evolution. The appearance of time "lags" between predicted and simulated storm peaks is another type of error that, together with the under-and over-predictions mentioned in the paper may lead to important differences when simulating erosion or flooding processes. 20 The spatial distribution of sediment fluxes introduces further uncertainties in the simulations, whose errors may go from 20 % for hydrodynamics to more than 50 % for the computed morphodynamics. This is due to the integrative nature of the calculated sea bed evolution (responding to wave and surge conditions plus the effects of sediment characteristics, the presence of barriers etc.). The result evolution comes from mutually interacting long-shore and cross-shore transports that must appear explicitly in the computations (e.g. the discussed offshore transport during storms that is then captured by the alongshore flux and leads to enhanced by-pass around the tip of breakwaters or groins). The river and outfall solid discharges should be also considered in the analyses since their signature is apparent in some of the data. The lack of quantitative measurements, however, has precluded considering it in the paper, adding another source of uncertainty. The morphodynamic computations also "integrate" the sediment 5 fluxes at the lateral and shoreline boundaries, all of which are difficult to quantify. This explains the large discrepancies of simulations and observations (up to 1.5 in bed level in the worst cases) and justifies BSS index values around 0.30. The quality of simulated results gets worse in the presence of structures (groins, revetments, etc.), where the modified sediment fluxes are even more uncertain (more complex hydro and mor-10 phodynamics, partial barrier effects, etc.).
The quality of the predictions may be improved by semi-empirical coefficients (e.g. a 10 % increase in the driving wind velocity) or by improved physical formulations (swash zone parameterized fluxes for the shore boundary condition). More energetic conditions may lead to larger errors if the input or boundary conditions are perturbed 15 but they may also converge more quickly towards the "true" solution if these conditions are properly specified.
This illustrates the complex nature and propagation of errors in hydromorphodynamic modeling suites and the need to combine measurements with detailed "physics" to improve the accuracy and robustness of the calculations. This would be the 20 way to achieve more informed coastal decisions. The Mataró BSS is in purple; the Mataró Bias in green. The x axis indicated the "unperturbed" simulation (base value) and the numerical results correspond to increments (or decrements) in the wind field input.