High-resolution marine flood modelling with coupled overflow and overtopping processes : framing the hazard based on historical and statistical approaches

A modelling chain was implemented in order to propose a realistic appraisal of the risk in coastal areas affected as well by overflowing as overtopping processes. Simulations are performed through a nested downscaling strategy from regional to 15 local scale at high spatial resolution with explicit buildings, urban structures such as sea front walls and hydraulic structures liable to affect the propagation of water in urban areas. Validation of the model performance is based on hard and soft available data analysis and conversion of qualitative to quantitative information to reconstruct the area affected by flooding and the succession of events during two recent storms. Two joint probability statistical approaches (joint exceedance probability and environmental contour) are used to define 100 years off-shore conditions scenarios and to investigate the 20 flood response to each scenario in term of: (1) maximum spatial extent of flooded areas, (2) volumes of water propagation inland and (3) water level in flooded areas. Scenarios of sea level rise are also considerate in order to evaluate the potential hazard evolution. Our simulations show that for a maximising 100-year hazard scenario, for the municipality as a whole, 38% of the zones are prone to overflow flooding and 62% to flooding by propagation of overtopping water volume along the seafront. Results also reveal that for the two kind of statistic scenarios a difference of about 5% in the forcing conditions 25 (water level, wave height and period) can produce significant differences responses in terms of flooding like +13.5% of water volumes propagating inland or +11.3% of affected surfaces. In some areas, flood response appears to be very sensible to the scenario chosen with differences of 0.3 to 0.5 m in water level. The approach developed enable to bracket the 100-year hazard and to characterise spatially the robustness or the uncertainty over the results. Considering a 100-year scenario with mean sea level rise (0.6 m), hazard characteristics are dramatically changed with an evolution of the 30 overtopping/overflowing process ratio and an increase of 384% in volumes of water propagating inland and 247% in flooded surfaces. Key word : flood hazard, numerical modelling, joint probability, sea level rise, Mediterranean sea Nat. Hazards Earth Syst. Sci. Discuss., doi:10.5194/nhess-2017-147, 2017 Manuscript under review for journal Nat. Hazards Earth Syst. Sci. Discussion started: 2 May 2017 c © Author(s) 2017. CC-BY 3.0 License.


Introduction
Awareness of the increasing vulnerability of coastal cities to storms and expected effects of global warming lead to more and more studies focusing on the risks of coastal flooding in low-lying coastal areas.These studies often conclude that even a relatively slight rise of mean sea level will, in areas that are not actually exposed or where the hazard is currently manageable, trigger more frequent hazardous and potentially disastrous consequences (Hunter, 2012;Tebaldi et al., 2012).On many low-lying coastlines, a "tipping point" is likely to be reached with a mean rise in sea level of 0.5 m (Sweet and Park, 2014).
Published by Copernicus Publications on behalf of the European Geosciences Union.

208
A. Nicolae Lerma et al.: High-resolution marine flood modelling coupling overflow and overtopping processes Apart from a failure in flood defences, coastal flooding is mainly triggered in two ways.Overflow flooding occurs when the static sea level rises above the level of the natural terrain or flood defence.Overtopping occurs when a combination of high sea level and breaking waves causes successive sheets of seawater to sweep over the seafront.
Coastal flooding hazards are usually defined by the intensity of flooding (spatial extent, water height, flow speed, etc. or a combination of these parameters) associated with the probability of occurrence, usually defined as the "return period".
Low-lying areas exposed to waves can be flooded successively or simultaneously by overflowing and overtopping along the same coastline.In these conditions, risk mapping using simple methods (cross-referencing topography and sea level), decametric digital terrestrial model (DTM) resolutions or without taking buildings into account will not produce adequate or realistic details of the risks to urban areas.High-resolution numerical modelling has therefore become the preferred approach to characterize flooding hazards in the most exposed and vulnerable sites (e.g.Guimarães et al., 2015;Le Roy et al., 2015;Gallien, 2016).
Models of overflow flooding are now relatively accurate and usually based on well-proven physical and numerical methods that have been applied to river, coastal and estuarine contexts and are capable of representing well the extent of flooding.These include the semi-static method (e.g.Breilh et al., 2013), cellular automata (e.g.Dearing et al., 2005;Hawick, 2014) and hydrodynamic modelling (Martinelli et al., 2010;Gallien et al., 2011;Smith et al., 2012;Fortunato et al., 2013).
Models simulating overtopping are much more recent and still require substantial research developments (Hubbard and Dodd, 2002;Gallien, 2016).In the last few years, several process-based models have been developed and applied to address coastal flooding risks: VOF (volume of fluid) model (Tomás et al., 2016), Boussinesq model (Lynett et al., 2010;Andrade et al., 2013), nonhydrostatic phase-averaged model (Smith et al., 2012;Gallien, 2016), and NLSW (nonlinear shallow water) model (Suzuki et al., 2012;Guimarães et al., 2015;Le Roy et al., 2015).These, and especially the SWASH model (Simulating WAves till SHore), are able to reproduce the dynamics of wave surges and overtopping to an appropriate degree of reliability for coastal flooding studies (Suzuki et al., 2012).However, questions remain as to the order of magnitude of overtopping volumes, whether estimated empirically (Laudier et al., 2011;Gallien et al., 2014) or by digital modelling (Smith et al., 2012;Gallien, 2016).In some cases, the estimated value can vary by as much as an order of magnitude (Lynett et al., 2010).These estimated uncertainties as to the reproduction of overtopping volumes can also be attributed to the inadequacy of validation data, which are often qualitative and partial (Battjes and Gerritsen, 2002;Poulter and Halpin, 2008;Reeve et al., 2008;Anselme et al., 2011;Gallien et al., 2012).
As yet, only a few studies have attempted to couple flooding by overflowing and overtopping (Gallien et al., 2014;Stansby, 2013;Le Roy et al., 2015;Gallien, 2016).Le Roy et al. (2015) have attempted to integrate the spatial and temporal variability of overtopping by simulating overtopping in 2-D.This type of model requires a high spatial resolution (less than 2 m).Computing resources and time required to cover sites over several kilometers in extent are still prohibitive.The solution used in our study was to link several models into a chain in order to reproduce, on the one hand, variations in mean sea level, including tides, storm surges and wave setup.This work was realized by coupling a hydrodynamic model (MARS, hydrodynamical Model for Applications at Regional Scale; Lazure and Dumas, 2008) with a spectral wave model (SWAN, Simulating WAve at Nearshore; Booij et al., 1999).On the other hand, to assess runup and overtopping volumes at the seafront, we use an NLSW model (SWASH; Zijlema et al., 2011).The chained modelling enabled us to model the different coastal flooding process (overflowing or overtopping) and consequences at a high resolution over a spatial extent of several kilometers.Overflowing and overtopping processes are characterized by very different flow velocity dynamic and can cause different impacts on structures and building.Using this modelling approach, we aim identify areas prone to one or another kind of flooding and analyse the evolutions of these two kinds of flood hazards related to local mean sea level rise (SLR).Due to the specificities of the two kind of hazards, results can be useful for vulnerability studies, to adapt public safety measures, to elaborate evacuation plans or for risk management actions.
Coastal flooding hazards are also usually associated with a return period (i.e.probability of occurrence).The classic approach recommended in several countries in the EU, for example in the Water Framework Directive, or in the areas of potentially significant flood risk and national directives like the Plan de Prevention des Risques Littoraux in France involves running several scenarios with different return periods (10, 50 or 100 years) plus several scenarios for the same return period.Numerous studies have focused on multivariate extreme value analyses of interdependent meteorological and marine variables (for a review, see Jonathan and Ewans, 2013;Monbet et al., 2007).The complexity of a multivariate extreme value analysis is due to the inadequacy of current knowledge on the interdependence of variables in the tail of the multivariate distribution.An estimation of the tail behaviour is therefore required.Among the existing statistical models used to represent dependence in the tail of the distribution, the semi-parametric model for conditional extreme values first derived by Heffernan and Tawn (2004) is increasingly used in hydrological, coastal and ocean engineering applications due to its great flexibility and applicability (see, among others, Zheng et al., 2014;Wyncoll and Gouldby, 2015;Gouldby et al., 2014).In particular, it does not require any assumption about the depen- dence structure and can be easily extended to larger dimensions (typically larger than 2-3).While the return period in a univariate case is clear and unambiguous (i.e.related to the exceedance probability of the variable), it is much less so when two variables or more are considered together (Salvadori et al., 2011).For example, in the bivariate case, one can consider the "AND" return period (i.e.related to the joint exceedance probability of the variables), the "OR" return period (i.e.related to the exceedance probability of one or the other variable) or other definitions of multivariate return period (see e.g.Salvadori et al., 2016).Moreover, in most risk studies involving several variables, it is considered in the selected scenarios that the return period defined on the basis of the input variables corresponds to the return period of the system response, i.e. flooding in our case.However, the system response is often complex, and whenever a problem addressed has more than one dimension, the one-to-one relationship is rarely valid (Idier et al., 2013).For example, let us imagine we are interested in the 100-year return value of the inundated area.We perform a bivariate extreme value analysis based on wave heights and water level and select scenarios based on the "AND" 100-year return period of input variables.There is no guarantee (and generally it is not the case) that the inundated area resulting from the propagation of such scenarios will be the 100-year return value we want to assess.In the related field of structural engineering (design of coastal defences, offshore renewable energy systems, etc.), it is usual to refer to the environmental contour concept in conjunction with the inverse first-order reliability method (iFORM) (Winterstein et al., 1993;Jonathan and Ewans, 2013;Huseby et al., 2013Huseby et al., , 2015)).Such an approach focuses on extreme system responses rather than on the combinations of extreme environmental loads.The idea is to identify a set of design environmental loads (e.g.contours in 2-D) within which extreme responses with a given return period should lie.In other words, given the failure probability (or return period) of the system response, the objective is to identify what kind of restrictions this imposes on possible designs.This approach is rarely used to assess coastal flooding risks.In this paper, it will be compared with the more classic method, where in choosing the scenarios it is assumed that the return period defined from the input variables corresponds to the return period of the system response, in order to analyse the differences that arise from the methods used to define the scenarios.
This study therefore has two aims: the first is to use tools able to produce realistic representations of flooding by comparing the simulations to existing qualitative and quantitative data from past events and differentiating between the processes causing the flooding.The second aim is to assess, using the same tools, the risk of coastal flooding with a low probability of occurrence according to different statistical methods for defining scenarios with the current mean sea level and with the mean SLR expected during the 21st century.
In the first section, we describe the study site and the methods implemented.The second section presents the results.In the third section, we discuss the methods used and the results, before presenting our conclusions.
2 Flood modelling: data and methods

Study site
Like many other beach resorts in the Languedoc Roussillon area (SW France), our study site is highly exposed to coastal flooding hazards.The municipality of Leucate lies on the western side of the Gulf of Lion, with the Mediterranean on its eastern side and the Salses-Leucate lagoon to the west (Fig. 1).The coast has a microtidal regime (0.2 to 0.4 m) with low-energy mean wave conditions at significant wave height (H s ) of 0.67 m and peak period (T p ) of 4.6 s (observation period: December 2006-March 2013) and prevailing winds from the northeast (Fig. 1).The circulation pattern of winter storms, characterized by significant storm surges (0.6 to 1 m) and very intense wave conditions from the east-southeast (over 6 m in height with peak periods of 10 to 12 s), is damaging the seafront and causing recurrent flooding in different parts of the district (seafront, harbour and lagoon passes, Fig. 2).
On the site the coastal flooding hazard is mainly due to two aspects: the hazard is related to a general low-lying topography, particularly in the inner part of the lido were exchanges of water between sea and lagoon are constrained.Second, the vulnerability is high due to a massive urbanization and the fact that some neighbourhoods have been built directly onto the foredune (Fig. 2).For analysis, three areas were distinguished in the study site: Leucate Plage (Zone A), the naturist village (Zone B) and Port Leucate (Zone C).We also used several beach profiles along the coastline illustrating the spatial variability of sea front in topography and main structures (Wall, Built Sea Front, Back Beach Lows) (Fig. 2).

Forcing data
Different sources of wave data were used for the study (Table 1): (i) observation data from the Candhis 01101 buoy (hourly intervals over a discontinuous period from 2007 to 2015) for local wave parameters used as benchmarks for sea-state modelling and for the statistical analysis; (ii) data from IFREMER MEDNORD, code WaveWatch III, 0.5 • × 0.5 • resolution (IOWAGA project, IOWAGA project, http://wwz.ifremer.fr/iowaga),used as forcing data to model past events; (iii) a time series extracted from retrospective simulations (NOAA-CFSR-med_10 m forcing) with the SWAN model (Booij et al., 1999) on a Mediterranean grid (42-44 • N, 2-8 • E) with a resolution of about 1 km (Stepanian et al., 2014).This last source of data (abbreviated here as GuLWa for Gulf of Lion wave database), covering a 31-year period  at hourly intervals, was collected at the Candhis 01101 buoy location and used both for the statistical analysis and especially to adjust the marginal distribution of H s peaks (see Sect. 2.4).The tide gauge closest to the study site is located at Port la Nouvelle (SHOM/CR LRO), about 15 km to the north of the site (Fig. 1), but could only provide recent data (2013)(2014)(2015).These data were used to reproduce an event in the recent past (November 2014).For our analysis of earlier events and to conduct the statistical analysis, the SHOM/CR LRO tide gauge offshore from Sète, 80 km to the north (Fig. 1), was the only one able to provide sufficient data.The wind data used for the study are from the Leucate semaphore (Météo France data).

Topobathymetric data, built structures, surface roughness
High-resolution modelling of coastal flooding hazards requires a finely detailed representation of the bathymetry and topography.Significant data-collecting efforts were needed to produce an accurate representation of the study sector, including the land-sea continuum, the land areas, the lagoon and passes, the harbour and the nearshore and offshore areas (Fig. 2, Table 2).All data are presented in French national topographic reference (NGF).Numerous studies of the land area have shown that urban structures such as walls and banks can have a determining role in the dynamics of flooding (water flow and extent) and therefore need to be included in the representations of urban environments produced by digital models (Bernatchez et al., 2011;Brown et al., 2007;Fewtrell et al., 2008;Gallegos et al., 2009;Gallien et al., 2011;Mignot et al., 2006;Poulter and Halpin, 2008;Néelz et al., 2006).To represent these structures, altimetric data from lidar grids (DEM, DTM at 1 m resolution) are essential core data.To represent buildings, the necessary data were extracted from the Litto3D DEM via cross-referencing with the "built-up" layer (undifferentiated, industrial and outstanding buildings) from the IGN Topo database (only areas > 20 m 2 were taken into account).This "built-up" layer was then draped over the Litto3D DTM.This enabled us to include only buildings likely to obstruct water flow and to filter out any vegetation or noise in the raw model.
The horizontal resolution (1 m) of the core data and their degree of vertical accuracy, usually ±20 cm, were not sufficient to represent some structural elements that are fundamental in constraining and reproducing inland flows propagation.Some localized retouching should therefore be considered (Poulter and Halpin, 2008;Smith et al., 2012) to incorporate these details into the model.
A ground survey was carried out in June 2015 to set up control points for the different data sources and to make an inventory of elements that were not detectable or only represented discontinuously in the available lidar dataset.The topographic elevations and functional hydraulic characteristics of coastal retaining walls and hydraulic structures liable to affect the propagation of water masses were measured and incorporated so that the DTM grid cells concerned are automatically enhanced by the D-GPS survey values.
Based on these data, two topobathymetric models at different spatial resolutions were built up (Fig. 2), one at 2 m resolution (rank 0: 825 × 827 grid cells) covering the entire stretch of the Salses-Leucate lagoon, and one at 5 m resolution (rank 1: 606 × 1576 grid cells), covering the Leucate municipality (Port-Leucate and Leucate-Plage).Cross-shore profiles at 1 m resolution were also used to model overtopping along the sea front of the study area.
To ensure that flows are properly represented, it is necessary to consider the land use is represented in models.Land use is incorporated into the models through a variable friction coefficient that depends on the soil type and the type of urbanisation according to density (Le Roy et al., 2015;Bunya et al., 2010, see Table 3).
In this study, a spatialized representation of terrain roughness was obtained from a synthetic land use classification based on 2006 Corine Land Cover data.However, the Corine Land Cover data are to the scale of 1 : 100 000, which is not suited to the scale of our study.The data were therefore resampled and their footprint modified from orthophotographs to generate a suitable 20 m resolution roughness map.New urbanizations or land use changes from 2006 were also updated using orthophotographs and field observations.The values used to characterize roughness are those recommended by different sources as applicable to studies conducted in the marine and coastal domains (Bunya et al., 2010;Goutx and Ladreyt, 2001).Specific processing was carried out to represent roads, which are zones where water circulates easily due to the lack of obstacles and the nature of road materials (concrete and tarmac).

Flood modelling chain
Modelling of coastal flooding involved running several chained models: the MARS-2DH hydrodynamic model (Lazure and Dumas, 2008), the SWAN spectral wave model (Booij et al., 1999) and the SWASH NLSW (Zijlema et al., 2011) (Fig. 3).We used the MARS computing code (Lazure and Dumas, 2008) to assess the regional hydrodynamics based on tidal components and meteorological data.The processes represented by the model are associated with long wavelengths only (tides and storm surges).We used the 2DH version of the model, which resolves the Saint Venant equations that govern horizontal free-surface flows in two dimensions, after vertical integration of the Navier-Stokes equations.
When linked to the SWAN wave model, the MARS-2DH model includes short wavelength interactions between waves, sea level and currents (swells and wind sea), mainly in the coastal zone, and can thus calculate the additional water height resulting from wave setup.MARS-2DH thus calculates the speed and direction of currents, averaged to the vertical, and water heights, according to the limit conditions imposed at the edges of the computed domain (boundaries) and the meteorological forcing applied at each node in the model.
In sectors prone to coastal flooding by overtopping across the seafront, the overtopping volumes are estimated via 1-D modelling with the SWASH model (Zijlema et al., 2011).The SWASH model is a time domain model for simulating nonhydrostatic, free-surface and rotational flows.The governing equations are the shallow water equations including a nonhydrostatic pressure term.This model, whose performance in reproducing overtopping volumes was demonstrated by Suzuki et al. (2011), is used here to estimate runup and water volumes likely to overtop seafront walls according to their geometry.The water volumes along the length of the zone concerned are reinjected into the calculation for flooding behind the seafront and seawalls, in order to reproduce the inland propagation of overtopping volumes.
After completing the simulations, the coastal flooding hazard is defined by the intensity of submersion, described here by three types of information: the maximum spatial extent of flooded areas (written as S flood ), the volumes of water reaching inland (written as V flood ) and the spatially variable height of the floodwater (written as H flood ).
2.5 Exploiting historical data: storm conditions and flooded area

Water level and wave conditions
Two storm events in March 2013 and November 2014 were analysed in order to assess the performance of the linked models in reproducing the observed flooding events.These two events were characterized by different marine conditions and consequences in terms of coastal flooding.The data available to characterize the storm conditions are described in Table 1.For the November 2014 storm, the sea level data are from the Port-la-Nouvelle tide gauge.Because no data from this station were available for the March 2013 storm, the water level forcing data are from the Sète tide gauge.However, our analysis of the periods common to both tide gauges covering this stretch of the Gulf of Lion coast (Sète, Port la Nouvelle, Banyuls (Fig. 1) shows that for the events studied, the associated storm surges were fairly uniform along the Languedoc-Roussillon coast.The peak water level of the March 2013 storm was a 0.15 m difference at the Sète (0. 97 m/NGF) and Port Vendres (1.12 m/NGF) tide gauge, located, respectively, at 80 km northeast and 40 km south from the study site.
The wave data used to reproduce these events were extracted from the IOWAGA MEDNORD model at the limits of the domain investigated (rank 0) (Fig. 4).The quality of the reproduction of wave conditions in the domain studied was cross-checked with data from the Leucate buoy, with very good fitting (not shown).

Flood observations and field measurements
To help characterize the quality of coastal flooding modelling, we use several source of available information, from "hard" to "soft" data (Smith et al., 2012).Information was compiled from a wide range of sources: "hard" data from photographs, reports from technical departments and "soft" data from press, interviews and eyewitness accounts.This material enabled us to reconstruct the zones affected by flooding and the succession of events during the storms.Although often qualitative, these observations allowed us to estimate water levels reached locally, based on urban landmarks (pavements, walls, jetties, etc.).Each observation point was then cross-checked against lidar and/or DGPS measurements to produce quantitative information "hard data" from the qualitative validation "soft data" (Table 4).The limits reached by floodwaters in the worst-affected sec-  tors were also mapped with the help of the municipal agents who worked on the ground during the storms.
As no local tide gauge data on water levels in the harbour were available, we were constrained to extract validating material from these documents to assess the quality of the model's reproduction of water levels and flooding at different points across the study area (Fig. 4, Table 4).
During the 2013 storm, a breach observed in the seawall to the north of the municipality caused flooding in a large area of the village.The breach, 15 m in length, occurred because the seawall had not been designed to withstand the full weight of the water accumulating through wave action.Based on the limits and heights of the floodwater described by the Leucate municipal agents and inhabitants, and observed from photographs, we were able to reconstruct the extent and height of the flooding in the village of Leucate-Plage (Fig. 4).The volume of water that flooded the village as a result of the breach was estimated (by cross-referencing topographic and water level data) at a minimum of 37 000 m 3 (this figure is taken as a minimum because several instances of overtopping were observed in the non-urbanized zone to the south, which are hard to quantify) (Fig. 4, Table 4).
Simulations were run to reproduce the breach-induced floods extends and the water heights that have impacted the sector.Two methods using a predetermined breaching scenario (i.e.location, section, time and duration of digging predetermined at the start of the simulation) were used to quantify the water volumes and the flood extend.
The first involved the flooding model only (MARS-flood).Locally, the breach was simulated by applying the laws of hydraulic thresholds (flooded and dewatered conditions) to calculate the upstream to downstream flows from the breach.This incorporated the breach into the grid as a hydraulic singularity without modifying the topography in the model.The geometry of the breach was simplified into a rectangle with a fixed width and a variable threshold over time.
The second method involved simulating the breaking waves by running the SWASH model with a profile facing the breach zone.During the simulation, the seawall was erased at the point in time when the breach occurred.The water volumes coming through the breach were then injected into the flooding model.With this approach, the accelerating speeds and water volumes likely to flow through the breach are not taken into account.
The results obtained in terms of flooding were compared to available information from the ground on the extent of flooded areas (S flood ) and water volumes inland (V flood ) as estimated via GIS methods.

Statistical approach
The aim here was to produce scenarios for offshore marine conditions with a low probability of occurrence that propagate to the shoreline and then inland.To do so, a multivariate extreme value analysis (waves and water levels) was conducted to artificially enlarge the dataset so that scenarios could be selected from the results of two different methods, one based on the return period of offshore marine conditions ("joint exceedance contours") and the other on the return period of the hazard ("environmental contours").

Multivariate extreme values method
The interdependence of offshore forcing variables is modelled here using the semi-parametric approach developed by Heffernan and Tawn (2004), referenced hereinafter as H&T04 method.This approach extrapolates the joint probability density of the offshore marine variables (H s , still water level or SWL) in the extreme values domain by considering the structure of dependency between the variables.For detailed description of the method, readers may consult Heffernan and Tawn (2004), Gouldby et al. (2014) and Wyncoll and Gouldby (2015).Here, we provide an outline of the main steps followed to implement our case study.

Data preparation
The available data from Sète make up a continuous series for the 1996-2015 period, corresponding to 16.4 actual years Concerning storm dynamics in the Gulf of Lion, focusing, respectively, on surges and waves, Ullmann et al. (2008) showed that marine storms do not last longer than 3 days.We therefore decided to select the maximum H s values per 3-day block, with a minimum of 1.5 days between peaks to make sure of their independence.For each peak H s value, the SWL maximum was then sought within a 12 h window with the H s peak at its centre.Each H s value was associated with the corresponding peak period T p and peak direction D p .Several quadruplets (H s , T p , D p , SWL) were thus selected, corresponding to about 6 years of common data covering 111 events per year on average.Given the exposure of the coastline and the wave direction during storms, only waves from the 60-210 • sector were kept for the analysis.
The T p and D p variables are treated as covariables of H s : as the peak period is highly dependent on H s and is not an amplitude variable like H s and SWL, it is considered here as a covariable, as is the peak direction D p .

Marginal distributions
Adjustment of marginal probability distributions F i for each variable X i : when a properly selected high threshold u i is exceeded, this is modelled via a generalized Pareto distribution (GPD).Below this threshold, the empirical distribution F i of each variable is used: where ξ i and σ i > 0, respectively, are the GPD form and scale parameters, and z + for z ∈ R is defined as z + = max(z, 0).
The Languedoc coastline has a microtidal regime that does not warrant the use of indirect methods, i.e. separating the deterministic signal (tide) from the random signal (storm surge), to calculate extreme water levels (Haigh et al., 2010).A direct method was therefore employed to analyse the extreme signal values.
The marginal SWL distribution was calculated from the truncated Sète tide gauge series (see above), i.e. covering about 16.4 years.The time series was first resampled in the same way as to make up the sample of (H s , SWL) pairs over the common time span, i.e. by taking the maximum water level per 3-day block, then a statistical threshold was chosen beyond which the GPD is adjusted to the data.The threshold was chosen by applying several techniques based on a visual appreciation of quantile-quantile graphs, "mean residual life plots", "modified scale and shape parameters plots" and statistical tests such as the χ 2 test and the Kolmogorov-Smirnov test (Coles et al., 2001;Nicolae Lerma et al., 2015).The best fit among three methods for estimating GPD parameters (ξ and σ ), namely maximum likelihood, method of moments or probability weighted moments, was then chosen on the basis of visual and statistical tests (Nicolae Lerma et al., 2015).For the SWL variable, the best fit was achieved with the method of moments beyond the 0.96 m Z.H. threshold (Figs.5a and 6; p value of KS test = 0.98; p value of χ 2 test (10 classes) = 0.82).
The wave observation data (Candhis) cover only 7 years, discontinuously, which is too short a period to extrapolate the distribution of probability in the extreme range and to consider long return periods (typically 100 years).This is a classic problem for any analysis of extreme values from observation data.In order to extrapolate probability distributions in the extreme range, the amount of data must be sufficient to reduce statistical uncertainties to a reasonable level and thus produce meaningful results.When observation data cannot be used or are unavailable, a possible alternative is to use model output (reanalyses).However, in this case, errors attributable to the model (e.g.lack of precision in spa- Similarly to the treatment of water levels, the time series (observation data and model output) were first resampled taking the maximum H s per 3-day block, then a statistical u s threshold was chosen (based on observation data only) beyond which the GPD is adjusted to the data using the HIBEVA method.This also requires a "historical perception threshold".In this case, the threshold was set at u s + 0.15 m so that the lower limit of the interval I would be equal to at least u s .
Figure 5b shows the results of applying the HIBEVA method for H s .The u s threshold is set at 2 m.The chosen GPD parameters (solid red curve) correspond to the mode of the a posteriori joint probability distribution of GPD parameters (see Bulteau et al., 2015, for details).
By combining GuLWa and Candhis data, the actual duration of observations for statistical wave analysis can be extended from 7 years (Candhis data only) to 35 years.With Candhis data only, the maximum return period that could be considered was around 30 years (about 4 times the duration of observations; Pugh, 2004).The maximum return period now is around 140 years.

Fitting the dependency model in the Gumbel space
Original variable X i are transformed into common standard Gumbel margins Y i using the standard probability integral transform.Then, if Y −i is the vector for all variables except Y i , the nonlinear multivariate regression model is as follows: where a and b are parameters vectors (one value per parameter for each pair of variables), ν a threshold to be defined and W a vector of residuals.The model is adjusted using the maximum likelihood method on the assumption that the residuals W are Gaussian with a mean and variance to be calculated.
For our case study, the threshold selected for ν Eq. (3) was 0.75 (expressed as a probability of non-exceedance) using the diagnostic tools described in Heffernan and Tawn (2004).

Monte Carlo simulation
The next step was a Monte Carlo simulation to artificially generate Y , keeping to the original proportion of events where each Y i is a maximum.
For our case study, we simulated 1 110 000 events, representing a fictitious 10 000-year period.These 10 000 years should not be construed as a prediction or forecast for the future, but they are representative of currently available data. Figure 7 shows the results of the simulation.
Finally, the Gumbel variables Y i are transformed back into the original space.The final output is a large sample of artificial offshore sea conditions where at least one variable is extreme (exceeding a defined threshold) and which respects both the individual marginal distributions and the structure of dependence between variables.3.2 Defining the multivariate scenarios

Joint exceedance contour
Once the sample of offshore marine data has been artificially enlarged, scenarios for the return period T considered (here, T = 100 years) were selected for propagation.A commonly used approach in the field of coastal risks involves choosing combinations of forcing factors with a joint exceedance return period equal to T .The idea is then to calculate the joint exceedance contour (written here as "jec"), i.e. the contour (x, y) within the space (SWL, H s ) whereby the joint exceedance probability P (SWL > x, H s > y) is constant (and equal to the probability associated with T ) at every point around the contour (see Fig. 7): where λ is the average number of events per year (111 in our case).
We then need to find the maximum response Z (e.g.flooded area, maximum water height inland) along the contour (Hawkes et al., 2002).Practically speaking, this means separating the contour into a number of discrete combinations (SWL, H s ) that will all propagate inland (Fig. 7).The maximum response from these propagations is then associated with a return period T and written as z jec T .As emphasized previously, this approach rests on the assumption that the return period of the response is equal to the return period of joint exceedance of the input variables.In reality, the joint exceedance probability of the input variables is an underestimation of the true exceedance probability of the response (Hawkes et al., 2002;Idier et al., 2013;Bulteau et al., 2018).The reason for this is simply that combinations which do not belong to the space (SWL > x H s > y) can still produce values for the response variable Z in excess of z T .
Figure 8. 100-year JEC (blue) and ENC (red).B + is the surface delineated by the tangent to the contour which does not contain the convex surface B. The tangent is a linear approximation of the true limit state function (Huseby et al., 2013).

Environmental contour
A second approach involves using environmental contours (written here as "enc"), which are commonly used in offshore structural engineering (e.g.Huseby et al., 2013Huseby et al., , 2015;;Jonathan and Ewans, 2013).These contours are defined within the input variables space but are based on the probability of exceedance of the response variable.These methods rest on an approximation of the limit state curve and are independent from the model.A classical way of defining such environmental contours is to use the inverse firstorder reliability method (iFORM) (Winterstein et al., 1993).
Here we preferred to use the approach developed by Huseby et al. (2013Huseby et al. ( , 2015) ) as it overcomes some limitations of the iFORM and it is especially suited to Monte Carlo simulated datasets.An environmental contour defined in this way is an (x, y) contour in the space (SWL, H s ), outlining a convex inner surface.The probability for the space outlined by the tangent to the contour and not containing the convex surface is constant (and equal to the probability for T ) at every point along the contour (see Fig. 8): A. Nicolae Lerma et al.: High-resolution marine flood modelling coupling overflow and overtopping processes We then need to find the maximum response Z along the contour, and this step is done identically to "jec" approach.The maximum response is then associated with a return period T and written as z enc T .Here, as in Bulteau et al. (2018), in considering these two methods to define scenarios it is assumed that in normal conditions, the two approaches ("jec" and "enc") will calculate upper and lower boundaries of the true response z T and thus delineate the hazard resulting from the propagation of forcing conditions from the open sea to the coast: (5)

Covariates
Once the (H s , SWL) combinations are identified for "enc" or "jec", each H s must be associated with a value for peak period and peak direction.In this study, only waves from the 60-210 • sector were retained (see above).
The normalized frequency of peak directions observed per H s segment in the time series of peak H s from the Candhis buoy (i.e. the sample of systematic data that was used to adjust the GPD law to the H s with the HIVEBA method) shows that as from H s > 2.75 m, the most probable peak direction is between 100 and 110 • (Fig. 1).The value D p = 105 • was therefore retained for future simulations.
To model the peak wave period, we used an approach identical to that of Gouldby et al. (2014): the data for the peak period are first transformed into wave steepness by means of the equation: Next, a conditional regression model in H s , taking into account the heteroscedastic relationship between H s and St whereby the wave steepness tends towards a constant as H s increases, is adjusted to the data (wave characteristics from the sample used to apply the H&T04 method).In the Monte Carlo simulation, a value for the wave steepness (and therefore the peak period) was thus associated with each simulated value for H s (Fig. 7b).
Based on the data from the Monte Carlo simulation and given that the pattern of change in the T p vs. H s relationship tends towards a deterministic law, it was decided to attach a single T p value to each H s produced by the combinations selected from "enc" and "jec", taking the median of the periods simulated for each significant wave height considered (Fig. 7b).

Selecting multivariate scenarios
Table 5 shows the characteristics of the quadruplets (H s , T p , D p , SWL) selected for "jec" and "enc" scenario, respectively (see also Fig. 8).These scenarios were propagated via For each scenario, a 24 h period of evolutionary conditions (water level, waves, overtopping and propagation of inland flooding) was taken to simulate the storm conditions (including a 2 h spin-up period for water level and wave conditions).This simulation time corresponds to the duration of the peak of the storm conditions regularly observed at the study site.For each scenario, the mean water level and wave dynamics at the rank 0 limits are modelled following the shape of the 2013 storm, with concomitant water level and wave peaks at t + 12 h.
To analyse how flood hazards would evolve with the mean SLR anticipated as a result of climate change, the scenarios were run with a uniform mean SLR, with SLR = 0 corresponding to current mean sea level conditions, SLR0.2 = SLR + 0.2 m and SLR0.6 = SLR + 0.6 m.The 0.2 m value of SLR was chosen in order to estimate the impact of a slight SLR (corresponding to a median scenario for 2046-2065 compared to the 1986-2005 global average; source: IPCC WG1 Ch13; Church et al., 2013).The 0.6 m value corresponds to the mean SLR in the Mediterranean forecast by Slangen et al. ( 2014) and Kopp et al. (2014) for 2100 (RCP8.5).It should be stressed here that these are values chosen solely in order to demonstrate changing patterns of hazard in scenarios for a gradual SLR and that considerable uncertainties remain over the values for SLR in the Mediterranean, particularly because of the complex ocean processes taking place in the Gibraltar Straits.

Simulating past events
Reproducing two different flood events makes it possible to assess modelling performance for water levels, overtopping volumes and the reproduction of water flows in the zones most affected during the events.

Simulating flood water levels
The water levels obtained by the simulations were compared with those deduced from the analysis of topographic landmarks photographed during the storms (on jetties, roads, etc.), whether affected or not (Fig. 4).From different landmarks across the entire harbour zone, we were able to determine the mean water level in the harbour during the peak of the storm at 0.85 m NGF ± 5 cm in 2013 and at 1.05 m NGF ± 5 cm in 2014 (Table 4).
The water heights in the harbour obtained by simulation are of a similar order to the water levels estimated from photographs (with a difference of less than 5 cm for 2013 and an underestimation of about 10 cm for 2014).The wind action (maximum in 2013: 102 km h −1 , direction 90 • N; in 2014: 89 km h −1 , direction 115 • N) on the water height was slight, raising the water level by less than 5 cm in the harbour and 3 cm on the seafront during both events.However, the contribution of wave setup (maximum 27 cm in 2013 and 9 cm in 2014) appears to be a determining factor in reproducing water levels observed in the harbour.On the beaches, wave setup contributed up to 50 cm, although we do not have the measurements needed to assess the quality of reproduction of wave setup and runup on the beaches.Nevertheless, photographs taken during the storms show that wave runup regularly overtopped the berm on Port Leucate beach causing accumulation of water in back beach lows but did not produce overtopping at the sea front.Results of the SWASH model simulations concur with these observations.They show that with the given water level and wave conditions, wave runup overtops the first row of discontinuous dunes and fulfils back beach lows zone but without reaching the seafront (Fig. 9a).
The qualitative relevance of the reproduction of wave runup and overtopping during the 2013 storm is also supported by the results obtained for Zone B. In this sector, the first row of buildings sits directly on the upper beach, so that the seafront is affected by wave action during storms.The simulations produced results that concur with the observations of large overtopping volumes along the seafront (Fig. 9b).

Simulating breach flooding
The breach in the seawall during the 2013 storm was caused by both wave action and pressure due to the accumulation of swash water on the structure.
We were not able to reproduce the consequences of the breach with the first method used, because the static water level reached at the height of the storm was below the level of the terrain where the breach occurred.In contrast, the second method, including wave dynamics simulated by the SWASH model, produced V flood levels that were quite close to the V flood levels deduced from compiled information and GIS treatment.
The propagation of the nonstationary water volumes obtained with Method 2 shows the water height varying from 10 to 40 cm and locally exceeding 50 cm (Fig. 10a).The results of the simulation show that the extent of the flooded zones is consistent with observations, although slightly less extensive towards the south.These results also show an underestimation of water heights in the same zone.This is because regular overtopping by sheets of seawater in this sector is not taken into account.Simulations of overtopping only in this sector show considerable amounts of water entering the southern part of the neighbourhood (Fig. 10b).When the propagation of water volumes flooding through the breach are coupled to simulation of overtopping volumes across the seafront, the water heights in the southern zone are better reproduced (Fig. 10c, Table 6).Water volumes from the breach in the seawall and from overtopping propagate through the urbanized area controlled mostly by the topography.Consistent with observations, the combined water volumes flow towards the low-lying parts of the urbanized zone and then towards the natural area to the south, which is lower still.In   light of the information available, it appears that the method used overestimates the extent of the flooded areas (southern part of the neighbourhood).Although this zone is known to have been flooded during the event but being a non-urbanized area, no information exists against which the degree of overestimation can be assessed.

Simulating 100-year return events
All of the joint 100-year scenarios (combined sea level-wave characteristics) were simulated in order to determine the sce-nario with the greatest impact for each of the two statistical methods used.The results were analysed in terms of flooded surfaces (S flood ), associated water volumes inland (V flood ) and water height (H flood ) in each of the three zones of the municipality (see Fig. 2).The scenario with the greatest impact in terms of S flood and/or V flood is the environmental contour 1 scenario (ENC1).The scenario with the greatest impact using the "jec" method is JEC2 (Table 7).
For both types of scenario ("jec" and "enc"), although the processes causing flooding in each zone are different in nature (overflowing and/or overtopping) and patterns (timing and coastal flooding patterns), the maximising scenario in both S flood and V flood response is the same (ENC1) for the three zones affected (A-C).
In the northern zone (Zone A), ENC1 triggers major flooding (in terms of flooded area) across the entire village with water depths generally below 0.50 m but exceeding 1 m locally (Fig. 11).Flooding in this sector is caused exclusively by overtopping across the seafront and affecting the area several hours at a stretch.The water then floods the entire neighbourhood, which is more low-lying than the seafront itself.The floodwater circulates through the whole southern part of the neighbourhood and overflows into the natural area to the south, filling the hollows.
With the ENC1-SLR0.2and ENC1-SLR0.6scenarios, the flooded areas (S flood ) reach further inland, especially towards the north.With the ENC1-SLR0.6scenario, sheets of water sweep across almost the entire seafront of the neighbourhood, and overflowing occurs at the southern extremity.A change can be observed here in the nature of the processes causing flooding, which in turn significantly increases both water volumes (V flood ) and heights (H flood ).The H flood values also almost systematically exceed 1 m.
In Zone B, flooding in the ENC1 scenario is also mainly associated with overtopping.As shown by historical observations, wave action rather than water height is liable to cause the most damage to the seafront.The zone behind it is flooded by the accumulation of overtopping water volumes.The inner edges of the neighbourhood, along the first line of buildings, are also affected by overflow flooding.The ENC1-SLR0.2scenario shows that overtopping volumes are much larger along the seafront and also affect the southern part of the urbanized area.With the ENC1-SLR0.6scenario, the situation is especially critical because, except for the southernmost part of the urbanized area, all buildings are affected by floodwater and access roads are submerged.
In Zone C, S flood values produced by the ENC1 scenario are fairly close to those observed during events in the recent past.Only the harbour zone is affected on the quays by H flood values of around 20-30 cm.Given the width and morphology of the beach (see Fig. 2, profiles 8 to 11), the overtopping sheets of seawater do not reach the seafront buildings (or only slightly in the southern part).With ENC1-SLR0.2, the 0.2 m rise in sea level increases the number of sectors submerged by overflowing floodwater and also accentuates potential overtopping along the seafront.
With a SLR of 0.6 m (ENC1-SLR0.6), S flood extends over a much larger area in zones along the harbour (overflow flooding), with several sectors under more than 0.5 m of water.Floodwater from overtopping propagates from the seafront to the lower parts of the lido, then into numerous areas in the north and centre of Zone C. In the southern part, the overtopping water volumes accumulate in the natural area, submerging roads under several decimetres of water (see Fig. 11).
The simulations with a mean rise in sea level show the extent to which the site is affected in general by threshold effects: with ENC1-SLR0.2,V flood increases by 188 % and S flood by 160 %.With ENC1-SLR0.6, the situation is critical, with a 384 % increase in V flood and 247 % in S flood .

Overflowing vs. overtopping flooded area
Identifying the zones affected by each the two types of flooding shows why both types, overflow (O flow ) and overtopping (O topp ), have to be taken into account to show and characterize the exposure of the Leucate municipality to flood hazards (Figs. 12 and 13).
Zone A is affected exclusively by O topp with the scenarios where mean SLR is less than 20 cm.When only the maximising scenarios for each statistical method (ENC1 and JEC2) are considered, overflow flooding occurs only in scenarios with a mean SLR of +0.60 m and represents only 11 and 15 % of inland V flood .
In Zone C, the O topp /O flow ratio changes considerably with the different scenarios.With the ENC1 and JEC2 scenario, for example, almost all V flood (97 and 98 %) is caused by overflow.The ratio between V flood triggered by O flow and O topp is significantly different in the scenarios with mean SLR of +0.20 and +0.60 m: O flow is still the main process associated with flooding although O topp accounts for about one-third of V flood with the SLR0.2 scenarios and a little under one-quarter of V flood with the SLR0.6 scenario.These differences show, on the one hand, that the characteristics of flooding process are significantly changing with SLR scenario and, on the other hand, that respective contributions do not change linearly and that they depend on topographic particularities and threshold effects.
The simulation runs for scenarios with no mean SLR show that at present, the majority of coastal flooding in the municipality is due to overtopping (O topp ).The low-lying areas affected directly by O topp and indirectly by the propagation of the resulting water volumes account for 62 % of S flood (38 % of these sectors are flooded by O flow ).A moderate SLR of less than +0.2 m does not affect this distribution of flooding patterns (O topp = 63 % as against O flow = 37 %).However, a larger rise in mean sea level of +0.60 m (by 2060-2080 in the IPCC's BAU scenario) significantly affects the ratio between sectors flooded by O topp and O flow , which for the municipality as a whole tends to equalize, with a ratio for S flood of O topp = 54 % and O flow = 46 %.

Determining the 100-year uncertain flooded area
The two statistical methods used to build up the scenarios, i.e. different combinations of offshore marine forcing condi- tions with a given return period, can -once propagation has taken place -produce significantly different results.
In this section, we will therefore analyse the differences in S flood and H flood obtained after simulating the scenario with the greatest impact defined with each of the statistical methods used and on the assumption of a mean SLR of +0.6 m (ENC1-SLR0.6 and JEC2-SLR0.6scenarios) (Fig. 14).The illustration proposed here focuses on the central part of Zone C, because in the built-up sectors in the other zones the differences in extent and water height are relatively slight (mostly less than 0.1 m with both scenarios considered, ENC1 and JEC2).Indeed, most of the differences across the municipality are of less than 0.1 m, which may be considered as not very significant.This order of uncertainty is identical or below that obtained when comparing levels produced by modelling and actually observed during recent events.Furthermore, lidar topographic data are usually characterized by errors below 0.2 m.We have therefore considered that the uncertainty associated with the statistical method chosen is not significant for the zones shown in blue (Fig. 14).
However, as Fig. 14a shows, both S flood and H flood can differ significantly in Zone C depending on the scenario.For example, differences in S flood can be observed that are related to the statistical method used (zones in red in Fig. 14a).Here, the zones in red are considered to be zones of "uncertainty" regarding characterisation of the hazard.These sectors are not greatly flood-prone, if at all, with a JEC2 scenario but may be subject to H flood of 0.1 m to more than 0.5 m with ENC1.
These differences may be considered as moderate (green and yellow from 0.1 to 0.3 m) to large (red zone from 0.3 to 0.5 m) and show that the hazard intensifies considerably in the zones subject to threshold effects (topographic hollows).Significant differences were observed between the JEC2 and ENC1 scenarios in Zone B and in Zone C with the JEC2-SLR0.6and ENC1-SLR0.6scenarios.
Looking now at the marine forcing used for the two types of scenarios, a difference of 0.04 m in the offshore sea level and of 0.4 m and 0.3 s, respectively, for H s and T p (i.e. a difference of about 5 % in the forcing conditions) produces differences in H flood > 0.3 m in some streets in the town centre subject to O flow and O topp hazards.In other words, the response in term of flooding is highly sensitive to variations in the parameters chosen, especially when a rise in mean sea level is considered.The differences for total V flood and S flood show that a variation of about 5 % in the forcing parameters results in V flood = +13.5 % and S flood = +11.3%.With the SLR0.6 scenario, the relative differences (with V flood = +8.5 % and S flood = +5.3%) become smaller because zones A-C are affected by flooding.
Without making an analysis of the sensitivity of the linked models to forcing parameters, which was not the object of this study, our interpretation is as follows: given the statistical approaches used to determine the forcing scenarios to be propagated, one considered to be minimising (jec) and the other maximising (enc), we can consider that if S flood_jec = S flood_enc , the zone is very likely to be subject to a 100-year flood hazard (zones in blue, Fig. 14a).Given the generally small differences in these zones, with H flood_jec = H flood_enc (±0.1 m), we can also consider that the assessment of water heights is satisfactory.However, the zone in red can be considered as a zone of uncertainty in defining the 100-year hazard.
Considering the hazard characterisation for the whole study area, the H flood uncertainty arising from the statistical method used translates into a moderate impact on the spatial extent of flooding.However, the differences locally can be considerable, radically changing the nature of the hazard.
These differences are due above all to threshold effects, when a small change in water height exceeds a topographic threshold and allows propagating a great deal of water inland, which accumulate in topographic hollows.In our case study, these zones are mainly located in zones B and C. In the latter, they only become very evident with the SLR0.6 scenarios.

Discussion
The work undertaken to characterize the flood hazard at the Leucate site is the outcome of a succession of approaches.The first was to apply the recommendations of the French Risk Prevention Plan using a fixed elevation and available observations from tide gauges along the French Mediterranean coast (DREAL, 2008).Subsequently, Anselme et al. (2011) showed that the additional water height caused by wave setup and runup has to be taken into account to approach the values observed during past storms and to characterize the hazards to the seafront.However, the parametric method applied cannot be used to consider the hazard in zones not directly exposed to waves, such as harbour zones where the flooding pattern is different (O flow ).Our study shows that to map flood hazards, it is just as important to consider the overflow (O flow ) hazard as potential overtopping water volumes (O topp ).
The method applied in this study allowed the O flow hazard to be addressed by adding the wave setup contribution into the mean water level reach during the storm.The contribution at the storm surge due to wave setup can reach 50 cm on the beaches and 25-30 cm in the harbour, making it a decisive factor to address flooding along the inner part of the lido (up to one-third of the total rise).
The simulations to reproduce two events in the recent past produced a satisfactory representation of water levels in the harbour (average underestimation of 5 cm for 2013 and 10 cm for 2014).Besides the errors inherent to the simu-  lation method, there may be several reasons for the differences of a few centimetres that appeared between observations and modelling results.These include a lack of precise forcing data, the used of fixed bathymetry and potential reso-nance effects in the harbour that are not reproduced by the models used.Furthermore, sea levels in the northernmost pass are substantially underestimated (by 0.25 to 0.30 m).The underestimation is mainly due to the narrowness of the pass (15 to 20 m) and the potentially highly changeable bathymetry.These characteristics are the reason for the poor reproduction of water flows and levels in this sector but do not appear to alter the results for the other sectors in the studied area.
In contrast, the chain of models was able to handle zones potentially affected by overtopping by estimating the water volumes liable to overtop the seafront.As in the referenced studies, the information from recent events against which the reproduction of overtopping volumes was assessed for accuracy is less detailed for natural zones (few observers) and not easily quantified (overtopping simultaneous with overflow or taken together with rainwater flows).The simulations run did, however, indicate where overtopping occurred (to the south of Leucate Plage, north of the naturist village) or did not occur (Port Leucate beach), concurring with the available qualitative information (wave damage to the seafront, eyewitness accounts).However, this information is not sufficient to assess whether the reproduction of the overtopping volumes is accurate.This highlights the need to produce accurate validation data (see Gallien, 2016) to assess O topp on the field.It would be necessary, during future storms, to establish measuring protocols based on video data and topobathymetric monitoring data before and after the storm, in order to collect more precise data that would help to identify sectors subject to O topp .
The extreme value analysis undertaken in this study to define scenarios for propagation is innovative in two respects.First, by using a Bayesian approach (HIBEVA method), we were able to combine data of different types and different levels of accuracy and thus to calculate the marginal probability distribution for H s and consider long return periods.This would not have been possible by using only Candhis observation data, as the uncertainties over the estimated values would have been too great for return periods of more than 30 years.Second, the definition of offshore forcing scenarios to estimate 100-year coastal flooding hazards was based on two different statistical methods, one producing joint exceedance contours and the other environmental contours.The advantage of using the two methods is that while it is not possible to make a precise assessment of the 100-year flood hazard (there are not enough data on flooding available to analyse the extreme values of response variables directly), it is possible to frame the 100-year flood hazard between the values for the response variables (V flood , S flood and H flood ) that result from propagating the scenarios chosen with the two methods (see Eq. 8).This also gives an indication of the robustness of the result.For example, in our case study, the built-up areas in zones A and B are not very sensitive to the statistical method chosen, which indicates a sufficiently high level of confidence in the estimation of the 100-year hazard in these zones.For Zone C, however, there are notable differences depending on the statistical method applied, reflecting a greater uncertainty in the estimation of the 100year hazard for several neighbourhoods.To overcome this uncertainty arising from the choice of scenarios for propagation, one possible solution is to use a meta-model, which is, in essence, a mathematical approximation of a hydrodynamic model that predicts the modelled responses at a negligible cost in computing time (Idier et al., 2013).In this way, it becomes possible to estimate the response variables directly by "propagating" all the simulated combinations of forcing conditions obtained from the Monte Carlo simulation (see Sect. 2.4.1).This type of approach has been used in the coastal engineering field for regular and continuous modelling (Camus et al., 2011;Idier et al., 2013;Gouldby et al., 2014;Rueda et al., 2015).Unfortunately, in our case study, the complexity of the modelling chain prevents the use of classic meta-modelling techniques, and developing new alternatives is beyond the scope of this study.Additionally, the statistical model contains uncertainties that need to be outlined.In the GPD model, a main source of uncertainties is the choice of the statistical threshold above which the distribution is fitted to the data.Estimated quantiles are indeed highly dependent on the threshold, the selection of which is sometimes difficult and often subjective despite existing statistical tools to help threshold selection (Li et al., 2012).A second source of uncertainty comes from the potential nonstationarity of the environmental variables under study.Stationarity is a fundamental characteristic of variables required by classic extreme value analysis.Here, we assumed stationarity in the marginal distribution parameters as well as in the dependence structure of the variables H s and SWL.The long-term trend from the SWL time series was removed before conducting the analysis but seasonal and interannual variability of SWL and H s have not been dealt with, although this can lead to significant variations of extreme values in time (see e.g.Menéndez et al., 2009a, b).However, deriving timedependent ENC and JEC (see e.g.Bender et al., 2014) was beyond the scope of this study.To go one step further, this issue of potential nonstationarity of the environmental variables questions the relevance of the classic concept of average return period to characterize the risk of coastal flooding.Indeed, the average return period provides information about the probability of exceeding a threshold in any given year.It does not inform about the cumulative risk over a given period of time, which is of interest when it comes to the design of coastal defences for example (see the design life period of the structure and the concept of reliability; e.g.Read and Vogel, 2015;Rootzén and Katz, 2013).This discussion also leads us to question the general framework one uses to assess the risk of coastal flooding.A risk-based approach, starting from the end users needs rather than from a fully scientific analysis only based on physical and statistical considerations, might be better suited to take into account the planning horizon of the study (the design life period in the case of structure design) as well as various aspects such as risk perception (Idier et al., 2013) or economic factors (Rosner et al., 2014).
The differences in H flood between scenarios JEC2 and ENC1 show that threshold effects are liable to notably change the nature of the hazard, with sectors where small differences in forcing (around 5 %) can cause differences in water levels of 30 to 50 cm.It should be remembered here that modelling the inland propagation of coastal flooding is based on significant efforts to integrate terrain roughness, buildings, obstacles and flows and, conversely, on controlling the continuity of flows along the main traffic routes.However, as the propagation models are set at a spatial resolution of 5 m, they may trigger a threshold effect in some sectors (narrow street, topographic irregularity, etc.).

Conclusion
Using a modelling method based on a chain of several MARS-SWAN-SWASH models, we were able to reproduce water levels, O flow and O topp for two recent events consistently with the quantitative and qualitative information available for the site.
Scenarios for the forcing conditions of 100-year joint return period were determined by means of two different statistical methods (joint exceedance contours and environmental contours) in order to analyse the differences arising from the method used to define the scenarios.Simulations of the different 100-year scenarios show that the choice of statistical method used to define the forcing conditions for the scenarios produces notable differences in the response variables considered (V flood , S flood and H flood ).The largest differences are in Zone B with a sea level scenario based on the current mean sea level and in Zone C with a mean SLR of +0.6 m.
Because the "jec" method is minimising and "enc" maximising, using the two types of scenarios enabled us to calculate minimum and maximum values for the spatial extent and height of floodwater, thus framing the 100-year hazard.This also enabled us to characterize the uncertainty over the results that arises from the type of scenario chosen: whereas the results are robust when the S flood_jec response = S flood_enc and H flood_jec = H flood_enc (±0.1 m), the uncertainty is greater when these conditions are not met.In some sectors, this uncertainty can translate into differences of 0.3 to 0.5 m.The simulations of the different scenarios also bring out two major characteristics of the flood hazard in the Leucate municipality.
The first is that the types of flooding that affect the municipality are spatially different.This means that a realistic appraisal of the risk requires joint simulations of flooding by overflow and overtopping.With a maximising 100-year hazard scenario, for the municipality as a whole, 38 % of the zones are prone to overflow flooding and 62 % to flooding by propagation of overtopping water volume along the seafront.
The second is that the nature and scale of the hazard is likely to evolve drastically as the mean SLR.For a 100year event, our results show that overflow flooding affecting built-up zones is limited in extent.The hazard mainly arises from overtopping along the seafront, which is likely to cause significant flooding in the northern part of the municipality (Zone A).Although the hazard increases with a scenario based on a +0.20 m mean sea level (SLR0.2), the newly affected zones are mainly natural areas or roads, with little change in the characteristics of the hazard (ratio between zones affected by overflow flooding and overtopping).In contrast, the SLR0.6 scenarios illustrate what is meant by a tipping point (Sweet and Park, 2014), since they produce a 250 % increase in flooded areas in a 100-year hazard situation, with flooding across the entire municipality, builtup sectors severely affected by overflow flooding (zones A and C) and traffic and evacuation roads becoming almost impassable.
A further point to be made here is that this study focused only on the consequences of climate change under different assumptions of mean SLR.It did not address the consequences of potential changes in marine conditions (waves) or of an intensification of weather conditions during storms.Given the current exposure of the study site to wave overtopping, scenarios assuming an increase in storm intensity (atmospheric surge or wave conditions) would most certainly lead to more intense flooding by overtopping waves and exacerbate the flood hazard in general.
These changes in the flood hazard, and especially in the ratio between zones subject to flooding by overflow and/or overtopping, will not only alter the structural vulnerability of urban areas but also require changes in the messages to be communicated to the public on flood risk awareness and steps to be taken for crisis management in case of a flooding event.
Data availability.Observations and simulation data used in the study can be requested to referents organisms as cited in text.Additional field data can be made available on request to the corresponding author or are consultable on the project website (http: //crissis2015.free.fr/).
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Location map (a) and magnification (b).Circles are for tide gauges locations, crosses for buoys locations and diamond for Leucate meteorological station.The red rectangle delimits the domain used for simulation (R 0 ).

Figure 2 .
Figure 2. Study site and simulation domains: R 0 (extension 16.5 × 16.5 km) resolution 20 m, R 1 (extension 3 × 8 km) resolution 5 m and 11 topobathymetric profiles over three studied zones (A: Leucate Plage; B: the naturist village; C: Port Leucate).Main sea front characteristics are presented as W for sea wall, BSF for built sea front, D for dune, BBL for back beach low and R for road.

Figure 4 .
Figure 4. Example of "hard" information relative to water level during the 2013 and 2014 storm events.Reached water levels on pictures were measured on field using D-GPS in order to estimate quantitative water levels.Red points are related to 2013 event and blue points to 2014 event information (photographs source: Leucate municipal agents).

Figure 5 .
Figure 5. On the left, (a) GPD adjusted to the Sète tide gauge data.Threshold of the law is fixed at 0.96 m Z.H. Parameters of the law are estimated using the method of moments.Confidence intervals are calculated by parametric bootstrap (Mazas and Hamm, 2011).On the right, (b) GPD law adjusted to H s data (Candhis et SCOT) by the HIBEVA method.The threshold is fixed at 2 m.For illustration purposes, the SCOT data are presented by the central values of each interval.

Figure 6 .
Figure 6.PP plot (a) and QQ plot (b) for GPD of SWL (threshold at 0.96 m Z.H.).Hazen plotting position is used for empirical distribution.

Figure 7 .
Figure 7. (a) On the left are the results of Monte Carlo simulation for variables H s and SWL based on 6 common years between SWL data and H s data (Candhis).Declustered data are in black, simulated data (10 000 years) in red.(b) On the right are the results of Monte Carlo simulation for variables H s and T p .Black dots are declustered data.Grey dots are simulated data.In red is the median of the periods simulated given H s .

Figure 9 .
Figure 9. Overtopping observed during the March 2013 event at zone C (a) and zone B (b) sea front.

Figure 10 .
Figure 10.Flood simulation results for the zone A (Leucate Plage) after sea front wall breach (the yellow star indicate the location of the breach).The blue dotted line represents the reconstructed flood extension, and the red line the simulated flood extension.(a) Propagation of the water volume passing through the breach, (b) propagation of overtopping water volume and (c) propagation of the two sources of water.

Figure 14 .
Figure 14.Differences between ENC1-SLR0.6(environmental contour method -scenario 1 with mean sea level rise of 60 cm) and JEC2-SLR0.6(joint exceedance contour method -scenario 2 with mean sea level rise of 60 cm) in the C zone.(a) Maximal extension of flooded area, where blue surfaces are flooded for both scenarios and red surfaces are only flooded for ENC1-SLR0.6scenario.(b) Differences (in metres) in water level.

Table 1 .
Observed and simulated database.

Table 3 .
Land types and Stickler coefficients used.

Table 4 .
Observed vs. simulated water level, qualitative and deduced quantitative information.

Table 6 .
Observed vs. simulated water level reproducing flood after breaching.