# Modeling dependence and coincidence of storm surges and high tide: methodology, discussion and recommendations based on a simplified case study in Le Havre (France)

Amine Ben Daoued, Nassima Mouhous-Voyneau, Yasser Hamdi, Philippe Sergent

► To cite this version:

Amine Ben Daoued, Nassima Mouhous-Voyneau, Yasser Hamdi, Philippe Sergent. Modeling dependence and coincidence of storm surges and high tide: methodology, discussion and recommendations based on a simplified case study in Le Havre (France). Natural Hazards and Earth System Sciences, 2020, 20 (12), pp.3387-3398. 10.5194/nhess-20-3387-2020. hal-03539665

HAL Id: hal-03539665

https://hal.science/hal-03539665v1

Submitted on 21 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Modelling coincidence and dependence of flood hazard phenomena in a Probabilistic Flood Hazard Assessment (PFHA) framework: case study in Le Havre

Amine Ben Daoued<fnref n="1" />☐, Nassima Mouhous-Voyneau<fnref n="1" />, Yasser Hamdi<fnref n="2" />, Claire-Marie Duluc<fnref n="2" />, Philippe Sergent<fnref n="3" />

Received: 8 April 2019 / Accepted: 17 December 2019 / Published online: 8 January 2020
© Springer Nature B.V. 2020

## Abstract

Many coastal urban areas and many coastal facilities must be protected against pluvial and marine floods, as their location near the sea is necessary. As part of the development of a Probabilistic Flood Hazard Approach (PFHA), several flood phenomena have to be modelled at the same time (or with an offset time) to estimate the contribution of each one. Modelling the combination and the dependence of several flooding sources is a key issue in the context of a PFHA. As coastal zones in France are densely populated, marine flooding represents a natural hazard threatening the coastal populations and facilities in several areas along the shore. Indeed, marine flooding is the most important source of coastal lowlands inundations. It is mainly generated by storm action that makes sea level rise above the tide. Furthermore, when combined with rainfall, coastal flooding can be more consequent. While there are several approaches to analyse and characterize marine flooding hazard with either extreme sea levels or intense rainfall, only few studies combine these two phenomena in a PFHA framework. Thus this study aims to develop a method for the analysis of a combined action of rainfall and sea level. This analysis is performed on the city of Le Havre, a French urban city on the English Channel coast, as a case study. In this work, we have used deterministic materials for rainfall and sea level modelling and proposed a new approach for estimating the probabilities of flooding.

**Keywords** PFHA · Joint probability · Dependence · Coincidence · Hazard curve · Aggregation

<figure><img src="image_2.png" /></figure>

<footnote>Amine Ben Daoued
amine.ben-daoued@utc.fr</footnote>

<footnote n="1">¹ Sorbonne University, Université de Technologie de Compiègne, EA7284 Avenues, 60203 Compiègne, France</footnote>

<footnote n="2">² Institute for Radiological Protection and Nuclear Safety, Fontenay-aux-Roses, France</footnote>

<footnote n="3">³ Centre d'études et d'expertise sur les risques, l'environnement, la mobilité et l'aménagement (Cerema), Compiègne, France</footnote>
# 1 Introduction

Many exceptional flood events have occurred in the past in several countries. The flood events of 1984, 2000 and 2003 that impacted France’s Channel coast and particularly the lowlands of the city of Le Havre led to a debate about taking into account the combination of different flood phenomena in flood risk assessment. Partial flooding of the Blayais site in the Gironde estuary in France (1999) triggered several statistical studies on flood hazard [e.g. IPSN (2000)]. The aim was to complete the actions to be carried out on the sites in terms of protection. These studies were considerably accentuated following the Fukushima event (Japan, 2011) (ASN 2013).

The probability of exceedence and/or the frequency of occurrence of an extreme event are key parameters in the design of protection structures of urban sites against floods. Indeed, the approaches traditionally used are purely deterministic, basically because they involve a single (or a few) event scenario (even though this scenario is based on an extreme value analysis). Deterministic models are based on scenarios that are independent of the site specificities and are conservative. Protection must be improved to ensure compliance with the increased safety requirements. To complement knowledge gained from traditional deterministic analysis, the probabilistic approach has been identified as an effective tool for assessing risk associated with hazards as well as for estimating uncertainties. It is in this background that a probabilistic approach must be introduced to complement the deterministic models.

Probabilistic hazard approaches have already been developed in fields other than flooding, such as for probabilistic modelling of earthquakes (Probabilistic Seismic Hazard Assessment or “PSHA”) and tsunamis (Probabilistic Tsunami Hazard Assessment or “PTHA”). References in literature are multiple for these hazards. A typical PSHA method is presented in McGuire and Arabasz (1990). Many research programs use this method, in particular the International Atomic Energy Agency (IAEA)’s work on taking into account the seismic hazard in Probabilistic Safety Assessment (PSA) (IAEA 1993); the thesis (Beauval 2003) carried out in the French Institute for Radiological Protection and Nuclear Safety (IRSN) on the “analysis of the uncertainties in a probabilistic estimate of the seismic hazard, example of France” and more work in the same background can be found in Gupta (2002, 2007), and Bozzoni et al. (2011). Some studies have dealt with PTHA, such as the Tsunamis en Atlantique et MaNche (TANDEM) project (Hebert et al. 2014).

However, experience with Probabilistic Flood Hazard Assessment (PFHA) is limited. In 2013, the U.S. Nuclear Regulatory Commission (NRC) considered that a probabilistic approach is applicable to the flood hazard and established its development in a short term (U.S. Nuclear Regulatory Commission 2014, 2015). In a paper published in 2015 (Bensi and Kanney 2015), Bensi and Kanney presented a framework of a PFHA performed specifically for the storm surge hazard, referred to as Probabilistic Storm Surge Hazard Assessment (PSSHA) for United States Nuclear Power Plants. The joint probabilities method (JPM) was discussed among other existing methods for a probabilistic characterization of storms. Other research offices also refer to the JPM for the probabilistic characterization of storms, including USACE (Norberto et al. 2015). The European project ASAMPSA-E (2013–2016) (Rebour et al. 2016) aimed to define good practices for the development and application of PSA approaches. The flood hazard assessment and its implementation in a risk assessment framework were among the activities of this project. The project also introduces some notions related to the dependence between phenomena. Klugel (2013) described a methodology for the characterization of the external flood hazard (PFHA) in
the case of a riverine nuclear site as well as the linkage of this hazard study with a PSA and its application on a Swiss case study. Some authors use a “Source-Pathway-Receptor-Consequence” method (SPRC) to deal with the evaluation of the flood risk [e.g. Narayan et al. (2011), Horrillo-Caraballo et al. (2013)]. This SPRC model focuses on the linkage between hazard and risk. A similar method called MADS (Méthode d’Analyse du Dysfonctionnement du Système) was used in the “SAO POLO” project to evaluate the submersion risk in Le Havre with a climate change and a multiscale background (Sergent et al. 2012). A research work was carried out at IRSN (2014–2016) and allowed to explore probabilistic approaches and their applicability to flood hazard (Ben Daoued et al. 2016). Indeed, work in the literature is rare in terms of well-established and widely accepted probabilistic studies to characterize the flood hazard by integrating all flood phenomena (local rainfall, river flood, storm surge, tide, etc.) and by considering the combination and the dependence of these phenomena. A Ph.D. research is ongoing in the University of Technology of Compiègne (UTC), with a collaboration with IRSN, to develop a PFHA method that takes into account the combination and the dependence of flood phenomena (Ben Daoued et al. 2018). The PFHA should provide, for high return periods (1000 years for example), a total probability (or frequency) of exceeding a certain water depth at the site regardless the flood initiating phenomenon.

The flood phenomena can be characterized by several random variables, some of which can be correlated. For example, intense rainfall is described by its intensity and duration, the correlation of which is not usually negligible. In addition, numerous studies have shown that univariate frequency analysis does not allow to estimate in a complete way the probability of occurrence of the extreme values in the case of a hazard characterized by several variables (Chebana and Ouarda 2011). According to Salvadori and De Michele (2004), the modelling of the dependence allows a better understanding of the hazard and avoids under/over-estimating the risk.

On the other hand, several pairs of phenomena, such as rain and storm surge, can be independent but occur at the same time or with an offset time. This background of combining flood phenomena is rarely studied in the literature. Only a few references use the term “coincidence” [e.g. Apel et al. (2016)] to denote the occurrence probability of two phenomena at the same time. The notion of the offset time between the two phenomena is not considered. In this paper, we introduce a new method of taking into account the coincidence while studying flood phenomena. This new method consists in taking into account the offset time of the events of the flood phenomena in the calculation of the final hazard curves. This offset time can have a significant impact on both the flood magnitude and its probability. Indeed, the combined occurrence of natural hazards can be very destructive and can have a significant impact sometimes even when the two phenomena are of low intensity, e.g. the flooding of the site of Blayais (France) was the consequence of a conjunction of a high tide and a high storm surge caused by extreme winds (IPSN 2000), the littoral districts of France were hit in 2010 by deadly floods following the Xynthia event which also caused a conjunction of violent winds with high tides. Details on the coincidence as well as the state of the art are drawn further in this paper.

This work aims to develop and implement a probabilistic approach to evaluate the flood hazard (Probabilistic Flood Hazard Assessment or PFHA). This issue is topical for the scientific community in the fields of urban and nuclear because it has been shown that the deterministic approach exploring several scenarios has certain limits. Indeed, univariate and multivariate frequency analyses are methods that are used to estimate probabilities (or return periods) of exceeding (or not exceeding) return levels (rain intensity, sea level, etc.). These approaches are very often deployed in deterministic contexts by designers and
decision makers, usually by taking into account only one scenario and one phenomenon (e.g. designing a sewerage network with a 10-year rainfall return period). However, the work carried out in this paper will make it possible to exploit the know-how used in the deterministic approaches for implementing a novel probabilistic approach. This probabilistic approach aims to identify and combine all possible hazard scenarios to cover all possible sources of the flood risk. These scenarios are then propagated using a hydraulic model and then their contributions are aggregated, with the aim of obtaining hazard curves. The PFHA allows to characterize a variable(s) of interest (maximum water depth, volume, duration, etc.) at different points of a site based on the distributions of the different phenomena of the flood hazard as well as the characteristics of the site.

The paper is organized as follows. Section 2 presents the method of the PFHA in terms of identifying phenomena of flooding and association of random variables, probabilization, scenario construction and propagation and aggregation of hazard curves. In Sect. 3, the case study of Le Havre in France is presented. Observed sea level and rainfall data are analysed in Sect. 4. In Sect. 5, the PFHA method is applied on Le Havre. Finally, we discuss the results and conclude in Sects. 6 and 7.

# 2 Methods

## 2.1 Flood phenomena and random variables

The knowledge of the initiating phenomena allows to better understand the physics of flood phenomena. An initiating phenomenon may be one (or more) of the meteorological conditions (atmospheric pressure, wind, temperature, etc.) underlying the flood phenomena (sea level, rain, etc.). In this paper, we are interested in flood phenomena whose variables are measurable (e.g. rainfall height, sea level). The selection of these flood phenomena and their number depend on the geographical and morphological characteristics of the site as well as the risks to which it is exposed.

A preliminary step in the process consists in identifying and selecting candidate phenomena for the study of the flood hazard. This is based on a “screening” process which aims to identify the phenomena that have a non-negligible contribution to the hazard level. For generalization of this method, we will adopt two phenomena denoted $\phi_1$ and $\phi_2$. In the application of this method, $\phi_1$ and $\phi_2$ represent rainfall and sea level, respectively (as shown further in this paper).

Three sets of hypotheses of increasing complexity were considered to evaluate the combined effects of the two selected explanatory variables (i.e. rainfall intensity and maximum sea level) on the statistical distribution of flood depths:

(1) The phenomena $\phi_1$ and $\phi_2$ are non- coinciding (do not occur together);

(2) The phenomena $\phi_1$ and $\phi_2$ are coinciding (occur simultaneously or with an offset time) and independent;

(3) The phenomena $\phi_1$ and $\phi_2$ are coinciding (occur simultaneously or with an offset time) and dependent.

The flooding of the site can be caused by one of the two phenomena $\phi_1$ or $\phi_2$ independently and separately (case 1) or by these two phenomena when they are independent but occur at the same time (or with a small temporal offset time) (case 2) or when they are dependent
and occur together (case 3). The approach considering all the three cases is an important contribution of this paper, because this aims to take into account all the origins of flood while considering the combination and the possible dependence to estimate total probabilities of flooding. Each of these 3 cases requires a different sampling step.

A random variable (RV) is defined as a quantity describing a random phenomenon. A RV $X$ is a function, often with real values, defined on a sample space. This space corresponds to the set of possible values of a random phenomenon. In hydrology, the RV $X$ can represent the intensity of a rainfall, the duration of a rainfall, the maximum sea level, the water depth produced at a point of interest on a site, the duration of flooding of this point of interest, etc. By convention, the RV’s are designated by capital letters ($X, Y, Z, \text{etc.}$) and the realizations of these RV’s are indicated by lowercase letters ($x, y, z, \text{etc.}$). We denote by $F_X$ and $f_X$ the cumulative distribution function (CDF) and the probability density function (PDF) of a RV $X$, respectively.

$$ 
\begin{aligned}
F_X(x) &= \Pr(X \le x) \\
&= \Pr(\text{the RV } X \text{ takes a value less than or equal to } x) \\
&= \int_{-\infty}^{x} f_X(y)dy
\end{aligned} 
\quad (1) $$

In the following, we consider that the flood phenomena $\phi_1$ and $\phi_2$ are each described by a RV of intensity or amplitude ($V1$ and $V2$) and a duration RV ($D1$ and $D2$). The RV’s $X1$, $X2$ and $X3$ describe the maximum water depths obtained by propagating the RV’s of intensity $V1$ and $V2$ of durations $D1$ and $D2$ corresponding to flood phenomena $\phi_1$ and $\phi_2$ and the coincidence of both $V1$ and $V2$. This joins the event-based concept published in Mazas (2017, 2019), i.e. the RV’s ($V1$, $D1$) and ($V2$, $D2$) represent the source events and the RV’s $X1$, $X2$ and $X3$ represent response events. This event-based framework is already used in the univariate case for extreme sea levels (Bernardara et al. 2014).

A data analysis step on the sequential data of $V1$ and $V2$ (e.g. hourly observations of sea levels and rainfall intensity) is necessary to extract independent samples and to study correlation of these samples, which will allow to validate the chosen case. Frequency analyses (rather univariate or multivariate) are conducted then to obtain statistical models for the RV’s $V1$ and $V2$. These statistical models will be used to build flood scenarios, which represent the event-describing variables (e.g. peak sea level/rainfall intensity and the associated sinusoidal and triangular hydrographs).

## 2.2 Coincidence of flood phenomena

Coincidence is defined by the chance of occurrence of two phenomena at the same time or with an offset time $\delta$. This coincidence does not imply any dependence between the phenomena. Work on this subject in literature is rare.

Apel et al. (2016) distinguish between dependence and coincidence in the combination of fluvial and rainfall floods. Dependence implies a functional or causal relationship between the two phenomena, whereas coincidence does not involve any dependency relationship. If the independence between fluvial and rainy floods is assumed, the coincidence of both phenomena is taken into account by calculating the probabilities of coincidence. The joint probability of occurrence of fluvial and rain events is calculated by the product of the individual probabilities (coming from the rain and fluvial CDF’s) multiplied by the probability of coincidence, i.e. in this context, the probability
of coincidence is the chance that two events of fluvial and rainfall flood occur at the same season, e.g. in Apel et al. (2016) the probability of coincidence of fluvial and rain events is calculated from the length of the flood season and the duration of the fluvial flood events.

Some works in the literature have dealt with this question of coincidence in the case of dependence and in a context of selection of combinations (couples) in particular. Chebana and Ouarda (2011) proposed the decomposition of a contour of iso-probability in a “naive” part (tail of the distribution) and a “proper” part (central), then they proposed to select the couples forming the ends of the clean part. Salvadori et al. (2011) proposed to introduce a weight function and to select a single pair (combination) that minimizes this function. Dung et al. (2015) used the concept of coincidence density in the case of dependence introduced by Volpi and Fiori (2012). The pairs (discharges, volumes) of the bivariate contours have the same exceedence (non-exceedence) probability *p*. But, the probability of coincidence of the pairs on the same contour is not the same, i.e. the likelihood that an extreme discharge value occurs at the same time as a small volume is low (vice-versa). The method of estimating this coincidence density proposed by Volpi and Fiori (2012) consists in defining the occurrence probabilities of the bivariate pairs in curvilinear coordinates. The points of intersection of the bivariate contours with the total dependence line (first bisectrix), called the vertex, define the origin of this coordinate system and form the most probable combinations of the iso-probability contour. The computation of the probabilities of coincidence is realized by the numerical resolution of implicit equations describing the contour of iso-probabilities. In other papers, the distinction between the joint probability of exceedence (or non-exceedence) and the probability of coincidence has not been made, e.g. Yan and Chen (2013) define the coincidence probability of precipitation as the probability of combination calculated through a trivariate frequency analysis making it possible to estimate the distribution of joint probabilities of the studied variables.

In this paper, we propose a new method to take into account the coincidence of independent and dependent flood phenomena in a PFHA context. In a first step, the correlation of the RV's $V1$ and $V2$ is evaluated. This dependence can be quantified by a Chi-plot technique (Abberger 2005) or nonparametric estimators (e.g. upper tail dependence) (Serinaldi et al. 2015). This will allow us then to decide modelling the dependence of the RV's $V1$ and $V2$ using the copula theory (Sklar 1959) (if $V1$ and $V2$ are dependent) and to only consider the univariate CDF's in case of independence. Indeed, we do not aim to show details on how evaluating the dependence or how to model the RV's with a copula. We rather focus on how to deal with the problem of coincidence (in a PFHA context) while the RV's are dependent or not. This will be shown progressively through the next sections of this methodology part.

The next step in dealing with the coincidence is to consider *n* equiprobable offsets $\delta$ between the peaks of the RV's $V1$ and $V2$. This assumes considering one of the RV's as a reference, say $V2$. For each offset time, a set of combinations of couples ($V1, V2$) is considered. The joint exceedence probabilities of these combinations are then affected regarding the dependence or the independence of the RV's. This will be detailed in the following sections.

The non-coincidence case is presented in this paper to serve as a benchmark background that will help interpreting the results of the coincidence case. As it will be shown further, the non-coincidence case is already involved in the coincidence case.
## 2.3 Scenario construction

A scenario is a temporal signal that best represents the physics observed in the rainfall and sea level data. This signal is first defined by a shape (e.g. triangular, sinusoidal, etc.). The maximum amplitude of the scenario (or its intensity) is extracted from the statistical fit of the RV's $V1$ and $V2$, its duration is given by the discrete RV's $D1$ and $D2$ and its probability is provided by the CDF's of $V1$ and $V2$ denoted $F_{V1}$ and $F_{V2}$, respectively. A scenario is therefore a temporal process whose signal form depends on the phenomenon.

We although distinguish between simple scenarios that only concern a flood phenomena and combined scenarios which are composed by combinations of couples ($V1, V2$). In the case of a combined scenario, we consider the combination of two signals (one of $V1$ and one of $V2$). The offset time $\delta$ is defined between the peaks of the two signals. The probability of a combined scenario is given by $Pr\{V1, V2\} = (1 - F_{V1}) \times (1 - F_{V2})$ if the RV's $V1$ and $V2$ are independent and by $Pr\{V1, V2\} = 1 + C(V1, V2) - F_{V1} - F_{V2}$ if $V1$ and $V2$ are dependent, where $C$ is a bivariate copula that links the CDF's $F_{V1}$ and $F_{V2}$.

In this paper, $n$ equiprobable offsets $\delta$ are considered. For each offset time a set of $14(\text{rain}) \times 13(\text{sea level}) = 182$ combined scenarios is considered.

## 2.4 Propagation and hazard curves

This step allows to join deterministic knowledge to the probabilistic concept. It makes it possible to propagate all the scenarios built from the statistical tools from the source of the phenomena $\phi_1$ and $\phi_2$ (measurement station) to the point(s) of interest in the site. These scenarios are injected in the boundary conditions of the model and a water depth variation is obtained at the output, as presented with the example of Fig. 1.

One of the strengths of the proposed method is that it takes into account a set of scenarios. The scenarios are constructed not only from extreme values but also from frequent values of the CDF's.

A variable of interest is chosen according to needs and scientific interests thereafter. In this paper, the variable of interest is the maximum water depth $X$ at a point of interest.

* $X1$ is a RV describing the maximum water depth coming from $\phi_1$;

* $X2$ is a RV describing the maximum water depth coming from $\phi_2$;

* $X3_{\delta_i}$ is a RV describing the maximum water depth coming from a combination of $\phi_1$ and $\phi_2$ with an offset time $\delta_i$;

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Example of a scenario propagation (here a sea level scenario is presented)</figcaption></figure>
Let $h$ be a water depth. We actually assume that $\text{Pr}\{X1 > h\} = \text{Pr}\{V1 > v1\}$. (We suppose that a rain of a probability $p$ gives rise to a maximum water depth of a same probability $p$.) Similarly, $\text{Pr}\{X2 > h\} = \text{Pr}\{V2 > v2\}$. The curves $X1 = f(\text{Pr}\{X1 > h\})$ and $X2 = f(\text{Pr}\{X2 > h\})$ are called hazard curves.

For the combined scenarios, the situation is different, as many combinations with different joint probabilities may give rise to a maximum water depth $X3_{\delta_i}$ that exceeds a same water depth $h$. But, what is the total probability of exceeding $h$? In fact, the water depth $h$ can be exceeded by $X3_{\delta_i}^{C_1}$ OR $X3_{\delta_i}^{C_2}$ OR ... OR $X3_{\delta_i}^{C_n}$, where $X3_{\delta_i}^{C_j}$ is the maximum water depth given by the combination $C_j$ of $V1$ and $V2$ with an offset $\delta_i$. All these combinations give rise to an iso-water depth contour ($h$), as shown in Fig. 2a (only two contours are shown). Theoretically, the total probability of exceeding a water depth $h$ (for a given offset $\delta_i$) is given by the Poincaré formula:

$$ 
\begin{aligned}
\text{Pr}\{X3 > h\} &= \text{Pr}\{(X3_{\delta_i}^{C_1} > h) \text{OR} (X3_{\delta_i}^{C_2} > h) \text{OR} \dots \text{OR} (X3_{\delta_i}^{C_n} > h)\} \\
&= \text{Pr}\{(X3_{\delta_i}^{C_1} > h) \cup (X3_{\delta_i}^{C_2} > h) \cup \dots \cup (X3_{\delta_i}^{C_n} > h)\} \\
&= \sum_{k=1}^{n} (-1)^{k+1} \sum_{1 \le u_1 < u_2 < \dots < u_k \le n} \text{Pr}\left\{\bigcap_{r=1}^{k} (X3_{\delta_i}^{C_{ur}} > h)\right\}
\end{aligned} 
\quad (2) $$

A summary of methods similar to that presented through Eq. 2 and Fig. 2 can be found in DEFRA (2005), and Hawkes et al. (2002).

Figure 2b presents the resolution of Eq. 2 for the water depth $h$. The intersection terms in Eq. 2 corresponds to parts that are common between the combinations (of rain intensity and sea level) in the iso-water depth contour (orange cross-shaped points). First, the joint probabilities of all combinations on the iso-water depth contour (red dots in Fig. 2b)

<figure><img src="image_4.png" /><figcaption>**Fig. 2** a Construction of iso-water depth contours: red dots are propagated scenarios that give rise to a same water depth (only two contours are shown); **b** illustrating the method of calculating a total probability of exceeding $h$ using Eq. 2; ASL = NGF + 4.38 m</figcaption></figure>
are summed. Then, the joint probabilities of the intersection common parts (orange cross-shaped points in Fig. 2b) are subtracted from the total probability.

By applying Eq. 2 for the different iso-water depth contours of Fig. 2b, we obtain a hazard curve $X3_{\delta_i} = f(\text{Pr}\{X3_{\delta_i} > h\})$ for the offset $\delta_i$. We get as many hazard curves as offsets.

As it was shown in the last section, the joint probabilities $\text{Pr}\{X3_{\delta_i} > h\}$ are given by $(1 - F_{V1}) \times (1 - F_{V2})$ if the RV's $V1$ and $V2$ are independent and they are equal to $1 + C(V1, V2) - F_{V1} - F_{V2}$, otherwise.

We currently have three situations. The first one concerns the case of non-coinciding phenomena, we have two hazard curves each giving the maximum water depths from $\phi_1$ and $\phi_2$ according to their exceedence probabilities. The point of interest can be flooded with only $\phi_1$ OR with only $\phi_2$. Furthermore, in the second and the third situations concerning coinciding-independent/dependent phenomena, we have 13 hazard curves $X3_{\delta_i} = f(\text{Pr}\{X3_{\delta_i} > h\})$ for each case. But in all these three situations, what is the total probability that a water depth $h$ is exceeded? The answer to this question will be detailed in the next section devoted to the aggregation of the hazard curves.

## 2.5 Aggregation and final hazard curve

Aggregation provides decision makers with a more complete information about return periods of exceeding water depths (or other variables of interest) at a point of interest in the site. This is a key step and is one of the contributions of this paper. This step makes it possible to obtain a single (total) statistical probability that is aggregated from several probabilities coming from several phenomena.

### 2.5.1 Non-coincidence case

We aim to calculate a total probability of exceeding any water height $h$. There will be a part of this total probability that will come from each of the two hazard curves, as the water height $h$ can be exceeded either by only rain or by only sea level. Theoretically, the following total probability should be calculated:

$$ 
\begin{aligned}
P_{\text{total}} &= \text{Pr}\{X > h\} = \text{Pr}\{(X1 > h) \text{ or } (X2 > h)\} \\
&= \text{Pr}\{(X1 > h) \cup (X2 > h)\} \\
&= \text{Pr}\{X1 > h\} + \text{Pr}\{X2 > h\} - \text{Pr}\{(X1 > h) \cap (X2 > h)\}
\end{aligned} 
\quad (3) $$

By writing $P_{\text{total}} = \text{Pr}\{(X1 > h) \text{ or } (X2 > h)\}$ we already suppose that the data sets of $V1$ and $V2$ are non-coinciding, which means that they are disjoint. (This notion will be presented in the next section.) Equation 3 is then written as following:

$$ 
\begin{aligned}
P_{\text{total}} &= \text{Pr}\{X > h\} = \text{Pr}\{(X1 > h) \text{ or } (X2 > h)\} \\
&= \text{Pr}\{(X1 > h) \cup (X2 > h)\} \\
&= \text{Pr}\{X1 > h\} + \text{Pr}\{X2 > h\}
\end{aligned} 
\quad (4) $$

By iterating Eq. 4 for several water depths $h_i$, we obtain the final hazard curve: (Water depth) = $f(P_{\text{total}})$. Where $P_{\text{total}}$ are total probabilities of exceeding water depths.

The probabilities of Eq. 4 come from the rainfall and sea level CDF's. In order to obtain annual exceedence rates, all probability terms in Eq. 4 are multiplied by their annual event rate $\lambda$. This rate is greater than or equal to unity if the sample comes from a
“Peaks-Over-Threshold” (POT) declustering process and it is equal to 1 if the Annual Maxima (AM) declustering process is used (only one event per year is retained). Equation 4 then becomes:

$$ \lambda\{X > h\} = \lambda_1 \times \Pr\{X1 > h\} + \lambda_2 \times \Pr\{X2 > h\} \quad (5) $$

where

* $\lambda_1$: is the annual rate of V1 events;

* $\lambda_2$: is the annual rate of V2 events.

Some studies in the literature use an equation similar to Eq. 5. This is particularly the case for Probabilistic Seismic Hazard Assessment (PSHA) and Probabilistic Tsunami Hazard Assessment (PTHA). In PSHA and PTHA, Eq. 5 is classically taken in its integral form in which the PDF’s of seismic sources and other seismic parameters as well as the occurrence rates of events are involved (McGuire and Arabasz 1990; Beauval 2003). Bensi and Kanney (2015) presented a similar equation (in an integral form) in a framework of a Probabilistic Storm Surge Hazard Assessment (PSSHA). Other work is also carried out by the USACE (Norberto et al. 2015), using an equation similar to Eq. 5 in a context of aggregation of tropical and extra-tropical storm parameters.

The annual exceedance probabilities (AEP) can be calculated from the annual exceedance rate. This calculation is based on the Poisson process and is already done in a similar way by Lang et al. (1997) to get annual exceedence return periods from a rate-based return period (a return period with many over-threshold exceedences). This process assumes that the events are independent, an assumption that is already assumed in this paper. The Poisson’s law is written as follows:

$$ P(X > h, T) = 1 - e^{-\lambda \times T} \quad (6) $$

where $\lambda$ is the annual rate of flood events, whose maximum depth $X$ exceeds $h$, during a duration of exposure $T$. For floods (as for tsunamis), $T$ is taken equal to 1 year. Thus the resulting probability is annual. The probability $P(X > h, T)$ of Eq. 6 can be approximated to $\lambda$ (for $T = 1$ year).

In order to obtain the AEP’s, as accurate as possible, we transform the exceedence rates of Eq. 5 into AEP’s using the Poisson process given by Eq. 6. The following equation is obtained:

$$ \begin{aligned} \text{AEP} &= 1 - e^{-\lambda\{X>h\}\times T} \\ &= 1 - e^{-[\lambda_1\times\Pr\{X1>h\}+\lambda_2\times\Pr\{X2>h\}]} \end{aligned} \quad (7) $$

The exceedence return periods of the aggregated hazard curve are calculated using the following equation:

$$ \begin{aligned} T_{\text{total}} &= \frac{1}{\text{AEP}} \\ &= \frac{1}{1 - e^{-[\lambda_1 \times \Pr\{X1 > h\} + \lambda_2 \times \Pr\{X2 > h\}]}} \end{aligned} \quad (8) $$
### 2.5.2 Coincidence case

In the coincidence case, there are n hazard curves, each corresponds to a different offset time $\delta_i$ between the RV's V1 and V2. We therefore aim to calculate a total probability of exceeding any water depth h.

We consider that all offsets are impacting, a part of the total probability will then come from each of the n hazard curves.

Theoretically, the total probability is given by the Poincaré formula:

$$ 
\begin{aligned}
P_{\text{total}} &= \Pr\{X3 > h\} \\
&= \Pr\{(X3_{\delta_1} > h \text{ AND } \delta = \delta_1) \text{ OR } (X3_{\delta_2} > h \text{ AND } \delta \\
&= \delta_2) \text{ OR } \dots \text{ OR } (X3_{\delta_n} > h \text{ AND } \delta = \delta_n)\} \\
&= \Pr\{(X3_{\delta_1} > h \cap \delta = \delta_1) \cup (X3_{\delta_2} > h \cap \delta = \delta_2) \cup \dots \cup (X3_{\delta_n} > h \cap \delta = \delta_n)\} \\
&= \sum_{k=1}^{n} (-1)^{k+1} \sum_{1 \le i_1 < i_2 < \dots < i_k \le n} \Pr\left\{\bigcap_{j=1}^{k} (X3_{\delta_{ij}} > h \cap \delta = \delta_{ij})\right\}
\end{aligned} 
\quad (9) $$

We actually should distinguish between two “disjoint” and “independent” events:

* V1 and V2 are disjoint (non-coinciding): no common events $\rightarrow \Pr\{V1 \cap V2\} = 0$;

* V1 and V2 are independent: presence of common events and these common events are independent $\rightarrow \Pr\{V1 \cap V2\} = \Pr\{V1\} \times \Pr\{V2\}$;

The offset variable $\delta$ cannot be equal to $\delta_i$ AND $\delta_j$ at the same time, the elements are then disjoint. Equation 9 can then be simplified as following:

$$ 
\begin{aligned}
P_{\text{total}} &= \Pr\{X3 > h\} = \Pr\{X3_{\delta_1} > h \cap \delta = \delta_1\} + \Pr\{X3_{\delta_2} > h \cap \delta = \delta_1\} \\
&\quad + \dots + \Pr\{X3_{\delta_n} > h \cap \delta = \delta_1\} \\
&= \Pr\{\delta = \delta_1\} \times \Pr\{X3_{\delta_1} > h\} + \Pr\{\delta = \delta_2\} \times \Pr\{X3_{\delta_2} > h\} \\
&\quad + \dots + \Pr\{\delta = \delta_n\} \times \Pr\{X3_{\delta_n} > h\}
\end{aligned} 
\quad (10) $$

The offsets are equiprobable, their weight equal $\frac{1}{n}$. Moreover, if we suppose that the offset RV $\delta$ and the maximum water depth RV X3 are independent, then Eq. 9 becomes:

$$ 
\begin{aligned}
P_{\text{total}} &= \Pr\{X3 > h\} = \Pr\{X3_{\delta_1} > h \cap \delta = \delta_1\} + \Pr\{X3_{\delta_2} > h \cap \delta = \delta_1\} \\
&\quad + \dots + \Pr\{X3_{\delta_n} > h \cap \delta = \delta_1\} \\
&= \Pr\{\delta = \delta_1\} \times \Pr\{X3_{\delta_1} > h\} + \Pr\{\delta = \delta_2\} \times \Pr\{X3_{\delta_2} > h\} \\
&\quad + \dots + \Pr\{\delta = \delta_n\} \times \Pr\{X3_{\delta_n} > h\} \\
&= \frac{1}{n} \times \Pr\{X3_{\delta_1} > h\} + \frac{1}{n} \times \Pr\{X3_{\delta_2} > h\} + \dots + \frac{1}{n} \times \Pr\{X3_{\delta_n} > h\}
\end{aligned} 
\quad (11) $$

The annual exceedence rates are given by:
$$ \lambda\{X3 > h\} = \lambda_{m_{\delta_1}} \times \frac{1}{n} \times \Pr\{X3_{\delta_1} > h\} + \lambda_{m_{\delta_2}} \times \frac{1}{n} \times \Pr\{X3_{\delta_2} > h\} + \dots + \lambda_{m_{\delta_n}} \times \frac{1}{n} \times \Pr\{X3_{\delta_n} > h\} \quad (12) $$

where $\lambda_{m_{\delta_i}}$ is the mean annual rate of (V1, V2) events whose offset time is $\delta_i$.

The AEP's are calculated in the same way as in the previous section, using the Poisson process given by Eq. 6. We obtain the following equation of the AEP's:

$$ \begin{aligned} \text{AEP} &= 1 - e^{-\lambda\{X3>h\}\times T} \\ &= 1 - e^{-[\lambda_{m_{\delta_1}} \times \frac{1}{n} \times \Pr\{X3_{\delta_1}>h\} + \lambda_{m_{\delta_2}} \times \frac{1}{n} \times \Pr\{X3_{\delta_2}>h\} + \dots + \lambda_{m_{\delta_n}} \times \frac{1}{n} \times \Pr\{X3_{\delta_n}>h\}]} \end{aligned} \quad (13) $$

The exceedence return periods of the aggregated hazard curve are calculated using the following equation:

$$ \begin{aligned} T_{\text{totale}} &= \frac{1}{\text{AEP}} \\ &= \frac{1}{1 - e^{-[\lambda_{m_{\delta_1}} \times \frac{1}{n} \times \Pr\{X3_{\delta_1}>h\} + \lambda_{m_{\delta_2}} \times \frac{1}{n} \times \Pr\{X3_{\delta_2}>h\} + \dots + \lambda_{m_{\delta_n}} \times \frac{1}{n} \times \Pr\{X3_{\delta_n}>h\}]} } \end{aligned} \quad (14) $$

## 3 Case study: Le Havre

The urban city of Le Havre is located in the Seine-Maritime department, on the English Channel coast in Normandy (France). The map of Fig. 3 shows the situation of the city of Le Havre in the northwestern of France. The name of the city of Le Havre means “the port” (or the harbour), which is among the largest ports in France.

Due to its location on the coast of the Channel, Le Havre's climate is temperate oceanic. Days without wind are rare. There are maritime influences throughout the year. According to the meteorological records, precipitation is distributed throughout the year, with a maximum in autumn and winter linked to North-Atlantic extra-tropical storms. The months of June and July are marked by some relatively extreme thunderstorms on average 2 days per month. One of the characteristics of the region is the great variability of temperature, even

<figure><img src="image_5.png" /><figcaption>Fig. 3 Site of Le Havre</figcaption></figure>
during the day. The prevailing winds are from north-northeast for the breezes and, from the southwest for strong winds.

The city of Le Havre is crossed from east to west by a “dead” cliff which marks an old separator between the upper town and the lower town. The lower city is therefore built in a former intertidal space, that is to say the space of propagation of tidal waves, where the sea level evolves between high tide and low tide. This so-called “low” city has been the site of various and varied developments: channelization of the watercourse, embankments, construction of basins, etc. However, this lower city and the port remain naturally subject to the risk of marine submersion (Deglaire 2010; Sergent et al. 2012). Several relatively large flood events have affected the Channel coast of France. Some of these events are indicated in the municipal archive of the city of Le Havre. Major floods in the lowland areas of the city of Le Havre took place in June 1982 following torrential rains, causing serious losses in several neighbourhoods and a major landslide on a cliff in the city. Major floods were also caused by the torrential rains of a storm that devastated northern France in December 1999. Other flood events have occurred following storms, including Xynthia in 2010 and Eleanor in 2018 which caused a lot of damage in the city of Le Havre. Indeed, one of the most important characteristics of Le Havre, as a case study, is the fact that it is a city subject to marine submersions and instabilities of coastal cliffs (Elineau et al. 2013, 2010; Maspataud et al. 2016). In particular, the lowland area of the city (Saint-François district) is likely to be flooded by sea level and rain floods.

# 4 Data analysis

In the application of this method, $\phi_1$ and $\phi_2$ represent rain and sea level, respectively. The RV's V1 and V2 represent the rain intensity and the maximum sea level (see Fig. 4), respectively. While, D1 and D2 represent the durations of these two phenomena (rain and sea level, respectively). D1 and D2 will be considered fixed ($D1 = 1h$ and $D2 = 12h$) and will not be probabilized.

The RV's X1 and X2 describe the maximum water depths obtained by propagating the RV's rain intensity of duration 1 h and maximum sea level of duration 12 h.

The observed hourly sea levels recorded between 1971 and 2015 at the port of Le Havre are provided by the French Oceanographic Service (SHOM—Service Hydrographique et

<figure><img src="image_6.png" /><figcaption>**Fig. 4 Random variables of rainfall and sea level (V1 and V2)**</figcaption></figure>
Océanographique de la Marine). The different hourly data sets of rainfall are provided by Météo-France. The characteristics of these data sets are presented in Fig. 5.

The hourly rainfall data come from three different stations (Fig. 5). A cross-validation step is carried out with the observations of the three stations. This cross-validation consists in checking out, through certain previously identified common events, the homogeneity and consistency of these chronicles. A mixed data set with observations from the three stations can therefore be constructed to obtain a more complete data set.

A method based on a pre-threshold and a time criterion ($\Delta t$) allows to build independent samples for the RV's V1 and V2. The principle of this method is to retain only values above a threshold (high enough to guarantee the independence) and the offset time between two successive values must be greater than the time criterion $\Delta t$ (chosen according to the phenomenon).

## 4.1 Non-coincidence case

All observations of the sea level (04/01/1938–17/01/2017) and rainfall (01/01/1994–01/01/2005) data sets are used to construct the samples. 13 data sets of rainfall are extracted from the entire data period. Each data set corresponds to a different offset time. (Offset time $\delta = 0h$ corresponds to the maximum sea level moment.) The purpose of these data sets is to take into account the coincidence of rainfall and sea level, which will be detailed further in this paper.

This independent sample of sea level events (V2) is obtained by applying a pre-threshold of independence equal to 7.4 m ASL (Above Sea Level) and a time criterion $\Delta t = 1$ day. The sample obtained is composed of 1720 independent events. Applying a pre-threshold of 3 mm/h and a time criterion $\Delta t = 1$ day to rain intensity observations of duration $D1 = 1h$ (V1) gives rise to a sample of 331 independent events. The basic assumptions of independence, homogeneity and stationarity of these samples are checked using the Wald-Wolfowitz, Wilcoxon and Mann-Kendall statistical tests, respectively.

It was identified through the literature that one of the advantages of the “Peaks-Over-Threshold” (POT) method compared to the Annual Maxima (AM) method is that it improves the sample (eliminating non-extreme values) and increases its size (Coles

<figure><img src="image_7.png" /><figcaption>**Fig. 5** Rain and sea level data sets</figcaption></figure>

<figure><img src="image_8.png" /></figure>
<table><caption>**Table 1** GPD's parameters for sea level and rainfall intensity in the non-coincidence case</caption><thead><tr><th>Sample</th><th>POT threshold</th><th>σ</th><th>ξ</th><th>λ (event per year)</th><th>Sample size (years)</th></tr></thead><tbody><tr><td>Maximum sea level</td><td>8.50 m ASL</td><td>0.14</td><td>-0.03</td><td>1.47</td><td>46.80</td></tr><tr><td>Rain intensity D = 1 h</td><td>8 mm/h</td><td>3.49</td><td>0.27</td><td>5.25</td><td>8.95</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 6** GPD fit to the rain intensity $D = 1 h$ (a) and maximum sea level (b)</figcaption></figure>

2001; Hamdi et al. 2014). In the rest of this paper, only the POT method is used. This approach is based on the distribution of excesses beyond a fixed threshold and is commonly referred to as the POT approach. In the context of extreme value theory, the POT frequency model is based on the Generalized Pareto (GP) distribution. The parameters of the CDF may be estimated using the maximum likelihood method (Coles 2001).

The choice of the POT threshold is based on several criteria (Coles 2001). The stability of the quantiles and the linearity of the parameters of the distribution function as well as a sufficient number of observations above the POT threshold are the main criteria for making this choice.

The thresholds of rain intensity and sea level as well as the parameters of the GPD are given in Table 1.

It is noteworthy that the shape parameter $ξ$ of General Pareto Distribution (GPD) for the rain intensity fit ($D1 = 1 h$) is positive, while it is slightly negative for the sea level fit. This parameter governs the tail behaviour of the GPD. The right tail of the distribution is much heavier for the sea level than the rainfall (Fig. 6).

The bounds of the 70% confidence interval of the rain intensity fit are quite wide at extreme values (e.g. a 100-year return period) compared to those of the sea level fit (Fig. 6). This means that there is more inherent random uncertainty associated with rainfall data than with sea level data. Indeed, the sea level sample size is larger (47 years) than the rain sample (9 years). The presence of gaps and outliers among observations is also a source of uncertainty. The adequacy of these two fits is then checked through the Kolmogorov–Smirnov statistical test (Wang et al. 2003). It would be also interesting to take into account the seasonal variation of rainfall observations.
## 4.2 Coincidence case

Only the common observation period of the sea level and rainfall series (01/01/1994–01/01/2005) is used to construct the samples. Moreover, 13 subsets of rainfall are extracted from the raw data set. Each subset corresponds to a different offset. (Offset $δ = 0 h$ corresponds to the maximum sea level moment.) The purpose of these subsets is to take into account the coincidence of rainfall and sea level.

The sample of independent sea level maximum events (V2) is obtained by applying a pre-threshold of independence equal to 7.4 m ASL and a time criterion $Δt = 1$ day. The obtained sample is composed of 375 independent events. The same procedure is applied with the rainfall subsets. A total of 13 rain intensity samples (corresponding to 13 different offsets) are extracted in the same way as the sea level. The basic assumptions of independence, homogeneity and stationarity of these samples are checked using the same statistical tests as in precedent section.

The same GPD fitting work is done with these independent samples in order to get the marginal distributions of rainfall and sea level, whose parameters are in Table 2.

The thresholds and event rates $λ$ (Table 2) are less than those presented in Table 1, because the data sets used in the coincidence case are shorter. Moreover, the $ξ$ parameter is almost null for the majority of the samples. An exponential distribution may better fit these samples. The adequacy of all these fits (Table 2) is also checked through the Kolmogorov–Smirnov statistical test.

The correlation of the maximum sea level and the rain intensity (for each offset time) is checked, and we noticed that it is very weak. However, in this paper, we aim to apply the method with the dependence case too. Then, the results of the dependence case will be compared to those of the independence one. Many copulas are fitted and according to the AIC (Akaike Information Criterion) and the BIC (Bayesian Information Criterion); we

<table><caption>**Table 2** GP marginal distribution's parameters for sea level and rainfall intensity in the coincidence case</caption><thead><tr><th>Sample</th><th>POT threshold</th><th>σ</th><th>ξ</th><th>λ (event per year)</th><th>Effective duration (years)</th></tr></thead><tbody><tr><td>Maximum sea level</td><td>8.45 m ASL</td><td>0.11</td><td>0.03</td><td>2.62</td><td>10.67</td></tr><tr><td>Rainfall intensity</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>δ = −6 h</td><td>2.90 mm/h</td><td>2.02</td><td>−0.04</td><td>3.51</td><td>10.54</td></tr><tr><td>δ = −5 h</td><td>3.00 mm/h</td><td>2.25</td><td>0</td><td>3.79</td><td></td></tr><tr><td>δ = −4 h</td><td>3.00 mm/h</td><td>1.98</td><td>0.03</td><td>3.32</td><td></td></tr><tr><td>δ = −3 h</td><td>2.60 mm/h</td><td>1.59</td><td>0.48</td><td>4.17</td><td></td></tr><tr><td>δ = −2 h</td><td>3.00 mm/h</td><td>2.37</td><td>0.02</td><td>3.03</td><td></td></tr><tr><td>δ = −1 h</td><td>3.00 mm/h</td><td>2.12</td><td>−0.02</td><td>2.47</td><td></td></tr><tr><td>δ = 0 h</td><td>3.00 mm/h</td><td>2.00</td><td>0.015</td><td>2.37</td><td></td></tr><tr><td>δ = +1 h</td><td>2.60 mm/h</td><td>1.68</td><td>0.30</td><td>3.60</td><td></td></tr><tr><td>δ = +2 h</td><td>3.00 mm/h</td><td>3.47</td><td>−0.02</td><td>2.84</td><td></td></tr><tr><td>δ = +3 h</td><td>3.00 mm/h</td><td>2.04</td><td>0.26</td><td>3.98</td><td></td></tr><tr><td>δ = +4 h</td><td>3.00 mm/h</td><td>2.78</td><td>0.03</td><td>2.84</td><td></td></tr><tr><td>δ = +5 h</td><td>2.90 mm/h</td><td>1.80</td><td>−0.01</td><td>3.51</td><td></td></tr><tr><td>δ = +6 h</td><td>3.00 mm/h</td><td>2.29</td><td>0.03</td><td>2.37</td><td></td></tr></tbody></table>
chose a Clayton copula for the maximum sea and each rain intensity sample (as it mini-
mizes the AIC and BIC criteria for the most offsets). The fitted Clayton copula parameter
oscillates around 0.07 with a minimum value of 0.05 (for $\delta = +2h$) and a maximum value
of 0.12 (for $\delta = -4h$).

5 Results

**5.1 Scenario construction**

In this paper, a simple triangular shape is used for rain scenarios and a sinusoidal shape is
used for sea level scenarios. These forms represent more or less the observed events of rain
and sea level.

Among all the functions developed in this work, two functions allow, by giving them the
values of the intensity $V1 = i$ of rain (of duration $D1 = 1$ h) and the amplitude $V2 = NM$ of
the sea level (of duration $D2 = 12$ h), to construct triangular and sinusoidal temporal pro-
cesses of the rainfall and sea level scenarios, respectively (Fig. 7). The principle consists in
solving first-order equations:

* a linear function for rain $f(t) = at + b$, $f$ being the intensity of rain at a moment $t$, $a$ and $b$ are parameters determined by knowing the values of $f$ at instants $t = 0$ (0 mm/h), $t = D1/2$ ($I = 2i$ in Fig. 7a) and $t = D1 = 1$ h (0 mm/h);

* a cosine function for the sea level $g(t) = A \times \cos(\phi \times t)$, $g$ being the sea level at a time $t$, $A$ and $\phi$ are parameters determined by knowing the values of $g$ at instants $t = 0$ ($NM$) and $t = D2 = 12$ h ($NM$, Fig. 7b).

5.1.1 Non-coincidence case

Drawings in the CDF’s allow to get $V1$ and $V2$ for several simple scenarios of rain and sea
levels. Table 3 summarizes all the maximum sea levels (of duration $D2 = 12$ h) and the
rainfall intensities (of duration $D1 = 1$ h) of these scenarios as well as their exceedence
probabilities which come from their CDF’s $1 - F_{V2}$ and $1 - F_{V1}$, respectively. Applying the

<figure><img src="image_10.png" /><figcaption>**Fig. 7** Scenarios of rainfall (a) and sea level (b)</figcaption></figure>
<table><caption>**Table 3** Scenarios of rainfall intensity (simple triangular of duration $D1 = 1$ h) and sea level (sinusoidal of duration $D2 = 12$ h) in the non-coincidence case</caption><thead><tr><th colspan="2">Rainfall</th><th colspan="2">Sea level</th></tr><tr><th>$V1 = i$ (mm/h)</th><th>$Pr\{V1 > v1\}$</th><th>$V2 = NM$ (m ASL)</th><th>$Pr\{V2 > v2\}$</th></tr></thead><tbody><tr><td>15.27</td><td>$1.91 \times 10^{-1}$</td><td>8.55</td><td>$6.78 \times 10^{-1}$</td></tr><tr><td>26.17</td><td>$3.81 \times 10^{-2}$</td><td>8.77</td><td>$1.36 \times 10^{-1}$</td></tr><tr><td>32.38</td><td>$1.93 \times 10^{-2}$</td><td>8.86</td><td>$6.78 \times 10^{-2}$</td></tr><tr><td>39.47</td><td>$1.00 \times 10^{-2}$</td><td>8.94</td><td>$3.39 \times 10^{-2}$</td></tr><tr><td>50.71</td><td>$4.32 \times 10^{-3}$</td><td>9.06</td><td>$1.36 \times 10^{-2}$</td></tr><tr><td>60.88</td><td>$2.30 \times 10^{-3}$</td><td>9.14</td><td>$6.78 \times 10^{-3}$</td></tr><tr><td>72.76</td><td>$1.24 \times 10^{-3}$</td><td>9.22</td><td>$3.39 \times 10^{-3}$</td></tr><tr><td>80.61</td><td>$8.62 \times 10^{-4}$</td><td>9.27</td><td>$2.26 \times 10^{-3}$</td></tr><tr><td>86.63</td><td>$6.68 \times 10^{-4}$</td><td>9.31</td><td>$1.70 \times 10^{-3}$</td></tr><tr><td>91.58</td><td>$5.48 \times 10^{-4}$</td><td>9.33</td><td>$1.36 \times 10^{-3}$</td></tr><tr><td>95.81</td><td>$4.67 \times 10^{-4}$</td><td>9.35</td><td>$1.13 \times 10^{-3}$</td></tr><tr><td>99.52</td><td>$4.08 \times 10^{-4}$</td><td>9.37</td><td>$9.69 \times 10^{-4}$</td></tr><tr><td>102.84</td><td>$3.62 \times 10^{-4}$</td><td>9.39</td><td>$8.48 \times 10^{-4}$</td></tr><tr><td>105.85</td><td>$3.27 \times 10^{-4}$</td><td>9.40</td><td>$7.54 \times 10^{-4}$</td></tr><tr><td>108.61</td><td>$2.98 \times 10^{-4}$</td><td>9.41</td><td>$6.78 \times 10^{-4}$</td></tr><tr><td>119.87</td><td>$2.09 \times 10^{-4}$</td><td>9.46</td><td>$4.52 \times 10^{-4}$</td></tr><tr><td>128.51</td><td>$1.63 \times 10^{-4}$</td><td>9.49</td><td>$3.39 \times 10^{-4}$</td></tr><tr><td>141.66</td><td>$1.14 \times 10^{-4}$</td><td>9.53</td><td>$2.26 \times 10^{-4}$</td></tr><tr><td>151.74</td><td>$8.90 \times 10^{-5}$</td><td>9.57</td><td>$1.70 \times 10^{-4}$</td></tr><tr><td>160.02</td><td>$7.34 \times 10^{-5}$</td><td>9.59</td><td>$1.36 \times 10^{-4}$</td></tr><tr><td>188.55</td><td>$4.03 \times 10^{-5}$</td><td>9.67</td><td>$6.78 \times 10^{-5}$</td></tr></tbody></table>

scenario construction method as shown earlier, we obtain the scenarios presented in Fig. 7. These scenarios will then be injected one by one in the boundary conditions of the hydraulic model. This step will be shown in the next section of propagation.

### 5.1.2 Coincidence case

Combinations of 13 sea levels (sinusoidal signal form) with 14 rainfall intensities (triangular signal form) gave rise to $14 \times 13 = 182$ scenarios. Moreover, 13 equiprobable offsets $\delta = \{-6h, -5h, -4h, -3h, -2h, -1h, 0h, +1h, +2h, +3h, +4h, +5h, +6h\}$ are considered between the peak of the sea level signal ($V2$) and that of the rainfall signal ($V1$). We therefore have 182 scenarios per offset and a total of 2366 combined scenarios.

The joint exceedence probabilities of these scenarios are given by $Pr\{V1, V2\} = (1 - F_{V1}) \times (1 - F_{V2})$ if the RV's $V1$ and $V2$ are independent and by $Pr\{V1, V2\} = 1 + C(V1, V2) - F_{V1} - F_{V2}$ if $V1$ and $V2$ are dependent, where $C$ is the Clayton copula that links the CDF's $F_{V1}$ and $F_{V2}$.

Table 4 shows some combined scenarios with their joint exceedence probabilities in case of independence and dependence for the offset time $\delta = -6h$.

It is noteworthy that the case of non-coincidence is implicitly involved in the coincidence case, by considering the scenarios of rain with a null sea level and vice-versa (Table 4). In
<table><caption>**Table 4** Combined scenarios of rainfall (of duration 1 h and a triangular signal) and sea level (of duration 12 h and a sinusoidal signal) in the case of coincidence (non-exhaustive list, $\delta = -6h$)</caption><thead><tr><th colspan="2">Rainfall</th><th colspan="2">Sea level</th><th colspan="2">Joint probs</th></tr><tr><th>V1 = i (mm/h)</th><th>Pr{V1 > v1}</th><th>V2 = NM (m ASL)</th><th>Pr{V2 > v2}</th><th>Pr<sub>indep</sub></th><th>Pr<sub>dep</sub></th></tr></thead><tbody><tr><td colspan="2">δ = -6h</td><td></td><td></td><td></td><td></td></tr><tr><td>22.29</td><td>6.26 × 10<sup>-2</sup></td><td>0</td><td>1</td><td>6.26 × 10<sup>-2</sup></td><td>6.26 × 10<sup>-2</sup></td></tr><tr><td>0</td><td>1</td><td>9.14</td><td>6.65 × 10<sup>-3</sup></td><td>6.65 × 10<sup>-3</sup></td><td>6.65 × 10<sup>-3</sup></td></tr><tr><td>22.29</td><td>3.5 × 10<sup>-6</sup></td><td>9.14</td><td>3.81 × 10<sup>-3</sup></td><td>1.35 × 10<sup>-8</sup></td><td>4.69 × 10<sup>-7</sup></td></tr></tbody></table>

those scenarios when the rain (sea level) is considered alone, the probabilities are given by the non-coincidence CDF's (fitted with the complete-period samples).

## 5.2 Hydraulic model

There are a variety of hydraulic models. The simplest models are composed by a few cells linked by weir relationships and more complex 2D mesh-node models are based on the numerical resolution of second-order mathematical equations (e.g. equations of Navier Stokes, Barre de Saint-Venant).

The model of Le Havre was provided by the COmmunauté de l'Agglomération Havraise (CODAH). This 1D-2D model is built with the Infoworks-ICM software and concerns only the lower part of the city of Le Havre.

The model is composed of 76, 758 triangles and 75, 637 elements. A roughness coefficient (Manning's $n$) is equal to 0.02, triangles' minimum area is 25 m² and maximum area is 1000 m². The total area of the modelled zone is 1621.17 ha and the boundaries of this modelled zone are considered in the model as vertical walls.

Rain is taken into account by small sub-catchments defined around manholes.

The initial conditions define the initial state of the components of the hydraulic model. A wastewater profile with a daily flow of 130l/day is used. An initial elevation of 2.62 m NGF is defined for the water level in the different basins of the modelled area.

## 5.3 Point of interest

Tests on the model made it possible to note the floodplains of the lower part of the city of Le Havre. It has been found that the St-François district is among the first spots to be flooded. This district is like an island, surrounded by sea water from all its coasts. Due to the historic floods of this district (1981, 1983, 1984, 1990, 1999), the protective structures of the district were reinforced by walls (Fig. 8). For all these reasons, a point of interest was chosen in the Saint-François district (Fig. 8). This point of interest is among the lowest points of the district (altitude = 3.87 m NGF).
<figure><img src="image_11.png" /><figcaption>**Fig. 8** Point of interest in the Saint-François district (red asterisk)</figcaption></figure>

## 5.4 Propagation of scenarios

The rain and sea level scenarios presented in Table 3 (Fig. 7) are propagated. Figure 9 shows the model outputs obtained by propagating all these scenarios.

The curves in Fig. 9a show the response of the hydraulic model coming from the propagation of the rain scenarios.

As it can be noticed, water begins to accumulate at the point of interest after 30 min, that is, around the peaks of intensity of the scenarios. This means that the network's concentration time (which is the longest time a drop of water takes on the ground to reach the network outlet) is about 30 min. It should be noticed that the rain duration must be greater than or equal to this time of concentration to get overflows of the network. However, a 10-year rainfall lasting 1 h produces a nonzero response. The rain duration of 1 h is then adequate for this sewerage system.

<figure><img src="image_12.png" /><figcaption>**Fig. 9** Water depths obtained of the propagation of rain (a) and sea level (b) scenarios</figcaption></figure>

<figure><img src="image_13.png" /></figure>
Furthermore, it can be seen that for increasingly rare scenarios (Fig. 9a), the sewerage network takes a longer time to evacuate rainwater. Indeed, a 10-year rainfall (blue curve at the bottom) takes about 3 h to evacuate. This confirms that this sewerage network is correctly designed. It is however noteworthy that the form of the rainfall signal counts a lot in the response of the network; a double-triangle rainfall (conventionally used) could have taken less time to evacuate.

Beyond the 50-year rainfall scenario (3rd curve from bottom to top in Fig. 9a), water accumulates for more than 6 h at the point of interest, the sewerage network is then almost dysfunctional.

The curves in Fig. 9b show the response of the hydraulic model coming from the propagation of the sea level scenarios. The first scenario that floods the point of interest is of 700-year return period. This is due to the protection structures surrounding the Saint-François district (see Fig. 8), whose average coastline is about 5.2 m NGF.

Moreover, it can be seen in Fig. 10 that the part of the sea level which overtops the coastline of the Saint-François district only lasts about an hour and a half for the 10,000-year return period scenario and only about 40 min for the 700-year return period scenario. Indeed, the volumes of water coming from the sea levels of small return periods are low and therefore relatively quickly evacuated by the sewerage network. The scenarios with high return periods (from 2000 years up to 3000 years) cause a saturation of the network; the water thus takes much more time to evacuate (which explains the presence of bearings). The maximum water depth is reached almost at the same time as the sea level scenarios peaks (at 6 o’clock). This is due to the fact that the point of interest is close to the sea. These peaks of the water depths are closer and closer to 6 o’clock for scenarios of high return periods (10,000 years for example). Indeed, the volumes of overflows are high for these rare scenarios; this causes notably high flows which makes the sewerage network to more and more quickly saturate.

The 2366 combined scenarios are propagated in the same way as before.

## 5.5 Hazard curves

In the case of rain and sea level, the variables of interest are:

* X1 is a RV describing the maximum water depth coming from rainfall scenarios;

<figure><img src="image_14.png" /><figcaption>**Fig. 10** Saint-François district
coast line: sea level scenarios</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 11** Hazard curves of rainfall $X1 = f(\Pr\{X1 > h\})$ (a) and sea level $X2 = f(\Pr\{X2 > h\})$ (b)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 12** Water drainage direction in the 1D network: Saint-François district</figcaption></figure>

* $X2$ is a RV describing the maximum water depth coming from sea level scenarios;

* $X3$ is a RV describing the maximum water depth coming from combined scenarios (rainfall and sea level);

### 5.5.1 Non-coincidence case

The curves $X1 = f(\Pr\{X1 > h\})$ and $X2 = f(\Pr\{X2 > h\})$ are the rain and sea level hazard curves (Fig. 11).

It can be noticed that the rainfall hazard curve shape (Fig. 11a) is close to that of the GPD distribution that is used to extract the sea level scenarios (Fig. 6a). Indeed, this linearity which allowed to keep the shape of the distribution, is potentially due to the fact that the rain scenarios give rise to relatively small volumes in the Saint-François district. These volumes of water are assimilated by the sewerage system that evacuates the water towards the northern part of the district as shown in Fig. 12. This linearity is therefore due to the well functioning of the drainage network (1D). This linearity also reflects that the hypothesis made on the probabilities of rain is verified. Indeed, a rain of return period $T$ gives rise to a maximum water depth of the same return period $T$.

Unlike the rainfall, the sea level hazard curve (Fig. 11b) shape looks more or less different from that of the GPD distribution (Fig. 6b). This is due to the malfunctioning of the drainage system which is unable to properly assimilate rare sea level scenarios. This network is out of capacity because of the high volumes of water discharged into it (Fig. 9b).
<figure><img src="image_17.png" /><figcaption>**Fig. 14** Hazard curves in case of independence and dependence (the same curves) for all the offsets</figcaption></figure>

<figure><img src="image_18.png" /></figure>

Under these conditions, only the 2D part of the hydraulic model works, thus inducing this effect of nonlinearity. The assumption of conservation of the distribution for sea level scenarios is then more or less weak.

### 5.5.2 Coincidence case

The situation is different for the combined scenarios, as many combinations with different joint probabilities can give rise to a maximum water depth $X3_{\delta_i}$ exceeding the same water depth $h$. We must therefore apply Eq. 2 to sum the joint exceedence probabilities of all combinations whose maximum water depths $X3_{\delta_i}$ exceed $h$, as shown in Fig. 2. Hazard curves for an offset $\delta = -6h$ are shown in Fig. 13.

The exceedence return periods can be over-estimated when one does not consider the coincidence between the rain and the sea level. Without coincidence, a water depth of 0.6 m is exceeded on average once every 400 years by only rain (Fig. 11a) and once every 1500 years by only sea level (Fig. 11b), whereas this water depth is exceeded each 100 years in case of coincidence of these two phenomena (Fig. 13).

The same hazard curve is obtained in both cases of dependence and independence (Fig. 13). This confirms that the dependence between the two phenomena is negligible in Le Havre.

Similar contours to those of Fig. 2 are obtained for all the offsets. The same method used to construct the hazard curves of Fig. 13 is applied for the rest of the offsets. We thus obtain a total of 13 hazard curves for the independence case (Fig. 14) and 13 hazard curves for the dependence case (the same curves in Fig. 14).

It can be noticed in Fig. 14 that the hazard curve linked to $\delta = 0h$ (discontinuous black curve) is below the other curves. This shows that the offset time $\delta = 0h$, classically considered in copula approaches, can underestimate the hazard level.
The return periods are the same in both cases of independence and dependence (same curves in Fig. 14). Indeed, the dependence between the rain and the sea level is negligible in Le Havre.

There is a change in the shape of the hazard curves for the different hazard curves for a water depth of 0.70 m corresponding to a return period of about 500 years. This could be due to sea level contribution becoming considerable at this return period. Indeed, as seen in Fig. 11b corresponding to the action of only sea level, the overflows start to take place at the point of interest from a return period of 700 years. However, the change of shape on Fig. 14 is observable at return periods of less than 700 years for most offsets (at 500 per year). This means that the coincidence of a 500-year sea level (no overflows when acting alone) with a rainfall (even a frequent one) gives rise to a change in the behaviour of the drainage network. The behaviour of the network then depends on the coincidence of the flood phenomena.

The curves of Fig. 14 can be approximately enveloped by two hazard curves corresponding to $\delta = -3h$ (upper bound) and $\delta = 0h$ (lower bound). This means that the most impacting rainfall (that induces the highest water depths) occurs around $3h$ before the maximum sea level.

## 5.6 Aggregation: final hazard curve

### 5.6.1 Non-coincidence case

As it was mentioned in the last section the hazard curves presented in Fig. 11 represent $X1 = f(\text{Pr}\{X1 > h\})$ and $X2 = f(\text{Pr}\{X2 > h\})$.

The total AEP's are given by Eq. 7. The total exceedence return period is calculated by Eq. 8.

In the application of this method, $\lambda_1$ and $\lambda_2$ come from the univariate frequency analyses done in Sect. 4 (Table 1). We thus note: $\lambda_1 = 5.25$ events per year and $\lambda_2 = 1.47$ events per year.

The calculations of the total AEP's and the total exceedence return periods carried out with Eqs. 7 and 8 allowed to obtain the aggregated hazard curve $h = f(T)$ presented in black in Fig. 15.

It is noteworthy that Fig. 15 only takes into account non-coinciding rainfall and sea level events.

For return periods smaller than 2300 years, the water depths are more frequently exceeded by the rain (blue curve) than by the sea level (green curve). The aggregated hazard curve (black curve) starts superimposed and remains close to the rain curve (blue

<figure><img src="image_19.png" /><figcaption>**Fig. 15** Hazard curves of rain (blue), sea level (green) and aggregated (black) in the non-coincidence case</figcaption></figure>
curve). Indeed, a water depth of 32 cm is exceeded more frequently by the rain (every 10 years) than by the sea level (every 855 years). The total exceedence return period of a 32 cm water depth is about 10 years.

Furthermore, the return period of exceeding a 75 cm water depth is the same for both sea level and rainfall ($T = 2300$ years, overlapping point of the blue and green curves). The total return period is approximately half of the rain (or sea level) return period.

Beyond the point A (Fig. 15, beyond the return period $T = 2300$ years), water depths are exceeded much more frequently by the sea level than by the rain, the sea level hazard curve passes above that of the rainfall. Indeed, the total exceedence return period becomes closer and closer to that of the sea level.

We also can notice that the total exceedence return period (of the aggregated curve) is always close to (but still less than) the smallest return period of the two hazard curves for a same water depth.

Point A is very important in understanding the behaviour of the network. In fact, as it was shown in the Sect. 5.4, the sewerage network behaves as if it were plugged for all rainfall scenarios of return periods greater than 50 years, and this network assimilates the water coming from sea level scenarios of return periods less than a value between 2000 and 3000 years. In the area 1 in Fig. 15, water depths are more frequently exceeded by rain because the drainage network is clogged. At a 2300-year return period (point A), the behaviour of the hydraulic model is identical for both rain and sea level. However, in the area 2, the 1D network becomes clogged for the sea level scenarios (this is remarkable in Fig. 9), the water depths are then exceeded more frequently by the sea level in this area. This point A is also remarkable in Fig. 11b. This confirms the change in the network behaviour.

It is noteworthy that for a same exceedence return period, the water depth given by aggregated hazard curve (black curve) is always greater than those given by only rainfall or only sea level (e.g. for a 1000-year return period, we notice 40 cm from sea level, 70 cm from rainfall, and 80 cm from the aggregated hazard curve). Considering only one phenomenon therefore underestimates the hazard level. The same reasoning could be made on return periods. Indeed, for a same water depth, the total exceedence return period is smaller than both return periods given by only rainfall or only sea level.

### 5.6.2 Coincidence case

The total AEPs are given by Eq. 13 and the total exceedence return periods are calculated using Eq. 14.

In the application of this method, the mean annual rates of events $\lambda_{m_{\delta}}$ are calculated on the basis of the POT frequency analysis of the univariate marginals over a common observation period for the different offsets. The values are shown in Table 5.

The calculations of the total AEP's and the total exceedence return periods gave rise to the final hazard curve $h = f(T)$ in both independence and dependence cases (black curve in Fig. 16).

The result of Fig. 16 is identical for both cases of independence and dependence. This confirms once again that the dependence between the rainfall and the sea level in Le Havre can be neglected.

The aggregated hazard curve (in black in Fig. 16) is like a mean of the other hazard curves. The offset $\delta = 0h$, classically used in the literature to deal with multivariate studies, can be not satisfactory. Indeed, taking into account many offsets revealed that an offset $\delta = -3h$ can be more impacting than that of $\delta = 0h$. The aggregated hazard curve
<table><caption>**Table 5** Mean annual rates of combined events (rain, sea level) for the different offsets</caption><thead><tr><th>Rate</th><th>Value (event per year)</th></tr></thead><tbody><tr><td>λ<sub>m<sub>δ</sub>=6h</sub></td><td>3.07</td></tr><tr><td>λ<sub>m<sub>δ</sub>=5h</sub></td><td>3.21</td></tr><tr><td>λ<sub>m<sub>δ</sub>=4h</sub></td><td>2.97</td></tr><tr><td>λ<sub>m<sub>δ</sub>=3h</sub></td><td>3.40</td></tr><tr><td>λ<sub>m<sub>δ</sub>=2h</sub></td><td>2.83</td></tr><tr><td>λ<sub>m<sub>δ</sub>=1h</sub></td><td>2.54</td></tr><tr><td>λ<sub>m<sub>δ</sub>=0h</sub></td><td>2.50</td></tr><tr><td>λ<sub>m<sub>δ</sub>=1h</sub></td><td>3.11</td></tr><tr><td>λ<sub>m<sub>δ</sub>=2h</sub></td><td>2.73</td></tr><tr><td>λ<sub>m<sub>δ</sub>=3h</sub></td><td>3.30</td></tr><tr><td>λ<sub>m<sub>δ</sub>=4h</sub></td><td>2.73</td></tr><tr><td>λ<sub>m<sub>δ</sub>=5h</sub></td><td>3.07</td></tr><tr><td>λ<sub>m<sub>δ</sub>=6h</sub></td><td>2.50</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**Fig. 16** Final hazard curve (in black) of combined scenarios (same curve for both independence and dependence cases)</figcaption></figure>

combined all these offsets with a same weight (1/13) for all of them. Supposing the RV $\delta$ to be equiprobable is the most neutral way to probabilize it. In some future work, it would be interesting to associate probabilities the RV $\delta$, based on the observed offsets (i.e. empirical method). This would allow to give a more realistic and significant importance to the most impacting offset.

# 6 Discussion

Overall, the coincidence can be neglected for low return periods (< 500 years). Indeed, for these low return periods, the blue curve (Fig. 17) translates the action of only rain.

The coincidence becomes non-negligible from about a 500-year return period (Fig. 17). Beyond this return period, the contribution of the sea level becomes non-negligible. Then, the results given by the non-coincidence case can over-estimate the return periods. The coincidence of rain and sea level is thus important for these high return periods.

Interactions between rainfall and sea level are observable beyond a 500-year return period in the coincidence case (Fig. 17). In fact, a frequent rainfall can cause the saturation of the evacuation network when it is coinciding with a moderate sea level. The contribution
<figure><img src="image_21.png" /><figcaption>**Fig. 17** Comparison of final hazard curves in cases of non-coincidence (black curve) and coincidence (same curve for both independence and dependence cases)</figcaption></figure>

of sea level then becomes non-negligible; it intervenes at the outlets of the 1D network by preventing the drainage of rainwater towards the sea. Given these interactions between the two phenomena, the sea level then becomes dominant beyond a 500-year return period [while in the case of non-coincidence, rain could still be dominant up to a 2300-year return period (Fig. 15)].

The results of Fig. 17 may be biased at low return periods (< 100 years). This is because of using an extreme value distribution function to extract frequent scenarios. One can get around this limit by choosing a new distribution function or by constructing a mixed one (composed of the GPD for the extreme values and some adapted distribution function for the body of the distribution).

The final hazard curves shown in Fig. 17 depend on the POT thresholds (mentioned in Tables 1 and 2). The choice of these thresholds is a well-known difficult task in the literature because one should make the compromise between the sufficiency of the number of the exceedences (annual rates), the criteria of stability of the parameters of the GPD and on the other hand the adequacy of the GPD fit with the chosen threshold. An additional difficulty in this work is that several samples, whose size is not the same, are used. Figure 18 shows the sensitivity of the results to the POT threshold of the rainfall (non-coincidence sample). It is noteworthy that the distance between the curves becomes greater for the higher threshold values. The results are sensitive to the changes in threshold value.

# 7 Conclusion

In the present paper, we provided detailed reasoning for the need, in a PFHA framework, to combine flood phenomena at the forefront of the state of the art. Few ideas have been proposed in the literature to tackle the coincidence of flood phenomena. The present work proposes a new concept to take into account the coincidence in a dependence background and in a PFHA framework. Two cases were investigated. The first one (non-coincidence)

<figure><img src="image_22.png" /><figcaption>**Fig. 18** Sensitivity of the final hazard curves (coincidence and non-coincidence cases) to the POT threshold of the rainfall</figcaption></figure>
showed up, in a non-combinatorial way, a final (aggregated) hazard curve that synthesized the classical univariate hazard curves. The second one dealt with a more complex concept where the flood phenomena can occur separately, simultaneously, or with an offset time.

Another consideration in this paper was applying and illustrating these approaches on the example of rainfall and sea level in Le Havre, northwestern France. It may be noted that the methodology is not exemplary developed for this case study; it applies to any site likely to experience a multi-phenomena flooding.

Results in terms of hazard curves plots in both cases of coincidence and non-coincidence are examined. Overall, the application has shown that the rainfall is a dominant phenomenon at the low return periods, in Le Havre. The sea level becomes dominant at high return periods. Moreover, as shown all over the results (e.g. Fig. 13) the dependence between rain and sea level is negligible in Le Havre. On the other hand, the coincidence of rain and sea level is negligible at low return periods (< 500 years). This is due to the fact that rain is a dominant phenomenon at these low return periods. However, the interaction between the rain and the sea level (even frequent values) makes the sea level “roughly” become dominant beyond a 500-year return period (saturation of the evacuation network related to the coincidence of frequent rains with moderate sea levels). The coincidence therefore becomes important at return periods greater than 500 years.

Furthermore, the uncertainties are present in all the steps of this work. Indeed, the size of the data sets (shown in Tables 1 and 2) is not enough big to extrapolate to rare scenarios of rainfall and sea level. Future work will allow to identify and quantify the possible sources of epistemic and aleatory uncertainties. Adding to that, global warming has a direct effect on flood phenomena. It induces overall a strong rise in sea level and a change in the characteristics of storms (change in the intensity and frequency of storms). This impact will make extreme sea level and rainfall events more and more frequent. The coincidence of these two phenomena and their dependence are also likely to evolve under the impact of climate change. It is thus suggested that the question of climate change is essential to be taken into account in our proposed work. Ongoing work will allow to bring more details on the impact of climate change on coincidence and dependence of flood phenomena.

Applying the developed method on other sites could be useful to thoroughly improve the proposed method.

**Acknowledgements** The authors are grateful to the SHOM (Service Hydrographique et Océanographique de la Marine) and Météo-France for providing data. Thanks to the CODAH for their collaboration and for providing the hydraulic model of Le Havre.

## References

* Abberger K (2005) A simple graphical method to explore tail-dependence in stock-return pairs. Appl Financ Econ 15(1):43–51. https://doi.org/10.1080/0960310042000280429

* Apel H, Martinez-Trepat O, Hung NN, Chinh DT, Merz B, Dung NV (2016) Combined fluvial and pluvial urban flood hazard analysis: concept development and application to Can Tho city, Mekong Delta. Vietnam. Nat Hazards Earth Syst Sci 16(4):941–961. https://doi.org/10.5194/nhess-16-941-2016

* ASN (2013) Protection des installations nucleaires de base contre les inondations externes, guide no. 13, 40 pp. (in french). Technical report, ASN (Nuclear Safety Authority)

* Beauval C (2003) Analyse des incertitudes dans une estimation probabiliste de l’aléa sismique, exemple de la france. Ph.D. thesis, Université Joseph-Fourier - Grenoble I

* Ben Daoued A, Guimier L, Hamdi Y, Duluc CM, Rebour V (2016) Development of a probabilistic flood hazard assessment (PFHA) for the nuclear safety. In: EGU general assembly
Ben Daoued A, Hamdi Y, Mouhous-Voyneau N, Sergent P (2018) Modelling dependence and coincidence of marine flooding phenomena: methodology and simplified case study in Le Havre in France. In: EGU general assembly

Bensi M, Kanney J (2015) Development of a framework for probabilistic storm surge hazard assessment for united states nuclear power plants. In: SMiRT-23

Bernardara P, Mazas F, Kergadallan X, Hamm L (2014) A two-step framework for over-threshold modelling of environmental extremes. Nat Hazards Earth Syst Sci. https://doi.org/10.5194/nhess-14-635-2014

Bozzoni F, Corigliano M, Lai C, Salazar W, Scandella L, Zuccolo E, Latchman J, Lynch L, Robertson R (2011) Probabilistic seismic hazard assessment at the Eastern Caribbean Islands. Bull Seismol Soc Am 101:2499–2521. https://doi.org/10.1785/0120100208

Chebana F, Ouarda T (2011) Multivariate quantiles in hydrological frequency analysis. Environmetrics 22(1):63–78

Coles S (2001) An introduction to statistical modeling of extreme values. Springer series in statistics. Springer, London

DEFRA (2005) Use of joint probability methods in flood management—a guide to best practice. Technical report, DEFRA/Environment Agency

Deglaire M (2010) Le havre—etude sur le risque de submersion marine. Technical report, Université Paris 1 Panthéon-Sorbonne/DIRM-CODAH

Dung N, Merz B, Bardossy A, Apel H (2015) Handling uncertainty in bivariate quantile estimation—an application to flood hazard analysis in the Mekong Delta. J Hydrol 527:704–717

Elineau S, Duperret A, Mallet P, Caspar R (2010) Le havre : Une ville cotière soumise aux submersions marines et aux instabilites de falaises littorales. In: Journées Impacts du changement climatique sur les risques cotiers

Elineau S, Duperret A, Mallet P (2013) Coastal floods along the english channel: the study case of Le Havre town (NW France). Caribbean Waves. (poster)

Gupta I (2002) The state of the art in seismic hazard analysis. ISET J Earthq Technol 9(4):311–346

Gupta I (2007) Probabilistic seismic hazard analysis method for mapping of spectral amplitudes and other design-specific quantities to estimate the earthquake effects on man-made structures. ISET J Earthq Technol 44(1):127–167

Hamdi Y, Bardet L, Duluc CM, Rebour V (2014) Extreme storm surges: a comparative study of frequency analysis approaches. Nat Hazards Earth Syst Sci. https://doi.org/10.5194/nhess-14-2053-2014

Hawkes PJ, Gouldby BP, Tawn JA, Owen MW (2002) The joint probability of waves and water levels in coastal engineering design. J Hydraul Res 40(3):241–251. https://doi.org/10.1080/00221680209499940

Hebert H, Abadie S, Benoit M, Creach R, Duluc CM, Gailler A, Garziglia S, Lemoine A, Loevenbruck A, Macary O, Maspataud A, Mercer R, Morichon D, Pedreros R, Rebour V, Ricchiuto M, Silva Jacinto R, Terrier M, Toucanne S, Hayashi Y (2014) Le projet tandem (tsunamis en atlantique et manche : définition des effets par modelisation) (2014–2017) : enjeux pour les vulnerabilites littorales aux tsunamis

Horrillo-Caraballo J, Reeve D, Simmonds D, Pan S, Fox A, Thompson R, Hoggart S, Kwan SS, Greaves D (2013) Application of a source–pathway–receptor–consequence (SPRC) methodology to the Teign Estuary, UK. J Coastal Res 65:1939–1944. https://doi.org/10.2112/SI65-328.1

IAEA (1993) Probabilistic safety assessment for seismic events. Technical report, International Atomic Energy Agency (IAEA)

IPSN (2000) Rapport sur l'inondation du site du blayais survenue le 27 décembre 1999. Technical report, Institut de Protection et Sûreté Nucléaire (IPSN)

Klugel JU (2013) Probabilistic safety analysis of external floods—method and application. Kerntechnik 78(2):127–136. https://doi.org/10.3139/124.110335

Lang M, Rasmussen P, Oberlin G, Bobée B (1997) Echantillonnage par valeurs supérieures à un seuil: modélisation des occurrences par la méthode du renouvellement. Revue des sciences de l'eau 10(3):279–320. https://doi.org/10.7202/705281ar

Maspataud A, Elineau S, Duperret A, Ruz MH, Mallet P (2016) Impacts de niveaux d'eau extremes sur deux villes portuaires de la manche et mer du nord : Le havre et dunkerque. Journées REFMAR

Mazas F (2017) Extreme meteo-oceanic events. Theses, Université Paris-Est https://pastel.archives-ouvertes.fr/tel-01743728

Mazas F (2019) Extreme events: a framework for assessing natural hazards. Nat Hazards. https://doi.org/10.1007/s11069-019-03581-9

McGuire RK, Arabasz WJ (1990) An introduction to probabilistic seismic hazard analysis. In: Geotechnical and environmental geophysics. Society of Exploration Geophysicists, pp 333–354. https://doi.org/10.1190/1.9781560802785.ch12
Narayan S, Hanson S, Nicholls R, Clarke D (2011) Use of the source–pathway–receptor–consequence model in coastal flood risk assessment. EGU 2011, Vienna, Austria

Norberto CNC, Jeffrey AM, Victor MG, Andrew TC (2015) Coastal storm hazards from Virginia to Maine. Technical report, US Army Corps of Engineers (USACE)

Rebour V, Georgescu G, Leteinturier D, Raimond E, La Rovere S, Bernadara P, Vasseur D, Brinkman H, Groudev P, Ivanov I, Turschmann M, Sperbeck S, Potempski S, Hirata K, Kumar M (2016) Report 2: guidance document on practices to model and implement external flooding hazards in extended PSA. Technical report, IRSN-PSN-RES-SAG-PSN-RES-SAG

Salvadori G, De Michele C (2004) Frequency analysis via copulas: theoretical aspects and applications to hydrological events. Water Resour Res 40

Salvadori G, de Michele C, Durante F (2011) On the return period and design in a multivariate framework. Hydrol Earth Syst Sci 15:3293–3305. https://doi.org/10.5194/hess-15-3293-2011

Sergent P, Prevot G, Mattarolo G, Luck M, Brossard J, Nguyen T, Morel G, Mar NF, Benoit M, Ropert F, Guillou N, Bouttes F, Kergadallan X, Trichet JJ, Delisle JR, Menon JM, Mallet P, Voyneau N, Lam M, Le Banner G (2012) Projet sao polo adaptation des structures cotieres au changement climatique

Serinaldi F, Bardossy A, Kilsby CG (2015) Upper tail dependence in rainfall extremes: would we know it if we saw it? Stoch Environ Res Risk Assess 29(4):1211–1233. https://doi.org/10.1007/s0047 7-014-0946-8

Sklar M (1959) Fonctions de répartition à n dimensions et leurs marges. Publ Inst Stat Univ Paris 8:229–231

U.S. Nuclear Regulatory Commission (2014) Probabilistic flood hazard assessment research plan, ADAMS accession no. ML14296A442

U.S. Nuclear Regulatory Commission (2015) Workshop on probabilistic flood hazard assessment. Accessed 24 May 2015

Volpi E, Fiori A (2012) Design event selection in bivariate hydrological frequency analysis. Hydrol Sci J 57(8):1506–1515. https://doi.org/10.1080/02626667.2012.726357

Wang J, Wan Tsang W, Marsaglia G (2003) Evaluating kolmogorov’s distribution. J Stat Softw 8:1–4. https://doi.org/10.18637/jss.v008.i18

Yan B, Chen L (2013) Coincidence probability of precipitation for the middle route of South-to-North water transfer project in China. J Hydrol 499:19–26. https://doi.org/10.1016/j.jhydrol.2013.06.040

**Publisher’s Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.