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 by overflowing as well as overtopping processes. Simulations are performed through a nested downscaling strategy from regional to local 15 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 approaches (joint exceedance contour and environmental contour) are used to define 100-year offshore conditions scenarios and to investigate the flood response to each scenario in terms of: 20 (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 considered 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 affected 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 (water level, wave height and period) can 25 produce significant differences 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 sensitive to the chosen scenario with differences of 0.3 to 0.5 m in water level. The developed approach enables one to frame the 100-year hazard and to characterize 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 overtopping/overflowing process ratio and an increase of a factor 4.84 in 30 volumes of water propagating inland and 3.47 in flooded surfaces. Key word : flood hazard, numerical modelling, joint probability, sea level rise, Mediterranean sea


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 hazard and potential disastrous consequences (Hunter, 2012;Tebaldi et al., 2012). On many low-lying 5 coastlines, a "tipping point" is likely to be reached with a mean rise in sea level of 0.5 m (Sweet and Park, 2014).
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 cause 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 10 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 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 15 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., 2006;Hawick, 2014), and hydrodynamic modelling (Martinelli et al., 2010;Gallien et al., 2011;Smith et al., 2012;Fortunato et al., 2013). 20 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 (Tomas et al., 2014), Boussinesq model (Lynett et al., 2010;Andrade et al., 2013), NLSW (Non Linear Shallow Water) model (Suzuki et al., 2011;Guimarães et al., 2015;Leroy et al. 2015), and Non-Hydrostatic Phase-Averaged Model (Smith et al., 2012;Gallien, 2016). These, and especially the SWASH model 25 (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., 2011). 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 30 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 et al., 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 2D. 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, 5 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, 2007) 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 a 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. 10 Overflowing and overtopping process are characterized by very different flow velocity dynamic and can cause differents impacts on structures and building. Using this modelling approach, we aim identify areas prone to one or another kind of flooding and analyze the evolutions of these two kinds of flood hazards related to local mean sea level rise. Due to the specificities of the two kind of hazards, results can be useful to vulnerability studies, to adapt people safety measure, elaborate evacuation plans, or also for risk management actions . 15 Coastal flooding hazards are also usually associated with a return period (i.e. probability of occurrence). The classic approach recommended in several EU for example in the Water Framework Directive (WFD), or in the Areas of Potentially Significant Flood Risk (APSFR) and national directives like the Plan de Prevention des Risques Littoraux (PPRL) 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 20 (for a review, see Jonathan andEwans, 2013 andMonbet, 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 25 applicability (see, among others, Zheng et al., 2013;Wyncoll and Gouldby, 2015;Gouldby et al., 2014). In particular, it does not require any assumption about the dependence 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), 30 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 5 energy systems, etc.), it is usual to refer to reliability methods (Jonathan and Ewans, 2013;Winterstein et al., 1993;Huseby et al., 2013;Huseby et al., 2015). Reliability methods focus on estimating extreme system responses rather than on the combinations of extreme events that produce these responses. In this case, the scenarios are chosen to obtain a system response with the return period under consideration. Such a method 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 10 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 15 according to different statistical methods for defining scenarios with the current mean sea level and with the mean sea level rise 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 20

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 lowenergy mean wave conditions at significant wave height (Hs) = 0.67 m; peak period (Tp) = 4.6 s (observation period: 12/2006 25 -03/2013) and prevailing winds from north-east (Fig 1).

Figure1
The circulation pattern of winter storms, characterized by significant storm surges (0.6 m to 1 m) and very intense wave 30 conditions from the east-south-east (over 6 m in height with peak periods of 10 to 12 seconds), is damaging the seafront and causing recurrent flooding in different parts of the district (seafront, harbour and lagoon passes, Fig 2).

Figure2
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 5 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 Low) (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°x 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 15 from retrospective simulations (NOAA-CFSR-med_10 m forcing) with the SWAN model (Booij et al., 1999) on a Mediterranean grid (42°N-44°N/2°E-8°E) with a resolution of about 1 km (Stépanian et al., 2014). This last source of data (abbreviated here as GuLWa for Gulf of Lion Wave data base), covering a 31-year period  at hourly intervals, was collected at the Candhis 01101 buoy location and used for the statistical analysis, and especially to adjust the marginal distribution of Hs peaks (cf. section 2.4). 20 Table 1 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 25 (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).

Topo-bathymetric data, built structures, surface roughness 30
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).
5 Table 2 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 10 al., 2011;Mignotet al., 2006;Poulter and Halpin, 2008;Néelz et al., 2006). To represent these structures, altimetric data from LIDAR grids (DEM, DTM at 1m 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² 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 15 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 localised retouching should therefore be considered (Poulter and Halpin, 2008;Smith et al., 2012) to incorporate these details into the model. 20 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. 25 Based on these data, two topo-bathymetric models at different spatial resolutions were built up (Fig. 2), one at 20 m resolution (Rank 0: 825 x 827 grid cells) covering the entire stretch of the Salses-Leucate lagoon, and one at 5 m resolution (Rank 1: 606 x 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 30 incorporated into the models through a variable friction coefficient that depends on the soil type and the type of urbanisation according to density (Leroy et al., 2015;Bunya et al., 2010, see Table 3). Table 3 In this study, a spatialised 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 re-sampled and their footprint modified from ortho-photographs to generate a suitable 20 m-resolution roughness map. New urbanizations or land use changes from 2006 were also updated using ortho-5 photographs 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).

Figure 3 15
We used the MARS computing code (Lazure and Dumas, 2007) 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. 20 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. 25 In sectors prone to coastal flooding by overtopping across the seafront, the overtopping volumes are estimated via 1D 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. 30 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 Sflood), the volumes of water reaching inland (written as Vflood) and the spatially variable height of the floodwater (written as Hflood).

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 10 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 north-east and 40 km south from the study site. 15 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 crosschecked with data from the Leucate buoy, with very good fitting (not show).

Flood observations and field measurements 20
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 to estimate water levels reached locally, based on urban landmarks (pavements, walls, jetties, etc). Each 25 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 sectors 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 30 the study area (Fig. 4, Table 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 5 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, 10 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 15 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 20 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 (written as Sflood) and water volumes inland (written as Vflood) 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 values 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 30 ("environmental contours").

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

Data preparation
The available data from Sète make up a continuous series for the 1996-2015 period, corresponding to 16.4 actual years due to record data interruptions. For the statistical analysis, the long-term linear trend in sea-level rise was eliminated and the values Concerning storm dynamics in the Gulf of Lion, focusing respectively on surges and waves, Ullmann (2008) and Gervais (2012) showed that marine storms do not last longer than 3 days. We therefore decided to select the maximum Hs values per 3-day block, with a minimum of 1.5 days between peaks to make sure of their independence. For each peak Hs value, the SWL 15 maximum was then sought within a 12-hour window with the Hs peak at its centre. Each Hs value was associated with the corresponding peak period Tp and peak direction Dp. Several quadruplets (Hs, Tp, Dp, SWL) were thus selected, corresponding to about 6 years of common data covering 111 events/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 Tp and Dp variables are treated as covariables of Hs: as the peak period is highly dependent on Hs and is not an amplitude 20 variable like Hs and SWL, it is considered here as a covariable, as is the peak direction Dp.

Marginal distributions
Adjustment of marginal probability distributions for each variable : when a properly selected high threshold is exceeded, this is modelled via a Generalised Pareto Distribution (GPD). Below this threshold, the empirical distribution ̃ of 25 each variable is used: Where and > 0 respectively are the GPD form and scale parameters and + for ∈ ℝ is defined as + = max( , 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 re-sampled in the same way as to make up the sample of (Hs,SWL) pairs over the common 5 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 χ² test and the Kolmogorov-Smirnov test (Coles, 2001;Nicolae Lerma et al., 2015). The best fit among 3 methods for estimating GPD parameters ( and ), namely maximum likelihood (ML), method of moments (MOM) or probability 10 weighted moments (PWM), was then chosen on the basis of visual and statistical tests (Nicolae Lerma et al., 2015).  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 20 distributions in the extreme range, the amount of data has to 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 (re-analyses). However, in this case, errors attributable to the model (e.g. lack of precision in spatiotemporal resolution, bathymetry or forcing data) are transferred to the statistical analysis and generate uncertainties as to the results (Caires and Sterl, 2005;Mínguez et al., 2012). 25 Bulteau et al. (2015) developed a method (called HIBEVA for Historical Information in Bayesian Extreme Value Analysis) for using historical data from archives to analyse extreme water level values. The flexibility and overall Bayesian framework of HIBEVA justify its use in this study to estimate the marginal probability distribution of significant wave heights via a combination of observation data and model output. The observation data (Candhis) are treated as systematic data and the modelled data (GuLWa) are treated as uncertain historical information. We therefore only used the GuLWa data for 1979-30 2006 (i.e. before the Candhis data came on line).
To estimate the uncertainties relating to the model output data, a comparison (not presented here) was made between the two datasets over the common period from 2007 to 2009. From this we deduced a working hypothesis: for the 1979-2006 period, the true Hs peak values fall within an interval I = [peak Hs from GuLWa -0.15m; peak Hs from GuLWa + 1.60m].
Similarly to the treatment of water levels, the time series (observation data and model output) were first re-sampled taking the maximum Hs per 3-day block, then a statistical us 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 us + 0.15m so that the lower limit of the interval I would be equal to at least us. Where a and b are parameters vectors (one value per parameter for each pair of variables), a threshold to be defined and a vector of residuals. The model is adjusted using the maximum likelihood method on the assumption that the residuals are Gaussian with a mean and variance to be calculated.
For our case study, the threshold selected for Eq3 was 0.75 (expressed as a probability of non-exceedance) using the diagnostic tools described in Heffernan and Tawn (2004). 20

Monte Carlo simulation
The next step was a Monte Carlo simulation to artificially generate Y, keeping to the original proportion of events where each 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 25 the results of the simulation.
Finally, the Gumbel variables Yi 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. Figure 7 3.2 Defining the multivariate scenarios 5

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,Hs) whereby the joint exceedance 10 probability ( > , > ) 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 15 et al., 2002). Practically speaking, this means separating the contour into a number of discrete combinations (SWL, Hs) that will all propagate inland (Fig. 7). The maximum response from these propagations is then associated with a return period T, and written as .
As underlined 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 20 underestimation of the true exceedance probability of the response (Hawkes et al., 2002;Idier et al., 2013;Bulteau et al., in prep). The reason for this is simply that combinations which do not belong to the space ( > , > ) can still produce values for the response variable Z in excess of .

Environmental contour
A second approach involves using environmental contours (written here as enc), which are commonly used in offshore 25 structural engineering (e.g. Huseby et al., 2013Huseby et al., , 2015Jonathan 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. The approach followed here was developed by Huseby et al. (2013Huseby et al. ( , 2015. An environmental contour defined in this way is an (x,y) contour in the space (SWL, Hs), 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 Figure 8): 5 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 .
Here, as in Bulteau et al., in prep, 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 and thus delineate the 10 hazard resulting from the propagation of forcing conditions from the open sea to the coast: Figure 8

Covariates 15
Once the (Hs,SWL) combinations are identified for enc or jec, each Hs must be associated with a value for peak period and peak direction.
In this study, only waves from the 60°-210° sector were retained (cf. above).
The normalised frequency of peak directions observed per Hs segment in the time series of peak Hs from the Candhis buoy (i.e. the sample of systematic data that was used to adjust the GPD law to the Hs with the HIVEBA method) shows that as 20 from Hs > 2.75 m, the most probable peak direction is between 100° and 110° (Figure 1). The value Dp = 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 Hs, taking into account the heteroscedastic relationship between Hs and St whereby the wave steepness tends towards a constant as Hs 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 Hs (Figure 7b).
Based on the data from the Monte Carlo simulation and given that the pattern of change in the Tp vs. Hs relationship tends towards a deterministic law, it was decided to attach a single Tp value to each Hs produced by the combinations selected from enc and jec, taking the median of the periods simulated for each significant wave height considered (Figure 7b). 5 Table 5 shows the characteristics of the quadruplets (Hs,Tp,Dp,SWL) selected for jec and enc scenario respectively (see also  Table 5  10 For each scenario, a 24-hour period of evolutionary conditions (water level, waves, overtopping and propagation of inland flooding) was taken to simulate the storm conditions (including a 2h spin-up period for water level and wave conditions).

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 (Figure 4). From different landmarks across 5 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 under-estimation of about 10 cm for 2014). The wind action (maximum in 2013: 102 km/h, direction 90°N; in 2014: 89 km/h, direction 115°N) on the water height was slight, raising the water level 10 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 15 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 (Figure 9, a). 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 20 produced results that concur with the observations of large overtopping volumes along the seafront (Figure 9 b).

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 25 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. At the contrary, the second method including wave dynamic simulated by the SWASH model produced Vflood levels that were quite close to the Vflood levels deduced from compiled information and GIS treatment. 30 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 (Figure 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 5 neighbourhood (Figure 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 (Figure 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 the light of the information available, it 10 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 exist 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 scenario with the greatest impact for each of the two statistical methods used. The results were analysed in terms of flooded surfaces (Sflood), associated water volumes inland (Vflood) and water height (Hflood) in each of the three zones of the municipality (cf. Fig.  20 2). The scenario with the greatest impact in terms of Sflood and/or Vflood is the Environmental Contour 1 scenario (ENC1). The scenario with the greatest impact using the jec method is JEC2 (Table 7). Table 7 25 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 Sflood and Vflood response is the same (ENC1) for the 3 zones affected (A,B,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 (Figure 11). Flooding in this sector is caused exclusively by 30 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. 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.2 scenario shows that overtopping volumes are much larger along the seafront and also affect the southern part of the urbanized area. With the ENC1-SLR0.6 scenario, the situation is especially critical 15 because, except for the southernmost part of the urbanized area, all buildings are affected by floodwater and access roads are submerged.
In Zone C, Sflood 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 Hflood values of around 20-30 cm. Given the width and morphology of the beach (cf. Fig. 2, profiles 8 to 11), the overtopping sheets of seawater do not reach the seafront buildings (or only slightly in 20 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 sea level rise of 0.6 m (ENC1-SLR0.6), Sflood 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 25 water volumes accumulate in the natural area, submerging roads under several decimetres of water (cf. Figure 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, Vflood increases by 188% and Sflood by 160%. With ENC1-SLR0.6, the situation is critical, with a 384% increase in Vflood and 247% in Sflood.

Overflowing vs overtopping flooded area
Identifying the zones affected by each the two types of flooding shows why both types, overflow (Oflow) and overtopping (Otopp), have to be taken into account to show and characterize the exposure of the Leucate municipality to flood hazards (Fig.   12, Fig. 13).

Figure 12
5 Zone A is affected exclusively by Otopp with the scenarios where mean sea level rise 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 sea level rise of +0.60 m and represents only 11% and 15% of inland Vflood.
Flooding in Zone B is also mainly by overtopping along the seafront, affecting 77% (JEC2) to 84% (ENC1-SLR0.6 and JEC2-SLR0.6) (cf. Table 7). The Otopp/Oflow ratio is fairly stable considering the different SLR scenarios.  The simulation runs for scenarios with no mean sea level rise show that at present, the majority of coastal flooding in the 20 municipality is due to overtopping (Otopp). The low-lying areas affected directly by Otopp and indirectly by the propagation of the resulting water volumes account for 62% of Sflood (38% of these sectors are flooded by Oflow). A moderate sea level rise of less than +0.2 m does not affect this distribution of flooding patterns (Otopp = 63% as against Oflow = 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 Otopp and Oflow, which for the municipality as a whole tends to equalise, with a ratio for Sflood of Otopp = 54% and 25 Oflow = 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 conditions with a given return period, can -once propagation has taken place -produce significantly different results. 30 In this section, we will therefore analyse the differences in Sflood and Hflood obtained after simulating the scenario with the greatest impact defined with each of the statistical methods used and on the assumption of a mean sea level rise of +0.6 m (ENC1-SLR0.6 and JEC2-SLR0.6 scenarios) (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 5 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 (Figure 14). Figure 14 10 However, as Fig. 14a shows, both Sflood and Hflood can differ significantly in Zone C depending on the scenario. For example, differences in Sflood can be observed that are related to the statistical method used (zones in red in Figure 14a). Here, the zones in red are considered to be zones of "uncertainty" as regards characterisation of the hazard. These sectors are not greatly floodprone, if at all, with a JEC2 scenario but may be subject to Hflood 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 m à 0.3 m) to large (red zone from 0.3 m to 0. 5 15 m), and show that the hazard intensifies considerably in the zones subject to threshold effects (topographic hollows).
Significant differences between the JEC2 and ENC1 scenarios were observed in Zone B, and considerable differences in Zone C with the JEC2-SLR0.6 and ENC1-SLR0.6 scenarios.
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 20 0.4 m and 0.3 s respectively for Hs and Tp (i.e. a difference of about 5% in the forcing conditions) produces differences in Hflood > 0.3 m in some streets in the town centre subject to Oflow and Otopp 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 Vflood and Sflood show that a variation of about 5% in the forcing parameters results in Vflood =+13.5% and Sflood =+11.3%. With the SLR0.6 scenario, the relative differences (with Vflood = +8.5% and Sflood =+5.3%) become smaller 25 because zones A, B and 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 Sflood_jec = Sflood_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 30 zones, with Hflood_jec = Hflood_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 Hflood 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 5 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 (PPR) using a fixed elevation and available 10 observations from tide gauges along the French Mediterranean coast (DREAL LR, 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 (Oflow). Our study shows that to map flood hazards, it is just as important to consider the overflow (Oflow) hazard as 15 potential overtopping water volumes (Otopp).
The method applied in this study allowed the Oflow 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 1/3 of the total rise). 20 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 simulation 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 resonance effects in the harbour that are not reproduced by the models used. Furthermore, sea levels in the northernmost pass are substantially underestimated (by 25 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.
On the other hand, 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 studies in the bibliography, the information from recent events against which 30 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 (cf. Gallien et al., 2016) to assess Otopp on the field. It would be necessary, during future storms, to establish measuring protocols based on video data and topo-bathymetric 5 monitoring data before and after the storm, in order to collect more precise data that would help to identify sectors subject to Otopp.
The extreme values 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 Hs and consider long return periods. This would not 10 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. Secondly, 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 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 (cf. §2.4.1). This type of approach has been 25 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 30 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 non-stationarity of the environmental variables under study. Stationarity is a fundamental characteristic of variables required by classic extreme values analysis. Here, we assumed stationarity in the marginal distribution parameters as well as in the dependence structure of the variables Hs and SWL. Long-term trend from the SWL time series was removed before conducting the analysis but seasonal and interannual variability of SWL and Hs have not been dealt with, although this can lead to significant variations of extreme values in time (see e.g. Menéndez et al., 2009aMenéndez et al., , 2009b. However deriving time-dependent ENC and JEC (see e.g. Bender et al., 2014) was beyond the scope of this study.

5
The differences in Hflood 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 10 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, Oflow and Otopp for two recent events consistently with the quantitative and qualitative information available for the site. 15 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 (Vflood, Sflood and Hflood). The largest differences are in Zone B with a sea level scenario based on the current mean sea level, and in Zone C 20 with a mean sea level rise of +0.6m.
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 Sflood_jec response = Sflood_enc and Hflood_jec = Hflood_enc (± 0.1 m), the uncertainty is greater when these conditions are not met. 25 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 30 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 sea level rises. For a 100-year 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 5 flooding / overtopping). On the other hand, 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, built-up 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 10 assumptions of mean sea level rise. 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 15 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 flooding event.

Acknowledgements.
We would like to thank the Conseil Supérieur de la Formation et de la Recherche Stratégiques (CSFRS) for the funding of the CRISSIS research program that initiated this publication. We would like also to thank the municipality of Leucate, which 20 agreed to be a partner in the program. We thank particularly Nicolas Guilpain and Bruno Troqueraud from the capitainerie of Port Leucate, for their help in the field. We gratefully acknowledged all the work from members of the CRISSIS research program. Finally, thanks to Yves François Thomas who initiated this project and whose teachings allowed the original approach of the marine flood hazard developed in this work.