B Renard
email: renard@lyon.cemagref.fr

M Lang
email: lang@lyon.cemagref.fr

USE OF A GAUSSIAN COPULA FOR MULTIVARIATE EXTREME VALUE ANALYSIS: SOME CASE STUDIES IN HYDROLOGY

Keywords: multivariate analysis, extreme value analysis, extremes dependence, copula, field significance, regional frequency analysis, QdF models, design hydrographs, asymptotic properties, risk assessment

Risk assessment requires a description of the probabilistic properties of hydrological variables. In a number of cases, this description is made on a single variable, whereas most hydrological events are intrinsically multivariate. In this context, copulas have recently received attention in order to derive a multivariate frequency analysis. After a reminder of the general results in the field of multivariate extreme value theory, the paper gives a description of a very simple copula, the Gaussian copula. Four case studies demonstrate its usefulness in the contexts of field significance determination, regional risk analysis, Discharge-Duration-Frequency (QdF) models with design hydrograph derivation and regional frequency analysis. The limitations and potential errors related to this statistical tool are also highlighted.

I. Introduction

Extreme value theory (EVT) is widely used by hydrologists, especially for flood and drought mitigation. The severity of a hydrological extreme event is expressed as a non-exceedance probability, or equivalently, in terms of return period. This can be used to locate an observed event at a probabilistic scale and the theory can be used for extrapolation, i.e. for the determination of the probability of observing an event, even outside the range of observations. In some cases, hydrologists can be interested in several hydrological variables, which would lead to use of a multivariate statistical tool. Examples of such a situation in the flood mitigation field include the following:

• A flood event can be studied in different ways. The most commonly assessed trait is the peak flow of the event, but the volume or threshold exceedance duration can also be of interest. As an example, an event with a peak of a hundred-years return period could be less damaging than an event with a ten-years return period both in peak and volume. Unfortunately, these traits are not independent, thus preventing each variable from being studied separately. The multivariate distribution of the triplet (peak, volume, duration) is needed. Such an analysis has been undertaken by Adamson et al. [START_REF] Adamson | Bivariate extreme value distributions: an application of the Gibbs sampler to the analysis of floods[END_REF] or Grimaldi and Serinaldi [START_REF] Grimaldi | Asymmetric copula in multivariate flood frequency analysis[END_REF] as an example. • A flood event at a river confluence can be the result of high discharges in only one or both of the upstream flows. If the discharges of these two streams are independent, the probabilistic behavior of the downstream flow, equal to the sum of the two upstream ones, can be obtained by convolution. Nevertheless, both tributaries can tend to have simultaneous high flows. In this case, the preceding calculation can lead to a strong underestimation of the risk. Once again, the multivariate distribution of the two upstream flows is needed. Some examples are provided by Mousavi [START_REF] Mousavi | Composition des lois élémentaires en hydrologie régionale : application à l'étude des régimes de crue[END_REF], Le Clerc and Lang [START_REF] Clerc | Flood frequency analysis downstream confluences. Comparison between bivariate densities and experimental data[END_REF] and Favre et al. [START_REF] Favre | Multivariate hydrological frequency analysis using copulas[END_REF]. Coles [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF] describes the underlying methodology. • Water resources managers usually have to deal with several rivers in a given geographical area. A currently observed phenomenon is that a ten-year flood event is observed almost every year at a regional scale: the at-site non-exceedance probability is not suitable at the regional scale [START_REF] Troutman | Regional flood probabilities[END_REF]. This phenomenon can easily be quantified in the case of M independent sites: the probability of observing at least one event with non-exceedance probability p is equal to 1 1 M p p -≥-. In contrast, if all sites are perfectly correlated, then the regional and the local probabilities are equal. Between these two extreme cases, which are never encountered in practice, the regional behavior of flood events is related to the dependence between stations, which has to be taken into account through a multivariate distribution. Multivariate extreme value theory (MEVT) is studied by mathematicians (see Coles [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF] for a review of the topic). Unfortunately, theoretical results are often difficult to use in practice, compared to the classical EVT. Additional hypotheses therefore have to be made to take into account dependence between variables. In this context, copulas have recently received particular attention [START_REF] Caperaa | Bivariate distributions with given extreme value attractor[END_REF][START_REF] De Michele | Bivariate statistical approach to check adequacy of dam spillway[END_REF][START_REF] Favre | Multivariate hydrological frequency analysis using copulas[END_REF][START_REF] Genest | The joy of copulas: bivariate distributions with uniform marginals[END_REF][START_REF] Grimaldi | Asymmetric copula in multivariate flood frequency analysis[END_REF][START_REF] Salvadori | Frequency analysis via copulas: Theoretical aspects and applications to hydrological events[END_REF][START_REF] Song | Multivariate dispersion models generated from Gaussian copula[END_REF], because they are relatively easy to handle and they can be used in a wide range of situations. In counterpart, they are based on a model of the dependence structure, which has to be checked for adequacy. The aim of this article is to present some examples involving a particular copula, the Gaussian copula. In the first part, a review of some results of MEVT will be given, with a general description of copulas, and a more detailed description of the Gaussian copula (section II). The first application deals with the problem of the field significance determination in multi-testing problems (section III.1). The second example deals with regional risk estimation, as presented in the introduction (section III.2). In section III.3, the Discharge-Duration-Frequency (QdF) methodology will be studied, and a copula will be used to compute the return period of a design hydrograph. The problem of intersite dependence and its implications for regional frequency analysis will be explored in section III. [START_REF] Caperaa | Bivariate distributions with given extreme value attractor[END_REF]. Examples where a Gaussian copula is not able to model properly the dependence structure of observations will also be presented (section III.5). Finally, the results obtained in this paper will be discussed (section IV), before giving some conclusions about the advantages and drawbacks in using a Gaussian copula, and proposing some perspectives for improving multivariate hydrological extreme events analysis (section V).

II. Multivariate extreme value theory II.1. Some results of MEVT for componentwise maxima

Let 1 ( ,..., ) n X X be independent and identically distributed variables. Let M n denote the maximum of these variables, 1,..., max( )

n in i M X = = .
Under regularity conditions, it can be proven [START_REF] Fisher | Limiting forms of the frequency distribution of the largest or smallest member of a sample[END_REF] that the only possible limit distribution of M n , normalized by suitable values, is the general extreme value distribution (GEV), whose density is:

11 1 1( ) ( ) (; , , ) 1 e x p 1 xx fx ξ ξ ξβ ξβ αβξ αα α -  --    =- --          ≤ (1)
This theorem justifies the use of the GEV distribution to model annual maximum discharges, for example. In the multivariate case, such a theorem is also available. For clarity the bivariate case will be described. Let be independent and identically distributed vectors.

Let M

∫

and H is a distribution on [0;1] with a mean of ½.

(

) 2
Unfortunately, it appears that the limit distribution cannot be expressed in a classical parameterized way: an infinity of H distributions is able to validate the mean constraint. The choice of this function is the main difficulty of multivariate extreme value analysis. In practice, G(x,y) can be chosen in a parametric family of distributions, which will be able to model a wide range of dependence structures, ranging from independence to total dependence. An example of such a family is the logistic family [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF]. Up to this point, it is important to notice that this choice is a hypothesis which is not based on theoretical results, in contrast to the use of the GEV in the univariate case. Consequently, it is necessary to check the adequacy of the dependence model and to keep in mind that the overall uncertainty is made up of the uncertainty related to the estimation of the parameters, which can be quantified by standard statistical methods, and of the modeling uncertainty, which is more difficult to assess.

II.2. Copulas

Copulas are alternative tools for dealing with multivariate extremes, and have become very popular in recent years. Let (X (1) ,…,X (d) ) be a d-dimensional random vector with a probability distribution of:

(1) ( ) 11 ( ,..., ) Pr({ } ... { })

d dd Fx x X x X x =≤ ∩ ∩ ≤ (3) A copula is a function c verifying: 11 1 :[0;1] [0;1]
and ( ,..., ) ( ( ),..., ( ))

d dd c F x x c Fx Fx = a d (4)
The c function is used to model dependence between variables, while marginal distributions can be described in an usual way, with for example GEV distributions. Consequently, the copula function can be used in order to derive a multivariate distribution with the desired marginal distributions. Conversely, let F be a multivariate distribution with marginal distributions F 1 ,…,F d . It can then be proven that c exists [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF], which means that any multivariate distribution can be written in the form of equation ( 4). Moreover, if the marginal distributions F i are continuous, then c is unique. Nevertheless, as with the H distribution, no theoretical results are available to determine the copula: given marginal distributions, an infinity of multivariate joint distributions can be derived. Once again, a parametric family has to be chosen to model the dependence (see Favre et al. [START_REF] Favre | Multivariate hydrological frequency analysis using copulas[END_REF] for a review). As an illustration, Archimedean copulas can be constructed as follows. Let :

:[0;1] [0;1] (() ) i f (() ) ( 0 ) ( ,..., ) 0 otherwise d dd ii ii d ii c Fx Fx cx x ϕϕ ϕ ϕ - ==   ≤   =     ∑∑ a (5)
Classical choices for the ϕ generator function include ( )

1 tt α ϕ - = -(Clayton copula), 1 () l o g 1 t e t e α α ϕ  - =  -   (Frank copula) or ( ) ( log ) tt α ϕ =- (Gumbel-Hougaard copula).
The α parameter thus summarizes the dependence. Although such copulas are very useful in the bivariate case, they are more problematic with a high number of dimensions, because a single parameter is not sufficient to describe a random vector with contrasted levels of dependence between marginal components. It is thus necessary to generalize the previous generator functions, by including additional parameters. Grimaldi and Serinaldi [START_REF] Grimaldi | Asymmetric copula in multivariate flood frequency analysis[END_REF] provide an example of such a generalization. An alternative family may be easier to handle in highly dimensional cases and comprises the Elliptical copulas. These copulas use a symmetric and positive definite matrix in order to model dependence. The elements of this matrix can be interpreted as dependence measures between couples of variables, leading to an analogy with the correlations used in the case of a multivariate Gaussian distributions. The most well known elliptical copulas are the Gaussian copula, which will be described in details later, and the Student's copula. Whatever the family chosen, the number of parameters quickly grows with dimension. Estimation using standard methods such as maximum likelihood can be impossible because of numerical difficulties. A common estimation scheme consists in estimating marginal parameters separately as a first step, then estimating copula parameters knowing marginal estimations. On the other hand, the efficiency of estimators is not ensured. Although appealing, the copula theory also suffers from a number of drawbacks, essentially because of its lack of connection with standard multivariate extreme value theory. Although specific copulas have been proposed to account for the theoretical properties of extremes [START_REF] Caperaa | Bivariate distributions with given extreme value attractor[END_REF][START_REF] Genest | A characterization of Gumbel's family of extreme value distributions[END_REF], additional developments are still needed, especially in the field of spatial extremes. This problem may be a limit for extrapolation, i.e. for the computation of very low probabilities, corresponding to multivariate events lying outside the observations range. Mikosh [START_REF] Mikosch | How to model multivariate extremes if one must?[END_REF] thus lists a number of questions that remain unsolved, and provides a skeptical point of view about the use of copulas.

II.3. The Gaussian copula

The Gaussian copula is a member of the Elliptical copulas family. The dependence is modeled by means of a symmetric and definite positive matrix, whose elements are used to describe the dependence between couples of variables. This model is thus convenient when the number of dimensions is more than two or three, and has been widely studied in finance (see Cherubini et al. [START_REF] Cherubuni | Copula Methods in Finance[END_REF] and references therein). In the hydrology field, the meta-Gaussian density studied by Kelly and Krzysztofowicz [START_REF] Kelly | A bivariate meta-Gaussian density for use in hydrology[END_REF] and Herr and Krzysztofowicz [START_REF] Herr | Generic probability distribution of rainfall in space: the bivariate model[END_REF] can also be viewed as a Gaussian copula. The Gaussian copula is defined as follows:

( ) 11 11
( ,..., ) ( ),..., ( )

dd d cu u u u φφ -- =Φ ( 6
)
where φ is the cdf of the standard normal distribution N(0,1) and d Φ is the cdf of a multivariate normal distribution with mean 0 and covariance matrix Σ .

In other words, the multivariate cumulative distribution of the data will have the following form:

( )

11 1 11 11
( ,..., ) ( ( ),..., ( )) ( ( )),..., ( ( ))

dd d dd F x x c Fx Fx Fx Fx φ φ -- = =Φ d (7)
or equivalently, the multivariate density can be written as:

( ) ( ) ( ) ( ) 11
φ φ φ φ -- - -- - =× ××     Σ-     Σ-     ) (8)
The main advantage of this copula is its simplicity: once the data have been transformed by , the well known multivariate normal distribution is used to calculate probabilities. Figure 1 illustrates the principle of the copula, while Figure 2 shows the shape of the multivariate distribution obtained in a two-dimensional case. It also shows the limit of the Gaussian copula: with GEV marginal distributions, the distribution of the data has to be suitably described by this V-shaped dependence. In particular, a Gaussian copula will not be able to model properly more complex dependence structures. As an example, dependence between extreme variables may depend on marginal values. In other cases, it may be inappropriate to model the dependence of transformed variables by simple correlations. It is Mardia [40] for additional methods). The next step consists in estimating parameters. The estimation scheme presented in the previous section can be used to estimate first the marginal parameters, then the covariance matrix. Phoon et al. [START_REF] Phoon | Simulation of non-Gaussian processes using fractile correlation[END_REF] proposed an alternative method, which allows the marginal and dependence parameters to be estimated independently:

1. Estimation of Spearman's rank correlations r i,j between pairs of variables. These correlations are invariant by monotonic transformation.

Transformation toward Pearson correlation coefficients by ρ

Estimation of marginal parameters (Maximum likelihood for instance).

The drawback of this method is that the obtained Σ covariance matrix is not ensured to be non-negative definite. If this is not the case, Σ can be estimated by the classical observed covariance matrix between transformed data X % by 1 ˆT X X n Σ %% = .

To conclude with this presentation, notice that a Gaussian copula is an easy tool to simulate data with prescribed marginal distributions and correlations: as a first step, a multivariate normal data set is simulated from Φ ( ) d 0, Σ , where the elements of the diagonal of are equal to one. Marginal samples are then transformed by .

Σ 1 (( . ) ) i F φ -

III. Use of the Gaussian copula in Hydrology

III.1. Field significance

Field significance determination is a problem encountered in multi-testing studies. It has been studied by Livezey and Chen [START_REF] Livezey | Statistical field significance and its determination by Monte Carlo techniques[END_REF], Lettenmaier et al. [START_REF] Lettenmaier | Hydro-Climatological Trends in the Continental United-States, 1948-88[END_REF] , Douglas et al. [START_REF] Douglas | Trends in floods and low flows in the United States: impact of spatial correlation[END_REF], and Yue and Wang [START_REF] Yue | Regional streamflow trend detection with consideration of both temporal and spatial correlation[END_REF]. In the hydro-meteorological field, a major preoccupation is the determination of the impacts of climatic change on various variables. The method normally consists in testing a number of stations for stationarity. As an example, suppose that 100 stations are tested with a risk equal to 0.05. Under the hypothesis that all stations are stationary, approximately five stations should be detected as non-stationary, because of the 5% risk. But what would be the conclusion with 6 significant results? In other words, what is the minimum number of locally significant results to conclude, with a risk α', that all these results cannot be due to chance? In order to answer this question, the distribution of N, the number of locally significant tests under the hypothesis that all series are stationary, is needed. This distribution can be calculated theoretically as follows: let X (i) denote the data recorded at site i, i=1,…, p. Assume that all series are tested by comparing a test statistic S i =g(X (i) ) with a critical value c. Let k Ω denote the set of all possible combinations of k elements among p. The distribution of N can then be derived as:

( ) ( ) ( ) ( ) 11 1 1 1 1 () () ( ) () () () ( ) () 11
+ + + Ω Ω -- - - Ω Ω == > > ≤ ≤ => > ≤ ≤ =∈ + ∞ ∈ + ∞ ∈ -∞ ∈ -∞ =∈ ℜ ∑ ∑ ∑ XX X X XX X X X k Ω ∑ c ( 11
)
where ℜ is defined by:

k Ω ( ) { } 1 1 () () ( ) () () 1 1 1 1 ,..., (] ; [),..., (] ; [), (] ; ]),..., (] ; ]) p kk n nn n p g cg c g c g + -- - - ∈+ ∞ ∈+ ∞ ∈ -∞ ∈ -∞ (1) X= X X X X X X c (12)
This equation shows that the multivariate distribution of the tested series is needed, except in the case of independence. In this case:

( ) ( ) ( ) ( ) ( )
Ω Ω - Ω - == > > ≤ ≤ => > ≤ =- =- ∑ ∑ ∑ ≤ (13)
If all sites are independent, N follows a binomial distribution. With the preceding example, the minimum number of locally significant results to ensure a 5% field significance is thus equal to nine under the independence hypothesis, because ( 8) 0.937 pN≤ = and . ( 9) 0.972 pN≤= If the tested series are dependent, a Gaussian copula can be used to take into account spatial correlations in field significance determination, in the following way:

• parameters estimation • do i=1,…, M o simulation of a new multivariate data set thanks to the estimated copula o test of each series for stationarity o computation of N i , the number of significant results The simulated data series are all stationary, and agree with the observed correlations and marginal distributions. The sample of (N i ) can thus be used to approximate the distribution of N, and to compute a critical value. This method was applied to a set of 13 hydrometric stations in the North-East of France. Between 39 and 84 years of daily discharges are available. Annual maxima were extracted and a Mann-Kendall test [START_REF] Kendall | Rank correlation methods[END_REF][START_REF] Mann | Nonparametric tests against trend[END_REF] was applied to the 13 series obtained. Five stations showed a significant trend at risk 5% (Figure 3). The preceding procedure was then applied in order to evaluate the field significance of these five significant results. Marginal distributions are assumed to be Gumbel distributions, whose parameters are estimated with the method of moments on all available years. The Σ correlation matrix was directly estimated on transformed data by 1 ˆT n = Σ XX %% , because the Phoon method based on rank correlations led to a non-definite positive matrix. Only years shared by all stations were taken into account for this estimation (30 years). Estimated correlations ranged from 0.051 to 0.976. The adequacy of the Gaussian copula model was evaluated with the QQ-plot described in section II.3 (Figure 4), which shows no strong departure from the expected distribution. The distribution of N was finally approximated with M=1000 simulations, as shown by Figure 5. The decision depends on the combination of a local risk and a regional risk, which are not necessarily equal. Different combinations are summarized in Table 1. As an example, the five significant changes obtained with a 5% test were also regionally significant at a risk of 5%, but were not at risk 1%.

III.2. Regional risk estimation

In this section, the aim will be to estimate the probability of observing at least one event of given return period T, in a set of d stations, as explained in the introduction. The return period corresponding to this probability will be called the regional return period T R . A similar analysis can be found in the paper by Troutman and Karlinger [START_REF] Troutman | Regional flood probabilities[END_REF]. Let X (i) denote the annual maximum discharge of station i. If a Gaussian copula is suitable for modeling the multivariate data set, then this probability can be computed as follows:

( 1 )

(at least one event with return period T) 1 (no event with return period T)

1 ({ } ... { }) 1 ( ,..., ) 1(( ) , . . . ,( ) ) 1 (1 1/ ,...,1 1/ ) 1( 1 1 / ) , . dd TT d TT d Td T d PP PX q X q Fq q cF q F q cT T T φ - = - =- ≤ ∩ ∩ ≤ =- =- =-- - =-Φ - ( ) 1 .., (1 1/ ) 1/ R T T φ -- = (14)
This computation is only based on the dependence between stations, and does not depend on marginal estimations. More precisely, if the probability to compute was , then the uncertainty would have been greater, as it would include correlations and marginal parameter estimation uncertainties.

(1)

(1)

( ) ( ) ({ } ... { }) dd PX x X x ≤∩ ∩ ≤
The preceding equation was applied to the same data set as in section III.1, leading to the results presented in Figure 6. The 90% confidence interval was obtained by bootstrapping (i.e. years are bootstrapped, and the corresponding data of 13 sites are added to the new data set).

As an example, the regional return period of the event "observing at least one event of return period ten years" is approximately two years. In the same manner, the regional return period is approximately 15 years for a local return period of 100 years. It can be noticed that the confidence interval for the regional return period does not encompass the two extreme cases of independence (

13 1 1( 11 /) R T = -- T
) and total dependence (T R = T).

Moreover, from a manager's point of view, an interesting question may be: for the set of 13 stations, what are the discharges to be protected from in order to obtain a regional return period T R , if it is assumed that all local protection levels are identical? Figure 6 can be used to compute the local return period T, and the marginal estimations lead to the local quantiles with return period T. As an example, with a regional protection level T R = 10 years, the local protection level has to be around 70 years. This difference may appear quite large, but once the return period has been translated in terms of discharges thanks to the marginal estimations, this difference between regional and local risks leads to a difference ranging from 25 to 33% between the discharge with local return period 10 years and the discharge to be protected from in order to obtain a regional return period of 10 years. Finally, the preceding computations were processed by using point-estimates of the parameters. Consequently, sampling uncertainty has to be kept in mind when interpreting such results.

III.3. QdF Analysis

The Discharge-Duration-Frequency (QdF) model has been developed in order to generalize the well known Intensity-duration-Frequency (IdF) model used for precipitations analysis [19, 28-30, 44, 50]. It aims at describing the relationship between quantiles computed for mean discharges V(d,T) over a range of durations. More accurately, the model assumes that quantiles decrease as a hyperbolic function of duration: (0, ) (, )

1/

VTP VdT P d - =+ +∆ (15)
where P and are parameters which have to be estimated. ∆ One application of the QdF methodology deals with design hydrograph construction. In order to estimate the impact of floods on human activities, hydrologists are asked to provide a hydrograph with a given return period. Of course, this is an almost intractable problem, because a hydrograph can be described in a number of ways, such as peak flow, duration or volume. The QdF approach is used here to compute a Mono-Frequency Synthetic Hydrograph (MFSH), i.e. a hydrograph whose mean discharges over a given range of durations have the same return period. The method used to construct the MFSH can be found in Le Clerc [START_REF] Clerc | Revisiter la notion de scénario hydrologique de référence pour la caractérisation des inondations[END_REF]. It is based on the study of the hydrograph shapes. When the hypothesis of shape invariance can be accepted, a design hydrograph is computed using the mean of non-dimensional hydrographs, with the discharge divided by peak flow and a synchronization on peak flow. The recession is then corrected in order that each mean discharge over duration d matches with the various quantiles V(d,T) from a QdF analysis. When the various hydrographs have different shapes, Le Clerc [START_REF] Clerc | Revisiter la notion de scénario hydrologique de référence pour la caractérisation des inondations[END_REF] proposed dealing with different subsets of hydrographs, related to the various flood origins, e.g. snow melt, thunderstorms or frontal rainfall, or various spatial flood extents within the catchment. Another approach is based on the multivariate analysis of the mean discharges Vd over different durations d. The hydrograph shape invariance implies a high rank correlation between the mean discharges Vd (the event with strongest peak discharge has also the strongest discharges over durations d 1 , d 2 ,...). On the other hand, if it assumed that the discharges over different durations are not so dependent, the constructed MFSH could then have a smaller probability than that suggested by the return period used for its construction. A Gaussian copula can once again be used in order to estimate this dependency, and to derive what we will call the "multivariate MSFH return period" T M . Such an evaluation is described below by two case studies. The first deals with the Zorn river at Waltenheim (688 km²). Mean discharges were computed over durations d=1, 2, 3, 4 and 5 days, and annual maxima were then extracted on 80 years of data. Figure 7a shows the scatterplot matrix of the data: it appears that mean discharges are very strongly correlated over this duration range, which indicates that the shape invariance hypothesis is acceptable.

The Gaussian copula was then estimated with the Phoon method, and its adequacy was checked with the QQ-plot described in section II.3. The multivariate MFSH return period T M was finally calculated as follows:

(1)

(1) ( ) ( )

( 1 ) ( 1 ) () () 11 1( 1 )1 1( )1 1
Pr (all mean discharges have a return period greater than T) Pr({ } ... { }) , the notations c and d Φ have been introduced in equation ( 16). However, it should be noticed that there is no simple relationship similar to the equation 1 cc =which holds in the univariate case. Figure 7b illustrates the relationship between the return period used for MFSH construction and the multivariate MFSH return period. It shows that a hydrograph constructed with a 100-years return period has in fact a return period T M ≈ 160 years. Conversely, a hydrograph of return period T M = 100 years would be obtained from mean discharges quantiles V(d,T) with a return period T around 65 years. Translated in terms of discharge, this leads to a difference of about 7% between the mean discharge V(d,T=100 years) and the mean discharge which would lead to a hydrograph return period of 100 years, for the range of durations d considered. The constructed MFSH can thus be considered as suitable for flood scenarios. The same study was conducted on the Ubaye river at the Lauzet (946 km²), with 45 years of available data. Mean discharges are here computed over durations d = 1, 10, 20 and 30 days, because most floods are related to snow melt, leading to very slow events. The data present a weaker dependence than in the preceding case, especially between the daily duration and the others. Consequently, the ratio between marginal and multivariate return periods T and T M becomes larger: a hydrograph constructed with T=100 years has in fact a return period T M about 1000 years (Figure 8a). This result emphasizes the lack of shape invariance for snowrelated events: moderate rainfalls can be superimposed on the high baseflow created by snow melt, thus leading to a great variety of shapes at a daily resolution. By contrast, Figure 8b shows the results of the QdF analysis conducted on the same river, but only for events occurring between September and February (rainfall-related floods), with durations d=1, 2, 3 and 4 days. The results are now similar to those obtained with the Zorn River. The construction of a MFSH is thus possible for such flood events. These two examples provide an illustration of the link between the shape variability of the hydrographs and the dependence between the mean discharges Vd over different durations d.

Pr({ ( ) ( )} ... { ( ) ( )}) Pr({ ( ( )) (1 1/ )} ... { ( ( )) (1 1 dd TT dd Td d T d d Xq X q FX Fq FX Fq FX T FX φφ φφ -- - - => ∩ ∩ > => ∩ ∩ > => - ∩ ∩ > -
The Gaussian copula provides a quantification of the relationship between the return period T of each mean discharge and the multivariate return period T M, which can be very different from one catchment to another one. However, the meaning of the multivariate return period T M is far from obvious, as the sum of the probabilities (p 1 +p 2 ) can be far from one, with p 1 =Pr(hydrograph > MFSH) and p 2 =Pr(hydrograph < MFSH). A number of intermediate cases can occur, with a cross-over between an observed hydrograph and a MFSH: some mean discharges Vd have a larger return period, and others a smaller, than the prescribed return period used for the MFSH construction. Consequently, such multivariate analysis is more recommended for the study of the hydrograph shape (likelihood of various shape hypotheses) than for an exact assessment of the probability of a design hydrograph.

III.4. Regional Frequency Analysis

Regional frequency analysis is used by hydrologists in order to improve the estimation of high quantiles. The principle is to collect data originating from different locations and to derive a regional distribution of extreme streamflows or rainfalls. Since the precursory work of Dalrymple [START_REF] Dalrymple | Flood frequency analyses[END_REF], a number of methodological improvements have been proposed (e.g. [START_REF] Hosking | Regional Frequency Analysis: an approach based on L-Moments[END_REF]). However, a common drawback of almost all regional methods for extremes is that they ignore the spatial dependence of data. This problem has been addressed by Stedinger [START_REF] Stedinger | Estimating a regional flood frequency distribution[END_REF], Hosking and Wallis [START_REF] Hosking | The effect of intersite dependence on regional flood frequency analysis[END_REF][START_REF] Hosking | Regional Frequency Analysis: an approach based on L-Moments[END_REF] or Madsen and Rosbjerg [START_REF] Madsen | The partial duration series method in regional index-flood modeling[END_REF]. Roughly speaking, these authors found that ignoring intersite dependence led to underestimation of the variance of the estimates, but did not lead to any bias. In the following example, dependence will be explicitly taken into account in the probabilistic model, by means of a Gaussian copula. Six rainfall stations located around Paris, France, were used (Figure 9). 58 annual maximum values of daily rainfalls were extracted from the years between 1922 and 2003, with 1926, 1936, 1939-1944, 1948-1949, 1952, 1955-1958, 1981-1983, 1991 and 1997-2001 as missing years. The index-flood procedure of Dalrymple [START_REF] Dalrymple | Flood frequency analyses[END_REF] comprises an homogenization step, by dividing the at-site samples by the at-site means or medians. For the six study stations, the atsite median values ranged from 28.5 mm to 31.75 mm. Moreover, from a meteorological point of view, the study area can be considered as homogeneous, with little altitudinal range. Consequently, we assume that the six rainfall series arise from an identical GEV distribution. The parameters of this regional distribution are estimated using two models: In model 1 the intersite dependence is ignored. The multivariate density of a vector of six annual maxima (1) (6) ( ,..., ) tt x x at year t can be written as:

6 (1) (6) ( ) 1 1
( ,..., ) ( ; , , )

i tt t i fx x G E V xα βξ = = ∏ (17)
In model 2 the intersite dependence is modeled with a Gaussian copula, whose adequacy was checked with the Chi-square plot of section II.3. Using Equation ( 8), the multivariate density of a vector of six annual maxima (1) (6) ( ,..., ) tt x x at year t can be written as:

( )

6 1/ 2 (1) (6) ( ) 2 1 (1) (6) 1 (1) (6)
(1) 6) ;,,) and ( ; , , ) the cdf of the GEV distribution.

- = - --  =×        -      -    = ∏ ΣΣ V Σ IV V ( ) (
The dependence is thus summarized in a correlation matrix Σ. In order to obtain an acceptable number of parameters, the dependence , ij ρ between two stations is assumed to decrease as a function of the distance :

, ij d ( ) ,0 1 , , exp if 1 ij ij ii di ργ γ ρ  =-   =   j ≠ (19)
Such a parameterization can be compared to classical variograms used in geostatistics. Because of the complexity of this model, a Bayesian estimation scheme using MCMC algorithms is applied. Some details of these methods in a hydrological context can be found for instance in the papers by Perreault et al. [START_REF] Perreault | Bayesian change-point analysis in hydrometeorological time series. Part 1. The normal model revisited[END_REF][START_REF] Perreault | Bayesian change-point analysis in hydrometeorological time series. Part 2. Comparison of change-point models and forecasting[END_REF], Thyer et al. [START_REF] Thyer | Quantifying parameter uncertainty in stochastic models using the Box-Cox transformation[END_REF], Marshall et al. [START_REF] Marshall | A comparative study of Markov chain Monte Carlo methods for conceptual rainfall-runoff modeling[END_REF] and Renard et al. [START_REF] Renard | An application of Bayesian analysis and MCMC methods to the estimation of a regional trend in annual maxima[END_REF]. The likelihood of the data is simply computed from the product of the multivariate densities of observations, i.e. for one of the models k=1,2 described above: The prior distribution of the shape parameter ξ is chosen to be a Gaussian distribution with zero mean and standard deviation 0.3, which implies that the interval [-0.6;0.6] encompasses more than 95% of the density. This prior distribution can be compared with the Martins and Stedinger [START_REF] Martins | Generalized maximum-likelihood generalized extreme-value quantile estimators for hydrologic data[END_REF] geophysical prior, which is less variable and is entirely included in the interval [-0. Because almost non-informative priors are used, the posterior distribution will be mostly influenced by the likelihood, thus making the ML-estimates relevant. For the two remaining parameters of model 2, the exponential model described in equation ( 19) is fitted to the empirical correlations estimated by the method of Phoon, using a classical least-square approach. As a second step, these estimations are used as starting parameters of a Metropolis algorithm (see Renard et al. [START_REF] Renard | An application of Bayesian analysis and MCMC methods to the estimation of a regional trend in annual maxima[END_REF] for a detailed description), which was run for 100 000 iterations, and whose convergence was checked using the approach suggested by Gelman et al. [START_REF] Gelman | Bayesian data analysis[END_REF]. Finally, the last 50 000 iterations are used to perform the inference. Figure 10 describes the relationship between distance and dependence. The solid line represents the median exponential decrease described in Equation [START_REF] Galéa | Notions de base et concepts utiles pour la compréhension de la modélisation synthétique des régimes de crue des bassins versants au sens des modèles QdF[END_REF], and the dashed lines denote a 90% posterior confidence interval. The model of dependence/distance relationship seems acceptable for the six series studied. Figure 11 compares the posterior distributions of parameters α,β and ξ obtained with the two models. For the first two parameters, the results are almost identical for the mean posterior values, but the estimates of model 2, which takes into account the intersite dependence, have a larger variance. This is consistent with the conclusions of Stedinger [START_REF] Stedinger | Estimating a regional flood frequency distribution[END_REF], Hosking and Wallis [START_REF] Hosking | The effect of intersite dependence on regional flood frequency analysis[END_REF][START_REF] Hosking | Regional Frequency Analysis: an approach based on L-Moments[END_REF] and Madsen and Rosbjerg [START_REF] Madsen | The partial duration series method in regional index-flood modeling[END_REF] described earlier. Conversely, the posterior variances of the shape parameter are almost identical for the two models, with a slight shift between the two distributions. This result shows that the effect of intersite dependence is not necessarily the same for all parameters. The implications for quantiles estimates are shown in Figure 12: the difference between the two models, in terms of posterior variance, is stronger between moderate quantiles (i.e. 10years return period) than between high quantiles (i.e. 100-years return period). These results also have implications for regional methods of extrapolation. As an illustration, the FORGEX method [START_REF] Faulkner | The FORGEX method of rainfall growth estimation -III: Examples and confidence intervals[END_REF][START_REF] Reed | The FORGEX method of rainfall growth estimation -II: Description[END_REF][START_REF] Stewart | The FORGEX method of rainfall growth estimation -I: Review of requirement[END_REF] uses the series of the annual maximum of the standardized values observed over a measurement network (netmax series) in order to extrapolate the distribution of a target site. This method uses the following observation [START_REF] Reed | Inter-site and inter-duration dependence in rainfall extremes[END_REF]: if the N stations of the network are independent, then the netmax variable Y has the following distribution: In other terms:

( )

log( log( ( ))) log log ( ) log( ) log( log( ( ))) N YX X Fy Fy NF   -- = - -     =- - - , ( 23
)
so that in a Gumbel repair, the distribution of the netmax lies at a distance of log(N) to the left of the population growth curve (i.e. the distribution of the standardized values). In order to take into account the effect of intersite dependence, this distance is in fact replaced by log(N e ), where N e is an equivalent number of stations. The latter concept has been used for a long time [START_REF] Matalas | Information content on the mean[END_REF] in order to take into account the information redundancy caused by dependence. However, in the present situation, this approach suffers from two drawbacks, whatever the method used to estimate the equivalent number of stations. First of all, the preceding results show that the effect of intersite dependence is not identical for all parameters, or equivalently for the whole quantile curve. Consequently, using a unique value as N e to summarize the effect of dependence is questionable. Secondly, the parallelism between at-site and netmax growth curves is only true under the hypothesis of independence.

As an illustration, a simple bidimensional case is explored: for a network of two sites, whose dependence is described with a Gaussian copula with correlation 0.8, the at-site distribution (a GEV(1,0,-0.5)) is plotted in a Gumbel repair, together with the netmax distribution under the independence hypothesis and the "true" Gaussian copula hypothesis. Figure 13 (left panel)

shows that the at-site and the netmax distributions are parallel assuming independence, but not with the prescribed model of dependence. In the right panel, the gap between these two distributions is shown: in the independence case, the gap is constant and equal to log(2), but this not the case in the Gaussian copula case, where this gap depends on marginal values, and can therefore not be equal to log(N e ), whatever the N e used. The same computation is made for the preceding case study (Figure 14), with identical conclusions. In this case, the distributions are plotted using point-estimates of the parameters. Consequently, they are affected by sampling uncertainty, which is not shown in the figure for clarity, but which is far from negligible, thus complicating the conclusions of the comparison.

The case study could be enhanced in several ways. First of all, the assumption of a single shared distribution may be relaxed. Rather than dividing the at-site samples by an index variable, the spatial variability of the parameters could be modeled with covariables such as altitude or distance from the sea. Such an approach has been proposed by Diggle et al. [START_REF] Diggle | Model-based geostatistics[END_REF],

Cooley [START_REF] Cooley | Statistical Analysis of Extremes Motivated by Weather and Climate Studies: Applied and Theoretical Advances[END_REF] and Cooley et al. [START_REF] Cooley | A Spatial Bayesian Hierarchical Model for a Precipitation Return Levels Map[END_REF]. Moreover, from a Bayesian perspective, the use of a proper prior distribution could improve the accuracy of estimates. Bayesian model checking or Bayesian comparison of models could also be useful for deriving a more complete validation of the model assumptions. Finally, it is clear that the preceding conclusions are dependent on the Gaussian copula hypothesis, which remains questionable. Further investigations are needed to refine regional frequency analysis, using more complex dependence structures.

III.5. The Gaussian copula can fail

The need to check the adequacy of the copula before using it for any computation has been pointed out in section II. The following case study illustrates the fact that a Gaussian copula is far from being a universal tool. 55 flood events were selected from the daily discharge series of the Ubaye river at Barcelonette (549 km²;1904-2001). The scatterplot of the peak discharge of the hydrograph versus the volume computed above a threshold equal to the half of the peak discharge is shown in Figure 15a. At first sight, the dependence structure of the data seems complex. Nevertheless, the shape of the scatterplot may be due to the marginal structure of the observations: on this river, flood events can result from snow melt or heavy rainfalls. Both marginal distributions are therefore likely to derive from a population mixture model. In order to erase the influence of marginal distributions, a normal score transformation is applied to the data, leading to the scatterplot shown in Figure 15c. Snow-related events are denoted by crosses, and rainfall-related ones by triangles. Despite this transformation, the two types of events still lead to different dependence structures: for snow-related events, peak flows and volumes are more correlated than for rainfall-related ones. It is thus clear that modeling the dependence structure with a single correlation parameter will lead to very poor results. As an illustration, the contour of the estimated bivariate Gaussian density is also plotted in the Figure 15c. It clearly appears that the distribution of the points is not consistent with a Gaussian description. This departure from normality can also be viewed in the Chisquare plot of Figure 15d, where high quantiles have a tendency to deviate from the theoretical line y=x.

A more subtle problem arises from the asymptotic properties of the copula, also known as the tail dependence properties. More accurately, the Gaussian copula implies the asymptotic independence of marginal values, which means that ( )

p FX p FY p → >>   → Pr . 1 ( ) | ( ) 0 XY
In other words, the dependence weakens for very extreme events. This property has a strong influence for the computation of very low probabilities, concerning events lying outside the range of observations. More explicitly, the risk can be strongly underestimated if an asymptotically independent model is used with asymptotically dependent data. For this reason, empirical or physical evidence of tail independence is necessary before using such a copula for extrapolation. Coles et al. [START_REF] Coles | Dependence measures for extreme value analyses[END_REF] provide a number of tools for exploring the asymptotic properties of a multivariate data set.



This phenomenon will be illustrated with a simulated case study. The following variables are considered: 2 ~( 0 ; 1 ) exp( / 4) , where | ~0;

. 4

XN X YXZ Z XN  - =+   (24)
The idea behind this construction is to create a couple of variables whose dependence increases with marginal values. A sample of (X,Y) with size 10 6 was thus generated. In order to deal with known marginal distributions, data were finally transformed to fit a GEV distribution with parameters (1,0,-0.5):

{ } { } ( 11 ˆ(,) () , () XY ) X YG F X G F Y -- = %% ( 25
)
where G is the cdf of the GEV distribution and is the empirical cdf of X. ˆX F The first 100 couples were used for the copula estimation. The scatterplot of these data is shown in Figure 16a. In order to focus on the problems caused by the dependence structure, the marginal parameters were not estimated but fixed at the true value (1,0,-0.5). Figure 16b shows the scatterplot of the normalized data, and the contour of the bivariate density estimated with the Phoon method. At first sight, the Gaussian copula seems to be a reasonable description of the dependence of the data. The Gaussian copula is now used to extrapolate far outside the observations range. More accurately, the probability Pr( , ) p Xu Yu => > %% is computed, with u possibly lying outside the observations range. This probability is estimated with the Gaussian copula as ( )

) p Gu Gu φφ -- 11 ( ( )), ( ( d )

=Φ

, and is compared to the theoretical value, empirically estimated with the total sample of length 10 6 . The results are shown in Figure 17, where the probability described above is expressed in terms of return period (T=1/p). On the left-hand side, it can be observed that the Gaussian copula is adequate in the observation range. Conversely, the quantile curves differ strongly in the extrapolation range: the curve associated with the Gaussian copula grows faster than the theoretical one, which leads to a severe underestimation of risk. As an example, the quantile is under-estimated by about 38% for a return period of 10 3 years, and by about 58% for a return period of 10 4 years. This can be explained by the fact that the simulated data are not asymptotically independent. The probability of observing an extreme value in both components is thus higher than that suggested by the Gaussian copula model.

IV. Discussion

The case studies presented in sections III.1 to III.4 demonstrate the usefulness and the relative simplicity of the Gaussian copula for dealing with multivariate events. Conversely, section III.5 highlights its limitations and emphasizes the error risk related to extrapolation. As no theoretical argument can be employed to justify the use of this copula, this raises the question of how to decide, in practical cases, if this model is acceptable or should be avoided. Two different aspects have to be explored for this purpose. The first is to check that the Gaussian copula is consistent with the observed data. This implies that both the marginal distributions and the dependence structure have to be in agreement with the observations. For the marginal distributions, a number of standard univariate tools can be used (e.g. goodness-of-fit tests or QQ-plots). Moreover, univariate extreme value theory provides a strong theoretical justification for the use of extreme value distributions. After marginal transformations, the distribution is determined by a simple multivariate Gaussian model. Diagnoses for multivariate normality have also been extensively studied, thus providing a number of tools for model checking [START_REF] Mardia | Tests of univariate and multivariate normality[END_REF].

The second aspect is more difficult to handle, and is related to the asymptotic dependence properties of the data. From a theoretical point of view, it can be shown that the limit distribution of the componentwise maxima (Equation ( 2)) has the property of asymptotic dependence, except in the case of exact independence. Consequently, the Gaussian copula might be considered at first sight to be inadequate for describing these data, because it is asymptotically independent. Alternative multivariate distributions may be used [START_REF] Coles | Modelling extreme multivariate events[END_REF], e.g. using the logistic family. Such a development has been undertaken by Tawn [START_REF] Tawn | Modelling multivariate extreme value distributions[END_REF] and Schalter and Tawn [START_REF] Schalther | A dependence measure for multivariate and spatial extreme values: Properties and inference[END_REF] for an arbitrary number of dimensions, and may be applied with additional hypotheses in order to reduce the number of parameters (2 d -1). However, this limitation has to be mitigated, because the limit distribution of Equation ( 2) is valid when the block size tends toward infinity. In most hydrological applications, the block size is typically of 365 days, with non-independent daily data, and a number of zero values in the case of rainfall. It is thus possible that the asymptotic arguments used to derive the limit distribution of componentwise maxima do not hold in real-life hydrological studies. Although the latter argument can also be argued for univariate extreme value distributions, it seems that the problem is more critical in the multidimensional case. For instance, Bortot et al. [START_REF] Bortot | The multivariate Gaussian tail model: an application to oceanographic data[END_REF] reported empirical evidence of asymptotically independent data in the oceanographic field. Alternatively, Hosking and Wallis [START_REF] Hosking | The effect of intersite dependence on regional flood frequency analysis[END_REF] argued that a Gaussian copula model was fairly well supported by British annual flood series, although in some cases, stronger dependence seemed to occur at higher levels.

From a practical point of view, a given data set of componentwise maxima will generally not include enough data to provide empirical evidence of asymptotic dependence or independence. Physical arguments may be used to discriminate between these two hypotheses, by using for example meteorological model simulations, or by exploiting the knowledge of historical extreme events. Consequently, our point of view is that the Gaussian copula should not be used for the computation of very low probabilities, unless strong arguments can be derived to justify the asymptotic properties of the data. Conversely, the Gaussian copula can be very useful within the observations range. As an illustration, in the case study in section III.4, the copula can be used to take into account intersite dependence in regional frequency analysis, and to derive a more realistic quantification of the uncertainties than with the independence hypothesis, which is in any case a very crude model for describing dependent data.

V. Conclusion and perspectives

The aim of this article was to present some possible applications of the Gaussian copula in the flood mitigation field. The main advantage of this method is its simplicity, even with more than two or three dimensions. Rough estimates can thus be obtained in problems involving an evaluation of the dependence of variables. Moreover, the Gaussian copula can be used to simulate values with prescribed correlations and marginal distributions. Four case studies were used to demonstrate its usefulness in the contexts of field significance determination, regional risk analysis, QdF models with design hydrograph derivation and regional frequency analysis. Nevertheless, this tool is far from universal, as it is based on no theoretical justification. The suitability of the Gaussian dependence model therefore has to be properly checked. Moreover, even if the copula seems well fitted to observed values, caution is needed before using it for extrapolation.

Alternative copula families can be used in order to improve the dependence model arising from the Gaussian copula. For instance, the Student copula may be more appropriate if tail dependence is observed in the data. Nevertheless, the choice of the more suitable copula for modeling a data set is not straightforward. Various criteria may be used for this purpose, but they will only reflect a particular feature of the data. For instance, two models may lead to an adequate fit to the data on the basis of a given criterion, but to contrasted results in extrapolation. The choice of a copula in such a case is problematic, especially if the data set is not informative enough to provide relevant indications about the asymptotic dependence properties. Modeling uncertainty can thus be an important part of the overall uncertainty. Another problem concerns the parameter estimation uncertainties, because the maximum likelihood method can usually not be used, due to the great number of parameters to be estimated. Computations based on the Fisher information matrix are thus impossible. Finally, whatever the copula used, the lack of theoretical justification remains a limitation for extrapolation.

The preceding approaches only consist in modeling the extreme part of the variables distribution. An alternative consists in studying the whole distribution, and evaluating how dependence evolves for high values. For instance, Bortot et al. [START_REF] Bortot | The multivariate Gaussian tail model: an application to oceanographic data[END_REF] applied a model where the tail dependence is described using a multivariate Gaussian distribution, after suitable marginal transformations. In spirit, this model is similar to that of the Gaussian copula. The major difference is that the data are not block maxima, and that the Gaussian tail model is only applied above a multivariate threshold , which has to be high enough to ensure that the asymptotic assumptions hold, both for marginal and joint distributions. Moreover, given that the distribution belongs to the domain of attraction of a multivariate extreme value distribution, some theoretical results related to the theory of regular variation [START_REF] Bingham | Regular variation[END_REF] can be used in order to compute the probability of extreme sets. An example of such study is provided by De Haan and De Ronde [START_REF] De Haan | Sea and Wind: Multivariate Extremes at Work[END_REF]. Such an approach is appealing, because it is physically more convincing: couples of values are describing the same physical event, while componentwise maxima lead to the use of couples of values which may describe two distinct events. Nevertheless, some difficulties still hold. Firstly, obtaining a sample of independent values is not obvious. As an example, if the data consist of daily discharges at two locations, a sampling strategy has to be found which ensures the independence of successive values and preserve enough data for estimation. Secondly, this approach can only be used to compute the probability of sets for which both components are extreme. This can be a strong limitation: in the confluence problem, as an example, an extreme downstream discharge can result from a high value in only one of the two upstream flows. Heffernan and Tawn [START_REF] Heffernan | A conditional approach for multivariate extreme values[END_REF] recently proposed a semi-parametric conditional approach which overcomes this difficulty. More generally, we believe that multivariate extreme events analysis is likely to improve hydrological risks assessment, as emphasized by Katz et al. [START_REF] Katz | Statistics of extremes in hydrology[END_REF]. In a number of situations, multivariate events are used without a formal statistical model to account for dependence. As an example, the well-known index flood methodology [START_REF] Dalrymple | Flood frequency analyses[END_REF] is intended to improve quantiles estimations by using data from several sites. Unfortunately, among other problems, ignoring spatial dependence leads to an underestimate of the quantiles uncertainty. Alternatively, most hydrological extreme events (floods or droughts) are intrinsically multivariate, as noted by Adamson et al. [START_REF] Adamson | Bivariate extreme value distributions: an application of the Gibbs sampler to the analysis of floods[END_REF]. Potential damage is thus likely to be a function of several random variables. Hydrologists involved in risk assessment should thus be attentive with the progress achieved by statisticians in multivariate extremes theory.

LIST OF CAPTIONS TABLES Table 1: Comparison of observed and critical number of significant tests for various local and regional risks. A star denotes a regionally significant result.

FIGURES

Figure 1. Principle of the Gaussian copula Figure 2. Isolines of a bivariate distribution arising from a Gaussian copula with marginal distributions (1, 0, 0.2) GEV and correlation coefficient equal to 0.8. Figure 3. Annual maxima of 13 stations in the north-east of France. Stations with significant trend at risk 5% are denoted with a star in their figure box. Figure 4. Chi-square plot (see text for explanation) for the 13 stations studied. Figure 5. Estimated distribution of N, the number of locally significant results, under the assumption that all series are stationary. Figure 6. Local versus regional return period (see text for definition). 90% confidence intervals are obtained by Bootstrap. Figure 7. QdF Analysis of the Zorn River. (a) Scatterplot matrix of annual maxima, (b) return period of the mean discharges Vd versus multivariate hydrograph return period. Total dependence and independence cases are respectively denoted by dash and dash-dot lines. Figure 8. QdF Analysis of the Ubaye River at Lauzet. (a) all events; (b) events occurring between September and February. Total dependence and independence cases are respectively denoted by dash and dash-dot lines. Figure 9. Location of the six rainfall stations. Figure 10. Relationship between intersite distance and intersite dependence. Crosses denote Phoon's estimates of correlations, the median line with 90% confidence intervals are obtained by means of the posterior distribution of parameters. Figure 11. Posterior marginal distributions of scale, location and shape parameters. Dashed line: independence hypothesis, solid line: Gaussian copula. Figure 12: Posterior distribution of quantiles with probabilities 0.9 and 0.99. Dashed line: independence hypothesis, solid line: Gaussian copula. Figure 13. Bidimensional case: Comparison between at-site distribution (triangles) and netmax distribution, in the cases of independence (crosses) and with a Gaussian copula (circles). Left: distributions in a Gumbel repair, right: difference with at-site distribution. Figure 14. Network of six rainfall series. Comparison between at-site distribution (triangles) and netmax distribution, in the cases of independence (crosses) and with a Gaussian copula (circles).

  [0;1] [0; ] ϕ →+ ∞ be a continuous decreasing function such that (1+∞ and k=1,…,d, where d is the number of dimensions. The following function is then a copula:

k

  5;0.5]. For other parameters, almost non-informative priors are set, by using uniform distributions with large variances, namely ~U[0,1000] the independence between these prior marginal distributions is assumed, in order to derive the multivariate prior distribution of the parameters of the model k.The posterior distribution is finally obtained up to a constant of proportionality by: Samples arising from the posterior distribution of the parameters are generated with MCMC algorithms in the following way. First of all, at-site maximum likelihood estimations are used to derive a rough estimate of the posterior means and variances of the parameters α, β and ξ .

1 (

 1

Figure 15 .

 15 Figure1. Principle of the Gaussian copula Figure2. Isolines of a bivariate distribution arising from a Gaussian copula with marginal distributions (1, 0, 0.2) GEV and correlation coefficient equal to 0.8. Figure3. Annual maxima of 13 stations in the north-east of France. Stations with significant trend at risk 5% are denoted with a star in their figure box. Figure4. Chi-square plot (see text for explanation) for the 13 stations studied. Figure5. Estimated distribution of N, the number of locally significant results, under the assumption that all series are stationary. Figure6. Local versus regional return period (see text for definition). 90% confidence intervals are obtained by Bootstrap. Figure7. QdF Analysis of the Zorn River. (a) Scatterplot matrix of annual maxima, (b) return period of the mean discharges Vd versus multivariate hydrograph return period. Total dependence and independence cases are respectively denoted by dash and dash-dot lines. Figure8. QdF Analysis of the Ubaye River at Lauzet. (a) all events; (b) events occurring between September and February. Total dependence and independence cases are respectively denoted by dash and dash-dot lines. Figure9. Location of the six rainfall stations. Figure10. Relationship between intersite distance and intersite dependence. Crosses denote Phoon's estimates of correlations, the median line with 90% confidence intervals are obtained by means of the posterior distribution of parameters. Figure11. Posterior marginal distributions of scale, location and shape parameters. Dashed line: independence hypothesis, solid line: Gaussian copula. Figure12: Posterior distribution of quantiles with probabilities 0.9 and 0.99. Dashed line: independence hypothesis, solid line: Gaussian copula. Figure13. Bidimensional case: Comparison between at-site distribution (triangles) and netmax distribution, in the cases of independence (crosses) and with a Gaussian copula (circles). Left: distributions in a Gumbel repair, right: difference with at-site distribution. Figure14. Network of six rainfall series. Comparison between at-site distribution (triangles) and netmax distribution, in the cases of independence (crosses) and with a Gaussian copula (circles). Figure 15. (a) Scatterplot of peak daily discharge versus volume of selected hydrographs on the Ubaye River at Barcelonette. (b) Scatterplot of peak discharge versus volume, transformed

Acknowledgements

The financial support provided by Cemagref and EDF for the PhD research of B. Renard is gratefully acknowledged. Discharge data series are derived from the national HYDRO database from the French Ministry of Environment, and the rainfall series have been reviewed by Météo-France within the European IMFREX project. Both organisms are acknowledged. We also thank two anonymous reviewers for their helpful comments.


O Souhar
email: drsouhar@hotmail.com

J-B Faure
email: faure@lyon.cemagref.fr

A Paquier
email: paquier@lyon.cemagref.fr

Automatic sensitivity analysis of a finite volume model for two-dimensional shallow water flows

Keywords: Automatic differentiation, Hydraulic model, Sensitivity analysis, Uncertainty propagation, Steady flow. O. Souhar et al

Given a numerical model for solving two-dimensional Shallow Water Equations, we are interested in the robustness of the simulation by identifying the rate of change of the water depths and discharges with respect to a change in the bottom friction coefficients. Such a sensitivity analysis can be carried out by computing the corresponding derivatives. Automatic Differentiation (AD) is an efficient numerical method, free of approximation errors, to evaluate derivatives of the objective function specified by the computer program, Rubar20 for example. In this paper AD software tool Tapenade is used to compute forward derivatives. Numerical tests were done to show the robustness of the model and to demonstrate the efficiency of these AD-derivatives.

Introduction

The two-dimensional Shallow Water Equations are a relevant basis for modelling free surface flows, but their solution requires detailed information about the river configuration as the topographical data and the friction coefficient according to Manning law. It is also necessary to provide information bound to a particular event: these are the initial conditions of the flow as well as the upstream and downstream boundary conditions of the river. The model input data are submitted to uncertainties difficult to evaluate and which may lead to meaningful uncertainties on the model output. The output results may be a map showing a flooded urbanized area; consequently, the induced damages may be very high. In these situations, the river can go far away from their mean water channel, submerging large areas, eventually with high velocity. In order to evaluate the uncertainties of the water level and river flow width according to the input data of a one-dimensional hydraulic model (Giraud et al. [START_REF] Giraud | Hydrologic modeling of a complex wetland[END_REF]), Dhervillez ([5]) and Faure et al. ( [START_REF] Faure | Automatic evaluation of the uncertainty on a flooded area[END_REF]) tested several noising methods. Their study showed that if the bottom level is gaussian then consequently the water level is gaussian too in steady flow. On the other hand, they showed in the same paper that if the bottom level is gaussian then the river flow width is not gaussian. In such a case, sensitivity analysis can be lead by computing partial derivatives of the output functions with respect to the corresponding parameters. The most familiar schemes to approximate derivatives of functions specified by computer programs are the central differences and the so-called forward finite differences. So for a smooth function F : IR n -→ IR m and for some direction d ∈ IR n , the central approximation, for ε << 1, is defined as

∇F (x).d ≈ F (x + ε.d) -F (x -ε.d) 2.ε , (1.1)
which have a truncation error of O((ε.d) 2 ). Thus, using this central operator, the gradient requires at least (2n) evaluations of the function F . However, the forward operator defined below requires only (n + 1) evaluations of F , but the truncation error is only of

O(ε.d): ∇F (x).d ≈ F (x + ε.d) -F (x) ε . (1.2)
From a theoretical point of view, a good accuracy of the approximate (1.1) or (1.2) is obtained for very small values of ε. Nevertheless, ε may not be arbitrarily small, since for very small values, the difference between F (x + ε.d) and F (x) is of the discretization error order, thus losing its significance. Nowadays, methods and softwares for sensitivity analysis of several engineering problems (fluid dynamics, chemistry, meteorology or environment studies etc) have demonstrated that forward sensitivities can be computed reliably and efficiently via Automatic Differentiation (AD). AD is superior to divided differences because AD-generated derivative values are free of approximation errors, see Section 4 for more details. In this paper the Tapenade software tool, developed by Tropics team INRIA Sophia-Antipolis, is used to compute the desired derivatives of the objective functions specified by the finite volume two-dimensional shallow water model Rubar20 with respect to the roughness coefficients. The output results or the objective functions are the water depth and the discharge computed by Rubar20 (Paquier [START_REF] Paquier | Modélisation et simulation de la propagation de l'onde de rupture de barrage[END_REF]). The outline of this paper is as follows. The governing equations are first presented. Details of the numerical model Rubar20 are discussed. An overview of the AD strategy is given. The numerical results are presented and the computed results are verified by comparing them with other results obtained from another hydraulic model. Finally, the sensitivity of the model according to the roughness coefficient is studied.

Shallow Water Equations

The two-dimensional shallow water equations describe the unsteady free surface flows, which simply express the conservation of mass and momentum. Usually these equations are described in terms of water depth and discharge under the following assumptions (Lai [START_REF] Lai | Numerical modelling of unsteady open-channel flow[END_REF]): hydrostatic pressure distribution, very small bottom slopes, uniform velocity in the vertical direction and incompressible fluid. The resulting equations are

∂h ∂t + ∂q x ∂x + ∂q y ∂y = 0 (2.1) ∂qx ∂t + ∂( q 2 x h +g h 2 2 ) ∂x + ∂( qxqy h ) ∂y = -gh ∂Z ∂x -gq x √ q 2 x +q 2 y C 2 h 2 + D( ∂ ∂x (h ∂( qx h ) ∂x ) + ∂ ∂y (h ∂( qx h ) ∂y )) (2.2) ∂qy ∂t + ∂( qxqy h ) ∂x + ∂( q 2 y h +g h 2 2 ) ∂y = -gh ∂Z ∂y -gq y √ q 2 x +q 2 y C 2 h 2 + D( ∂ ∂x (h ∂( qy h ) ∂x ) + ∂ ∂y (h ∂( qy h ) ∂y )) (2.3)
where t = time; h = water depth; Z = bed elevation; q x and q y = 'discharge' (water depth multiplied by velocity) along respectively the x-axis and the y-axis; C = Chezy friction coefficient which is considered as a constant or equalled to ST K × h 1 6 where ST K is the Strickler's friction coefficient; g = acceleration due to gravity; D = diffusion coefficient which is estimated as a constant.

Numerical model

Based on a finite volume method, the model Rubar20 is designed to solve the two-dimensional shallow water equations on unstructured triangular and rectangular grids. The model, written in Fortran 77 has been developed at Cemagref since 1988 (Paquier [21]). It uses a Van Leer explicit projectioncorrection scheme with second-order spatial and temporal accuracy in order to avoid spurious oscillations. The two main features of the numerical scheme are described here below. In case of hydraulic structures, the classical relations for weirs, gates, ... are used instead of the shallow water equations. In fact, an hydraulic structure can be defined as a set of a few cells where fluxes through one edge are computed either from the relations for the hydraulic structure, either from an addition of a flow computed from shallow water equations and a flow coming through an hydraulic structure. This last possibility may be useful, for instance, in the case of an embankment which is submerged although some discharge crosses it through a culvert. The model has been applied to simulate different situations of flows (Mignot [17], Mignot et al. [START_REF] Mignot | Modeling floods in a dense urban area using 2D shallow water equations[END_REF], Mignot et al. [START_REF] Mignot | Comparison of numerical and experimental simulations of a flood in a dense urban area[END_REF], Paquier [START_REF] Paquier | 1-D and 2-D models for simulating dam-break waves and natural floods, Union Européenne CADAM meeting[END_REF], [START_REF] Paquier | Surface flows during high floods in towns[END_REF], Paquier and Farissier [START_REF] Paquier | Use of a 2-D model fo simulating the flooding of a plain[END_REF], Paquier and Mignot [START_REF] Paquier | Use of 2D models to calculate flood water levels: calibration and sensitivity analysis[END_REF], Paquier et al. [START_REF] Paquier | Estimation of the flood levels for a flash flood in urban area: comparison of two hydrodynamic models on the 1988's Nîmes flood[END_REF], Proust et al. [START_REF] Proust | Flow in compound channel with abrupt floodplain contraction[END_REF]). The scheme is stable for Courant number less than 1. The two main steps of the numerical scheme are described here below.

Riemann problems

In a finite volume approach, the updated values of the main variables (h,q x ,q y ) are calculated using the fluxes crossing the edges. To estimate these fluxes, the 2-D problem constituted by the first member of the equations (2.1) + (2.2) + (2.3) is simplified in a 1-D problem in which the left and right values (on each side of the edge) are constant. If x means the direction perpendicular to the edge, the simplified 1-D problem writes:

∂h ∂t + ∂q x ∂x = 0 (3.1)
∂q x ∂t + ∂( q 2 x h + g h 2 2 ) ∂x = 0 (3.2)
∂q y ∂t + ∂( qxqy h ) ∂x = 0 (3.3)
The equations (3.1) + (3.2) + (3.3) constitute a nonlinear hyperbolic system for the 3-D variable W = (h, q x , q y ). Eigenvalues are u -c, u and u + c where u is the velocity perpendicular to the edge and c = √ gh is the wave celerity. The second field is linearly degenerated whereas the first and the third ones are truly nonlinear. The Riemann invariants can be very simply computed as they respectively equal u + 2c and v, u + 2c and u -2c, u -2c and v where v is the velocity parallel to the edge. To solve this Riemann problem, it is possible to use a Roe type linearization (Roe [START_REF] Roe | Approximate Riemann solvers, parameter vectors and difference schemes[END_REF]) which directly gives an estimate of the flux but alternative methods keeping the properties of the fields (and thus rarefactions and shocks) exist. Generally, these latter methods have the advantage to provide also an estimate of the values of the variables on the edge, which can be used when computing second member of the system (2.1) + (2.2) + (2.3) (Paquier [START_REF] Paquier | Modélisation et simulation de la propagation de l'onde de rupture de barrage[END_REF]). In order to obtain a second-order scheme in space, the values considered in the Riemann problem are computed by using slopes in the x and y directions (for each of the 3 variables h, q x and q y independently). For every cell, these slopes are computed by the method of the least squares (on each of the 3 variables independently) from the basis of the values at the centre of the cells which have one common edge with the given cell (Naaim [START_REF] Naaim | Numerical modelling of hydrodynamic effects caused by a solid sliding in a reservoir[END_REF]). The slope for any variable is later limited if the value of the given variable A (which may be h + Z, q x or q y ) at the middle of one edge is not located between the value at the centre of the cell A c and αA n + (1 -α)A c where A n is the value at the centre of the corresponding adjacent cell. The Van Leer correction factor α should be selected between 0.5 and 1. In order to obtain second-order accuracy in time, the values of fluxes should be assessed at time t n+1/2 = t n + 0.5∆t instead of time t n (Paquier [START_REF] Paquier | Modélisation et simulation de la propagation de l'onde de rupture de barrage[END_REF]).

Treatment of second member

Using the fluxes across edges, it is possible to obtain provisional values of the variables at time t n+1 = t n + ∆t. Then we add the contribution of the second member which includes, at least, 2 terms: i. Gravity or slope terms (-gh ∂Z ∂x or -gh ∂Z ∂y ).

ii. Bed friction by the use of a Chezy's coefficient (small water depths) or a Strickler's coefficient.

Gravity terms are treated as fluxes in such a way as an horizontal water surface remains strictly horizontal (Paquier [21]). For instance, the integral on one cell -gh ∂Z ∂x dxdy becomes 0.5

edges α i L i S (h c + h i )(Z i -Z c ),
in which α i is the x-co-ordinate of the normal to the edge i, L i is the length of the edge, h c and Z c are the water depth and the bed elevation at the centre of the cell respectively, h i and Z i are the water depth and the bed elevation at the middle of the edge i respectively, S is the area of the cell. Bed friction terms are more simply assessed at the centre of the cell. Their computation uses an implicitation in time in order to avoid numerical instabilities when a rapid change in water depth or velocity occurs.

4 Automatic differentiation

Generality

Derivatives play a fundamental role in various areas of computational science including optimization, data assimilation, assessment of uncertainties, sensitivity analysis, etc (Buecker et al. [4], Gilbert and Lemaréchal [START_REF] Gilbert | Some numerical experiments with variable storage quasi-Newton algorithms[END_REF]). A classical way to obtain these derivatives (from computer models) is the divided differences method. The rounding errors and cancelation are important drawbacks of this approximation method; furthermore, these errors and complexity increase when higher derivatives are calculated. To overcome these problems we use AD tool (Souhar [START_REF] Souhar | Evaluation of uncertainties for 1-D fluvial model. Automatic differentiation of 2-D Shallow water flows by Tapenade[END_REF], Souhar and Faure [START_REF] Souhar | Two gradient techniques to study sensitivities applied to 1-D and 2-D hydraulic models[END_REF], [START_REF] Souhar | Uncertainty propagation and sensitivity derivatives of a hydraulic model via automatic differentiation[END_REF]).

The AD, just like the divided differences, requires only the original program.

In mathematics, it's a method to numerically compute the partial derivatives of a function F specified by a computer program ( Buecker et al. [START_REF] Buecker | Automatic Differentiation: Applications, theory and tools[END_REF], Carle and Fagan [START_REF] Carle | ADIFOR 3.0 Overview[END_REF], Giering [START_REF] Giering | Tangent linear and Adjoint model compiler[END_REF], Griewank [START_REF] Griewank | Evaluating Derivatives: Principales and Techniques of Algorithmic Differentiation[END_REF], Hascoet and Pascual [START_REF] Hascoet | Tapenade 2.1 user's guide, Rapport INRIA n 300[END_REF]). An AD tool (Carle and Fagan [START_REF] Carle | ADIFOR 3.0 Overview[END_REF] or Hascoet and Pascual [START_REF] Hascoet | Tapenade 2.1 user's guide, Rapport INRIA n 300[END_REF] for example) builds a new source program that computes the derivatives of F . The tool behaves like a compiler, which first understands the meaning of the original program, represents it in an internal form, performs analyses to get deepen understanding, and generates an object code. AD exploits the fact that any computer program that implements a vector function Y = F (X) can be decomposed into a sequence of elementary assignments, any one of which may be trivially differentiated by a simple table lookup. These elementary partial derivatives are combined to form some derivative information of F . This process leads to exact numerical derivatives, because derivatives are computed analytically.

In AD there are two methods to compute derivatives, the so-called tangent linear model that computes directional derivatives and is easy to implement on computer. The second method is the reverse model that computes the gradient of a function at one execution of the adjoint program. This latter program is relatively expensive and should therefore be optimized whenever possible.

Let M be a general nonlinear model, a mapping from IR n of the input variables X = (x 1 , ..., x n ) T (initial conditions, topographic data, roughness coefficient, boundary conditions ... in hydraulic model cases) to IR m of the output variables Y = (y 1 , ..., y m ) T (model state, objective function ...):

F : IR n -→ IR m X -→ Y = F (X)
The vectors X ∈ IR n and Y ∈ IR m may be represented with respect to the ordinary unit basis vectors { e i ; i = 1, ..., n} and { f k ; k = 1, ..., m} as

X = n i=1 x i e i , Y = m k=1 y k f k .
The entries x i and y k forming the vectors X and Y are called input independent and output dependent variables respectively. And all the other variables used in the program to evaluate F (X) are called intermediate. Consider a perturbation dX of the input variables. Its effect on the outputs may be obtained via the linear approximation of the model M in terms of its Jacobian matrix, evaluated in the point X according to dY = ∇F (X).dX (4.1) Equation (4.1) is called the tangent linear model (Hascoet and Pascual [START_REF] Hascoet | Tapenade 2.1 user's guide, Rapport INRIA n 300[END_REF], Inria [START_REF]TAPENADE tutorial[END_REF]).

δy k = n i=1 ∂F k (X) ∂x i δx i , (4.2)
where, dX = n i=1 δx i e i and dY = m k=1 δy k f k .

The cost (algorithmic complexity) for computing gradients of scalar-valued functions with either divided differences method or tangent linear model of AD grows linearly with the number of variables, whereas the so-called reverse model of AD, an extremenlly efficient way, can compute such gradients at constant cost (Griewank [START_REF] Griewank | Some bounds on the complexity of gradients, Jacobians, and Hessians[END_REF]). Thus, if there are many input variables (10 6 -10 8 , for large-scale shallow water applications), it becomes prohibitive to proceed directly as in (4.2), and for these large applications the reverse model is used, see (Castaings et al. [START_REF] Castaings | Automatic differentiation: a tool for variational data assimilation and adjoint sensitivity analysis for flood modelling[END_REF], Gilbert [START_REF] Gilbert | Automatic differentiation and iterative processes[END_REF], Griewank [START_REF] Griewank | Some bounds on the complexity of gradients, Jacobians, and Hessians[END_REF], Hascoet and Pascual [START_REF] Hascoet | Tapenade 2.1 user's guide, Rapport INRIA n 300[END_REF], Inria [START_REF]TAPENADE tutorial[END_REF], Rall [START_REF] Rall | Automatic Differentiation Techniques and Applications[END_REF]) for more details.

In this work we use Tapenade, an AD software tool, to compute automatic derivatives of a two-dimensional hydraulic model. This tool is developed at INRIA Sophia-Antipolis by Tropics team (Hascoet and Pascual [14], Inria [START_REF]TAPENADE tutorial[END_REF]).

Building and validation of AD models

Derivation of models is not so easy, indeed a lot of difficulties can occur if direct models (the source models which must be differentiated) have not been written in the goal to be differentiated. The first step in using Tapenade, in particular to derive Rubar20 model, is to load the input source program and to identify which routine implements the mathematical function that we want to differentiate. We should then share the model in two parts corresponding respectively to the reading of data and initialization of the model and to the calculation a long the time. Therefore only the part of the model that achieves the computations is going to be differentiated.

A pre-treatment step for AD strategy must be done. Tapenade may issue a number of error and warning messages when differentiating a file is done. For example, in the source of the direct model there exists some subroutines with subroutine as arguments. Tapenade doesn't accept this type of instructions and hence obliges the user to modify manually the source model. A post-treatment step is also required, it consists to insert read data and write results instructions in the differentiated subroutines and the necessary calls of subroutines like those of initializations. After this, a part of debugging compilation errors must be done (Souhar [START_REF] Souhar | Evaluation of uncertainties for 1-D fluvial model. Automatic differentiation of 2-D Shallow water flows by Tapenade[END_REF]).

Tapenade can be utilzed as a server, which runs at INRIA Sophia-Antipolis at http://tapenade.inria.fr.8080/tapenade/indes.jsp. It can also be installed locally, in that case it is run by a simple command line, which can be included into a Makefile. Figure 1 shows a typical Makefile built to make easy the derivation of the code after each modification. So the options are defined as follows.

• head allows to specify the top differentiation routine that implements the mathematical function.

• inputlanguage allows to select the source code language.

• outputlanguage allows to select the output differentiated code language.

• outvars allows to define the output variables and whose derivatives are required (here the outputs are the water depth and discharge).

• vars allows to choose the input variables with respect to which differentiated must be made (here the inputs are the friction coefficients).

• O point out the directory of the differentiated code.

The command make d (resp. make b) permit then to obtain the tangent linear model (resp. the reverse model) of the source program defined by its arborescence given in FILES option of Fig. 1.

Tapenade will probably not produce a perfect differentiated program at first try. After AD, a validation step is necessary to check the computed derivatives. It should be noted that a common technique to check validation of these derivatives is to compare them to finite differences approximations, for example:

∂F (x) ∂x i = lim ε→0 1 ε [F (x 1 , ..., x i + ε, ..., x n ) -F (x 1 , ..., x i , ..., x n )]
for i = 1, ..., n. In (Hascoet and Pascual [START_REF] Hascoet | Tapenade 2.1 user's guide, Rapport INRIA n 300[END_REF], Inria [START_REF]TAPENADE tutorial[END_REF]), the developers of Tapenade proposed this technique to validate tangent linear models and it is called "gradient test". The so-called "dot product test" is proposed in the same user's guide to validate reverse models. Figure 2 below shows the validation of derivatives of Rubar20 obtained by Tapenade. This figure shows well that the error function r(ε) -→ 0 when ε -→ 0, where r is defined as follows:

r(ε) = F (x + ε.d) -F (x) ε -∇F (x).d 2 , (4.3)
where d = (1, ..., 1) T is the unit vector of IR n and ∇F (x).d is computed by Tapenade. In all of our tests, the precision for water depths is 10 -4 m, the correction Van Leer's coefficient is 0.6, the Courant number is 0.5 and the duration of the simulation is 80 s.

Validity of Rubar20

We are interested here in a complex problem. The flow, a portion of a river of length 28.5 m and width 20 m, is perturbed by some obstacles. The portion is crossed by a bridge with two triangles representing the bridge abutments and two cylindric bridge piers schematized by two octagons (see Fig. 3). The cross sections of the main channel are trapezoidal of slope 2/1, a discharge of 62 m 3 /s constitutes the upstream boundary condition, the downstream condition is imposed such that Z + h = 0. The friction coefficient according to Manning law is homogeneous and uses a Strickler's coefficient equal to 40. We study now the recirculation behind the bridge piers and the bridge abutments for two different diffusion coefficients D. Fig. 4 and Fig. 5 give the velocity fields which permits to assess the recirculation's length. For D = 5.10 -2 (Fig. 4), the results are near to those computed by Télémac (Télémac is a numerical model developed by Electricité de France) (Galland et al. [START_REF] Galland | Code TELEMAC, Version 1.0, Note de validation[END_REF]). For D = 5.10 -3 , Télémac finds eddies behind the bridge piers. Rubar20 gives recirculations behind the bridge abutments (Fig. 5). This difference results mainly from a small difference of grids used for computations in the two codes.

Sensitivity analysis

Assessment and sensitivity analysis play an important role in decision making on flood control, so minimization of the risk should be the general objective.

In the following we are interested in the sensitivities of water depth and discharge with respect to the roughness coefficient. We compute via AD the corresponding derivatives.

The study has a homogeneous roughness, then it is important to define the subspace IK = span{ϕ} of IR n , where n = 548 is the number of computational points and the basis vector is ϕ = (1, ..., 1) T ∈ IR n , in this case the function F of Subsection 4.1 becomes:

F : IK -→ IR m (ST K) -→ F (ST K) = (F 1 (ST K), ..., F m (ST K)) T , (5.1)
where F j (ST K) represents the water depth or the discharge at section j and m = 548, and ST K is the friction coefficient.

According to Taylor's formula, which illustrates the physical definition of the gradient, the values of the function F at the points ST K and ST K + δϕ are related by:

F (ST K + δϕ) = F (ST K) + δ∇F (ST K).ϕ + O(δ 2 ϕ), (5.2)
with

∇F (ST K) =         dF 1 (ST K) dx . . . dFm(ST K) dx         (5.3)
Fig. 6 and Fig. 7 show gradients of water depths (m) and discharges (m 2 /s) related to the friction coefficient (defined by (5.3)) computed using the AD tool Tapenade. So the sensitivities of the water depths and discharges are also symmetric for both diffusion coefficient values. We remark that these sensitivities vary in the opposite of the diffusion coefficient; i.e., sensitivities are less high for D = 5.10 -2 than for D = 5.10 -3 . We remark also that discharge's sensitivities are high behind the bridge piers and the bridge abutments for both diffusion coefficient values. According to Fig. 4 and Fig. 6 (D = 5.10 -2 ), an increase of the friction coefficient may lead to find recirculations behind the bridge piers, whereas a decrease of the friction coefficient may lead to find recirculations behind the bridge abutments. From Fig. 5 and Fig. 7 (D = 5.10 -3 ), we remark that an increase of the friction coefficient may lead to find more important recirculations behind the bridge piers and to increase recirculations zones behind the bridge abutments.

Conclusion and remarks

The two-dimensional hydraulic model Rubar20 provides satisfactory results faced to a complex problem. The model can provide detailed velocity field and water levels. But a computation with fine grids should be used to detail and understand the recirculation problem more accurately. The computational technique can be easily extended to channels with arbitrary cross section, for computing the propagation of a flood surge on a dry bed, also for the computation of the flood wave resulting from the break of a dam, etc. We have shown how AD tool can aid in the sensitivity analysis of a two dimensional finite volume model of shallow water flows. So the numerical sensitivities of water depths and discharges related to the friction coefficient are also symmetric as the results provided by the model and because of the symmetry of the channel. Furthermore, these sensitivities become more important for a small diffusion coefficient. Such sensitivities employed the forward derivatives computed using the AD software Tapenade. But it is still necessary to manually modify the source codes before automatic derivative step because AD tools still remain unable to take into account all programming instructions.

[ 1 ]

 1 , Castaings et al. [3], Corliss et al.

Figure 1 :

 1 Figure 1: A typical Makefile of "Rubar20"

Figure 2 :

 2 Figure 2: Gradient test for a 2-D model: "Rubar20"

Figure 3 :

 3 Figure 3: Grid for computation of the study

Figure 4 :

 4 Figure 4: Velocity field in steady flow using Rubar20: D = 5.10 -2

Figure 5 :

 5 Figure 5: Velocity field in steady flow using Rubar20: D = 5.10 -3

Figure 6 :

 6 Figure 6: Sensitivities of water depths and discharges to bottom friction: D = 5.10 -2

Figure 7 :

 7 Figure 7: Sensitivities of water depths and discharges to bottom friction: D = 5.10 -3