Generation of Spatially Heterogeneous Flood Events in an Alpine Region — Adaptation and Application of a Multivariate Modelling Procedure

Flooding often has a negative impact on society. In particular, widespread flood events can cause a lot of damage. These events are often spatially and temporally heterogeneous and should be duly considered for an appropriate analysis of flooding. Therefore, a conditional multivariate approach is adapted and applied in order to (i) contribute to a better understanding of the spatial characteristics of fluvial floods and (ii) to deliver sets of synthetically generated flood events. The present paper focuses on a simulation procedure consisting of careful data preparation and selection and the application of a conditional multivariate approach. The conditional approach is adapted to account for the seasonality of runoff data. Model checks attuned to the model are presented to ensure the consistence of simulated and observed data. The Austrian Province Vorarlberg was chosen as the study area. A thorough data analysis of runoff time series showed that the hydrological behaviour is characterized by a strong seasonality that was considered within the applied modelling procedure. The analysis of the spatial dependence of high river flows identified regions where floods likely occur simultaneously and regions with low spatial dependence. The main result of the modelling procedure, a large set of widespread flood events, was successfully generated.


Introduction
Fluvial floods are a natural phenomenon that regularly cause significant losses to property and human life.Through efficient flood risk management, the adverse consequences of flooding can be mitigated or even avoided.Therefore, a thorough flood risk analysis provides the basis for decision making of involved parties, such as flood management agencies, public authorities, and insurance industries.
In this context, flood risk is usually understood as a combination of the probability of events and the potential adverse consequences [1][2][3].A traditional procedure to derive the probability of flooding is to apply extreme value statistics to observed runoff series at river gauging stations.By means of this flood frequency analysis, the return period (T) of flood events in a year is determined for a single point i.e., gauging station) and being assigned to an entire catchment.In the course of flood risk analyses, these homogenous flood scenarios, assuming constant return periods of discharges within a river reach or an entire region, are frequently used.
As flood events in most cases are heterogeneous in space and time [4,5], the use of homogenous flood scenarios may lead to inaccurate results [6,7].Therefore, heterogeneous flood scenarios, which take into account the spatial variability of flood events, should be considered in flood risk analysis [6,8,9].As the number of past observed flood events leading to flood damages is usually low, the corresponding number of heterogeneous flood scenarios is limited and a statistical evaluation of the flood risk is difficult to achieve.
A promising approach to overcome this shortcoming is the simulation of synthetic flood events considering a vast range of possible flooding situations.The flood characteristics within a region of interest is typically represented by several gauging stations, and therefore a multivariate method, which is also capable of reproducing the spatial dependence structure in the investigated region, has to be applied.
In the field of hydrology, a popular approach to combine different marginal distributions and dependence structures is the copula approach.Although, in recent years, copula-based models have been applied in hydrology and related fields, such as multivariate frequency analysis, risk assessment, and multivariate extreme value analysis [8,[10][11][12][13][14][15][16][17][18][19], most applications are limited to bivariate and trivariate cases.Jongman et al. [8] used a stepwise conditional copula approach to account for a higher dimensional case in flood risk analysis.Apart from this study, the semi-parametric conditional exceedance model introduced by Heffernan and Tawn [20] (henceforth referred to as the HT model) is able to capture the multivariate nature of widespread flooding and can be used to simulate flood events.
The HT model can be interpreted as a multisite peak-over-threshold approach [21] and is able to model the joint probability of large sets of variables (e.g., river flow or sea-level data).It is an appropriate model to describe the probability that one or multiple variables are extreme.Because of these properties, the HT model is well-suited (i) to show the dependence structure of runoff in a region and (ii) to describe and reproduce the joint occurrence of high river flows at different sites.The HT model, thus, delivers high and extreme runoff peaks at one or multiple sites representing heterogeneous flood situations.A key feature of the HT model is its flexibility in modeling different dependence structures of extreme events.It can be used for both asymptotically dependent and asymptotically independent data.None of these two dependence structures should be excluded by the applied model (as it is done e.g., by a copula model) because it is likely that, for example, two gauges along the same river are asymptotically dependent, whereas it is possible that not-nested sites are asymptotically independent.Furthermore, the HT model is able to describe the changing dependence structure when runoffs become extreme.
The general idea and methodology of this conditional model was introduced by Heffernan and Tawn in 2004.Since that time, the HT model has been applied successfully particularly in hydrology and related fields; for example, in order to analyse the spatial dependence of runoff and precipitation in Great Britain (GB) [22], to simulate various ocean variables for extreme conditions [23][24][25], and in the context of flood risk analysis [21,[26][27][28][29][30][31].These flood risk studies show that the HT model is an appropriate approach to generate widespread flood events and that it is well-suited to consider different physical sources of flooding such as high river flows and sea levels.To the authors' knowledge, except for one study [27], all previous flood risk studies in which the HT model was applied focused on flooding in GB.The HT model makes the assumption that flooding in different seasons of the year have the same spatial dependence [22].The seasonal features of the processes were not explicitly considered in the above-mentioned flood related studies, but are taken into account in the present work.
This study introduces a modelling procedure that comprises suitable methods for data review and preprocessing, the application of the HT model with properly preprocessed and deseasonalized input data, and model checks attuned to the conditional model.This approach is therefore suitable to be applied in any setting in which the spatial dependence between seasons differs.The modelling procedure was applied and tested in a mountainous study area in which the seasonal characteristics of runoff in summer are different from the characteristics in winter.The main results were synthetically generated flood events on point scale, or, in other words, a data set of simulated runoff at the spatial extent of river gauges.A comparison of observed and simulated stream flow data shows the plausibility of the simulation.
The Austrian Federal Province of Vorarlberg has been chosen as the study area.In the past decades, Vorarlberg has been affected by serious flood events (e.g., May 1999, August 2002, August 2005 and June 2013) where each event led to severe economic losses [32].Thus, a probabilistic flood risk analysis model was developed to analyse the risk of flooding.The model consists of three major modules: (i) a hazard module where a data set of synthetic flood scenarios is generated which serves as surrogate to determine the probability of widespread flooding; (ii) an impact module with results being loss-probability-relations on a community level representing the potential flood damages; and (iii) a risk assessment module in which the probability of flood risk in the study area is evaluated statistically [33].The present paper focuses on the generation of spatially heterogeneous flood events.
The present study is organized as follows.Section 2 consists of a short description of the study area and the data used in this study.The modelling procedure begins with Section 3, where the focus lies on the data selection and preprocessing, which is inevitable to account for the seasonality of runoff.As the second part of the modelling procedure, Section 4 provides a review of the HT model.In Section 5, results about the seasonal characteristics of runoff, the spatial dependence measure, the application of the HT model including model checks, and the simulation results are presented.Finally, in Section 6, conclusions are drawn and future steps are discussed.

Study Area and Data Used
The Austrian Federal Province of Vorarlberg has been chosen as the study area (Figure 1).Vorarlberg is located in the eastern Alps and is characterized by its complex mountainous topography where elevation ranges from 395 m (Lake Constance) to 3312 m above sea level (Piz Buin).The region under investigation covers 2601 km 2 , with 91% of the area belonging to the river Rhine basin and 9% draining into the river Danube.
Vorarlberg is located north of the Alpine divide and experiences a relatively high amount of mean annual precipitation; for example, the region Bregenzerwald (mostly drained by the river Bregenzerach) receives approximately 2000 mm year −1 [34].A high percentage of the precipitation in winter falls as snow.Therefore, the seasonality of river flow is characterized by low runoff in winter and high runoff in spring and summer, especially in the high mountain regions in the southern part of the study area [35].
Runoff time series from river gauging stations, which are a direct measure for river flooding, were used as the primary data source for driving the HT model.For the present study, 17 representative gauges were chosen where time series of daily maximum runoff from 1 January 1976 to 31 December 2013 with no missing data were available (Table 1).A preliminary data analysis helps to understand the hydrological characteristics of the investigated region.Firstly, daily maximum and mean runoff data were compared.The advantage of using daily maximum runoff data is that information about peak discharge induced by short rainfall events (e.g., synoptic weather patterns) is included, which can be relevant in small catchments.A comparison of annual maximum series (AMS) derived from daily mean data (AMS dm ) and from daily maximum data (AMS dx ) for coinciding time period shows that the use of daily mean data would lead to an underestimation of the flood situation in Vorarlberg.On average, the medians of AMS dx are 87% higher than the medians of AMS dm .The use of daily maximum data is inevitable in the presented study, although, in most previous flooding related studies [21,26,28], daily mean data were used to run the HT model.
Secondly, the flood process typology was reviewed that also gives some indication of the length of flood events.Merz and Blöschl [36] analysed runoff conditions in Austria and classified annual maximum peaks (time series 1971-1997) in categories such as long-rain floods, short-rain floods, flash floods, rain-on-snow floods, and snowmelt floods.According to this classification, the flood events in the 17 watersheds considered in this study were triggered mostly by long-rainfall events (45%), by short-rainfall events (32%), rain-on-snow events (14%), and others (9%) [36].

Data Review and Preparation
Severe flood events usually have a low probability of occurrence and thus their investigation is based on a small amount of data.Hawkes et al. [37] emphasized that the data selection and preparation is probably the most important component when analysing extremes.This section introduces and reviews some appropriate approaches and aspects of data preparation and selection attuned to the HT model.

Event Definition
A key aspect of this study is to analyse widespread flooding i.e., flood events where multiple stations are affected).A flood event is defined as the maximum runoff of each site that occurs in a certain time interval of length τ simultaneously-in other words, the maximum value within regularly spaced blocks.This block maxima are used for further analyses.The choice of length of the time interval τ rests on the following statistical and hydrological considerations and finally relies on expert judgment.
1.An analysis of runoff time series, where the consecutive days exceeding a specific threshold (defined by the pth quantile q p ) are counted, provides a first estimation of τ.Denote the daily maximum runoff at gauge i ∈ [1, . . ., n] and time t by D i,t , and then the probability of peaks of length L is defined as 2. The travel-or concentration time provides additional information about the hydrological response of a watershed and therefore contribute to a thorough event definition.In this context, the concentration time T c (in hours) is a possible measure to indicate the time interval in which a receptor (e.g., flood-prone area) can be affected by simultaneous flooding from two or multiple tributaries.Grimaldi et al. [38] summarized and presented selected formulas for estimating the concentration time, such as, the formulas of (a) 'Department of Public Works', (b) 'Giandotti', (c) 'Kirpich', and (d) 'Viparelli', where only few input parameters are necessary.Full details concerning these formulas of the concentration time and their specified restrictions (e.g., regarding the catchment size) can be found in Grimaldi et al. [38].Although the empirical formulas to estimate the concentration time are associated with large uncertainties [38], they still provide a valuable instrument to calculate the travel time in investigated watersheds.3. The event definition (i.e., selection of τ) may also depend on the intended application of the model, such as reconstruction of flood situations in the past lasting a certain number of days or application in a distinct (re)insurance context [21].

Event Categorization
Besides the definition of flood events, the categorization of their severity can be conducted as follows: a severe flood event (e.g., that causes damage) is defined as the simultaneous occurrence of high runoff at different sites where at least one site exceeds a certain threshold.The severity of events can be determined using the proxy called unit of flood hazard (UoFH) [33,39].The UoFH is defined as the number of sites that experience runoff that equals or exceeds a certain threshold i.e., a runoff that corresponds to aspecific T).In this study, a return period T = 30-year was selected as a threshold to define the UoFH.

Seasonality of Runoff
The HT model relies on the assumption that fluvial flooding has the same spatial characteristics in winter months and in summer months [22].Thus, the knowledge about the runoff seasonality may help to analyse if this assumption is fulfilled.Directional statistics [40] are a particularly useful approach to anaylse the seasonality of high flows.Therefore, the occurrence dates of flood peaks (exceeding a certain threshold) are translated into a location on a unit circle, where the start of the year is plotted on the easternmost point and further months are plotted in a counterclockwise sense.Generally, the starting point depends on the hydrological characteristics of the study area i.e., splitting of high flow month has to be avoided).The monthly frequency of flood peaks (above a threshold) can be represented by rose diagrams, which helps to identify high and low flow months [41].A second measure for analysing the seasonality is the 'Burn vector' [42,43], which consists of the pair θ and r. θi represents the mean occurrence date of flood peaks of a specific site i and ri indicates the variability of the occurrence date (details can be found in [42]).

Interpretation of Extremes
When analysing several sites with different properties (e.g., catchment size, mean annual runoff), the use of the probability scale is convenient.In the field of hydrology, the interpretation of probabilities is usually done using return periods of events instead of quantiles of the variables.For a given quantile p, the return period (T) is defined as where n py is the average number of events per year.It is calculated as n py = 365 r s τ , where r s is the ratio that is required in order to take the type of the data set into account (for an all-year (AY) series r s = 1, a half-year (HY) series r s = 0.5, and a seasonal series r s = 0.25) and τ is the length (in days) of the time interval, which is used for event definition (cf.Section 3.1).

Heffernan and Tawn Model
This section reviews the measures to characterize the spatial dependence of river runoff, the statistical model (HT model) used to generate simulated events, and the estimation and simulation strategy including some details about the simulation procedure for the nonparametric part of the HT model.
For a set of gauges I, let Q i be the variable representing maximum flow within a time interval τ at gauge i ∈ I (henceforth real scale).Denote the number of gauges by n and an observation at time t by Q i,t .For either a set of observations or a random variable X, denote its pth quantile by q p (X), meaning that q p (X) has a probability p of not being exceeded.

Spatial Dependence Measures
The first measure of spatial dependence can be interpreted as an exploratory summary measure of bivariate dependence and is defined as [22] so P i,j (p) is the probability that gauge i exceeds q p (Q i ) given gauge j is extreme i.e., exceeds q p (Q j )).
The second dependence measure describes the proportion of gauges that are extreme given that gauge j is extreme and is defined as An efficient method for analysing the spatial dependence has to be able to qualitatively describe the behaviour of measured data for moderate and high p-values.The calculation of spatial dependence measures P i,j (p) and N j (p) is based on approx.Forty-five and five values in the present study for moderate i.e., p = 0.99) and high i.e., p = 0.999) p-values, respectively.The information about P i,j (p) and N j (p) over a range of p helps to characterize the dependence structure of different sites for various levels of extremeness.

Statistical Model
The statistical model is built up of a model for the marginal distribution for the runoff series Q i and a conditional dependence model for the joint distribution of these transformed variables given one of them is extreme.The model for the marginal distribution is used to transform all data to the same scale and ensures a relatively simple semiparametric form of the conditional dependence model.

Marginal Model
For the standardisation of the original data, we used a transfomation to Laplace marginal distribution [44] instead of a Gumbel marginal distribution as it is recommended in the primary HT model publication [20].
We use F i to denote the marginal cumulative distribution function (CDF) of the runoff series Q i of gauge i.Then, the transformation is defined by The transformed random variable X i (henceforth transformed scale) has a Laplace distribution with both the parameters of location and scale equal to 1. Before applying this marginal model, an inference on the distribution of Q i has to be done (see Section 4.3.1).

Dependence Model
In order to present the key idea of the conditional dependence model [20], we split the vector of transformed river runoff series (X i ) i∈I to a one-dimensional conditioning component X = X j and the vector of dependent components Y = (X i ) i∈I\{j} .
In the following, vector algebra is to be interpreted as componentwise and (in)equalities involving vectors should hold for all components.The HT model relies on the assumption that there are normalising functions a : R + → R n−1 and b : R + → R + n−1 such that for some non-degenerate distribution function G (with no mass at +∞).
On the one hand, the assumption of the existence of normalising functions is quite common in extreme value modelling, for example when fitting a Generalized Pareto distribution to onedimensional data.On the other hand, this assumption means that, conditional on X > u, the random variables X − u and Z = Y−a(X) b(X) converge to independent random variables for u → ∞.For a wide range of copula dependence models, the normalising functions a and b belong to a simple class of functions.When using the transformation to Laplace marginal distribution, this class can be represented as [44] a(x) = αx and b(x) = x β (7) When applying the HT model, the asymptotic assumptions are supposed to hold exactly for X > q p sim , for a high threshold q p sim (which is defined as the empirical quantile of the data for some high non-exceedance probability p sim ).Then, the dependence of the components in Y on the conditioning variable X is given by the parametric model where X and Z are independent.This represents n − 1 univariate regression type models , where the dependence of Y i on X is modelled parametrically, but the dependence structure of the random vector Y is given non-parametrically by the distribution G of Z.

Estimation and Simulation
The aim of the simulation is to generate a set of synthetic events S such that at least one of the river flows is extreme-in other words, has an exceedance probability higher than a given p sim .Several steps are involved in fitting the HT model to the data and subsequently to use the fitted model to generate simulated extreme events.These steps are described in some detail below and examples i.e., of sites Kennelbach (j = 2) and Thal (i = 5)) are outlined in Figure 2a-d.) observed data (threshold for marginal distribution q p GP =0.9 ) in transformed scale, incl.q p sim =0.98 ; (c) observed data in transformed scale (values above q p sim ), incl.fitted model (red color indicates that the GP distribution was applied above the threshold q p GP =0.9 ; grey color indicates q p GP =0.98 ); and (d) observed and simulated data (n y = 38 year) on a real scale.

Estimation of Marginal Distribution
As suggested in [20], estimation of the marginal distribution Fi of the runoff variable Q i for 1 ≤ i ≤ n is done by using the empirical distribution ECDF i of the data below some threshold q p GP , where q p GP is the empirical quantile of the data for some high non-exceedance probability p GP .The Generalized Pareto distribution is fitted to the data above the threshold q p GP using a maximum likelihood method.The CDF of this distribution is denoted by GP σi , ξi .Here, σ i and ξ i are the scale and shape parameter of the Generalized Pareto distribution [45].Hence, Fi (q) = 1 − (1 − p GP )(1 − GP σi , ξi (q)) for q > q p GP ECDF i (q) for q ≤ q p GP (9)

Marginal Transformation
Transformation to Laplace margins is done using the estimated CDF Fi instead of F i for the marginal transformation using Equation (5).When simulation results X i are available on a transformed scale, and backwards transformation to real scale is done by where F−1 i (p) is the inverse of the empirical distribution for p ≤ p GP and the inverse of the Generalized Pareto distribution GP σi , ξi evaluated at (p − p GP )/(1 − p GP ) for p > p GP .

Parameter Estimation
To estimate the model parameters α i and β i in Y i = α i X + X β i Z i , we use a quasi maximum likelihood approach [20].Let µ i be the mean of Z i and σ i its standard deviation, and then For a pair of transformed observations (x t , y i,t ) t=1,...,N t , where x t is extreme (x t > q p sim for all t), we set and estimates of the parameters are given by (α i , βi , μi , σi ) = arg min The distribution of Z i is not estimated parametrically.Instead, simulation will be based directly on the estimates Ẑi,t = y i,t − αi x t This estimation procedure is performed for all possible pairs of data (x t , y i,t ), where x t is extreme.

Simulation of the Conditioning River Flow X
The HT model provides a method to describe the joint conditional distribution of all runoff variables conditional on one being extreme.Before the simulation of the conditional is done, the data to condition on has to be provided, as follows (similar to [28]): 1. Number of simulated extreme events per simulation period: The number of extreme events i.e., where at least one site exceeds q p sim ) per simulation period (e.g., all-year-, half-year series) is approximated by a negative binominal distribution.2. Selection of conditioning gauge: The set of extreme events S is partitioned to S = n i=1 S i , where S i is the set of extreme events where gauge i has the highest non-exceedance probability i.e., is the most extreme one).Then, for an event E ∈ S, the probability Pr(E ∈ S i |E ∈ S) that gauge i is the most extreme if at least one gauge is extreme has to be estimated using given data.The gauge that is selected as conditioning variable X is drawn from a multinomial distribution with these probabilities.Denote the index of the selected gauge by i .3. Drawing the conditioning x i : According to the transformation to standard Laplace margins and since p sim is close to one and hence q p sim > 0 set x i = q p sim + e, where e is drawn from an exponential distribution with mean one.

Simulation of Z
Since no parametric model for the distribution of Z is known, we apply the data based random vector generation method of Taylor and Thompson [46] to the set of vector valued estimates { Ẑt } t=1,...,N t = {( Ẑi,t ) i∈I\{i } } t=1,...,N t .This method uses m nearest neighbours of a randomly chosen element in { Ẑt } t=1,...,N t , where m is a smoothing parameter.A simulated vector Z is generated by the following steps: The smoothing parameter m is chosen and the limits of the uniform distribution from which the weights w i are drawn from are calculated such that the variance and covariance structure of the data { Ẑt } t=1,...,N t is unchanged, at least approximately.

Simulation of the Dependent Event
When the conditioning value x i and a sample Z of { Ẑt } t=1,...,N t are drawn, it is straightforward to compute a sample of the dependent gauges Y by Y = αx + x βZ (14) If all components of Y are smaller than x , a valid synthetic event is generated.Otherwise, the simulation of Z and computing Y has to be repeated.

HT Model and Seasonal Characteristics of Runoff
The HT model relies on the assumption that river flooding in the winter half-year has the same spatial characteristics as flooding in the summer half-year (cf.Section 3.3 and [22]).The spatial dependence measure N j (p) is an appropriate tool to investigate whether the spatial characteristics of different series (all-year (AY), May to October (M2O), and November to April (N2A) series) differs significantly.If the assumption does not hold, then the data has to be split into periods that reveal the same spatial characteristics and the HT model has to be applied to both data sets independently.

Results and Discussion
Runoff data from selected river gauging stations of the study area Vorarlberg were used (i) for event definition and seasonality analysis, (ii) to conduct a spatial dependence analysis of river flows, (iii) to apply the present multivariate conditional model, and thereby (iv) to generate synthetic (extreme) flood events.The results are based on daily maximum runoff data from 17 representative river gauges in the study area.
As a preliminary remark, note that the number of years n y simulated with the HT model can be arbitrarily high.This property is essential for the generation of synthetic flood events and for the subsequent evaluation of flood risk.Nevertheless, in this paper, a relatively small simulated data set i.e., n y = 38 year) was used to make the comparison with the observed data (38-year series) reasonable.

Event Definition and Seasonality of Runoff
Daily maximum runoff data (D i ) were analysed using the methods introduced in Section 3.1 in order to identify the most appropriate τ.Firstly, the results of the measure R(p, L) are summarized in Table 2.For all investigated quantile values p, the probability that an event lasts longer than three days is less than 5%.Secondly, the concentration time (cf.Section 3.1 (ii)) of the selected catchments is displayed in Table 3. Apart from the Rhine catchment at Lustenau (E) (ca.6000 km 2 ), all sites react relatively quickly (T c less than half a day).Thirdly, the flood process typology of annual peaks [36] (cf.Section 2) reveals different relevant flood controls, which hints that the time interval with length τ of three days is a good compromise.Fourthly, (cf.Section 3.1 (iii)), it should be noted that the presented work is embedded in a flood risk project with a (re)insurance background.In this industry, a time span of 72 h ('hours clause') is often used for event definition [47,48].Based on the three criteria, a time interval with length τ = 3 days was chosen.This time interval guarantees in most cases that consecutive independent flood peaks induced by short rainfall events were considered.On the other hand, τ = 3 days is long enough to ensure that flood peaks in neighboring catchments induced by a single weather pattern will be interpreted as one event.Besides the aspects mentioned above, the idea of the spatial dependence measure P i,j (p) can be applied to investigate the temporal dependence in individual series.Therefore, Equation (3) has to be adapted to account for daily maximum data D i,t .The dependence measure is then given by P (D,L) j,j (p) = Pr(D j,t+L > q p (D i )|D j,t > q p (D j )) [22].
For analysing the temporal dependence of all individual series, was calculated for different p-values.Thus, P (D,L) (p) shows averaged values of the dependence measure over all sites.For a lag of L = 2 days P (D,L) (0.99) = 0.11, P (D,L) (0.997) = 0.03, and P (D,L) (0.999) = 0.01.For high thresholds (corresponding to p = 0.997 or higher), the temporal dependence is less than 5%.Only for low thresholds (corresponding to p = 0.99) is the temporal independence is not fulfilled; however, this small level of extremeness is not relevant when applying the HT model.Generally, a time interval with a length τ = 3 days ensures temporal independence of events for high thresholds (corresponding to pth quantile).As regularly spaced blocks (i.e., block maxima) were used in this step, it might happen that a block limit falls within a runoff peak at an individual site that can cause the same event to be selected twice.However, as extreme peaks that are relevant in terms of model estimation (above u) are very short (cf.Table 2), a splitting of high peak is extremely rare and thus these cases are negligible.So far, daily maximum runoff data (D i ) were used for data analyses and event definition.The following results are based on maximum runoff in a time interval τ = 3 days i.e., Q i values).Note a quantile value p = 0.997 applied to daily data (D i ) equals T ≈ 0.82-year and therefore roughly corresponds to a quantile value p = 0.99 applied to Q i values (T ≈ 0.91-year).
The runoff in the investigated catchment is characterised by strong seasonality, meaning that high flow conditions mainly occur in summer.The mean occurrence date (represented by θi ) of AMS of all sites is from July to mid of August (Table 3).Following Merz and Blöschl [36], 94% of the catchments experience medium to strong and only 6% weak seasonality, which indicates that in most catchments annual peaks are concentrated around the mean occurrence date.Figure 3 shows monthly frequency of flood peaks above certain threshold values q p (Q i ).The results of all sites were combined and The pronounced seasonality suggests that there is a different spatial characteristic in M2Oand N2Aas well as in the all-year series.Figure 4 shows scatter plots of the spatial dependence measure N j (p) of 17 sites and for three different p-values for half-year and all-year series.The threshold q p used in Equation ( 4) to calculate N j (p) was derived from the all-year series.To make the results comparable i.e., data from half-year-and all-year series), the pth quantile of half-year series p HY was calculated using p HY = 1 − 1−p AY 0.5 , where p AY is the pth quantile of all-year series.In Figure 4a, the dependence measure of the all-year series N AY j (p AY ) versus the N M2O j (p HY ) (summer month) is plotted, whereby a strong correlation can be especially observed for high p-values.For p AY = 0.95, 0.99 and 0.997, the dependence measure for AY series is systematically lower than for M2O series.Figure 4b shows that the spatial dependence between AYand N2A series differs and that the dependence measure for all-year series is systematically higher than for the winter half-year.The distribution-free Wilcoxon test [49] was applied to assess whether N N2A j (p HY ) and N M2O j (p HY ) differ significantly.The null hypothesis that the spatial dependence of winter and summer half-year do not differ was rejected for a significance level of 1%.However, the HT model relies on the assumption that flooding in the winter half-year has the same spatial characteristics as flooding in summer half-year [22].Thus, half-year series were used for further analyses.

Spatial Dependence in River Flows
The analysis of the spatial dependence of river flows can help to understand the characteristics of widespread flooding.The spatial dependence measures P i,j (p) and N j (p) were calculated for M2O and N2A series, respectively.P i,j (p) explains the dependence of a pair of variables for a certain threshold q p , where P i,j (p) = 0.4 means that there is a 40% probability of Q i being large i.e., above q p ) if Q j is above the pth quantile.The spatial dependence measures P i,j (p) of four selected stations (Kennelbach, Thal, Gisingen, Schruns) over a range of p are depicted in Figure 5, where P M2O i,j (p HY ) (summer half-year) is shown in the upper triangle and P N2A i,j (p HY ) (winter half-year) is shown in the lower triangle.P i,j (p) of observed data (solid cyan lines) shows the range of dependence between the stations.There is a relatively strong spatial dependence for low p-values between nested catchments where the downward propagation of flood wave plays an important role (e.g., Kennelbach and Thal as well as Gisingen and Schruns), whereas P i,j (p) varies markedly for high p-values.The dependence of summer half-years between the stations Kennelbach and Gisingen is fairly equally spread over the entire range of p.These watersheds are relatively large and therefore more vulnerable to long-rainfall events triggered by advective storm systems that affect the entire study area.In Figure 5, the lowest dependence can be found between the stations Thal and Schruns.This is because the watersheds are relatively small and their hydrological characteristics differ significantly [35].Generally, the spatial dependence between stations in the summer half-year is stronger than in the winter half-year; in 65% of all possible combinations, i.e., 272 cases), P M2O  i,j (p HY ) is higher than P N2A i,j (p HY ) (averaged over p).The most obvious differences between P N2A i,j (p HY ) and P M2O i,j (p HY ) appear between station Kennelbach and Gisingen as well as Thal and Gisingen, where P N2A i,j (p HY ) is high for large p.The second spatial dependence measure N j (p) can be interpreted, as synopsis of P i,j (p) over a range of i.In contrast to Keef et al. [22], who investigated the spatial dependence of river flow and precipitation across GB by analysing N j (p) of stations within a defined radius (30 and 60 km), we analysed the spatial dependence over the entire study area i.e., consideration of all stations).This is reasonable since the entire study area is relatively small and widespread flooding across the study area is relevant for the Province of Vorarlberg.
Figure 6 shows the spatial dependence measure N j (p) for the summer half-year.Similar to the first dependence measure, a comparison of N M2O j (p HY ) and N N2A j (p HY ) reveals that the dependence in summer is higher than in winter i.e., in 77% of possible combinations, N M2O j (p HY ) is higher than N N2A j (p HY )). Figure 6a illustrates N M2O j (p HY ) over a range of p HY .Overall, the plot indicates that the estimated dependence between the sites decreases with increasing p HY .In general, the stations located at the same river (nested catchments) show high dependence (e.g., Kennelbach, Mellau, Hopfreben).Data from small headwater catchments and catchments that are situated on the edge of the study area often show low dependence (e.g., Thal, Unterhochsteg, Lochau).In particular, for these stations, N M2O j (p HY ) increases with higher p HY i.e., fewer data involved in the analysis), because these small catchments seem to be particularly vulnerable to synoptic weather patterns, which in turn are often responsible for the most extreme events in small catchments.Furthermore, the above-mentioned stations with low dependence lie within a radius of 5 km, which indicates that spatial proximity does not guarantee similar behaviour.As the entire study has ,a distinct flood risk background, spatial dependence of rare events was studied in detail.Therefore, we illustrate N M2O j (p HY ) for p HY equal to 0.99, 0.998, and 0.999, which correspond to T ≈ 1.6, 8.2, and 16.4-year in maps of Vorarlberg (Figure 6b-d).The gauges are coloured according to N M2O j (p HY ). Figure 6c is particularly suitable for detailed discussion as N M2O j (p HY ) is based on a reasonable number of events and p HY = 0.998 is high enough to analyse the dependence of severe flood events.Overall, the stations can be divided into three groups according to their spatial dependence: (i) low dependence N M2O j (0.998) < 0.3 (stations: Thal, Hoher Steg, Unterhochsteg, Schruns, Lustenau (E), Lochau); (ii) moderate dependence 0.3 ≥ N M2O j (0.998) < 0.4 (stations: Lingenau, Enz, Lustenau (H), Garsella, Gisingen, Lech); and (iii) high dependence N M2O j (0.998) > 0.4 (stations: Mellau, Kennelbach, Hopfreben, Schönenbach, Laterns).The sites with high dependence are all located in the region Bregenzerwald, which experiences a high amount of annual precipitation.

Estimation Results and Model Checks
The estimation and simulation procedures outlined in Section 4.3, which were implemented in MATLAB (The Mathworks, Natick, MA, USA), were carried out using the data from summer and winter half-year.The interim results of the bivariate case are illustrated in Figure 2 for data of stations Kennelbach (i = 2) and Thal (i = 5) for data of M2O. Figure 2a shows the observed data (Q 2 , Q 5 ) on a real scale.The same data are depicted in Figure 2b in a transformed scale, whereas the GP distribution was applied for data above the threshold q p HY GP =0.9 .Figure 2c illustrates X 2 and Y 5|2 in a transformed scale above the threshold q p sim .This figure includes data where the marginal GP distribution was applied above the threshold of q p HY GP =0.9 (red) and above q p HY GP =0.98 (grey), respectively.Furthermore, the fitted model (Y 5|2 = α2 X 2 + X β2 μ2 ) was applied for values with an exceedance probability higher than p HY sim = 0.97, 0.98, and 0.99. Figure 2c indicates the influence of the choice of the model parameter on the regression model.The red bold line illustrates the fitted model for the selected model parameters p HY GP = 0.9 and p HY sim = 0.98.The corresponding quantile values of an all-year series would be p AY GP = 0.95 and p AY sim = 0.99. Figure 7 shows data and fitted regression models for periods M2O and N2A for exemplarily chosen stations.Thus, the difference of the models for M2O and N2A periods are apparent, which confirms the necessity of using half-year series.After parameter estimation was done, the simulation procedure starts with the selection of the conditioning river flow.Therefore, the number of extreme events per year was simulated with a negative binomial distribution.Figure 8a shows that the distribution of simulated data is compatible with the empirical distribution of the observed data.The next step was the selection of the conditioning station i i.e., the station with the highest non-exceedance probability.The distribution of the conditioning station is shown in Figure 8b.The cyan crosses represent the probability that one of the 17 gauges is the most extreme based on observed data, and the percentile ranges result from drawing from the multinomial distribution (100 replications).The distribution of the sites with the largest probability of being extreme (Figure 8b) may appear counter-intuitive, as one could expect a uniform distribution.However, the results are in accordance with the findings of Keef et al. [28], who find similar properties of investigated runoff data in GB.Sites that lie on the same river exhibiting high dependences have low probability of being used as conditioning gauge (e.g., sites i = 1 and 2).Sites with low dependence are expected to have higher probability of being used as conditioning gauges; the four sites with the highest probability (i = 11, 12, 15 and 17) are categorized as stations with low spatial dependence (see Section 5.2).After the conditioning site i was selected, its value x i was drawn.Before simulating the dependent component (i.e., Y), Z has to be drawn as outlined in Section 4.3.5.The smoothing parameter was chosen with m = 5.Finally, a set of spatial dependent extreme events S was simulated.Figure 9d shows a scatterplot (Kennelbach and Thal) in real scale, where green triangles represent simulated events representing n y = 38 years of the simulation period.Note that, although the simulated events can be assigned to a selected time stamp (i.e., year), no continuous time series of river flow was made available.As no general Goodness-of-Fit test of the HT model exists, several aspects were tested to ensure that the model assumptions were not violated.The first assumption that needs to be checked is that the tails of the marginals are Generalized Pareto distributed.Asymptotically, this assumption is supported by the Pickands-Balkema-de Haan theorem [45].To assess whether the threshold q p GP was high enough, we used the Goodness-of-Fit test for GP distribution [50].For thresholds corresponding to p GP ≥ 0.9, the Generalized Pareto distribution of the tails was never rejected for the 17 investigated river gauges on the 5% significance level; on the 10% significance level, very few rejections appeared.Furthermore, stability of the estimated parameters (see e.g., chap.4.3.4 in Coles [51]) was checked by visual inspection of parameter estimate plots.This shows that a value of p GP of at least 0.9 is reasonable for all gauges.
Secondly, the assumption of independence of the conditioning variable X and the vector of normalized dependent variables Z = Y−αX X β , given that X is above some high threshold, is checked

Table 4.
Comparison of observed and simulated data for selected criteria (P i,j (p) and N j (p), and medians of yearly maximum runoff for summer and winter half-year).The results exhibit the proportion of possible combinations (P i,j (p), N j (p)) and sites (median(MR)), respectively, where the observed data lies within the percentile range of simulated data.

Simulated Extreme Events
The main output of the HT model is a set of spatially distributed flood events S on defined locations i.e., river gauges).As the runoff intensity of the summer half-year (see also Figure 3) is considerably higher than the intensity of the winter half-year, only events from M2O were considered here.In a first step, the information about the intensity of these events was made available in real scale i.e., m 3 s −1 ) for each individual event and every station.To enable a further usage of the synthetic events, it was necessary to define the so-called level of extremeness [21] for all events and every single station on a probability scale.This probability scale corresponds to the concept of return period of hydraulic loads.
The link between physical scale (m 3 s −1 ) and 'return period' scale, which corresponds to 'probability' scale, was constructed by applying the Generalized Extreme Value (GEV) distribution.
The method of L-moments [53] was used to fit the GEV distribution to the AMS of observed data of all stations and a χ 2 Goodness-of-Fit test was carried out to check the suitability of this distribution.The GEV distribution was never rejected at a 5% significance level for the 17 investigated stations.Furthermore, previous studies pointed out that the GEV distribution is an appropriate distribution function within the region [35,54,55].
Figure 11 shows five observed and five simulated flood events where the color of the triangles represents the return period that occurs on each station.The synthetic events reveal a reasonable spatial dependence structure that can also be found in observed flood events.For further analyses, only severe flood events are of interest, where one or multiple stations exceed a T ≥ 30-year.The severity of the illustrated flood events was categorized by the measure unit of flood hazard (see Section 3.2), which is indicated on the bottom of each map.For comparison, following the categorization of the UoFH, nine events were identified in the period 1976-2013, and the most severe observed flood event (August 2005) was categorized as a 12 UoFH-event i.e., the runoff was equal to or exceeded a threshold that corresponded to a return period T = 30-year at 12 river gauges).

Conclusions
In this paper, we have demonstrated that the HT model [20] can be successfully applied in a complex mountain topography with high runoff seasonality.However, before applying the conditional dependence model, a priori data review and preparation is essential in order to avoid violation of the model assumptions.A thorough data analysis of daily maximum data of 17 river gauging stations in the study area Vorarlberg shows that a time interval of three days is most suitable for event definition.Because of the pronounced seasonality of runoff and the different spatial characteristics in the half-year series, a separation of the all-year series was necessary.Thus, the conditional model was driven with time series from M2O and N2A.
The spatial characteristics of widespread and often spatially heterogeneous flood events were analysed.Therefore, the spatial dependence measures P i,j (p) and N j (p) were calculated in order to improve the understanding of the spatial characteristics of widespread river flood events.River gauging stations were identified where high runoff likely occurs simultaneously with other stations and, on the other hand, stations were found that show hardly any spatial dependence with other stations.
By applying the HT model, a set of synthetic flood events were generated.These extreme events represent a large range of possible flooding situations as they could occur in the investigated region.The set of synthetic flood events serves as a surrogate for probabilistic flood risk analysis.The severity of individual flood events was categorized by a simplistic hydrological proxy (units of flood hazard).The present modeling procedure can be applied wherever a dense network of river gauging measurements exist that capture the spatial characteristics of flood events in a specific region.
The analysis that is shown here does not include considerations about the potential adverse consequences of flooding.However, for a comprehensive analysis of flood risk, the following aspects are indispensable.Besides the generation of synthetic flood events, the adverse consequences of flooding have to be quantified (e.g., the monetary damage of buildings) and the overall flood risk within a region has to be assessed.Schneeberger et al. [33] introduces a probabilistic framework for risk analysis, which takes into account all relevant aspects and requires a set of spatially distributed flood events.

Figure 1 .
Figure 1.Study area Vorarlberg and the location of river runoff gauges.

Figure 2 .
Figure 2.Scatterplot of observed and simulated data from half-year series (May to October) of gauges Kennelbach (i = 2) and Thal (i = 5), statistical model, and model thresholds: (a) observed data in real scale; (b) observed data (threshold for marginal distribution q p GP =0.9 ) in transformed scale, incl.q p sim =0.98 ; (c) observed data in transformed scale (values above q p sim ), incl.fitted model (red color indicates that the GP distribution was applied above the threshold q p GP =0.9 ; grey color indicates q p GP =0.98 ); and (d) observed and simulated data (n y = 38 year) on a real scale.

Figure 3 .
Figure 3. Rose diagram of runoff (for all Q i ) above threshold values (q p (Q i )).

Figure 4 .
Figure 4. Comparison of spatial dependence measure N j (p) derived from all-year and half-year series: (a) N AY j (p) versus N M2O j (p) and (b) N AY j (p) versus N N2A j (p).

Figure 8 .
Figure 8. Model application: (a) empirical CDF of number of extreme events i.e., at least one site exceeds q p sim ) per year and (b) probability that gauge i is most extreme.

Figure 10 .
Figure 10.Comparison of observed and simulated data of period M2O for selected 17 sites.Medians of observed time series are depicted by blue crosses and medians of simulated data (100 × 38 year) are shown by selected percentiles.

Figure 11 .
Figure 11.Observed and examples of simulated flood events (T in year) including a categorization of the events severity by means of unit of flood hazard (UoFH).

Table 1 .
Selected River gauging stations, including the ID (station ID) according to the Austrian Hydrographic Service, and the catchment area (km 2 ).
1. Randomly choose one Z t 0 from { Ẑt } t=1,...,N t .2. Find the m nearest neighbours of Z t 0 denoted by Z t 1 , . . ., Z t m and compute the vector of componentwise means Z of Z t 0 , . . ., Z t m .3. Draw a random sample w 0 , . . ., w m from a uniform distribution with lower bound 1 m+1 −

Table 2 .
Probability R(p, L) that peaks (exceeding a threshold q p ) last L days.