Chi Cong Nguyen

Eric Gaume
email: eric.gaume@ifsttar.fr

O Payrastre

Regional flood frequency analyses involving extraordinary flood events at ungauged sites: further developments and validations

Keywords: Floods, regional analysis, statistics, Extremes, GEV

Regional flood frequency analyses involving extraordinary flood events at ungauged sites : further developments and validations

Introduction

Although a large number of statistical inference methods have been progressively developed, the question of estimating extreme design floods is still problematic due to the generally limited amount of available data. Continuous discharge series at gauged sites are generally too short to provide reliable estimates of extreme quantiles -typically the 100-year or higher return period quantiles (NERC, 2000). To cope with this difficulty, hydrologists have tried to complement the available data sets, either through a "temporal extension", incorporating data on historical and paleofloods (Hosking and Wallis, 1986a,b;[START_REF] Stedinger | Flood frequency analysis with historical and paleoflood information[END_REF][START_REF] Cohn | Use of historical information in a maximum likelihood framwork[END_REF][START_REF] Gary | Regional regression of flood characteristics employing historical information[END_REF][START_REF] Sutcliffe | The use of historical records in flood frequency analysis[END_REF][START_REF] Minghui | Flood frequency analysis with regional and historical information[END_REF][START_REF] Sheffer | Paleofloods and historical floods of the Ardche river, France[END_REF][START_REF] Reis | Bayesian GLS regression with application to LP3 regional skew estimation[END_REF][START_REF] Neppel | Flood frequency analysis using historical data: accounting for random and systematic errors[END_REF][START_REF] Payrastre | Usefulness of historical information for flood frequency analyses: Developments based on a case study[END_REF], or through a "spatial extension", merging data sets in regions considered as statistically homogeneous, "trading space for time" according to the words of Hosking & Wallis [START_REF] Hosking | Regional frequency analysis: An approach Based on L-Moments[END_REF][START_REF] Heinz | Using regional regression within index flood procedures and an empirical Bayesian estimator[END_REF][START_REF] Charles | Development of regional regression relationship with censored data[END_REF][START_REF] Ouarda | Regional flood frequency estimation with canonical correlation analysis[END_REF][START_REF] Kjeldsen | Regional flood frequency analysis in the KwaZulu-Natal province, South Africa, using the index-flood method[END_REF][START_REF] Merz | A process typology of regional floods[END_REF][START_REF] Seidou | A parametric Bayesian combination of local and regional information in flood frequency analysis[END_REF][START_REF] Ribatet | Usefulness of the reversible jump Markov Chain Monte Carlo in regional flood frequency analysis[END_REF][START_REF] Norbiato | Regional frequency analysis of extreme precipitation in the eastern italian Alps and the August 29, 2003 flash flood[END_REF][START_REF] Wallis | Regional precipitation-frequency analysis and spatial mapping for 24-hour and 2-Location River name Date Q (m 3 /s) S (km 2 ) Pont de Rolandy Ardèche 22/09[END_REF][START_REF] Kjeldsen | A formal statistical model for pooled analysis of extreme[END_REF][START_REF] Haddad | Regional flood frequency analysis in eastern Australia: Bayesian GLS regression-based methods within fixed region and ROI framework: Quantile Regression vs. Parameter Regression Technique[END_REF].

Recently, [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF] observed that estimated extreme peak discharges at ungauged sites are often available, but never really used in flood frequency studies and proposed a method to incorporate such information in regional flood frequency analyses.

The proposed approach is based on the index flood principle [START_REF] Dalrymple | Flood frequency analyses. Water Supply: Geological Survey[END_REF], assuming that, within a statistically homogeneous region, all local statistical distributions are identical apart from a site-specific scaling factor: the index flood. Usually, the index flood corresponds to the mean of the local series [START_REF] Hosking | Regional frequency analysis: An approach Based on L-Moments[END_REF]. The approach proposed by [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF] is based on the calibration of an index flood relation linking the characteristics of the watersheds and the index flood value. Although this relationship represents an additional homogeneity requirement that may limit the extent of the region used for the statistical analysis, it also enables to estimate the index flood at ungauged sites, and thus to incorporate the corresponding ungauged extremes in the regional sample.

Based on several case studies, [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF] showed the possible great value of such an approach, depending on the characteristics of available extreme flood inventories. The index flood relationship proposed, of the form S β (S being the area of the watershed and β a parameter to be calibrated), appeared satisfactory in the test regions. The presented inference results were based on a Bayesian MCMC framework [START_REF] Castellarin | Probabilistic behavior of a regional envelope curve[END_REF][START_REF] Reis | Bayesian GLS regression with application to LP3 regional skew estimation[END_REF][START_REF] Seidou | A parametric Bayesian combination of local and regional information in flood frequency analysis[END_REF][START_REF] Ribatet | Usefulness of the reversible jump Markov Chain Monte Carlo in regional flood frequency analysis[END_REF][START_REF] Castellarin | Probabilistic envelope curves for design flood estimation at ungauged sites[END_REF][START_REF] Payrastre | Use of historical data to assess the occurrence of floods in small watersheds in the French Mediterranean area[END_REF][START_REF] Payrastre | Usefulness of historical information for flood frequency analyses: Developments based on a case study[END_REF] to adjust the regional growth curve with associated 90% credibility intervals. The results showed that the incorporation of ungauged extremes could lead to a significant reduction of the width of the computed credibility intervals.

In the initial version of the method [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF], the index flood relation was adjusted prior to the calibration of the regional growth curve, and the uncertainties associated with its calibration were not taken into account. This led certainly to underestimate the credibility intervals and over-rate the added value of the ungauged extremes and of the proposed method. The effects of possible variations (heterogeneities) in the average relation calibrated in a given region should also be considered for a fair comparison with other statistical methods. This paper proposes an extension of the initial method to account for uncertainties in the calibrated index flood relation. It also tests the effect of possible regional variations in the average relation on the efficiency of the proposed statistical inference approach.

The performances (i.e. widths and correctness of computed credibility intervals) of the proposed approach and of the standard regional frequency approach proposed by [START_REF] Hosking | Regional frequency analysis: An approach Based on L-Moments[END_REF] are first compared in the case where gauged data only are considered. The comparison is based on samples generated through Monte Carlo simulations in order to be able to verify the accuracy of the calculated credibility intervals and to introduce controlled heterogeneities in the samples. In a second step, both approaches are applied to the statistical analysis of a data set from the Ardèche region in France composed of 168 records at 5 gauging stations and 18 estimated ungauged extremes.

The paper is organised as follows: section 2 presents the basics and adaptations of the two regional flood frequency methods: Hosking & Wallis and the proposed approach. The performances of the two approaches are compared based on simulated samples of random variables in section 3. In section 4, the methods are applied to the real-world case study. Conclusions are drawn in section 5.

Tested regional flood frequency analysis methods

The index flood hypothesis

The two approaches considered in this paper are based on the same fundamental simple scaling hypothesis or index-flood principle [START_REF] Dalrymple | Flood frequency analyses. Water Supply: Geological Survey[END_REF]: in a statistically homogeneous region, all the local annual maximum peak discharge distributions are supposed to be identical apart from a site-specific scaling factor. This hypothesis is summarized in equation 1:

Q i (F ) = µ i q (F ) (1)
Where F is the probability of non-exceedance, i is the index of the site (i = 1, ..., s), s the total number of sites in the homogeneous region, Q i (F ) is the discharge quantile, q (F ) is the regional dimensionless (i.e. reduced) quantile and µ i is the index flood (or scaling factor).

The index flood may be any constant value proportional to the expectancy of the local distribution. Usually, when only data from gauged sites are considered, the index flood is estimated by the at-site sample mean [START_REF] Hosking | Regional frequency analysis: An approach Based on L-Moments[END_REF][START_REF] Castellarin | Probabilistic behavior of a regional envelope curve[END_REF][START_REF] Castellarin | Probabilistic envelope curves for design flood estimation at ungauged sites[END_REF]. A regional flood frequency method where the index flood is computed as the average of the local series of annual maxima will be called hereafter method of Hosking & Wallis. This, even if a likelihood based Bayesian MCMC procedure rather than a L-moment based procedure, as suggested by [START_REF] Hosking | Regional frequency analysis: An approach Based on L-Moments[END_REF], is used to calibrate the parameters of the regional growth curve. [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF] suggested an alternative approach to account for extreme discharge estimates that may be available at ungauged sites. An inventory of ungauged extremes may includes h extreme peak discharges Q k (k = 1, .., h), each Q k corresponding to the largest flood at site k during a period of length n k . In order to include this additional information in the regional dataset, [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF] proposed to use an index flood relation linking the index flood value to the catchment area S, since an average annual peak discharge can obviously not be computed at ungauged sites:

µ i = S β i and µ k = S β k (2)
Where S i and S k are the catchment areas at the corresponding sites, and β a coefficient to be calibrated. More complex relations based on various climatic and physio-geographic characteristics may be tested in the future, but at the price of an increased number of parameters to be calibrated. In the initial version of the method, the value of β was adjusted through a regression between the log tranform of the average annual peak discharges and the watershed areas at gauged sites.

It is proposed here to calibrate β along with the parameters of the regional growth curves using a modified version of likelihood as described below.

Likelihood of the observed sample

The inference approach applied herein is directly derived from [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF] and inspired by numerous previous works [START_REF] Reis | Bayesian GLS regression with application to LP3 regional skew estimation[END_REF][START_REF] Renard | An application of Bayesian analysis and Markov Chain Monte Carlo methods to the estimation of a regional trend in annual maxima[END_REF][START_REF] Payrastre | Usefulness of historical information for flood frequency analyses: Developments based on a case study[END_REF]: i.e. based on the likelihood of the available data sets and a Bayesian MCMC algorithm for the estimation of the parameters of the growth curve and of their posterior distribution according to the observed data set.

Considering the regional sample D described above, including both (i) s series of gauged annual maximum discharges Q i,j , j = 1, • • • , n i the index of the year, and (ii) the h estimated largest peak discharges Q k over n k years at h ungauged sites, the standard expression of the likelihood of the regional sample D would be the following:

(D | θ) = s i=1 n i j=1 f θ Q i,j µ i h k=1 f θ Q k µ k h k=1 F θ Q k µ k (n k -1) (3)
Where f θ and F θ are respectively the probability density function and the cumulative probability function of the selected statistical distribution for the regional growth curve, and θ corresponds to the vector of parameters to be estimated. The GEV distribution, often used to describe peak discharge growth curves [START_REF] Lu | Sampling variance of normalized GEV/PWM quantile estimators and regional homogeneity test[END_REF][START_REF] Stedinger | Appraisal of regional and index flood quantile estimators[END_REF][START_REF] Coles | Bayesian methods in extreme value modelling. A review and new developments[END_REF][START_REF] Coles | A bayesian analysis of extreme rainfall data[END_REF][START_REF] Heinz | Using regional regression within index flood procedures and an empirical Bayesian estimator[END_REF][START_REF] Seidou | A parametric Bayesian combination of local and regional information in flood frequency analysis[END_REF], was selected here (Eq. 4 and 5). The vector θ comprises the position, scale and shape parameters (ξ, α, κ) of the GEV distribution.

F θ (Q) = exp -1 - κ (Q -ξ) α 1/κ α>0 (4) f θ (Q) = 1 α 1 - κ (Q -ξ) α 1/κ-1 exp -1 - κ (Q -ξ) α 1/κ α>0 (5)
In equation 3, the first term corresponds to the probability of the gauged series. It is the only necessary term if continuous series of measured annual maximum discharges are used. The second term is the probability of the ungauged extremes. The third complementary term is the probability associated to the fact that the ungauged extreme value has not been exceeded during the remaining (n k -1) years at each ungauged site. The index values µ can be estimated before the calibrating the regional growth curve : method of Hosking & Wallis and initial version of the proposed approach. Equation 2can be also directly introduced in the formulation of the likelihood (eq. 3)

adding one parameter, β, to be calibrated. This replacement necessitates a slight but determining modification of the expression of likelihood (section 2.4). Depending on the option, the set of parameters to be calibrated will count either 3 or 4 parameters.

Bayesian Monte Carlo Markov Chain algorithm

The Bayesian Monte Carlo Markov Chain procedure is now relatively common for hydrological applications [START_REF] Reis | Bayesian GLS regression with application to LP3 regional skew estimation[END_REF][START_REF] Renard | An application of Bayesian analysis and Markov Chain Monte Carlo methods to the estimation of a regional trend in annual maxima[END_REF][START_REF] Seidou | A parametric Bayesian combination of local and regional information in flood frequency analysis[END_REF][START_REF] Ribatet | Usefulness of the reversible jump Markov Chain Monte Carlo in regional flood frequency analysis[END_REF][START_REF] Neppel | Flood frequency analysis using historical data: accounting for random and systematic errors[END_REF][START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF][START_REF] Payrastre | Usefulness of historical information for flood frequency analyses: Developments based on a case study[END_REF][START_REF] Viglione | Flood frequency hydrology: 3.A Bayesian analysis[END_REF] and will only be briefly presented here. Part of the algorithms used is included in the R software library nsRFA [START_REF] Viglione | Non-supervised regional frequency analysis R library[END_REF]. Recalling the Bayes' theorem, the likelihood of the sample given the parameters of the statistical model (D | θ) can be related to the likelihood or density of probability of the parameters given the sample p (θ | D) (posterior distribution):

p (θ | D) = (D | θ) p (θ) p (D) (6)
Where p (θ) is the so called prior distribution of θ, which summarizes any prior or alternative knowledge on the distribution of θ, and p (D) is the probability of the sample D which is unknown. When prior information on the distribution of θ does not exist, then p (θ) is often taken as uniform.

It is the case here, which implies that p (θ

| D) is proportional to (D | θ).
The statistical model being chosen, it is possible to compute the probability density function of its parameters according to the observed data sample.

p (θ | D) being known, or more precisely (D | θ) which is proportional to it according to our hypotheses (uniform prior), parameter sets θ will be sampled according to p (θ | D) to build their posterior distributions and compute the corresponding credibility limits for the discharge quantiles. The MCMC algorithms, combining random walk Monte Carlo methods with Markov chains, are a class of algorithms for the efficient sampling from multivariate random distributions [START_REF] Tanner | Tools for statistical inference: Methods for the Exploration of Posterior, Distributions and Likelihood Functions[END_REF][START_REF] Robert | Monte Carlo statistical methods[END_REF].

Figure 1 illustrates the Bayes MCMC inference procedure: posterior distributions of two of the estimated parameters for 3 different Markov chains run in parallel.

Modification of the likelihood formulation

A modification of the values of the index floods µ i and µ k , corresponding to a change of scale of the analyzed discharges, will also affect the values of the optimal parameters ξ and α, the variable (Q/µ -ξ)/α remaining unchanged. The absolute value of the cumulative density function (eq. 4) will be unaffected by such a change of variable, but not the absolute value of the density function, due to the factor 1/α in front of the expression of the density (eq. 5). Any method based on the maximization of likelihood expression 3, will therefore tend to minimize the value of α, and hence to maximize the values of the indices µ i = S β i if they are calibrated. Therefore, the use of likelihood expression 3 in an inference procedure where the parameter β is adjusted will inevitably lead to biased estimates of β (over-estimations)

and of the parameters of the regional growth curve which are linked to β. This is confirmed by numerical tests as illustrated in Figure 1. To avoid this disruptive effect, one possibility consists in including only cumulative densities in the likelihood expression (eq. 7).

(D | θ) = s i=1 n i j=1 F θ Q U i,j S β i -F θ Q L i,j S β i h k=1 F θ Q U k S β k -F θ Q L k S β k h k=1 F θ Q U k S β k (n k -1) (7)
Where the exponents U and L indicate the upper and lower estimate for the considered discharge values (upper and lower guesses for an estimated ungauged extreme or known uncertainty range for measured discharges). It has been verified that equations 7 and 3 provide very close results if the difference between Q U i,j and Q L i,j is small. In the rest of the manuscript Q U i,j

will be set equal to 1.01 Q i,j and Q L i,j to 0.99 Q i,j .

Figure 1 illustrates the type of results obtained through the application of 10 the MCMC algorithm using either expressions 3 or 7 of the likelihood, when the index flood relation is directly introduced in the likelihood formulation and the parameter β is calibrated along with the parameters of the regional growth curve. The simulated regional data sample is described in section 3.

These results show the tremendous under-estimation of the parameter α of the regional growth curve when equation 3, containing density functions, is used. The assessment seems correct with equation 7, even if affected by uncertainties, due to the finite number of records (i.e. to the limited information content of the data set).

MCMC convergence diagnosis

Like for any other optimization method, convergence of the MCMC algorithm towards the posterior distribution of the model parameters may sometimes be difficult to reach, especially when the number of parameters to be calibrated is increased. For the inference on a single data set, convergence Among the large number of methods developed for convergence monitoring of MCMC algorithms [START_REF] Cowles | Markov Chain Monte Carlo convergence diagnostics: a comparative review[END_REF][START_REF] Salaheddine | Comparison of methodologies to assess the convergence of Markov Chain Monte Carlo methods[END_REF], the popular Gelman & Rubin test has been selected [START_REF] Gelman | Inference from iterative simulation using multiple sequences (with discussion)[END_REF].

It consists in running several Markov chains in parallel and in computing the evolution during the iterations of the square root of the ratio between average within-chain and average between-chains variance of the likelihood.

This ratio should ideally be equal to 1 if convergence is achieved.

The test was applied herein with two chains and a threshold value of 1.05. Figure 2 shows the evolution of the Gelman & Rubin criterion with the number of iterations when MCMC algorithms are run on synthetic data samples presented in section 3. To achieve fast convergence, a first MCMC chain was systematically run with a limited number of iterations (30000) to explore the shape of likelihood function in the vicinity of the parameter values maximizing the likelihood and evaluate the covariance matrix of the parameters (important controling factor of the MCMC chain). Figure 2 shows that the convergence is rapidly achieved with the reference approach, when it appears to be more difficult with the proposed approach. According to the conducted tests, the maximum iteration number was set to 100 000 and the various MCMC outputs computed on the basis on the last 30 000

iterations. This number of iterations was sufficient to ensure the convergence of the MCMC algorithm for all the simulations presented in this paper.

First test of the two methods

The two regionalisation approaches (standard and proposed) have been first tested for the statistical analysis of simulated random data sets to evaluate their performances, before their application to a real-world case study (statistical analysis of data samples available in the Ardèche region, South of France). A special emphasis has been put on (i) the validation the credibility intervals computed with the MCMC procedure and (ii) the analysis of the sensitivity to regional variability in the index flood relation.

Method

(i) Monte Carlo simulation procedure.

The test set is composed of 1000 simulated random samples drawn from a regional GEV distribution (ξ=3.34, α=2.24, κ=-0.16) and with an index flood relation (β=0.76) adjusted on the Ardèche set analysed in section 4.

The samples have the same characteristics as the Ardèche regional gauged sample: catchment areas and record lengths, see table 1.

Both regionalisation approaches are then applied to each simulated dataset, including the estimation of the posterior distribution of parameters thanks to the Bayesian MCMC procedure.

Results

(i) Case of homogeneous data sets Figure 4 presents the inference results based on 1000 simulated samples, without introduction of regional variability in the index flood relation. The fluctuations of the maximum likelihood estimates QML i (0.99) (respectively qML (0.99) for the regional growth curve), and associated 90% credibility upper and lower bounds Q5 i (0.99) , Q95 i (0.99) , are presented both, for the estimates corresponding to each gauged watershed area of the region and for the regional growth curve. In order to facilitate comparisons, all estimations have been divided by the real quantile value Q i (0.99) (or q (0.99)). Therefore the position of the real quantile is 1 for all the box plots. areas and explains the large uncertainties for the estimated quantiles for the regional growth curve in the proposed approach.

An in depth analysis of both figures reveals nevertheless some differences of the two regionalisation approaches: the median values of the maximum likelihood and 5% and 95% credibility bounds for the 100-year quantile are slightly lower when the standard Hosking & Wallis method is used, but the uncertainties affecting these estimates (widths of the box plots) appear larger except for the smallest considered watershed area (63 km 2 ). The median maximum likelihood estimated value lower than 1 for the standard method suggests also that it is slightly biased.

The reliability tests conducted on the computed credibility intervals indicate furthermore that, while the intervals seem correctly estimated with the proposed method (Fig. 5.e and f) with the exception of the smallest watershed area where a tendency to the overestimation of the quantile appears. The credibility intervals computed with the standard method appear too narrow (Fig. 5.a and b). A complementary test presented in appendix Appendix A indicates that a 20% to 40% underestimation affects the computed widths of credibility intervals in this case study: a correction factor m equal to 1.2 to 1.4 should be applied to the standard-deviation of the posterior distributions of Qi (0.99) to obtain a uniform distribution of F (Q i (0.99)) (Fig. 5.c and d). This under-estimation of the variance of the posterior distributions of Qi (0.99) for the standard method is logical since the uncertainty associated with the estimation of the local sample means is not accounted for in the expression of the likelihood. Its effect on the quantile estimations can therefore not be evaluated by the Bayesian MCMC procedure.

The origins of the observed limited biases -tendency to underestimate the quantile on average for the standard method and to overestimate the quantile for the smaller areas for the proposed approach -are difficult to depict.

These results are overall extremely satisfactory. When the regional sample corresponds to the proposed statistical model (regional growth curve of the GEV type and power-law linking the index flood and the area of the watershed), the proposed approach appears consistent and accurate. The parameters of the regional statistical model seem correctly retrieved from samples with only a slight bias in the estimation of the parameter β affecting the estimates of the discharge quantiles of the watersheds with the smallest areas. In the standard method, an additional source of uncertainty seems to be introduced by the computation of local sample means, estimates of the local sample expectancies, leading to an unexplained minor but detectable underestimation of the quantiles and to an increase of the estimation uncertainties of these quantiles if compared to the proposed method. The credibility intervals computed through the MCMC algorithm underestimate in this case the uncertainties associated to the calculated quantiles.

The additional parameter β appears to be a factor of stabilization rather than a source of further uncertainties in a flood frequency analysis. Fixing an index flood relation seems to help to filter out the variability of the local sample means due to sampling fluctuation. This is an unexpected result.

The proposed approach, initially designed to introduce ungauged extremes in flood frequency analysis, is as efficient, if not more efficient than the standard approach when applied to gauged series. To be really able to draw such a conclusion, it is necessary to consider that the regional data set may not follow perfectly the simple suggested index flood relation. How does the method resist to some regional variability in the index flood relation? It is the focus of the next section.

(ii) Impact of some regional variability in the index flood relation Without surprise, the introduction of variability in the index flood relation increases both, the estimated credibility intervals and the uncertainties affecting the estimates of the maximum likelihood, 5% and 95% credibility bound estimates (Fig. 6).

The effect remains nevertheless limited when the variability is relatively low (Fig. 6.a). The results obtained for δ = 0.1 do not differ much from the results obtained in the ideal case (δ = 0 in Fig. 4.b). The maximum likelihood estimate is still unbiased and it could be verified that the credibility intervals were also only slightly biased. Uncertainties increase drastically when the regional variability is increased (Fig. 6.b), but still with a limited bias, on average, on the maximum likelihood estimate of the quantile (overestimation)

and some bias appearing in the estimated credibility intervals. Of course, the biases on the individual local estimates can be significant as shown by the widths of the box plots.

These results finally suggest that the proposed approach may resist to a certain level of heterogeneity of the proposed index flood relationship in the considered region. Important heterogeneities should nevertheless be avoided through a careful selection of the regional data set based on the plot of the local sample means versus watershed areas [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF]. This poses an additional constraint to the delineation of homogeneous regions if compared to the method of Hosking & Wallis, limiting the size of the regional gauged sets, but opening the possibility to valuate the ungauged extremes.

This constraint should be kept in mind for a fair comparison of both methods.

4. Application to a real-world example with ungauged extremes:

the Ardèche area in France

The Ardèche data set

The Ardèche case study had already been selected by [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF] for the initial tests of the regionalisation method involving ungauged extremes. It has been again used here to illustrate the changes introduced by the new inference procedure including the calibration of the index flood relationship.

The Ardèche region is located in the southeast of Massif Central in southern France (Fig. 7). It is one of the areas in Europe exposed to the most frequent and severe flash floods [START_REF] Gaume | A compilation of data on European flash floods[END_REF] Using the large set of documentary sources available in the region, ungauged largest peak discharge estimates over the last 50 years could be retrieved at 14 additional sites (Fig. 8). Similarly, information on the 50 years preceding the gauging period could be found at 4 of the gauged sites. The set of ungauged extremes is finally composed of 18 maximum values over 50 years, that is equivalent to 900 censored records, six times more than the regional gauged data set (Table 2). Note that the selected set of ungauged extremes has been reduced if compared to the set used by [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF]. A minimum distance between sites on the same stream has for instance been imposed to avoid redundancy (see appendix Appendix B for a short discussion about the sampling of extremes).

The homogeneity of the growth curves within the region has been tested using the Hosking & Wallis heterogeneity measure [START_REF] Hosking | Regional frequency analysis: An approach Based on L-Moments[END_REF]. The H 1 value based on the 5 gauged series is equal to 2.7, which suggests that the region is possibly heterogeneous. Homogeneity has nevertheless been assumed to proceed with the computations. The influence of regional heterogeneities in the growth curves on the inference results will be studied in a future publication.

Provided the preceding assumption, this case study finally offers an interesting opportunity to apply and compare the approaches described above, in two different contexts: (i) regional dataset limited to gauged sites, and (ii) availability of information on the largest discharge at ungauged sites. The two methods -Hosking & Wallis and the proposed approach -lead to close results when applied to the gauged data set (Table 3, cases a and b). The widths of the estimated credibility intervals are slightly larger in the case of the proposed approach, but these intervals are underestimated by about 20% in the reference approach according the previous results. The underestimation of the quantiles is not systematic when the reference approach is used as illustrated by table 3. The inference results are depending on the characteristics of the sample.

Results and discussion

The main advantage of the proposed approach lies in the possibility to incorporate information on ungauged extremes in the analysis. In the Ardèche case, the ungauged extremes appear consistent with the gauge series and the GEV distribution (Fig. 8). Their incorporation in the statistical inference leads to a clear reduction of the credibility intervals that are divided by a factor of about 2. This confirms the possible great value of the proposed approach in regions for which information on extremes floods at ungauged sites is available. The assumption of homogeneity of the index flood relationship limits the number of gauged sites that can be included in the homogeneous region. But, the extremes at ungauged sites represent a potentially large ad-ditional source of information (equivalent to 900 records in the present case), that may reduce significantly the uncertainties attached to estimated high return period quantiles. the Eyrieux, appears to have relatively smaller empirical mean discharge and could be possibly regarded as non-consistent with the other sites as far as the index flood relation is considered. This low mean value may nevertheless be attributed to sampling fluctuation but may likewise be explained by the fact that this site is located further North and less influenced by the Mediterranean climate. This last comment illustrates that any regional flood frequency method is based on homogeneity hypotheses that can hardly be verified. It is therefore certainly sound to test the robustness of the obtained inference results to different hypotheses and region delineations. In this perspective, the proposed approach enriches the available panel of methods and possible hypotheses.

Finally, the same Monte Carlo simulations as the ones used in section 3

were conducted with generated samples including 18 ungauged sites where only the highest discharge over 50 years is known to verify if the conclusions drawn on the real-world example can be generalized (Fig. 10.d). In fact, the results of an anlysis conducted on one single data set may be attributed to sampling variability. On average, the credibility intervals appear to be divided by two on a sample of 1000 simulated datasets when the "ungauged" extremes are taken into account, which is the order of magnitude obtained for the Ardèche data set. 100-year quantile estimates based on the L-moments (standard method suggested by Hosking & Wallis) have also been computed for sake of verification (Fig. 10.a). L-moment and maximum likelihood based estimates appear to have similar performances on this data sample (figures 10.a and b). The L-moment estimates exhibit a slightly lower dispersion (L-moment are less sensitive to sampling variability) but a higher bias.

Conclusions

The approach initially proposed by [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF], aiming at incorporating information on ungauged extremes in regional flood frequency analyses, has been consolidated and thoroughly tested. The proposed method and algorithm, after some necessary adjustments, in particular of the formulation of the likelihood, appear to be able to estimate accurately quantiles and their uncertainty bounds on the basis of a regional sample combining continuous gauged series and ungauged censored data. The introduction of an index flood relationship and of an additional parameter to be calibrated, proved to be a stabilization factor rather than a factor of complexity and a source of uncertainty. The method shows a relative robustness to some regional variability in the index flood relation. The additional homogeneity constraint limits the possible extent of homogeneous regions if compared to the reference approach, but this is compensated by the additional information brought by the estimated ungauged extremes. In the case of the Ardèche region the ungauged data set represents 900 censored recording station-years:

six times the amount of the available regional gauged series. Its incorporation in the regional flood frequency analysis led to a division by two of the estimated 100-year quantile 90% credibility intervals.

Uncertainties on the extreme peak discharge estimates were not considered herein and the ungauged extreme dataset that could be retrieved in the et al. (2011). It remains to be confirmed that the same conclusion holds for ungauged extremes when included in regional flood frequency analyses.

In the same line of thoughts, the possible variability in the growth curves in a region considered as homogeneous has not been considered yet. Its influence could be evaluated in the same manner as the influence of variability in the index flood relation has been considered in section 3. It may modify partly the judgment on the relative efficiency of the various tested regional flood frequency analysis methods.

Finally, it is important to keep the conclusion of section 4 in mind. All the regional flood frequency methods are based on homogeneity hypotheses that can hardly be verified. It is therefore certainly sound to test the robustness of the obtained inference results to different hypotheses and region delineations: to combine various methods. In this perspective, the proposed approach does not replace but enriches the available panel of methods and possible hypotheses. It may lead, as shown herein, to quantile estimates that are theoretically equally or even more accurate than the estimates obtained with standard methods. Its main constraint may also be considered as an advantage: the two necessary homogeneity hypotheses limit the extent of homogeneous regions that may be, on the other hand, hydrologically more consistent and less questionable than in standard regional flood frequency analyses.

old value to ensure exhaustiveness, even if this limits the information content of the inventories of large or extreme floods. For this reason, the collation period has been limited in the Ardèche to the last 50 years (typical depth of local memory). Likewise, the inventory has been limited to the largest known flood event at each ungauged site. This last choice facilitates also the computation of the plotting positions for the ungauged extremes, which can become tricky: the return periods of the ungauged extremes reported on the cumulated distribution plots are estimated through Monte Carlo simulations [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF]. 47

  problems may be solved by tunning the MCMC algorithm: better selection of the parameter starting values, adjustment of the Monte Carlo Markov chain random search controling factors. The tests presented hereafter implied repeated applications of the MCMC algorithm on a large number of data samples. It was therefore necessary to strengthen the algorithm to ensure rapid convergence in most cases, to verify automatically the convergence of the MCMC algorithm for each run and to avoid introducing results corresponding to insufficient convergence or convergence failures in the analysis.

(

  ii) Reliability of the computed credibility intervals.Monte Carlo simulations enable to test the relevance of the posterior distributions of the discharge quantile Qi (0.99) obtained through the Bayesian MCMC procedure. The real quantile value Q i (0.99) is known, and its probability of non-exceedance F (Q i (0.99)) can be computed according to the computed posterior distribution of the discharge quantile Qi (0.99) for each simulated sample. If these posterior distributions of the estimated discharge quantiles Qi (0.99) are unbiased -i.e. reflect the correct density of probability of the real quantile Q i (0.99), the values of F (Q i (0.99)) estimated for the 1000 samples should be uniformly distributed over [0, 1]. Figure3illustrates how the distributions of F (Q i (0.99)) should be interpreted. For instance in case (c), the real quantile value has an estimated probability of non-exceedance close to 0 or 1 for too many cases : the real quantile lies too often at the margins of the computed posterior distribution for the estimated quantiles, sign that the variance of the computed distribution is too low and that the computed intervals are too narrow.Likewise, case (d) illustrates the case of an overestimation of the variance of the posterior distributions and resulting too large credibility intervals. Case (b) reveals an underestimation of the quantiles. (iii) Introduction of regional variability in the index flood relation. The Monte Carlo simulations have been repeated with the introduction of a regional variability in the index flood relation to test the robustness of the new proposed method. For each generated sample, the index flood value at site i, µ i = S β i was replaced by a random value drawn from a lognormal distribution with mean S β i and standard deviation δ * S β i , δ being successively set to 0.1 and 0.3. Statistical tests indicate that values of δ equal to 0 or 0.1 are compatible with the gauged sample available in the Ardèche region.

Figure 4 .

 4 Figure 4.a and b appear overall very similar except for the regional growth curves: the maximum likelihood estimates and the width of computed credibility intervals fluctuate in similar ranges, particularly for the larger catchment areas. The uncertainties increase logically in the proposed method for smaller catchment areas, because the sensitivity of the index flood to the value of the parameter β is increasing with the watershed area according to equation 2. Conversely, the smallest areas have a lower weight in the determination of the index flood relation and as a consequence the calibrated index flood values are less accurate for lower surfaces. The regional growth curves have different significations in both methods: distribution of the reduced discharges (values divided by the expectancy of the distribution) in the method of Hosking & Wallis and distribution of the discharges for a 1 km 2 watershed in the proposed approach. This is largely in the extrapolation range of the index flood relation for low values of watershed

Figure 8

 8 Figure 8 and table 3 summarize the results obtained for the Ardèche case study. They should be interpreted at the light of the conclusions drawn in the previous section. The presented results correspond to the two extreme gauged catchment surfaces in the region: 63 km 2 (Borne watershed at Puylaurens), and 2240 km 2 (close to the outlet of the entire Ardèche watershed at Saint Martin).

Figure 9

 9 Figure 9 compares the modelled relation between average discharges and watershed areas and the observed locally estimated mean discharges at the light of the proposed Lognormal model to reproduce regional variability in the index flood relationship. This illustrates the consistency of the proposed index flood relation with the observed data. The second smallest watershed,

  Ardèche region may appear as particularly rich if compared to what could reasonably be obtained in other regions. Previous works, conducted on the usefulness of historic records in flood frequency studies, indicated that the accuracy level of the estimated extreme discharges and the number of extremes for which an estimated discharge is available, have little influence on the inference results when the historic period of record if fixed Payrastre

  the computation of the credibility intervals are based on the assumption of independence between the various recorded maximum annual peak discharges. But, extreme peak discharges induced by the same rainfall event sometimes occur at the same date at nearby ungauged sites, suggesting at least weak depedencies in the inventory of ungauged extremes. To prevent any effect of these dependencies on the results and especially to limit the risk of underestimation of the credibility intervals due to an implicite overrating of the information content of the ungauged extreme inventory neglecting depedencies, the data set had to be refined. When the largest reported flood occurred at the same date at close ungauged sites, one randomly selected site has been kept in the final ungauged data set. Keeping the largest value would have introduced a sampling bias. The refinement of the Ardèche ungauged extreme set led to eliminate more than half of the records initially used by Gaume et al. (2010), hence resulting in a slight increase only of the computed credibility intervals. Posterior distributions of F (Q i (0.99)): (a) perfect credibility intervals, (b) tendency to over-estimate the quantile, (c) too narrow estimated credibility intervals, (d) too large estimated credibility intervals. . . . . . . . . . . . . . . . . . . . . . . . . 4 Estimates of the 100-year quantile through Bayesian MCMC regional flood frequency inferences conducted on 1000 different samples: (a) Method of Hosking & Wallis and (b) proposed method. Box plots of the 1000 maximum likelihood and 5% and 95% credibility bounds for each considered watershed area. All the values have been divided by the real 100-year quantiles Qi (0.99). . . . . . . . . . . . . . . . . . . . 5 Distributions of the values of F (Q i (0.99)) estimated for the 1000 simulated samples: (a,b) Method of Hosking & Wallis, (c,d) Hosking& Wallis with correction factor m for the posterior distributions and (e,f) proposed method. . . . . . . . . . . 6 Results of Bayesian MCMC regional flood frequency inferences conducted on 1000 samples with the proposed method when a regional variability in the index flood relation is introduced: (a) δ = 0.1, (b) δ = 0.3. . . . . . . . . . . . . . . . . . . . . . . 7 The Ardèche region, location of the 5 gauged sites (triangles) and the 18 sites where ungauged extremes were retrieved (dots). Results of the Bayesian MCMC regional flood frequency inferences conducted on the Ardèche dataset. Results provided for the Ardèche at Saint Martin (2240 km 2 ): (a) Method of Hosking & Wallis and (b) proposed approach including ungauged extremes. Statistical distributions corresponding to the maximum likelihood (continuous line) and 5% and 95% credibility bounds (dotted lines). Gauged values (crosses) and ungauged extremes (dots). . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Calibrated relation between average discharges and watershed areas, sample mean discharges and variation ranges corresponding to the tested lognormal random models. . . . . . . . 10 Maximum likelihood estimate and 5% and 95% bounds of the computed credibility intervals for 4 tested methods on 1000 simulated samples having the same characteristics as the Ardèche data sample (St Martin case): (a) method of Hosking and Wallis with L-moments, (b) method of Hosking and Wallis used in this paper, (c) proposed method applied on gauged series only, (d) proposed method with 18 ungauged extremes. .

Figure 1 :Figure 2 :

 12 Figure 1: Posterior distributions for the parameters α and β. Calibration of a regional growth curve and an index flood relation, based on a simulated sample of 1000 records: α = 2.24, ξ = 3.34, k = -0.16 and β = 0.76. (a) totally biased result with the standard likelihood expression 3 and (b) correct result with likelihood expression 7.

Figure 3 :Figure 4 :Figure 5 :

 345 Figure

Figure 7 :

 7 Figure 7: The Ardèche region, location of the 5 gauged sites (triangles) and the 18 sites where ungauged extremes were retrieved (dots).

Figure 9 :

 9 Figure 9: Calibrated relation between average discharges and watershed areas, sample mean discharges and variation ranges corresponding to the tested lognormal random models.

Figure 10 :

 10 Figure 10: Maximum likelihood estimate and 5% and 95% bounds of the computed credibility intervals for 4 tested methods on 1000 simulated samples having the same characteristics as the Ardèche data sample (St Martin case): (a) method of Hosking and Wallis with L-moments, (b) method of Hosking and Wallis used in this paper, (c) proposed method applied on gauged series only, (d) proposed method with 18 ungauged extremes.

  . The region includes six main rivers: the Ardèche River 2380 km 2 , Ouvèze River 120 km 2 , Eyrieux River 860 km 2 , Doux River 630 km 2 , Ay River 110 km 2 and Cance River 380 km 2 . Five stream gauges provide accurate, continuous and long series of discharges in the region. The gauged dataset of annual maximum peak discharges counts 168 records. The situation of the 5 gauging stations is reported in figure ?? and the characteristics of corresponding samples are presented in table 1.

Table 1 :

 1 The sample size and the surface of the gauging stations in the regional Ardèche.

	Station	River name	Period of data	N(year)	S (km 2 )
	Beauvene	Eyrieux	1969-1998	30	392
	Vogüe	Ardèche	1966-2003	38	636
	Chambonas	Chassezac	1971-1997	27	507
	St.Laurent	Borne	1968-1997	30	63
	St.Martin	Ardèche	1963-2005	43	2240

Table 3 :

 3 Comparison the credibility intervals (CI) at St. Laurent and St.Martin gauging stations to three cases corresponding to return periods T=100 years.

	Surface	N cont	N hist Case	Q5 i (0.99)	QML i	(0.99)	Q95 i (0.99)
	(km 2 ) (years) (years)		(m 3 /s)	(m 3 /s)	(m 3 /s)
		168	0	(a)	423.8		520.9	712.5
	63	168	0	(b)	382.3		479.5	682.1
		166	900	(c)	442.8		500.4	574.2
		168	0	(a)	5538.6	6806.8	9310.8
	2240	168	0	(b)	5892.4	7602.2	10747.6
		166	900	(c)	5477.3	6560.0	7847.1

Acknowledgements

This work has benefited from the support of the French Ministry of foreign affairs.

Appendix A. Method for the evaluation of the underestimation factor for the computed credibility intervals In many cases, the proposed reliability test for the computed credibility intervals presented in section 3.1 indicates that the intervals are too narrow (Fig. 5). This is due the fact that some sources of uncertainties are not accounted for in the bayesian MCMC framework used to compute these intervals. The following procedure is proposed to evaluate the underestimation factor m of the computed intervals. Based on the hypothesis that the shapes of the posterior distribution of the discharge quantiles Qi (0.99) are correct and that only their dispersion is underestimated, the procedure to evaluate the underestimation factor m is the following.

The distance between the real quantile value Q i (0.99) and the median value of the computed posterior distribution M e Qi (0.99) is artificialy reduced by the factor m : the probability of non-exceedance Fm (Q i (0.99)) (eq. A.1) is computed rather than the probability F (Q i (0.99)) based on the MCMC results for each simulated sample.

The factor m leading to the distribution of Fm (Q i (0.99)) closest to a uniform distribution for the 1000 simulated samples, gives an estimate of the underestimation factor of the credibility intervals.

Appendix B. Sampling ungauged extremes in a region

The regional samples of ungauged extremes must be carrefully built to avoid introducing bias into the statistical analysis. At least, the extreme data collation method must seek to meet three conditions: (i) be based on an explicit sampling strategy that can be translated into the likelihood formulation, (ii) be comprehensive fo each selected location, (iii) limit possible dependencies between records and samples. The next three sections explain how these three conditions were accounted for to build the Ardèche regional sample of ungauged extreme floods.

Appendix B.1. Explicit and systematic sampling method for regional extremes Isolated extreme discharge values are often reported but can not be directly integrated in regional flood frequency analyses if their representativeness is unknown. Over which period of time can it be asserted that a larger flood event did not occur at the same site ? And, even more crucial, these values have been reported because they appear as extreme, i.e. this means implicitly that they seem large if compared the known flood events observed in neighbouring comparable catchments. Over how many years and for how many similar catchments (station) could it be finally asserted that the considered extreme discharge value -or rescalled discharge value -has not been exceeded ? Large uncertainties concerning the record length (station.years) corresponding to each extreme may drastically reduce their added value in regional flood frequency analyses as illustrated by the Slovakian test case in [START_REF] Gaume | Bayesian MCMC approach to regional flood frequency analyses involving extraordinary flood events at ungauged sites[END_REF]. An inventory of regional extremes should therefore be based, if possible, on a systematic data mining procedure. The various sites in a region where information on extremes may be available should be determined a priori: typically significant towns along rivers where large flloods may have induced damages and may have been reported. An inventory of past large floods should be realized at each of these sites based on local memory and archives and it must at least be verified that some given threshold values have or have not been exceeded at the considered sites in the recent past: non-exceedance is a statistically valuable information.

Appendix B.2. Comprehensiveness

Collation of extreme events over a given period of time provide censored data sets: only discharges exceeding a given threshold, that depends on the site and the considered period, are documented. Flood events exceeding the threshold should be exhaustively listed at each site over the considered past period for a correct statitical inference. To be on the safe side, it is generally preferable to limit the considered past period of time and to raise the thresh-List of Tables 1 The sample size and the surface of the gauging stations in the


Jan Van Eijck

Elements of Epistemic Crypto Logic

The talk presents an extension of DEL (dynamic epistemic logic) intended for model checking of cryptographic protocols. Key elements are a feasible epistemic representation of knowledge of large integers, using register models, and exchange of such knowledge over a network. I will demonstrate how the approach can be used for model checking Diffie-Helman key exchange and similar protocols.

The Number Guessing Game, Naive Version Scene: Gaia and Rosa are with their mother Heleen, and Gaia and Rosa quarrel about who can play with a toy.

Heleen proposes to play the number guessing game. She says: "I have a number in mind, in the range from one to ten. You may take turns guessing. Whoever guesses the number first gets the toy. It is Gaia's turn to start, for last time we played this game Rosa started the guessing." They all agree, and after a number of rounds Heleen announces: "Rosa, you have won." The Number Guessing Game, Streetwise Version Scene: Gaia and Rosa are with their mother Heleen, and Gaia and Rosa quarrel about who can play with a toy.

Heleen proposes to play the number guessing game. She says: "I have a number in mind, in the range from one to to ten. You may take turns guessing. Whoever guesses the number first gets the toy. It is Rosa's turn to start, for last time we played this game Gaia started the guessing."

Gaia does not agree. "How can we know you are not cheating on us? Please write down the number, so you can show it to us afterwards as a proof."

Now they all agree, and after a number of rounds Heleen announces: "Rosa, you have won." She shows the paper with the proof.

Knowing a Number

At the point where either Gaia or Rosa demands that the secret number gets written down, and that Heleen shows it as a proof that the procedure was honest, the important notion of a register arises.

The register allows Heleen to prove that she knew (had fixed) the number beforehand.

So what is it that Heleen knew when we say that she knew the number?

What Does it Mean to Know a Number?

We are working in the context of epistemic logic.

To know an integer number n is to be able distinguish a world where n is written in some register from all possible worlds where an integer number n different from n is written in that register.

This can be made more precise, by clearly distinguishing between a register and its contents.

And the representation can be made more concise, by lumping together all situations where the register does not contain the number.

To know an integer number n is to know the contents ṅ of register n, not to know a number is not to have access to register n.

Heleen knows a number

Mother knows n, the two children do not know n. We use n for the register, ṅ for the value of that register, and n for the domain (set of possible values) of the register.

This leads to the following register model.

0 : n 1 : n n → 6 n → {1..5, 7..10}
0 is a world that has register n, with number ṅ = 6 stored in it, while 1 differs from 0 in that it also has this register, but with all the possible values in it that differ from ṅ.

Gaia (dots) and Rosa (dashes) do not have access to the register.

Gaia Makes a Guess

Use g for the guess of Gaia. Gaia knows her guess, but has not yet revealed it to the others: solid lines for Heleen, dotted lines for Gaia, dashed lines for Rosa.

0 : ng 1 : ng 2 : ng 3 : ng n → 6, n → {1..5, 7..10} g → 5, g → {1..4, 6..10}
Without registers, this blows up to a model with 100 worlds.

k registers, each of m bits, blow up to (2 m ) k possibilities.

After Gaia has Revealed the Contents of her Register

Gaia reveals the contents of her register: g = 5.

Gaia states her guess: n = g.

0 : ng 1 : ng n → 6, n → {1..5, 7..10} g → 5

After Heleen's Reply

Since ṅ = ġ Heleen replies "you got it wrong", and we get:

0 : ng 1 : ng n → 6, n → {1..4, 7..10} g → 5

Truth of Equality Statements in Register Models

• It seems reasonable to represent the announcement of Gaia's guess as the announcement of the equality statement g = ġ.

• Since g means that the register does contain a different value from ġ, the statement g = ġ is false in both g worlds, so these drop out of the picture.

• Heleen can respond to Gaia's guess with n = g (this is what "you guessed it" means: your guess equals my number).

• n = g is different from n = ṅ.

Truth of (In)equality Statements in Register Models

• The difference between n = g and n = ṅ can perhaps be seen more clearly in the negative case.

• "You have it wrong" is paraprased as n = g, not as n = ṅ.

• In a world where g and n are both true, the equality n = g can be checked by checking whether ġ = ṅ.

• In a world where g is true and n false, the equality n = g means that every value of n except for ṅ equals g, which is certainly false (except in a domain with just two numbers, in case ṅ and ġ are different).

• In a world where g and n are both false the equality n = g means that every value for g except ġ equals every value of n except ṅ, which is certainly false.

Register Language for Guessing Games

Let p range over P and let N range over Z. Let i range over the set of agents I.

ϕ

::= | p | L i | p = E | ¬ϕ | ϕ ∧ ϕ | K i ϕ | [A]ϕ E ::= p | N
A is some update operation (see below). L i expresses that agent i is listening to announcements.

Let G ⊆ I. L G expresses that G are the listeners:

L G = i∈G L i ∧ i / ∈G ¬L i .

Update Operations

Revealing positive or negative information about p:

! p = E "You guessed it" ! p = E
"Your guess was wrong" Interpretations of ! p = E and ! p = E are action models in the sense of [START_REF] Baltag | The logic of public announcements, common knowledge, and private suspicions[END_REF].

Register Models

Let I be a finite agent set, let P be a denumerable set of propositions.

Definition 1 A register model for I and P is a tuple

M = (W, R, V, L, V , V )
where

• (W, R, V ) is a multimodal S5 model for I and P ∪ {L i | i ∈ I}, • Q = ( w∈W V w ) ∩ P is finite, • L : W → P(I), • V : W → Q → Z, • V : W → Q → P fin (Z), subject to constraints (below).

Register Models (2)

Q is the part of the vocabulary that is currently used (compare [START_REF] Van Eijck | Composing models[END_REF]). Call this the active vocabulary of the model.

L : W → P(I) is a function that indicates which agents are listening to broadcasts in world w.

The function V maps each w ∈ W to a function Vw that maps all members of Q to integers.

The function V maps each w ∈ W to a function Vw that maps all members of Q to interpretation domains (finite sets of integers).

Constraints:

• for all w ∈ W, q ∈ Q: Vw (q) ∈ Vw (q).

• for all w ∈ W, q ∈ Q: 2 ≤ | Vw (q)| < ∞.

Awareness Constraint

Agents are aware of their state of attention ([DHLoS13]):

• If L i ∈ V w then (w, w ) ∈ R i implies L i ∈ V w . • If L i / ∈ V w then (w, w ) ∈ R i implies L i / ∈ V w .
In register models satisfying this constraint, the following formulas are true in every world:

L i ↔ K i L i . ¬L i ↔ K i ¬L i .

Interpretation in Register Models

Definition 2 (Register Assignments) Let a register model M with world set W , active vocabulary Q and register function V be given.

The set of M-register assignments H ⊆ Z Q is the set of all functions h : Q → Z such that h(q) ∈ w∈W Vw (q) for all q ∈ Q.

Definition 3 (Agreement) h ∈ H agrees with w (notation w h) if it holds for all q ∈ Q that h(q) ∈ Vw (q) and moreover, q ∈ V w iff h(q) = Vw (q).

For w ∈ W let [ [w] ] M be the set {h ∈ H | w h}.
We will use these notions to define M, w, h |= ϕ for h agreeing with w by recursion.

Examples of (Dis-)Agreement

0 : ng 1 : ng n → 6, n → {1..4, 7..10} g → 5
The function h = {n → 6, g → 5} agrees with 0, but not with 1.

The function h = {n → 3, g → 5} agrees with 1 but not with 0.

Interpretation in Register Models: M, w, h |= ϕ

Let h be a register assignment that agrees with w.

M, w, h |= always M, w, h |= p iff p ∈ V (w) M, w, h |= L i iff i ∈ L(w) M, w, h |= p 1 = p 2 iff h(p 1 ) = h(p 2 ) M, w, h |= p = N iff h(p) = N M, w, h |= ¬ϕ iff not M, w, h |= ϕ M, w, h |= ϕ 1 ∧ ϕ 2 iff M, w, h |= ϕ 1 and M, w, h |= ϕ 2 M, w, h |= K i ϕ iff (w, w ) ∈ R i and w h imply M, w , h |= ϕ M, w, h |= [! p = E]ϕ iff M, w, h |= p = E implies M , (w, e), h |= ϕ where h = λ(w, e) → h(w) M = M • A A = e : p = E, e : , e I-G ←→ e • is product update with adjustment of V G is such that M, w, h |= L G M, w, h |= [! p = E]ϕ 2 iff . . .

Truth and Falsity at a World

From truth at a register assignment to truth and falsity at a world:

M, w |= ϕ iff

Method:

• Randomly generate a list of n agreeing assignments h 1 , . . . , h n for w.

• Check M, w, h i |= ϕ for each h i in the list.

• The probability that there is disagreement among the outcomes can be made arbitrarily small.

Truth Value Gaps Almost Closed

0 : n → 6 1 : n → {M = -2 63 ..2 63 -1} -{6} 1 18446744073709551615 n = 7 is false in 0. n = 7 is almost false in 1. n = 7 is true in 0. n = 7 is almost true in 1.
Article of faith of cryptanalysis: "Brute force attacks on number secrets are impossible."

Register Creation

Fix a large register size. Let M..K be the bounds of the register (say M = -2 63 , K = 2 63 -1).

Add to the language:

ϕ ::= [p i ← N ]ϕ p i
← N is the command to link p to the integer number N with the link known only to agent i. It is assumed that N fits the register size. • ?K a ϕ tests whether K a ϕ is true in the actual world.

• Open b adds b to the set of agents that are listening.

• !ϕ is the announcement of ϕ.

• Close b removes b from the set of agents that are listening.

• As far as b is concerned, the information ϕ could come from anywhere: no authentication.

• Anyone who is listening receives the info: the channel is insecure.

Feasible Computation

• Fast algorithms for primality testing (e.g., the probabilistic Miller-Rabin test [START_REF] Gary | Riemann's hypothesis and tests for primality[END_REF][START_REF] Michael | Probabilistic algorithm for testing primality[END_REF]).

• Fast algorithms for co-primality testing (Euclid's GCD algorithm).

• Fast algorithms for finding modular inverses (Euclid's extended GCD algorithm).

• Fast algorithms for addition and multiplication modulo, and for exponentiation modulo (see, e.g., [START_REF] Delfs | Introduction To Cryptography -Principles and Applications[END_REF][START_REF] Paar | Understanding Cryptography -A Textbook for Students and Practitioners[END_REF]).

• + Articles of faith: factorisation, discrete logarithm, . . . , are not feasible.

Extending the Register Language with Protocols

Let a countable set P of basic proposition letters and a finite set I of agents be given. Let p range over P and i over I. Let N range over Z.

Let f ( p) express any feasible computation from natural numbers named by p.

ϕ ::= | p | L i | p = E | ¬ϕ | ϕ ∧ ϕ | K i ϕ | [A]ϕ E ::= p | N A ::= Open i | Close i | !ϕ | ?ϕ | p i ← N | p i ← f ( p) | A; A | if ϕ then A else A

Informal Explanation

• We will use a broadcast mechanism for communication. L i expresses that agent i is listening to broadcasts.

• Open i adds i to the set of agents who are listening to broadcasts, while Close i deletes agent i from that set.

• !ϕ is the command to broadcast ϕ, while ?ϕ is the command to check whether ϕ is true.

• p i ← N is the command to link p to number N , with the link known only to agent i. p i ← f ( p) is the command to link p to the number computed by f from p, with the link known only to agent i. It is assumed that i knows all of p. Example:

Prelude> [ aˆ28 'mod' 29 | a <-[1..28] ] [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]

Explanation

We have: • Formulate complete logics. A complete calculus for register logic has axioms for propositional logic, S5 axioms for the K i operators, K axioms for the A operators, Modus Ponens, necessitation for K i and A, and reduction axioms for the command operators, in the style of [START_REF] Van Benthem | Logics of communication and change[END_REF].

(g a ) (p-1)-b = g a((p-1)-b) = g a(p-1) × g -ab = (g p-1 ) a × g -ab F ermat = 1 a × g -ab = g -ab mod p,
• Build an efficient model checker for this language, using register models as defined in this talk. See [START_REF] Gattinger | Epistemic crypto logic -functional programming and model checking of cryptographic protocols[END_REF] for a first prototype.

• Use the model checker for modelling cryptographic protocols, plus attacks on them.

Related Work

[Gat13]: Implementation in Haskell of (a slight variation on) the model checking algorithm presented in this talk.

[MS04, EO05]: Epistemic logic approach to Dining Cryptographers protocol.

[ [START_REF] Dechesne | Dynamic epistemic verification of security protocols: framework and case study[END_REF][START_REF] Wang | Epistemic Modelling and Protocol Dynamics[END_REF] Comparison between approaches based on dynamic epistemic logic and temporal epistemic logic.

[DHLoS13] Earlier approach to model checking cryptographic protocols with DEMO.

[Kra07, CD07] Intricate logics for cryptography, with lots of primitive notions.

[MW07] Epistemic analysis of the notion of 'secure channel'.

Overview•

  Epistemic Representation of Number Guessing Games • Register Models for Number Guessing Games • Language for Number Guessing Games • Interpretation in Register Models • Truth Value Gaps and What to Do About Them • Probabilities and Monte Carlo Style Model Checking • Register Creation as an Update • Communication: "a sends ϕ to b" • Feasible Computation; Protocols; Secret Key Exchange • Program for Epistemic Crypto Logic

  Interpretation: [ [p i ← N ] ] as the action model (in the sense of [BMS98]): p := {M..K}, ṗ := N, M ≤ N ≤ K, e : p := , e : p := ⊥, e n → 6, n → {M..K} -{6} g → 5, g → {M..K} -{5} Communication: "a sends ϕ to b" Sequence of commands: ?K a ϕ; Open b; !ϕ; Close b

  Example: Diffie-Hellman Key Exchange[START_REF] Diffie | New directions in cryp-tography[END_REF] Whitfield Diffie -Martin Hellman Diffie-Hellman Key Exchange Protocol Key Exchange Over Insecure Channel 1. Alice and Bob agree on a large prime p and a base g < p such that g and p -1 are co-prime.2. Alice picks a secret a and sends g a mod p = A to Bob. 3. Bob picks a secret b and sends g b mod p = B to Alice. 4. Alice calculates k = B a mod p. 5. Bob calculates k = A b mod p. 6. They now have a shared key k. This is because k = (g a ) b = (g b ) a mod p. Use of a Shared Secret Key for Secure Communication Let p be the prime that Alice and Bob have agreed on, and let k be their shared key. Then message m is encoded as m × k mod p. Such messages can be decoded by both Alice and Bob. Alice knows p, k, g b and a. She decodes cipher c with c × (g b ) (p-1)-a mod p. Bob knows p, k, g a and b. He decodes cipher c with c × (g a ) (p-1)-b mod p. Behind this: Fermat's Little Theorem (1601-1665) Theorem 4 (Fermat) If p is prime, then for every 1 ≤ a < p: a p-1 ≡ 1 (mod p).

  and therefore: c × (g a ) (p-1)-b = (m × g ab ) × g -ab = m × (g ab × g -ab ) = m mod p.Diffie-Hellman Key Exchange as a Register Language Protocola i ← N ; A i ← g a mod p ; Open j ; !x = A ; Close j ; b j ← M ; B j ← g b mod p ; Open i ; !y = B ; Close i ; k i ← y a mod p ; k j ← x b mod p ; ?k = kModelling Attacks on ProtocolsThe Dolev-Yao model of attacks on public key security protocols[START_REF] Dolev | On the security of public key protocols[END_REF] assumes that the attacker can manipulate the network.We can model such attacks as interleavings of the security protocol with an attack protocol.Example attack on Diffie-Hellman key exchange:Open e ; a i ← N ; A i ← g a mod p ; Open j ; Close i ; Close j ; !x = A ; Close j ; c e ← K ; C e ← g c mod p ; Open j ; !x = C ; Close j ; b j ← M ; B j ← g b mod p ; Open i ; !y = B ; Close i ; k i ← y a mod p ; k j ← x b mod p ; ?k = k Program for Epistemic Crypto Logic

  ∀h with w h : M, w, h |= ϕ.

	Truth Value Gaps K i closes the truth value gap Update Example
	M, w =| K i ϕ iff
	∀h with w iff	0 : n h : M, w, h |= ¬K i ϕ 1 : n 0 : n, L g 1 : n, L g	n → 6 n → {1..5, 7..10} n → 6 n → {1..5, 7..10}
	∀h with w	h∃w , h with wR i w , w	h and M, w , h |= ¬ϕ
	n = 7 is false in 0. iff
	∃h with w	h∃w , h with wR i w , w	h and M, w , h |= ¬ϕ
	iff
	∃h with w	h and M, w, h |= K i ϕ
	iff
	M, w |= K i ϕ.

M, w =| ϕ iff ∀h with w h : M, w, h |= ¬ϕ. n = 7 is neither true nor false in 1. n = 7 is true in 0. n = 7 is neither true nor false in 1. Suppose n = 5 is announced. Since Gaia is listening, we get: (0, e) : n → 6, L g (1, e) : n → {1..4, 7..10}, L g (0, e ) : n → 6, L g (1, e ) : n → {1..5, 7..10}, L g M, w |= ϕ iff for "enough" h with w h : M, w, h |= ϕ.

Numbers of agreeing assignments for register size 10.

Numbers of agreeing assignments for register size of 32 bits.

With 64 bit size: (2 64 -1) 2 = 340282366920938463426481119284349108225 agreeing assignments for world 3.

Monte Carlo Checking of (In)Equality Statements

Suppose the register size is large enough, say 64 bits.

Then in a world where n is false, there are

possibilities for the value of n.

Thus, in a world where n is false, equality statements n = X will be false for almost all values X.

Equality and inequality statements can be checked by means an approximate (Monte Carlo) model checking algorithm.


Axel Müller
email: axel.mueller@supelec.fr

Emil Björnson
email: emil.bjornson@ee.kth.se

Romain Couillet
email: romain.couillet@supelec.fr

Mérouane Debbah
email: merouane.debbah@supelec.fr

Analysis and Management of Heterogeneous User Mobility in Large-scale Downlink Systems

Modern cellular networks need to serve user terminals with large disparities in mobility, which incurs different accuracy of the channel state information for each user. The impact of such heterogeneous mobility on the multi-cell downlink is analyzed in this paper. The base stations serve a multitude of users by coordinated beamforming. We derive deterministic equivalents for the user performance in a large scale system where the number of transmit antennas and user terminals grow large at a fixed ratio. We show that low and high mobility users can coexist and be served simultaneously, since the CSI imperfections of a user only harms the performance of this particular user. Simulations are used to verify the applicability of our large scale approximations for systems of practical dimensions. Furthermore, we show that the performance of high mobility users can be improved by explicitly managing the user priorities in the network.

I. INTRODUCTION

In order to keep up with the exponential traffic growth in cellular networks, the infrastructure is becoming increasingly dense [START_REF] Hoydis | Green small-cell networks[END_REF]; the conventional macro base stations (BSs) are supplemented with small-cell BSs that act as hotspots inside the macro cells. The terminology of heterogeneous networks (HetNets) is commonly used to refer to these systems [START_REF] Baldemair | Evolving wireless communications: Addressing the challenges and expectations of the future[END_REF] and some of the key practical challenges are interference management and load balancing across BSs.

It is not only the infrastructure that is heterogeneous in modern networks, but also the user conditions and particularly the user mobility. While the conventional cellular networks, also known as mobile networks, were mainly deployed for outdoor usage, we can already see a large disparity in use cases: low mobility indoor users coexist with medium and high mobility outdoor users. The mobility is closely related to the quality of the channel state information (CSI) that can be acquired at the BSs; high mobility implies low-quality CSI.

In conventional networks with orthogonal allocation of time/frequency resources, heterogeneous user mobility mainly affects the scheduling decisions. However, modern cellular networks support multi-antenna transmission where a multitude of users are served simultaneously using the same time/frequency resources; inter-user interference is then controlled by spatial precoding [START_REF] Gesbert | Shifting the MIMO paradigm[END_REF]- [START_REF] Björnson | Optimal resource allocation in coordinated multi-cell systems[END_REF]. Since CSI is a prerequisite for effective spatial interference suppression [START_REF] Gesbert | Shifting the MIMO paradigm[END_REF], this raises the question of how to take the user mobility into account in the precoding design.

The majority of the existing literature on precoding with imperfect CSI considers homogeneous CSI quality among the users, which is not the case in HetNets. A notable exception is [START_REF] Shirani-Mehr | MIMO downlink scheduling with non-perfect channel state knowledge[END_REF] where the authors considered joint scheduling and precoding in the downlink of a single-cell system. The authors concluded that only low mobility users can be served jointly, while orthogonal transmission resources should be allocated to high mobility users. A similar result was obtained in [START_REF] De Kerret | Degrees of freedom of the network MIMO channel with distributed CSI[END_REF], where a cellular network with joint transmission between the BSs was considered. The analysis revealed that the multiplexing gain is limited by the high mobility users, thus it appears as if all other users suffer when coexisting with high mobility users.

In this paper, we consider a cellular network where the users in each cell have different mobility, which is a relevant distinction in HetNets where the high mobility users are served by macro BSs and the low mobility users are associated with the small-cell BSs. Different from [START_REF] De Kerret | Degrees of freedom of the network MIMO channel with distributed CSI[END_REF], we consider coordinated beamforming where each user is served by one BS but the precoding is coordinated across the cells to suppress inter-cell interference. To obtain fundamental insights, we consider the large-system regime where the number of transmit antennas and user terminals go to infinity at a fixed ratio. As a main contribution, we derive deterministic expressions for the asymptotic user rates, which also serve as accurate approximations in practical non-asymptotic regimes. These novel expressions generalize the prior work in [START_REF] Wagner | Large system analysis of linear precoding in MISO broadcast channels with limited feedback[END_REF] for single-cell systems and in [START_REF] Hoydis | Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?[END_REF] for multi-cell systems where only deterministic statistical CSI is utilized for suppression of inter-cell interference. To account for the different user mobility, the precoding is a simplified version of the optimal precoding parametrization in [START_REF] Björnson | Pareto characterization of the multicell MIMO performance region with simple receivers[END_REF]; more precisely, different priority weights are assigned to cells with different mobility.

Our analysis reveals that the mobility of a user has a detrimental impact on its own achievable rate, but has no direct impact on the other users. In particular, the low mobility users will not suffer from coexisting with high mobility users under coordinated beamforming-which is fundamentally different from the case of joint multi-cell transmission in [START_REF] De Kerret | Degrees of freedom of the network MIMO channel with distributed CSI[END_REF]. Furthermore, the heterogeneity can be managed in our framework by explicitly increasing/decreasing the priority of high mobility users, which changes their rates at the cost/gain of other users.

II. SYSTEM MODEL

In the following, we analyze cellular downlink multi-user multiple input multiple output (MU-MIMO) systems of the kind illustrated in Fig. 1. Each of the L cells consists of one base station (BS) associated with a number of single antenna user terminals (UTs). In more detail, the lth BS is equipped with N l antennas and serves K l UTs. We generally set N l ≥ K l in order to avoid scheduling complications. We assume transmission on a single narrow-band carrier, usage of a perfectly synchronized TDD protocol, full transmit-buffers, and universal frequency reuse among the cells.

The lth BS transmits a data symbol vector x l = [x l,1 , . . . , x l,K l ] T intended for its K l uniquely associated UTs. This BS uses the linear precoding matrix G l l ∈ C N l ×K l , where the columns g l l,k ∈ C N l constitute the precoding vectors for each UT. 1 We note that BSs do not directly interact with each other and users from other cells are explicitly not served. Thus, the received signal y l,k ∈ C at the kth UT in cell l is

y l,k = χ l l,k (h l l,k ) H g l l,k x l,k + k =k χ l l,k (h l l,k ) H g l l,k x l,k + m =l χ m l,k (h m l,k ) H G m m x m + n l,k (1)
where n l,k ∼ CN (0, 1) is the independent circularly symmetric complex Gaussian noise. The transmission symbols are chosen from a Gaussian codebook, i.e., x l,k ∼ CN (0, 1). We assume block-wise small scale Rayleigh fading, thus the channel vectors are modeled as h m l,k ∼ CN (0, 1 Nm I Nm ). The path-loss and other large-scale fading effects are incorporated in the χ m l,k factors. The scaling factor 1 Nm in the fading variances is of technical nature and utilized in the asymptotic analysis. It is canceled for a given arbitrarily sized system by modifying the transmission power accordingly.

A. Channel State Information and Effect of User Mobility

The UTs are assumed to perfectly estimate the respective channels to their serving BS, which enables coherent reception. This is reasonable at most user speeds, if proper downlink reference signals are alternated with data symbols.

Concerning the downlink, the transmitters in our model can only obtain imperfect channel state information (CSIT), for example due to imperfect pilot-based channel estimation, delays in the acquisition protocols, and user mobility. High user mobility naturally implies large CSI imperfections, because the coherence time of the channel diminishes with the speed. If the system allocates a constant fraction of the coherence time for pilot signaling, the accuracy of the channel estimation reduces with a shorter coherence time. If the system, on the other hand, keeps the pilot resources fixed, then a higher mobility induces outdated CSI.

To model imperfect CSIT without making explicit assumptions on the acquisition protocol, we define the estimated channel vectors ĥm l,k ∈ C Nm to be

ĥm l,k = χ m l,k (1 -(τ m l ) 2 )h m l,k + τ m l hm l,k (2)
where hm l,k ∼ CN (0, 1 Nm I Nm ) is the normalized independent estimation error. Using this formulation, we can set the accuracy of the channel acquisition between the UTs of cell l and the BS of cell m by selecting τ m l ∈ [0, 1]; a small number on τ m l implies low mobility and vice versa. Furthermore, we remark that these choices imply ĥm

l,k ∼ CN (0, χ m l,k 1
Nm I Nm ). For convenience later on we define the aggregated estimated channel matrices as Ĥm

l = [ ĥm l,1 , . . . , ĥm l,K l ] ∈ C Nm×K l .

B. Precoding and Power Constraints

Non-linear precoding schemes are infeasible in large scale systems due to their high complexity [START_REF] Hoydis | Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?[END_REF]. It is therefore important to find linear precoding schemes that perform well under imperfect CSI and heterogeneous mobility. We consider precoding matrices G m m , m = 1, . . . , L that have a generalized regularized zero forcing (RZF) structure. The generalizations are selected to closely resemble an optimal linear precoding structure derived in [START_REF] Björnson | Pareto characterization of the multicell MIMO performance region with simple receivers[END_REF]

]: 2 G m m = L l=1 α m l Ĥm l ( Ĥm l ) H + γ m I Nm -1 Ĥm m β 1 2 m (3) ∆ = α m m Ĥm m ( Ĥm m ) H + Z m + γ m I Nm -1 Ĥm m β 1 2 m where Z m = l =m α m l Ĥm l ( Ĥm l ) H .
The α m l are optimization variables and the classic regularization parameter γ m is chosen to be the inverse of the transmit signal to noise ratio (SNR). In [START_REF] Hoydis | Making smart use of excess antennas: Massive mimo, small cells, and tdd[END_REF] it is argued that such additional weights can be used to balance the precoding orthogonality to different (interference-) subspaces. This offers an intuitive explanation for the effectiveness of the generalized precoding observed in Subsection IV-C.

The factors β m are used to fulfill the average (w.r.t. the UTs) transmit power constraint P m of each BS m:

1 K m tr G m m (G m m ) H = P m . (4)
2 The truly optimal precoding matrices under perfect CSI would be parallel

to G m m = L l=1 K l k=1 α m l,k h m l,k (h m l,k ) H + γmI Nm -1 H m m β 1 2
m [START_REF] Björnson | Optimal resource allocation in coordinated multi-cell systems[END_REF].

C. Performance Measure

Many performance measures in cellular systems are functions of the signal to interference and noise ratios (SINRs) at each UT; e.g., (weighted) sum rate, outage probability and user geometry. Under the system model, the received signal power at the kth UT of cell l (UT l,k ) is

Sig (l) l,k = χ l l,k (h l l,k ) H g l l,k (g l l,k ) H h l l,k (5)
Similarly, the interference power is

Int (l) l,k = m =l χ m l,k (h m l,k ) H G m m (G m m ) H h m l,k + χ l l,k (h l l,k ) H G l l[k] (G l l[k] ) H h l l,k . (6)
where

G l l[k] = α l l Ĥl l ( Ĥl l ) H + Z l + γ l I N l -1 Ĥl l[k] β 1 2 l ( 7
)
and Ĥl l[k] is Ĥl l with its kth column removed. Hence, the SINR at UT l,k can be expressed as

γ (l) l,k = Sig (l) l,k Int (l) l,k + 1 . ( 8
)
In the following we focus on the sum rate (SRs), which is a commonly used performance measure utilizing the SINR values and straightforward to interpret. Under the assumptions that interference is treated as noise, the SR is

R sum = l,k R l,k = l,k log(1 + γ (l) l,k ). (9)

III. DETERMINISTIC EQUIVALENT OF THE SINR

In order to obtain tractable and insightful expressions of the system performance, we propose a large scale approximation. This allows us to state the sum rate expression in a deterministic and compact form that can easily be interpreted and optimized. Also, the large system approximations accurate even in massive MIMO systems of modest size, as will be evidenced later via simulations (see Section IV).

We will derive a deterministic equivalent (DE) of the SINR values that allow for a large scale approximation to the sum rate expression. The DE is based on the following technical assumption. Introducing the ratio c i = N i /K i , we assume

A 1. N i , K i → ∞, such that for all i we have 0 < lim inf c i ≤ lim sup c i < ∞.
This asymptotic regime is denoted N → ∞ for brevity.

By adapting the analytical approach in [START_REF] Wagner | Large system analysis of linear precoding in MISO broadcast channels with limited feedback[END_REF] to the SINR expression in [START_REF] Wagner | Large system analysis of linear precoding in MISO broadcast channels with limited feedback[END_REF], we obtain a DE of the SINR, which is denoted γ (l) l,k in the following. Proposition 1 (Deterministic Equivalent of the SINR). Under A 1, we have

γ (l) l,k -γ (l) l,k a.s. ----→ N →∞ 0 ( 10
)
where

γ (l) l,k = Sig (l) l,k Int (l) l,k + 1 (11)
with

Sig (l) l,k = β l (χ l l,k ) 2 e 2 (l) 1 -(τ l l ) 2 (f l l,k ) 2 Int (l) l,k = L m=1 β m 1 + 2c m l,k e (m) + α m l χ m l,k c m l,k e 2 (m) χ m l,k g (m) • (f m l,k ) 2
where c m l,k = α m l χ m l,k (τ m l ) 2 . The parameter β m , the abbreviations g (m) and f m l,k , as well as the corresponding fixed-point equation e (m) and e (m) are given in the following.

First, we define e (m) to be the unique positive solution of the fixed-point equation

e (m) =   γ m + 1 N m Km j=1 α m m χ m m,j f m m,j + 1 N m l =m K l k=1 α m l χ m l,k f m l,k   -1 (12)
where

f m l,k = 1 + α m l χ m l,k e (m) -1
. We have

β m = P m Nm Km g (m) (13)
where

g (m) = - 1 N m Km j=1 χ m m,j e (m) (f m m,k ) 2 (14)
and e (m) can be found directly, once e (m) is known:

e (m) = 1 N m Km j=1 (α m m ) 2 (χ m m,j ) 2 f m m,j + 1 N m l =m K l k=1 (α m l ) 2 (χ m l,k ) 2 (f m l,k ) 2 -e -2 (m) -1 . ( 15
)
Sketch of proof: The start of the proof conditions that Z m is fixed and follows the steps given in [START_REF] Wagner | Large system analysis of linear precoding in MISO broadcast channels with limited feedback[END_REF]Appendix II] for each term from Subsection II-C and the power normalization β m . After applying [8, Theorem 1] we obtain the fundamental equations e (m) . We now allow Z m to be random and apply [START_REF] Couillet | Random Matrix Methods for Wireless Communications[END_REF]Theorem 3.13] to obtain [START_REF] Couillet | Random Matrix Methods for Wireless Communications[END_REF]. It is admissible to apply the two theorems one after the other, as Z m is a bounded sequence with probability one. The DEs of the treated terms are found by, again, following [8, Appendix II] and (for the interference term) [8, Lemma 7].

By employing dominated convergence arguments and the continuous mapping theorem (e.g., [START_REF] Couillet | Random Matrix Methods for Wireless Communications[END_REF]), we see that Proposition 1 implies

R l,k -log 2 (1 + γ (l) l,k ) a.s. ----→ N →∞ 0.
(16) a 1 Fig. 2. Simplified simulation scenario.

IV. SIMPLIFIED MODEL

A. Simplified Simulation Scenario

In the following we introduce a (geometrically symmetric) simplified simulation scenario (see Fig. 2). It enables us to demonstrate the possibilities of interference analysis and mitigation offered by the general DEs introduced in Section III. Simultaneously, the simple structure allows us to draw intuitively evident conclusions that are not obscured by complications from more realistic scenarios.

Fig. 2 illustrates the simplified simulation scenario. It is composed of one macro cell (MC) BS, BS1, in the center, that generally occupies itself serving high mobility (HM) UTs. Furthermore, we have four small cell (SC) BSs arranged equiangular on a circle of radius a. These SC BSs are serving low mobility (LM) UTs. BS1 is equipped with N antennas and all others utilize N/2 antennas. All BSs have the same number of associated UTs (K). To simplify analysis further, the UTs are equally distributed and co-located according to the markings in Fig. 2, i.e., at the edge of the cells defined by the SC BSs. The radius of the SC cells is assumed to be 1. SCs are often placed within buildings, whose walls shield UTs that are not in the immediate vicinity of the building from interference. Thus, we assume that each UT only receives signals from the MC BS and the closest SC BS. The SCs employ traditional zero-forcing precoding, which is achieved by choosing ∀l, m = 2 . . . , 5: α m l,k = 1, if m = l and α m l,k = 0 otherwise, in our precoder (3). Furthermore, the value of γ m , m = 2 . . . , 5 is set to be very small. The distancedepended pathloss factors χ m l,k are calculated via 1/d n , where d is the distance between BS m and UT l,k and n is the pathlosscoefficient. For intuitively clear reasons, we introduce the definitions ∀m : τ m

B. Numerical Analysis of Mobility

The first step of the mobility analysis is to verify the usability of the system model and the accuracy of the derived DE. For this reason, we compare in Fig. 3 the per UT rate, obtained with τ LM = 0.1 and τ HM ∈ {0.1, 0.4} (solid lines). Furthermore, α LM and α HM are set to 1. As expected, we see that the rates reduce with higher mobility. In order to judge the applicability of the DE for systems of the given size, the results of Monte-Carlo (MC) simulations with 1000 channel realizations are given as marker symbols. We remind that the approximation can be improved by increasing the numbers of UTs and BS antennas.

Fig. 4 introduces two additional performance indicators, namely the average rate of only the LM UTs and the average rate of only the HM UTs. This allows to asses the impact of mobility of each user group in a distinctive manner. For this figure we choose α HM = 3.09 and α LM = 1.43 for reasons that will become evident in Subsection IV-C. The dashed lines in Fig. 4 show the optimum sum rate performance given perfect CSIT, as a reference. One can observe that performance of the LM UTs is unaffected by the CSI quality of the HM UTs. In other words, precoding with the bad channel estimates for the HM UTs of one cell does not cause additional interference to the LM UTs in the other cells. This is a rather surprising result, as others have previously found the multiplexing gain of the whole system to be limited by the high mobility UTs [START_REF] De Kerret | Degrees of freedom of the network MIMO channel with distributed CSI[END_REF]. Yet the result makes sense if one remembers that random beamforming is optimal in large systems and also that our model of small scale fading causes all channels to be independent. The mathematical explanation of this effect is that the statistics of our estimated channels are equivalent to those of the true channels (see (2)).

C. Numerical Mobility Management

In this subsection we take full advantage of the closed form expression of the sum rate derived in Proposition 1. In order to investigate how the performance degradation caused by mobility can be managed, we employ numerical optimization of the sum rate. More specifically, we note that one can use a single α LM instead of four different weights for the SC UTs, since the symmetric nature of the simplified scenario makes them all equal. The numerical optimization thus corresponds to a grid search over the α HM and α LM weights included in the generalized RZF precoding of BS1.

Fig. 5 shows how the weights α HM and α LM can be used to obtain the optimal average UT rate (HM+LM), under degrading CSIT for the HM UTs. The optimization is carried out by a 2-dimensional grid search for α HM and α LM . Recall that only the precoding at BS1 is adjusted, the other precoding matrices stay as originally defined (i.e., classical zero-forcing). We see that the optimal weights for the HM UTs (α HM ) are decreasing, when its CSIT is getting worse. This indicates that less importance is placed on the corresponding estimated channels. This is reversed for the channel weights of the LM UTs (α LM ). Overall, the increasing importance placed on the estimated channels (i.e., shifting from maximum ratio transmission to zero-forcing/interference cancellation) indicates that interference is becoming the dominant issue when CSIT becomes increasingly unreliable. It is especially informative to compare Fig. 5 with Fig. 4, as the latter uses the same starting values for the weights as the first, but does not adapt them to the mobility information. We observe that taking into account this information for precoding significantly improves overall system performance.

V. CONCLUSION

In this paper we analyzed the impact of heterogeneous UT mobility on the sum rate of a downlink multi-cell system employing generalized RZF precoding. We modeled the UT mobility in terms of a varying quality of the available CSIT and we derived a corresponding large scale approximation of the SINR values at each UT. Increasing UT mobility was found to reduce the sum rate of the system. In a more detailed analysis, we found that when the mobility of UTs in one cell increases, then the performance of these UTs decreases. However, the performance of UTs in other cells stays unaffected. The deterministic expression of the sum rate was successfully used in a practical mobility management scheme. Significant sum rate improvements were achieved by optimizing the weights of the generalized RZF precoding to balance between privileging served UTs and suppressing interference to other UTs. In general, the optimization prioritizes the UTs with the best CSIT. For future work, we expect large gains to be found by combining the generalized RZF precoding and proper scheduling.

E.

  Björnson is funded by the International Postdoc Grant 2012-228 from The Swedish Research Council. This research has been supported by the ERC Starting Grant 305123 MORE (Advanced Mathematical Tools for Complex Network Engineering).

Fig. 1 .

 1 Fig. 1. Illustration of a general heterogeneous downlink system.

1 ∆∆=α LM , α 1 1 ∆

 11 =τ HM , ∀m, l = 1 : τ m l = τ LM for the mobility parameters and l = 1: α 1 l =α HM for the precoding weights.Unless otherwise noted, the following system parameters values are used in this section: L = 5, N = 64, K = 16, n= 2, a = 1.5, P 1 = 20 dB, l = 1 : P l = 3 dB, τ LM = 0.1.

1Fig. 3 .

 3 Fig. 3. Average UT rate w.r.t. average transmit power under different CSIT quality τ HM ∈ {0.1, 0.4}.

Fig. 4 .

 4 Fig. 4. Average UT rate w.r.t. changing mobility of the HM UTs. Optimization variables α HM/LM constant and taken from Subsection IV-C.

Fig. 5 .

 5 Fig. 5. Average UT rate w.r.t. changing mobility of the HM UTs. Optimization variables α HM/LM chosen to maximize sum-rate.

Throughout this paper, superscripts denote the origin cell (e.g., m) and subscripts denote the receiving cell (e.g., l) or UT (e.g., l, k).