{'fraction_non_alphanumeric': 0.027084073478088232, 'fraction_numerical': 0.01958161793901904, 'mean_word_length': 4.595416276894293, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Trivariate copula to design coastal structures

Some coastal structures must be redesigned in the future due to rising sea levels caused by global warming. The 5 design of structures subjected to the actions of waves requires an accurate estimate of the long return period of such parameters as wave height, wave period, storm surge and more specifically their joint exceedance probabilities. The Defra method that is currently used makes it possible to directly connect the joint exceedance probabilities to the product of the univariate probabilities by means of a simple factor. These schematic correlations do not, however, represent all the complexity of the reality and may lead to damaging errors in coastal structure design. 10 The aim of this paper is therefore to remedy the lack of accuracy of these current approaches. To this end, we use copula theory with a copula function that aggregates joint distribution function to its univariate margins. We select a bivariate copula that is adapted to our application by the likelihood method with a copula parameter that is obtained by the error method. In order to integrate extreme events, we also resort to the notion of tail dependence. We can select the copulas with the same tail dependence as data. In the event of an opposite tail dependence structure, we resort to the survival copula. The tail 15 dependence parameter makes it possible to estimate the optimal copula parameter. The most accurate copulas for our practical case with applications in Saint-Malo and Le Havre (France), are the Clayton normal copula and the Gumbel survival copula. The originality of this paper is the creation of a new and accurate trivariate copula. Firstly, we select the fittest bivariate copula with its parameter for the two most correlated univariate margins. Secondly, we build a trivariate function. For this purpose, we aggregate the bivariate function with the remaining univariate margin with its parameter. We show that this trivariate 20 function satisfies the mathematical properties of the copula. We can finally represent joint trivariate exceedance probabilities for a return period of 10, 100 and 1000 years.


Introduction
The design of coastal structures requires the multiplicity of variables and their degree of correlation to be taken into account. 25 We must therefore address the lack of robustness in the modelling procedure of the dependencies between the different variables characterizing the sea state (Sergent et al., 2014;Hawkes, 2005) such as wave height H, wave period T and storm surge S. The design of coastal structures is based in particular on the return periods of wave overtopping or of armour damage (Ciria et al., 2007). Since the applications on wave overtopping and armour damage depend on the parameters of the coastal structure, we will not deal with the return periods of these quantities. The aim of this paper is however to improve the methods 30 of estimating them in order to avoid costly and inappropriate decisions (Li et al., 2008). To this end, we provide accurate estimates of the correlations between the variables H, T and S and obtain reliable return periods. Currently, in reference manuals such as the Rock Manual (Ciria et al., 2007), it is recommended that a factor be applied to the product of univariate survival functions in order to determine the joint period. This method is named simplified Defra method.
Copulas are mathematical tools for modelling the dependence structure of several random variables. The theory of copulas 35 was developed by the mathematician Sklar (1959). The copula is a written form of the joint distribution function that provides all the information on the dependency structure. The recent interest in copulas started in financial risk management and insurance. Its use in environmental science especially concerns hydrology with the works for example of De Michele and Salvadori (2003), Favre et al. (2004), Grimaldi and Serinaldi (2006), Genest and Favre (2007), Zhang and Singh (2007), Aghakouchak et al. (2010), Lee et al. (2013), Chang et al. (2016). 40 In coastal engineering, in order to estimate the probability of failure of coastal or offshore structures caused in particular by the critical appearance of the combinations of parameters during a storm, Salvadori et al. (2007) use a copula in order to link the intensity of storm surge to its duration. Using the copula theory, Hawkes (2005) obtains, for example, the set of pairs of variables wave height H and surge S for a given return period. The bivariate return period can be generalized to the multivariate case (Charpentier, 2014). 45 In this paper we propose the use of copulas to take into account the dependence between three variables H, T and S. We want to show the relative importance of the choice of the copula family and of the copula construction. Tiloy et al. (2020) illustrated the importance of having a range of bivariate models when attempting to capture interrelations between pairs of hazards. In this paper we compare at the same time the choice of the copula family and the choice of the copula construction.
Copulas generally aggregate only two random variables. The purpose of this article is the creation of a new trivariate copula 50 and the evaluation of its robustness. Nelsen (1985) mentions that the construction of a trivariate copula requires a specific attention. In the literature the Chakak and Koehler (1995) method is commonly used and in particular by Joe (1997) and Salvadori et al. (2007). This method is based on bivariate conditional distributions and requires the use of three bivariate copulas. The method has a compatibility problem. There is indeed no guarantee that the method gives the same result when the order of variables is changed. Corbella and Strech (2013) study trivariate copula based on storm magnitude, storm duration 55 and wave height. They show that the fully nested method of creating hierarchical copulas provides the best results for their case study compared to Chakak and Koehler (1995) and conditional mixture. Aas and Berg (2009) propose copula construction with conditional sets : the pair copula construction (PCC). Gouldby et al. (2014) propose also a methodology for deriving extreme nearshore sea conditions for structural design with waves, winds and sea levels as offshore variables using also conditional distributions. This model is referred to as a conditional extreme model 60 in Tiloy et al. (2020). The pair copula construction (PCC) provides a powerful tool to construct flexible multivariate distributions which can be used to model complex dependencies. This method often performs better than other methods of construction of trivariate copula. Jane et al. (2020) show for example with multivariate statistics between rainfall, ocean-side water level and groundwater level, that PCC better captures the dependence than any of the five tested standard higher dimensional copulas. 65 Despite its performance, we did not use PCC for two reasons. Our objective is firstly the construction of a trivariate copula that can be easily used by coastal engineers. PCC requires complicated calculation of full conditional probabilities and the construction of vine trees. Secondly, the bivariate copulas that are selected as the most promising in our application are Archimedean copulas. If they satisfy some properties, these Archimedean copulas enable simpler methods of construction of multivariate copulas. 70 Opting to a balance between results and complexity, we propose to use a fully nested hierarchical trivariate copulas and to test the sensitivity of the results to the method of construction and to the choice of the copula. Showing that Archimedean copulas give the best results, we can indeed adopt a fully nested hierarchical copula.
In a first part, we define the theory by presenting, partly in appendix, the marginal distribution, the recommended method of the Rock Manual, the copula, the bivariate copula, the tail dependence, the survival copula, the trivariate copula and contours 75 of equal joint exceedance probability for different return periods. We obtain a bivariate copula and the copula parameter by the method of maximum likelihood and the method of the error. We show that the trivariate function that is obtained satisfies the mathematical properties of a copula.
In a second part, we present contours of equal joint exceedance probability for applications at the ports of Le Havre and Saint-Malo (Northern France) with bivariate copulas corresponding to different return periods. We select Clayton copula and survival Finally, in a third part, we apply trivariate copulas in Le Havre.

Theoretical approach
The notations and the main notions of copula for a bivariate distribution function are recalled in appendix A. In order to determine the return period of events that lead to wave overtopping or armour damages, we choose to use survival functions. 85 As mentioned by Serinaldi (2015), this option is not unique and will lead to a specific return period that he denotes TAND. For two random variables, TAND is directly related to the bivariate survival function ̅ ( , ) that is also noted ( > , > ) in appendix A.
We present here the sets of data on the sites, the selection of the best bivariate copula and the construction of trivariate copulas.

Sets of data 90
The approach is applied in two ports in Northern France, Saint-Malo and Le Havre that are presented in Figure 1. To characterize oceanic forcing, we introduce three random variables, wave height H, wave period T and storm surge S. The wave height is the significant wave height that is noted H in order to simply the notation. By convention, the random variables are written in capital letters and the realizations of these random variables are written in lowercase (h, t, s). The wave data are derived from the Anemoc numerical database, which reconstructs sea conditions over a period of nearly 25 years by hindcast.
Sea levels come from the tide gauges of various ports. The storm surges are deduced from these observations by removing 100 from the water level the value of the astronomical tide obtained using the Shom Predit software. This software also allows to obtain the density of high tides.
As the study focuses on the integration of tidal range in a macrotidal environment in the calculation of the probability of joint occurrence of waves and water levels, the used data are those of waves and surges taken at high tide. The sample is made of 105 706 events per year using the same definition as in the Rock Manual. In coastal engineering, it is customary to calculate the probability of occurrence of extreme tidal sea levels by making the product of convolution of high tide densities with the survival function of storm surges. As 706 annual tides occur, the annual probability of exceeding a given level is chosen in reference to this figure of 1/706. In addition, because coastal structures are located at shallow depths, the conditions that require the structures in terms of their stability or wave overtoppings correspond to waves occurring with high water levels. That is 110 why the privileged situations are high tides. The analysed samples are therefore pairs of storm surge associated with wave height at high tide. However, for safety, the wave height that is chosen is the maximum value over the 3-hour time interval on either side of high tide. Kergadallan (2015) recommends selecting the maximum H value within a time window centered on the time of high water. Using the same data as him, this recommendation is followed. Statistical samples therefore contain tuples of waves and storm surges values at high tides at the rate of 706 annual pairs. Separated by about 12 hours, the values 115 may not fully meet the independence test. However, this aspect is not considered here as it is often accepted in practice (see for example Hawkes (2005)).
Another approximation is the assumption of the presence of a unique wave population. This assumption is also not completely valid when we consider the wave direction of extreme events. The topic has already been discussed by Hawkes (2002) and Mazas (2017Mazas ( , 2019, among others. The treatment of wave direction can also be considered as a fourth random variable of the 120 oceanic forcing but has not been included in this work. For low and moderate values the cumulative distribution functions (CDF) FH, FT and FS are the empirical functions. For the strongest and extreme values, these cumulative distribution functions result from an adjustment of the exponential law.
Survival functions, whether for storm surge or wave height, are therefore adjusted by piece by exponentially-type analytical functions as close as possible to empirical frequencies whereas extrapolations for extreme values use exponential laws. 125 2.2 Selection of the best bivariate copula by two methods

The error method
We illustrate the method for the random variables wave height H and storm surge S. This method consists in determining the mean error e between the calculated joint cumulative distribution function (ℎ, , ) with the copula C and its parameter θ 130 and the observed joint cumulative distribution function (ℎ, ).
with n the number of pairs of values (ℎ , ).
For each copula, we first determine the parameter θ that minimizes the error e. We then select the copula with the lowest minimum mean error. 135

The maximum likelihood method
Let us call X the sample of measures (x1, x2, ..., xn) with bivariate xi = (hi, si), i = 1, .., n. The likelihood function L is defined by equation (2) : where is the probability density function of the bivariate cumulative distribution function . θ is the parameter of the copula. 140 The maximum likelihood method consists in finding the parameter θ, which maximizes the probability of obtaining the sample (Tassi, 2004). Since likelihood is a product of density we take its log-likelihood in order to facilitate calculations. We can thus work with the sum and derive it with respect to θ.
The best copula is the copula with the largest likelihood.

Construction of a trivariate copula 145
For more than two variables, C is not generally a copula (impossibility theorem of Genest (1995)). According to Nelsen (2006), it is difficult to construct n-order copulas from n-1 copulas. We present two methods for the construction of trivariate copulas.
In the first method, a trivariate copula generalizes the bivariate copula with three random variables and one parameter. In the second method, a trivariate copula associates two bivariate copulas with their two respective parameters.

Definition of a copula in dimension d > 2 150
A copula in dimension d is a distribution function on [0,1] d whose marginal laws are uniform on [0,1].

Trivariate copula with one parameter : a multi-level Archimedean trivariate 155
Since we are looking for the correlation between three variables, the first idea is to generalize the bivariate copula ( 1 , 2 ) to obtain ( 1 , 2 , 3 ). We must check that ( 1 , 2 , 3 ) is a copula, which is difficult. However Archimedean copulas like Gumbel and Clayton can be extended to an order greater than 2 using the property of Archimedean copulas (see appendix A).
For a Clayton copula of order n, this gives : For Clayton copula of order 3, it gives : 160 For Gumbel copula of order n, it gives : For Gumbel copula of order 3, it gives : By taking a single copula parameter for the three variables, we do not differentiate the two-to-two correlations of the variables even though some variables may be more correlated than others.

Trivariate copula with two parameters : a fully nested hierarchical copula 165
To better take into account the correlations of variables two by two, one option is to build trivariate functions from bivariate copulas as a fully nested hierarchical copula: Corbella (2013) tests a fully nested hierarchical copula but he uses a unique bivariate copula and does not distinguish the two bivariate copulas 1 and 2 . 1 is a bivariate copula with 1 as copula parameter. 2 is a bivariate copula with 2 as copula parameter. We must check that this function (10) is a copula and satisfies the properties of equations (4). We first aggregate 170 the two most correlated variables with the copula 2 and its copula parameter. We then add the third random variable with the copula 1 and its copula parameter. We will show later that this order provides the most robust copula.

Validity of copula properties for 2.3.3
We do not know any general methods to build high order copulas from low order copulas (Durrleman, 2010). Generally is not a copula. To prove that ( 1 , 2 , 3 ) is a copula, we must check that ( 1 , 2 , 3 ) 175 satisfies the three properties of equation (4) with d = 3, which is difficult. However Charpentier (2014) points out that C is a copula if it satisfies i) or ii). i) 1 and 2 are both Clayton or Gumbel copulas with parameters 1 for 1 and 2 for 2 positive and growing.
ii) 1 and 2 are both Archimedean copulas of respective generator 1 , 2 with 2 o 1 −1 being the inverse of a Laplace

transform. 180
For Gumbel and Clayton copulas 1 and 2 that are Archimedean copulas we check the condition (ii) that there is a function f for which the inverse Laplace transform −1 satisfies : For 1 and 2 Clayton copulas we have as the generator of 2 and as the inverse generator of 1 : This gives : 185 We can find that : With Γ( , ) the incomplete Gamma function set by for a complex with real part(a) > 0 : We conclude that there is a function f such that ϕ 2 oϕ 1 −1 = −1 [ ] : For 1 and 2 Gumbel copulas we have as generator of 2 and as the inverse generator of 1 : This gives : 190 We can find that : With Γ Gamma function, defined by :

Determination of the contour of equal joint exceedance probability
The determination of the contour of equal joint exceedance probability ( > ℎ, > , > ) consists in obtaining all the variables (H, T, S) associated with different return periods : T10 (10-year event), T100 (100-year event) and T1000 (1000-year 200 event).

Bivariate probability without tide
We deal with a set of pairs of values (h, s) that satisfy : ̅ [ ̅ , ̅ ] = 10 , 100 or 1000 ̅ is the selected bivariate survival copula. ̅ , ̅ are survival functions associated with the variables. The values f10, f100 or f1000 are the frequencies corresponding to the ten-year, hundred-year and thousand-year periods i.e 1/7060, 1/70600 and 1/706000. 205

Bivariate probability with tide
The bivariate probability with tide requires the development of the copula connecting wave height and storm surge. We can then define the joint survival function of the wave height and the storm surge. The chosen calculation method favors high tide.
The sea levels considered are therefore the sums of the astronomical high tide (generated by the attraction of the moon and the sun without weather disturbance) and the storm surges raised at the time of these astronomical high tides. This method is of 210 course valid only for macrotidal seas. The equation (23) established by Simon (1994) gives the probability that the sea level at high tide N exceeds a given value n : z is the height of the high tide, between the minimum and maximum values and respectively at high tide.
( ) is the probability that the high tide is between z and z + dz. 215 ̅ (s) is the probability of observing a storm surge S larger than s, thus ̅ (s) = P(S > s).
The bivariate survival function for wave height H and sea level N is therefore written as follows : This can be written by introducing the survival copula ̅ : The set of pairs (h, n) corresponding to the different return periods, the ten-year, hundred-year and thousand-year periods, satisfies : 220 ∫ ( ) ̅ ( ̅ (ℎ), ̅ ( − )) = 10 , 100 or 1000 It is thus possible to represent the contour of equal joint exceedance probability associated with the variables wave height and sea level .

Trivariate probability without tide
Here we have chosen the method of construction of a trivariate copula with two parameters known as fully nested hierarchical copula. We have : 225 with ̅ 1 and ̅ 2 the selected bivariate survival copula. From equations (27) and (28) we therefore obtain the equation (29) : The triplets of values (h, t, s) corresponding to the different return periods, T10 (10-year event), T100 (100-year event) and T1000 (1000-year event) satisfy : 230 ̅ 2 ( ̅ 1 ( ̅ (ℎ), ̅ ( )), ̅ ( )) = 10 , 100 or 1000 It is thus possible to represent the contours of equal joint exceedance probability associated with the variables wave height, wave period and sea level.

Trivariate joint exceedance probability with tide
The trivariate survival function for wave height H, wave period T and sea level N is written as follows: 235 This can be written by introducing the selected survival copula ̅ 2 : This expression can be written by introducing the survival copula ̅ 1 connecting F ̅ and ̅ .
The triplets of values (h, t, n) corresponding to the different return periods, T10 (10-year event), T100 (100-year event) and T1000 (1000-year event) satisfy : It is thus possible to represent the contours of equal joint exceedance probability associated with the variables wave height, 240 wave period and sea level with tide.

Tail dependence of the sample
It is necessary to treat the extreme events that are characterized by a very low occurrence. The difficulty of taking them into account is of a statistical nature: the scarcity of observations. In order to take the extreme events into account, we introduce the concept of tail dependence. For a bivariate copula, it measures the probability of simultaneous extreme realizations (Clauss, 245 2009). It describes the dependences of distribution tails for the simultaneous occurrence of extreme values. It is a highly relevant tool for the study of extreme values. We distinguish lower and upper tail dependences. They are characterized by their lower and upper tail dependence coefficients that are deduced from the following conditional probabilities, whose value is given by equations (35) and (36) that are given by (Clauss, 2009) for cumulative distribution functions 1 and 2 .
Since we fix the lower tail dependence coefficient and upper tail dependence coefficient by equations (37) and (38) : We deduce the definitions of tail dependence coefficients.

Definition:
The lower tail dependence coefficient is defined by : The copula C has a lower tail dependence if exists with Є ]0,1]. 255 If = 0 then it does not have a lower tail dependence.

Definition :
The upper tail dependence coefficient is defined by : The copula C has an upper tail dependence if exists with Є ]0,1].

If
= 0 then it does not have an upper tail dependence.
In the following section, we use survival copula ̅ and survival function ̅. Upper tail dependence and lower tail dependence 260 will be inverted.
The tail dependences of the different copulas are determined in (Nelsen, 2006) and (Roncalli, 2002) from their tail dependence coefficients. They are expressed in Table 1. Ali-Mikhail-Haq copula will be noted AMH in what follows. 265 We find that some copulas do not have lower and upper tail dependence coefficients. They cannot deal with extreme dependence. Some copulas have a lower tail dependence, others have an upper tail dependence.
The tail dependence of the sample is firstly checked. For this we graphically represent the evolution of ̅ (̅,̅) ̅ and determine its limit when ̅ tends to 0. Since ̅ is used, we note the tail dependence coefficients ̅ and ̅ . We can therefore decide whether the sample has or has not a lower or upper tail dependence. 270 In choosing the copula, it is essential to satisfy the tail dependence of the sample.
If the sample does not have a tail dependence, then the use of Gaussian copula or other copula with the same tail dependence characteristics is recommended.
If the sample has a lower tail dependence, the use of a copula with a lower tail dependence or the survival copula with an upper tail dependence is recommended. 275 If the sample has an upper tail dependence, the use of a copula with an upper tail dependence or the survival copula with a lower tail dependence is recommended.
We can also deduce the parameter of the copula from the tail dependence coefficient given by the sample.
The method that is proposed here for assessing the sample dependence refers to lower tail dependence. Other methods exist such as the chi-plot proposed by Switzer (1985, 2001) and used in coastal analyses by Mazas (2017) for instance. 280

Results for bivariate copulas
We select the most appropriate copulas at both the Le Havre and Saint-Malo (Northern France) sites using two methods. We analyze the tail dependence of the two samples. We represent the contour of equal joint exceedance probability with the selected copulas for three return periods in order to assess the relevance of the copulas.

Statistical law for adjusting wave height, wave period and storm surge 285
The representation of the contours requires knowledge of the statistical laws of adjustment of the different parameters. We therefore present these laws. For the two sites of Saint Malo and Le Havre we have used data files that provide the values for wave height, wave period and storm surge at high tide over a time period of about twenty years. The file for Le Havre site includes, for example, around 15.000 values. The wave data are extracted from the Anemoc numerical database. Sea levels at high tide are extracted from tide gauge measurements. The astronomical tide is obtained from the Shom Predit software. 290 Adjustments of the statistical laws are made according to the POT method on the basis of the exponential law. The copula parameters are calibrated from samples where wave height values less than one meter are excluded (see Figure 2), thus reducing the sample size to about 3.000 values. The copulas are fitted to all pairs/triplets of observations where the wave 295 height exceeds one meter. This threshold of one meter that is used for filtering wave height excludes the swells and leaves only a very homogeneous population of pure wind waves. This treatment removes long wave periods and increases the dependence between wave height and wave period.

Current pratice : Defra method 300
The use of the simplified Defra method in Ciria et al. (2007) is common among European coastal engineers for the study of wave overtopping or armor damages in coastal structures. It refers to the Defra method presented for example by Hawkes (2005) that is based on the Gauss copula. The simplified Defra method refers to univariate survival functions ̅ and ̅ rather than cumulative distribution functions of wave height and storm surge as coastal engineers usually work with exceedance probability rather than with non-exceedance probability. In this simplified method, the bivariate survival function is related to 305 univariate survival functions by expression (41). In France, the order of magnitude for the FD coefficient is about 20. Kergadallan (2013) recommends however a minimum value of 25.
The bivariate survival functions ̅ of table 4.15 of Rock Manual (Ciria et al., 2007) are determined with equation (41). Figure   3 shows the differences between observed bivariate survival functions and calculated bivariate survival functions using the simplified Defra method. The points of calculations in blue lie far from the first bisector in black in the figure. This shows that 310 the use of the Defra simplified method is inappropriate. This is due mostly to the use of the simplified Defra method of equation (41) but the complete Defra method with Gauss copula would not represent also perfectly the extreme events because Gauss copula has not tail dependence as we will see later. In order to improve the results we now introduce the copula theory.

Analysis of the tail dependence
The sample is analyzed in order to determine its tail dependence. This will affect the choice of copula. Since the sample has a tail dependence, it should be known whether it has a lower tail dependence or an upper tail dependence. Indeed, the result will 320 condition the choice of the copula depending on whether the sample has the same tail dependence as the copula or not. In equations (22), (26), (30), (34) the survival copula ̅ is used with survival functions. In the rest of the work we focus on this survival copula ̅ . If we kept the standard copula C, the chosen survival copulas (for example Clayton and survival Gumbel) would be said survival Clayton and Gumbel. In the two approaches, we are interested in the extreme events with large wave heights and water levels. For survival copula ̅ , we determine below its limit for survival function ̅ tending to 0.  Using the survival function ̅ , these two samples have a lower tail dependence which justifies the use of the Clayton copula.
we can determine the value of its copula parameter in Saint-Malo and Le Havre with equation (42). This copula parameter is 0.43 and 0.76 respectively.

= − ln 2 ln ̅ (42)
Note : as the Gumbel copula has an upper tail dependence, the use of its survival copula is recommended. This analysis of the sample makes it possible to understand why the Gumbel survival copula gives a minimum of error close to the minimum error 335 of the Clayton copula. We can therefore expect Gumbel survival copula results to be close to the results obtained by Clayton copula. For the set of survival copulas we determine their maximum likelihood with their parameter. We will select the survival copula that has the same tail dependence as the sample with the largest likelihood. In bold in Table 2 are presented the survival copulas with a lower tail dependence: Clayton, survival Gumbel, Student and AMH when copula parameter is close to 1. We will come back later to this special property of AMH copula. The Gauss copula has a relatively large likelihood. However, it does not 345 have a correct tail dependence and cannot therefore correctly represent the tail dependence.

Selection of the best bivariate copula for Le
For the Saint-Malo sample, we choose the Clayton copula, which has the same tail dependence as the sample, with a loglikelihood of 291 in table 2. For the Le Havre sample, we also choose the Clayton copula, which has the same tail dependence as the sample, with a log-likelihood of 387.
The Clayton copula parameters obtained by the tail dependence coefficients come close to those obtained by the log-likelihood 350 method for the Le Havre sample (3.040 values) and the Saint-Malo sample (5.888 values).
For Saint-Malo, we obtain as 0.38 the parameter of the Clayton copula using the method of maximum likelihood and 0.43 with the tail dependence coefficient.
For Le Havre, we obtain 0.74 as the parameter of the Clayton copula using the method of maximum likelihood and 0.76 with the tail dependence coefficient. 355 Even if this comparison is satisfactory, the method can be sensitive to the data and the way to determine the limit. Caillault and Guegan (2005) propose for example two methods which allow to estimate the copula characterizing the bivariate distribution function of a pair of markets. One method privileges the extreme behavior of the bivariate distribution function of the pair and the second one is based on the estimates of the copulas' parameter using a pseudo log-likelihood method. They conclude that the two approaches give different estimates of the tail dependence.
The value of the log-likelihood of the Gumbel survival copula is as large as the log-likelihood of the Clayton copula. In addition, the Gumbel survival copula has the same tail dependence as the Clayton. It is therefore as suitable as the Clayton copula.

The error method for the Clayton, Gumbel and survival Gumbel Copula 365
In order to select the most relevant copula, we represent the mean error e between the calculated survival function (ℎ, , ) with the copula C and its parameter and the measured (ℎ, ).  close to that obtained with the Clayton copula. The curve of the error obtained by the Gumbel copula survival however has a very acute minimum. Obtaining the parameter of this copula will therefore be very sensitive to the value of its minimum error.
It will therefore be necessary to determine it very precisely.
We note Emin the minimum of the mean error e and Error rate = exp (Emin) -1.  Table 3 : Emin, error rate and copula parameter for the for the different survival copulas in the ports of Le Havre and Saint Malo : Clayton, Gumbel and survival Gumbel. Table 3 is used to verify that Clayton copula is the most robust copula. It also appears that Gumbel survival copula is also an 380 appropriate option.
We have therefore shown by two methods that the Clayton copula is the most relevant for the Saint-Malo and Le Havre sites.
The parameters of the copula obtained by the error method are close to those obtained by the method of maximum likelihood for the Clayton copula.

Comparison of observed and calculated joint frequencies
In order to assess the robustness of the copulas, we show in Figure 6 the observed and calculated joint frequencies for the Le Havre sample (3.040 pairs of values). The copula represents reality more closely as the points approach the bisector y = x.
The simplified Defra method currently in use does not give a good representation of the reality of the joint frequencies for wave height and storm surge. The points obtained by this simplified Defra method are very far from the bisector.
The Clayton copula provides a good representation of the reality of joint frequencies for wave height and storm surge. The 390 points obtained by the Clayton copula come close to the bisector.
In contrast, the Gumbel copula does not give a good representation of the reality of the joint frequencies for wave height and storm surge. The points obtained by the Gumbel copula move away from the bisector. The explanation is therefore in the analysis of the sample carried out in section 3.3: we showed that the sample had a lower tail dependence whereas the Gumbel copula has an upper tail dependence. 395 The Gumbel survival copula provides a good representation of the reality of joint frequencies for wave height and storm surge.
The points obtained by the Gumbel survival copula come close to the bisector. The explanation lies in the fact of introducing the survival copula. The tail dependence of the Gumbel survival copula is opposite to the tail dependence of the Gumbel copula. We therefore reestablish a right tail dependence which gives correct results.
The results obtained by the AMH are surprisingly correct. Kumar (2010) shows that the AMH copula does not have tail 400 dependence except if the copula parameter is equal to 1. In our case, the copula parameter is close to 1. The copula seems therefore to behave like a copula with a lower tail dependence.
We show the utility of the Clayton copula in comparison with the Gumbel copula and the Defra method that is currently in use.
The results highlight the importance in copula selection of the tail dependence analysis of the sample. If the sample has a tail 405 dependence it is necessary to select a copula with the same tail dependence. The Clayton copula that has the same tail dependence as the sample gives a calculated joint frequency close to the observed joint frequency. Conversely the Gumbel copula does not correctly represent the observed joint frequency: it moves away from the bisector for the extreme points. This is because the sample has a tail dependence opposite to that of the Gumbel copula. In order to restore the proper tail dependence, we resort to the survival copula. The latter comes close the bisector but is slightly less robust than the Clayton copula. It should 410 be noted that calibration is performed on the entire sample. By truncating the sample for joint frequency values below 0.01, we would have obtained a much larger parameter for the Gumbel copula with results that are closer to measurements.  On the contrary, the joint exceedance curves obtained using the Gumbel survival copula are very similar to those obtained with Clayton. Results are therefore very sensitive to the choice of copula. A poor choice may lead to undersizing and may have economic consequences. 425 With tide the effect of storm surge on the sea level is small. The tidal range, which has an amplitude much larger than the 430 storm surge especially for the port of Saint Malo, mitigates the variations due to the storm surge. In particular, for the port of Saint-Malo, it can be seen that sea level is less sensitive to variations in the return periods than storm surge (cf. Figure 8).

Conclusion on selecting of the best bivariate copula
We selected the Clayton copula for the ports of Le Havre and Saint-Malo using three methods. In order to validate the Clayton 435 copula, we analyzed samples from 19 sites of the French coast along the Atlantic and English Channel with the maximum likelihood method. We always obtained the greatest maximum likelihood with the Clayton copula or the AMH copula (see appendix C). The sample always has a lower tail dependence (see appendix B). We can therefore conclude that the Clayton copula is the most appropriate copula for our application. For this purpose, the Table 4 gives the parameters of the different sites. Even though in some sites the AMH copula provides a larger likelihood than the Clayton copula, it should not be chosen 440 because it has a particular kind of behavior. It has a lower tail dependence if the copula parameter is 1 (or close to 1 in practice).
If the parameter is not 1, the AMH copula does not have tail dependence and its interests disappears. Since the robustness depends on the copula parameter and on the site, it cannot be recommended for a general use.  We show in Figure 9 that there is a coastal area with a maximal dependence from Concarneau to Port-Bloc (in grey in the 445 figure). There are areas that are the most exposed to wind that comes mainly from the West direction along the French Atlantic Coast.

Methodology
We have tested hierarchical construction using a fully nested hierarchical Archimedean copula. In this type of construction, we build a bivariate copula between two parameters, then we create a trivariate copula with the previous copula and the third parameter. Unlike Corbella (2013) we introduce two parameters

Construction of the best trivariate copula for the port of Le Havre 455
With the fully nested hierarchical copula method, we first determine the most appropriate copula for two parameters: (T, S), (H, T) and then (H, S). We construct the bivariate distribution function using the selected copula for the two most correlated variables. We determine the most relevant copula between the function obtained with the two most correlated variables and the third variable.

Bivariate copula for the three random variables 460
To determine the best bivariate copula we assess the maximum likelihood between ( ̅ , ̅ ), ( ̅ , ̅ ) and ( ̅ , ̅ ) with the different copulas in Table 5. For all three combinations, the Clayton copula still has the largest maximum likelihood value. In addition, we find that for the combination (H, T) the log-likelihood is significantly higher. The parameters (H, T) are therefore the most correlated parameters. It is expected because we deal with a very homogeneous population of pure wind waves as it is noticed in section 3.  Table 5 : Log-likelihood and copula parameter for the different bivariate copulas between the parameters H and S, T and S then H and T. 470

Determination of the best trivariate copula
We determine the maximum likelihood between ̅ , and ̅ with the different copulas in Table 6.
We obtain the largest log-likelihood for Clayton copula, with a parameter of 0.56, which gives: In conclusion, we have thus aggregated the most correlated H and T parameters with the best performing Clayton copula. We also used Clayton copula to aggregate ̅ , and ̅ . The aggregation requires two different parameters.  Table 6 : Log-likelihood and copula parameter for different bivariate copulas between ̅ , and ̅ . 480

Contours of equal joint exceedance probability with a trivariate copula
We represent in Figure 10 trivariate joint exceedance probability for return periods of 10, 100 and 1.000 years. The trivariate copula used is therefore constructed from a Clayton copula parameter 2.37 connecting H and T and a copula parameter 0.56 connecting FHT and FS.
In order to better visualize the incidence of return periods on trivariate joint exceedance probability, cross-sections along (H, In Figure 10a, a constant wave period is fixed corresponding to an annual return period. We show the joint exceedance probability of wave height and storm surge for three return periods of 10, 100 and 1.000 years.
In Figure 10b, a constant wave height is fixed corresponding to an annual return period. We show the joint exceedance 490 probability of the storm surge and the wave period for three return periods of 10, 100 and 1.000 years.
In Figure 10c, a constant storm surge is fixed corresponding to an annual return period. We show the joint exceedance probability of the wave height and the wave period for three return periods of 10, 100 and 1.000 years.
In the three latter figures we recognize the usual pattern and the characteristics of a strong correlation for (H, T). In Figure 10c we recognize indeed the classic pattern of contours for very dependent variables. Wave height and wave period are the most 495 correlated variables. This result is not surprising all the more since we deal with pure wind waves after we have removed the swell.
In Figure 10d, a relationship between H and T is obtained with a trivariate copula with (H,S) satisfying a joint exceedance probability of 1.000 years and with T which maximizes the trivariate joint probability density function. This relationship enables us to obtain the wave period from the wave height and the storm surge. 500

Error rate and goodness of fit for trivariate copulas
In order to show the utility of the constructed trivariate copula, we determine the error rate of the different copulas in the Le Havre area using the formula of the error given by equation (1) and the definition of the error rate given by exp(e) -1 (see Table 7).  The results obtained by the trivariate copula constructed by two bivariate copulas and two parameters are generally good.
However, by aggregating the most correlated variables first, the robustness improves as was stated by Charpentier (2014).
As expected, with one parameter Archimedean copula is less robust than fully nested hierarchical copula with two parameters.
It can also be seen that by associating the most correlated variables (H, T), the Clayton copula gives better results than the 510 Gumbel copula. For a single parameter the trivariate copula constructed with the Clayton copula is significantly more accurate than the Gumbel copula. Table 7 shows finally that the choice of the copula is much more important than the choice of the trivariate method of construction. This result validates our choice of a simple method of construction that can even lead to the most robust results according to Corbella (2013). 515  Table 8 : Goodness of fit of the different trivariate copulas for the port of Le Havre.

KHI-2 KS
The best results are obtained with two parameters. With one parameter Archimedean copula and fully nested hierarchical copula are exactly the same copula as shown in Table 8.
The results highlight the contribution of trivariate copulas constructed as a fully nested hierarchical copula with the help of 520 two Clayton bivariate copulas and two parameters by first aggregating the two most correlated parameters.

Conclusion
Wave structure designers must accurately estimate return periods of parameters such as storm surge, wave height and wave period, and more specifically, their joint probabilities of exceedance. In present practice, this joint probability of exceedance 525 is related to the product of univariate probabilities by means of a simple factor. This method can cause damaging design errors.
After highlighting the limitation of the current simplified Defra method, the theory of copula is introduced. Copulas make it possible to couple the marginal laws in order to obtain a multivariate law.
Analysis of the tail dependence of the sample is used to make an initial selection of the copulas. This is because if the sample has lower tail dependence (upper tail dependence, respectively), the copula with the same tail dependence or an inverse tail 530 dependence is chosen by taking the survival copula. The correlation between the storm surge and wave height is modelled using the Clayton copula and the survival Gumbel copula.
In order to take into account the three variables (wave height, wave period, and storm surge), we show that a fully nested hierarchical trivariate copula with two parameters is the best construction technique. This function satisfies the mathematical properties of the copulas. The error rate of 3.8 % is lower than the trivariate copula obtained by generalizing the Clayton copula 535 with a single parameter (error rate of 8.8 %). We confirm that the best results are obtained by first aggregating the most correlated variables that are here wave height and wave period. Nevertheless, the choice of method of aggregation is much less important than the choice of the copula.

Appendices
Appendix A: Outlines of copula theory A.1 Bivariate cumulative distribution function 635 We denote by the cumulative distribution function (CDF) of a random variable defined by : where P is the probability.
We must model the correlation between, for example, wave heights H and storm surges S by proposing a relation defining the 645 joint cumulative distribution function from the univariate cumulative distribution functions. We thus seek to obtain a function C which links the bivariate cumulative distribution frequency ( , ) to the univariate cumulative distribution frequencies ( ) and ( ) by integrating a correlation parameter.

A.2 Current practice in coastal engineering
The simplified Defra method that is presented for example in Ciria et al. (2007)  The dependence factor FD depends on the correlation coefficient ρ obtained from the Gaussian copula (see definition in section A.3.2). The variables X and Y for the bivariate analysis are generally wave height H and storm surge S. The dependence factor is site specific and results from the analysis of the local correlation between wave heights and storm surges.
The correspondence table between the correlation coefficient ρ and the dependence factor FD is given by Kergadallan (2013). 660 This table recommends, for example, for the North Sea, Channel and Atlantic coast the use of a minimum dependence factor FD of 25 that is a weak dependence.

A.3 Copulas
The copula is a statistical tool to characterize the dependence between several random variables where linear correlations are generally not able to represent them accurately (Charpentier, 2014). According to the latter, copulas have become an important 665 tool for modelling a multivariate law that "couples" univariate cumulative distribution functions, hence the Latin name "copula" name chosen by Sklar (1959).
If C is the copula associated with a random variable vector (X, Y) then the copula C couples the univariate cumulative distribution functions ( ) and ( ) using (A.4).
They are symmetrical copulas. They are widely used in finance. They are implicit and therefore do not have an explicit analytical form.
For more information, refer to reference books such as Joe (1997) and Nelsen (1999). The reader may also refer to Clayton (1978).