Spatial Dependence Modeling of Flood Risk Using Max-Stable Processes: The Example of Austria

We propose a new approach to model the dependence structure for aggregating the risk of flood damages from a local level to larger areas, which is based on the structure of the river network of a country and can be calibrated with publicly available data of river discharges. Building upon a suitable adaptation of max-stable processes for a flood-relevant geometry as recently introduced in the literature, this enables the assessment of flood risk without the need for a hydrological model, and can easily be adapted for different countries. We illustrate its use for the particular case of Austria. We first develop marginal flood models for individual municipalities by intertwining available HORA risk maps with the actual location of buildings. As a second alternative for the marginal modeling, we advocate an approach based on suitably normalized historical damage data of municipalities together with techniques from extreme value statistics. We implement and compare the two alternatives and apply the calibrated dependence structure to each of them, leading to estimates for average flood damage as well as its extreme quantiles on the municipality, state, and country level. This also allows us to quantify the diversification potential for flood risk on each of these levels, a topic of considerable importance in view of the natural and strong spatial dependence of this particular natural peril.


Introduction
Natural hazards are a challenging topic for the modeling and managing of the associated risks, mainly due to the low frequency and high impact of their occurrences. As a prime example of such high impact low frequency hazards, we consider (fluvial) flood risk in this paper. We will present two different modeling approaches to estimate marginal flood risk at the municipality level. The first one is based on risk maps to estimate the flood damage distribution in a given municipality, whereas the second approach uses real damage data for the municipalities and extreme value statistics techniques to fit distributions into ranges where previous observations are scarce or not present at all. In order to aggregate the flood damage distribution of individual municipalities to larger areas, we will use a recently proposed dependence structure based on a Brown-Resnick max-stable process for joint discharges at river gauges (see [1]), which we will adapt to the respective given river network. While the proposed approaches are generally applicable for arbitrary countries and hierarchy levels, throughout this paper we will work with the concrete case of Austria, and the levels of municipality, states ("Bundesländer" in German), and the entire country. The raw loss data on the municipality level will be taken from the Austrian disaster relief fund and are available for four of the nine Austrian states. We will suitably normalize these data to the respective amount of building value and the appropriate building cost index to make the data comparable and usable.
While the use of real damage data is quite scarce in the literature, the first approach for marginal modeling using flood hazard maps is somewhat similar to methods that have already been used in the flood risk literature. Concretely, in a first step, an approximation of the inundation area for a flood of a given return period is evaluated. Subsequently, impact functions are employed that estimate the damage associated with the inundated area. In the following, we briefly summarize some related previous literature following this approach. [2][3][4] estimate the inundation area for different return periods. In [2] the return period of the yearly maximal discharge (in terms of the present climate) is estimated for present and future climate data runs. The authors then use the interpolation of inundation areas of floods with given return periods (the inundation areas are taken from [5]) to get the inundation area for a return period of the considered year. In [4], first, the economic impact of floods with given return periods of 2-5-10-20-50-100-250-500 years are calculated for a set of considered areas. To that end, for a given area, the inundation area of a flood with a given return period is calculated, and then impact functions depending on different land use categories and the inundation are used to estimate the impact of the flood. To get an approximation of the distribution of the flood impact aggregated over the considered areas, a Monte Carlo method is used. That is, the authors generate realizations of the aggregated flood damage for a year. The impact for each area is calculated as the interpolation of the impact of the previously calculated impacts of floods with given return periods. The realization of the aggregated impact is then the sum of the realizations of the individual impacts. In order to combine the results over the considered areas, a copula approach is used for dependence modeling. Using discharges taken from the LISFLOOD model, the concrete choice of copula is estimated within a general flexible family of copulas. The consideration of the dependence structure allows to also estimate quantiles of the damage distribution. Similar to [4], also in [3] at first the hazard maps and corresponding impacts are calculated for a set of given return periods for present and future climate and the yearly average damage is calculated by integrating over the damages from the impact curves. Since flood protection measures are not incorporated into the inundation areas in the studies [2][3][4], only damages from floods that exceed given thresholds are considered. The studies [6,7] use similar methods as in [2] to get the inundation area. However, instead of using risk maps for given return periods, the inundation areas are calculated for each year individually from the estimated discharge values. This has the advantage that no interpolation method is needed. The aforementioned studies are all global or continental, meaning that they consider flood risk for larger areas. Consequently, the estimated inundation areas are quite rough, as detailed regional modeling is not feasible.
Flood impact studies that use more detailed models, are for example provided in [8,9]. In these studies, besides the inundation areas for a single discharge event with a given return period, a continuous simulation of several years is used to get the river discharge and flood damage from the precipitation. These models are restricted to individual flood river catchments in Germany, although [9] indicates that the method can in principle work for larger catchments. In fact, [8] stands out as there a weather generator is used to generate a long time series of flood damages which allows for the calculations of higher quantiles in the flood damage distribution.
In the literature, there are different choices for the units in which the impact of a flood event is expressed. The impact can be measured by the number of affected people (e.g., [10]); damages to the residential building stock (e.g., [8,9]) or damage functions depending on the land use class (e.g., [4]). In addition to the direct loss estimated by impact functions, some studies consider the indirect impact in terms of welfare losses (e.g., [6]), or the effects on the finances of countries and implications for risk management (e.g., [11]).
In this paper, we measure the impact as the damages to the residential building stock. Since the risk maps used in this paper only provide the flood extent and not further variables like inundation depth, we cannot use the advanced damage functions provided in [12][13][14][15]. Instead, we will use simpler damage functions that rely on two different distributions for the damage, depending on the return period of the flood and the location of the buildings. Further, we will use actual flood damage data on municipality level to calibrate the model. We would like to note that there are considerable uncertainties in the calculation of flood damages [16,17] and hence results should not be interpreted as exact values.
We use the HORA risk maps ('Hochwasserrisikozonierung Austria', see [18]) as the hazard map for the estimation of the inundation area. This has the advantage that it takes more details into consideration (e.g., existing concrete flood protection measures in certain areas) than models for larger areas could provide. The downside of this approach is that one can only apply this method for areas where flood risk maps already exist, and that one is restricted to using the available return periods in that map. This is a little critical, since, usually, the smallest considered return period is 30 years. Further, to aggregate the damages of the different inundation areas, a dependence structure has to be used for the calculation of quantiles of the damage distribution (see for instance [4,11]). Using such statistical methods for the aggregation of the damages in the individual areas has the disadvantage that physically implausible results may arise from the simulation, see [9,19] for a discussion of this problem. However, the principle of statistically combining marginal information with particular dependence shapes is a standard approach in many areas of risk management (see e.g., [20]), see also [21] and the two studies [22], [23] for the particular flood modeling context. In this paper, we improve on the existing approach from [21] in two aspects. First, as in [4] we use the dependence of discharges as a proxy for the dependence of flood losses, but instead of using the flexible yet somewhat arbitrary copula of [4], we adapt a model for dependence that was recently published in [1] and that accounts for the river network. The second improvement to the approach in [21] is the use of satellite data combined with building stock data to intersect building stock data and inundation areas.
As a second approach, we will also use actual loss damage data of municipalities of Austria. In particular, we will utilize Extreme Value Analysis (EVA) to generate a damage model for each municipality and then aggregate the data using the same dependence structure as described above. For the estimation, we will use damage data from the payments of the catastrophe fund for four Austrian federal states, where the timespan of the available data ranges from 26 years for Styria to only 7 years for Upper Austria. The latter is obviously short for statistical purposes. In addition, one may argue that flood risk changes over time. For example, [14] showed that the main driver of changes in flood risk is urban sprawl, as with increased population cities develop into more flood-prone areas. At the same time, flood risk can also change due to other changes in human behavior. [24][25][26] show how human behavior can change the effects of floods, and that behavior will often depend on recent flood experiences. We do correct the data for the non-stationarity to some extent by normalizing not only for building cost changes (which is a better measure than inflation for our purposes), but also considering the changes in building stock in each area. This second data-driven approach provides an interesting supplement and additional benchmark to the (more common) first method, as it includes the historical evidence in the projections to the extent possible. We would like to point out that the approach followed in this paper can also be interpreted as a counterpart of the analysis in [27,28] for storm risk modeling in Austria, where the dependence structure bestowed on a detailed analysis of marginal storm risk was motivated and calibrated with wind speed fields. In the context of flood modeling, we also do not have sufficient historical data to convincingly fit a dependence model based on the data alone, so that we suggest here the use of a dependence model 'borrowed' and adapted from one of the river discharges for which sufficient historical data are available.
The paper is structured as follows. In Section 2 we provide details on the used data and data preparation. In Section 3 we present the used methodology in detail. We start by giving the rationale of identifying clusters that form the core of the hierarchical dependence structure. Subsequently, we adapt the max-stable processes as suggested in [1] for the dependence modeling of river discharges in river networks to the situation of Austria. This in particular asks for a further branching of hierarchical elements due to the complexity of its topography. We then show how the spatial dependence model that is calibrated on particular measurement stations in the river network can suitably be extended to the level of municipalities. After that, the marginal modeling with both the HORA and the EVA approach is described in detail and compared. Section 4 gives the results for the estimation of the parameters and discusses their implications on the concrete flood risk quantification for all the levels (municipality, state, and country), including the impact on the resulting flood risk diversification potential across these levels. Section 5 then provides some final remarks. Throughout the sections, some technical details as well as the statistical analysis underlying, accompanying, and justifying certain assumptions are delegated to the Appendix A.

Damage Data and Normalization
For assessing the magnitude of flood damages, we use the payments of the Austrian catastrophe fund. The fund can compensate private persons for the damage of catastrophic events. We collected the data of the payments for the provinces of Upper Austria (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013), Styria (1989Styria ( -2013, Lower Austria (2005-2013), and of Salzburg . Data for later years were partially available, but not used here since most of them were not reliable final numbers due to still ongoing settlement procedures. Although the time series of Lower Austria and Upper Austria are relatively short, they at least include the relatively big flood events of 2013 and, for Lower Austria, also 2005. We extracted the data for buildings (including inventory) and for the category "flood" (note, however, that some flood-relevant data points may be missing, as they are possibly stored in the categories "mudflow" (Vermurung) and "landslides" (Hangrutschung), respectively).
To make the damage data usable, we have to normalize them first. In [14] it was shown that the main driver of flood risk change is the change in urban sprawl. To correct for the latter, we correct the data for building stock. The Catastrophe Fund reimburses part of the costs for the rebuilding of the damaged residence buildings, but keeps track of the total costs, an aspect that we use in our analysis. For the normalization, we use the building value of all buildings within a municipality, defined by living space in m 2 times the construction cost per m 2 . To construct a time series of the building stock in a municipality, we use the data on the living space and living buildings in each municipality (data for the years 1971, 1981, 1991, 2001 were available on a spatial resolution of 5 km (number of flats in different categories) and on 250 m grid (number of buildings). For the period 1966 to 2001 these data were already gathered in the course of the former project [22]. For 2011 we collected the number of residence buildings in each municipality. To get an estimate of the living space in each municipality for 2011, we estimate the average living space per residence building and assume that this quantity stays constant for the time between 2001 and 2011. Between 1971 and 2011 we use cubic spline interpolation to generate the living space and beyond 2011 we extrapolate the function with constant slope, i.e., we assume the relative increase in buildings to be constant (since this concerns only the years 2012 and 2013, this assumption, in any case, has a minor impact on the results). Finally, we normalize the data with the building cost index to account for the aspects of inflation that are relevant to the considered costs.

Risk Maps
We use the risk maps provided by the HORA system, which provides the flooded areas for floods with return periods of 30, 100, and 200 years, respectively (see [21] for more details). These areas are also called the HQ zones and are available for the entire territory of Austria.

Discharge Data
For the modeling of the dependence structure, we use river discharges at river gauges and corresponding catchment basins, as well as a river network for Austria. Discharge data for gauges of Austria are obtained from "Hydrographischer Dienst Österreich" [29]. We only use stations with at least 30 years of data, which amounts to 431 stations. Data for the river network and the catchment areas are obtained from Copernicus [30]. Since the resolution of these two datasets does not coincide, we identify each discharge station with the point of the river network that is closest with respect to Euclidean distance. To obtain the basin of a discharge station, we first calculate all river segments upstream of the station and use the associated catchment areas from the Copernicus data set. Furthermore, we remove stations from the dataset that have the same catchment area (i.e., that are on the same river section), after which 425 measurement stations remain for our analysis. With the above-described data, the elevation weighted center of the basin for each discharge station is calculated for later use (as a digital elevation model we use a Digital Elevation Model (DEM) for Austria (at 10 m resolution inside Austria from [31]) and the EU-DEM v1.0 (25 m from Copernicus [32])).

European Settlement Map
For the distribution of buildings inside each municipality, we use the European Settlement Map (Release 2016) [33].

Dependence Modeling
We split the modeling of flood risk into two parts. The first part consists of modeling the marginal distribution of flood risk for each municipality. The second part is to use a copula to combine the flood risk of the different municipalities in an appropriate fashion. The choice of copula determines the dependence structure and is hence a crucial step in the modeling. Due to the small number of years with reported damage figures, this is a particular challenge here and cannot reasonably be based on the historical data alone. We, therefore, use the dependence structure inherent in discharge data (for which a much longer and more consistent time series across all of Austria is available) and use it as a proxy for the dependence structure of the flood damages. Such an assumption makes sense, as significant damage can only occur when an extreme discharge is present. At the same time, one has to keep in mind that the link between discharge amounts and flood damage sizes is not deterministic and prone to some inaccuracy. For the dependence structure of discharge levels, we focus on the dependence of extreme discharge only. The latter is the basis for the method proposed in [1], where an analysis for the upper Danube basin is provided. We will mimic this approach but apply it to the entire river network of Austria, for the purpose of which it needs to be extended. This involves the identification of clusters of discharge measurement stations and their combination. Altogether, the dependence modeling will be split into three steps. First, we apply a cluster analysis for regions with similar flood risk. Second, a calibration of the river discharge model will be performed, and finally, the model for measurement stations will be extended to a dependence model among municipalities. We describe these three steps now in more detail.

Cluster Analysis
A particular motivation to look for clusters comes from [34,35], where regions are defined for the additional consideration of the weather system. Our intention is to find clusters in summer floods, hence we will focus on May to October (cf. [36]). Note that this period is different from the one chosen in [1].
The first step in the cluster analysis is to de-cluster the time series of daily mean discharge. We consider all discharge events that exceed the 90% quantile and apply a 9-day rule, in the sense that events that occur within 9 days are not counted as independent events (such a rule was also used in [1]). This procedure provides us with 848 remaining events, where, in at least one catchment, a peak discharge level above the 90% quantile appeared. Note that an alternative declustering method was used in [37], where two peaks were considered independent if the number of days between them exceeds 4.0483 + log (A), where A is the size of the catchment area and the minimal discharge between the peaks is less than three quarters of the smaller discharge of the two peaks. Since in our case the median of the logarithm of the catchment area of the considered measurement stations is 5.1, our 9-day rule is in fact reasonably in line with that alternative de-clustering method.
We want to cluster measurement stations with respect to the likelihood that extreme discharge levels appear together. As a distance measure, we use the probability that two peaks occur at the same time. We approximate this by the number of common peaks of the 10 largest discharges across the entire period. For the calculation of the clusters, we use hierarchical clustering (the hclust-algorithm in R) with the method ward. D to identify the clusters. This results in 7 clusters in Austria which are depicted in Figure 1. Concretely, this map depicts the altitude-averaged centers of each river basin colored according to the cluster it belongs to. For more details and comments on the clustering, we refer to Appendix A.

The Dependence Model for River Discharge
We need to model extremal dependence in two-dimensional space and for that purpose, max-stable processes (which are an extension of extreme value distributions to the functional setting) have successfully be applied in the extreme value statistics literature (see e.g., [38][39][40]). Whereas the classical approaches were based on assuming that the degree of dependence between two locations depends on the usual Euclidean distance, it was recently proposed in [1] to replace the Euclidean distance by a distance concept that considers flow-connectedness in a river network and the possible exposure to common meteorological events. We will follow that latter approach and adapt it to our setting. In [1] it is assumed that the marginal distribution of the discharge at a point of a river follows a generalized Pareto distribution and that the discharges normalized to unit Fréchet margins follow a (Brown-Resnick) max-stable process defined by Here t represents the location on the river, U i are the realizations of a Poisson point process with intensity u −2 du, and W i (t) is a Brownian motion with mean 0 and variance function σ 2 (t). It was shown in [41] that the distribution of such a process depends only on the covariance function. For any two locations s and t on the river network, the latter can be defined through a kernel function Γ(s, t) by A useful measure for the extremal dependence between two random variables X and Y with distribution functions F X and F Y is the so-called extremal coefficient see for instance [42]. Note that 1 ≤ θ ≤ 2, where θ = 1 corresponds to complete dependence and θ = 2 to independence. It has been shown in [41] that where Φ(x) is the cumulative distribution function of the standard normal distribution, so that the kernel function Γ(s, t) directly 'encodes' and parametrizes the extremal dependence between two locations s and t. Small values of Γ(s, t) correspond to strong dependence and large values correspond to weak dependence. In order to incorporate an appropriate distance concept for the geometry of river discharges, [1] suggested a function Γ(s, t) of the form The first term in this sum considers hydrological distance in the river network and the second term takes into account dependence resulting from the geographical structure of the river basin and spatially distributed meteorological variables. The two components are combined with weights λ Riv and λ Euc . The first term is only different from its maximal value 1 if t is downstream of s or s is downstream of t. Here π s,t is the fraction of the integrated altitude of the catchment basin of the upstream location (which can be seen approximately proportional to the average runoff accumulated in this catchment), divided by the integrated altitude of the (larger) catchment basin of the downstream location. H(t) is the elevation-weighted center of the basin for discharge station t, and its use is motivated by the assumption that precipitation is higher for areas with higher elevation, and hence the elevation-weighted center is an approximation for the precipitation-weighted center. The matrix takes into account that distance in different directions might not have the same influence on dependence, due to the topography and typical tracks of weather events. For example, if typical events have a track from west to east, then two locations of given distance that lie in east-west direction will be stronger dependent than for a north-west direction. The norm ||.|| is the Euclidean norm, and the power α is another constant to be calibrated from the data. Finally, d(s, t) is the distance of locations s and t along the river network and the constant τ > 0 is the maximal distance between two locations until which Γ(s, t) is influenced by changes in d(s, t), and is estimated from data. For a more detailed motivation for this particular choice of Γ(s, t) and further background of its justification we refer to [1] and the references therein. Whereas [1] only considered the upper catchment area of the Danube, our cluster analysis of the previous section indicates that there might be different processes driving extreme discharge in areas north and south of the Alps. We will, therefore, modify the above process for our purpose to consider the entire region of Austria. Concretely, we will use various distance measures and combine them. The cluster analysis of the previous section is used to group individual stations together into m groups, and we assume that there are n underlying distance kernels Γ i -s that lead to extreme discharges. The motivation to use different kernels is that some processes leading to extreme discharge are more local in nature, while others affect a large number of municipalities. Further, we assume that there exist weights λ i,j i = 1, . . . , n, j = 1, . . . , m such that n i=1 λ 2 i,j = 1 and where j(s) is the cluster to which the location s belongs. Note that when s and t are in the same group, then For the calibration of the resulting dependence model and a discussion and justification for our final choice of m = 7 and n = 4, we refer to Appendix A.

The Dependence Model for Municipalities
The approach in the previous subsection provides an approach to model the spatial dependence between individual locations, which we now want to extend to relate the municipalities in Austria. For that, we have to associate the municipalities with river segments. To each municipality, we assign the river segment that is associated to the maximum number of built-up area, using the European Settlement Map in a HQ zone. River distances are then estimated from the end node of a river segment to the start node of another river segment. For the river distance, we additionally assume that if two segments have a common predecessor, then we use the average distance to this common predecessor as the distance between the river segments. This provides us with the river distance between municipalities. The center of the catchment basin associated with a municipality is the weighted average of the center of catchment areas of river segments in the municipality, where the weight is the percentage of built-up area in the HQ-zones associated to the river network. Finally, we assign the municipality to the group for which the mean distance to the next three stations is smallest (distance is between catchment areas).

Application of the Dependence Model
The resulting dependence model is based on max-stable processes. Since explicit calculations for such processes are not available, we will use Monte-Carlo simulations to generate outcomes and estimate resulting distribution functions and quantiles. The simulation of instances of the underlying max-stable processes is time-consuming; we therefore decide for 100,000 realizations as a reasonable trade-off between accuracy and computation time and power (as we will see later, this leads to confidence intervals for the quantiles of the yearly aggregated damage between 1.5% and 5%, which is quite acceptable). For the calculation of yearly average damage, the error will be higher, but for the latter one can apply the results for independent municipalities (as the average values do not depend on the dependence structure), which is what we will do. Finally, for the simulation of the max-stable process, we use an approximative algorithm, where values below 1/log(2) are not necessarily accurate, but since our focus is the tail, this inaccuracy is not significant for our purposes.

Risk Map-Based Model
We will improve a model that was used in [21,22]. The HORA system provides the flood areas of floods with return periods of 30, 100 respectively 200 years. These areas are called the HQ zones and are available for the entire territory of Austria. The main idea in [21] was to simulate dependent return periods (which marginally can be seen as realizations of a Pareto random variable with tail parameter 1), and then for each return period in a municipality to estimate the caused damage to the flood with that return period. For that purpose, we first intersect the HQ zones with data on the building stock to get the number of affected buildings in a flood with recurrence level 30, 100, respectively 200 years.
For each municipality, these data are then interpolated to get the affected number of buildings for floods with other return periods. Finally, a damage function provides the estimated damage, given the number of affected houses and the return period of the considered flood. These calculations lead to an estimate of the marginal distribution of each municipality. To get an estimate of the distribution of the aggregate (countrywide) damage of a flood event, we then apply in a second step the dependence model described above. We now list some further details on the concrete implementation.

Interpolation of HQ-Zones
We use the interpolation method of [22]. That is, in the HQ30 zone we use an interpolation according to a beta distribution with parameters α = 4, β = 1.5, and in all other HQ zones we use a linear interpolation of buildings (For Vienna, we assume that all buildings are in an 'artificial' HQ1000 zone and use linear interpolation to get the number of affected buildings for a flood with smaller return period. Note that the flood risk for Vienna has a very different nature in view of the regulated Danube (and its magnitude is almost negligible when compared to the magnitude of flood risk for the rest of Austria), so that we will include it in this study only in a very limited form).

Damage Function
The damage of a building is modeled by a lognormal distribution. As in [22,23], we assume that there exists a parameter set for small damages and one for large damages. We model the damage of a building with the small damage distribution (which can be interpreted that only the cellar of the building is affected), unless the flood is unusually large. The latter event is here defined to be one where the flood covers the complete next HQ zone (i.e., a building in the HQ30 zone is hit by a flood with return period larger than 100, or a building in the HQ100 area is hit by a flood with return period larger than 200). In that case, one assumes that the damage is considerably larger since the entire building is damaged, and we then use a different loss distribution.
We use a virtual overbuilt area of 100 m 2 to characterize a building (note that this value is only important for the number of Lognormal random variables used, and, as discussed in [22], this assumption is not essential for the magnitude of the results). We then approximate the expected damage for each municipality with the available damage data. For that purpose, we first approximate for each municipality the loss degree of damage per overbuilt area by the ratio of the respective maximal yearly damage (in Euro 2015) among the n years of data available in that state and the overbuilt area in the HQ30 zones of that municipality. Note that for n years of damage data, the probability that the maximal damage within these n years is at least as large as a marginal event with a return period of 30 years is given by qm 30 = 1 − 1 − 1 30 n . We hence take all these ratios across municipalities in the state and pick the value according to the (empirical) qm 30 -quantile, which we assign as an approximation of the expected damage per overbuilt area of a 30-year event. Finally, we take the average of this quantity of all four states for which data are available, and fix this as the expected damage per overbuilt area for small damages across all of Austria. For the final specification of the parameters of the Lognormal distribution, we then choose the variance such that the coefficient of variation is the one for the choice in [21] that was based on further empirical estimates of [16]. Furthermore, we use insight from [16,21] to choose the distribution parameters for the large damages (where the entire building is damaged) by a proportional up-scaling of the distribution for small claims. This entire procedure may be considered somewhat ad-hoc, but is still a rather reasonable approach to use all available damage data as well as insight from other studies to estimate one underlying damage function that is then applied across Austria (one should also keep in mind that in case of many damages, the law of large numbers suggests that finally mainly the expected values matter, which is also observed in the analysis). After normalizing with the estimated damage for 100 m 2 , this eventually leads to the following parameters for the Lognormal distribution of the damage per 100 m 2 : µ s = 2460, σ s = 1480, µ l = 5676, σ l = 3517.

Extreme Value Model
While the method of the previous subsection used the claim data only for modeling the average damage per building, not for its entire distribution, we now proceed to an extreme value approach that models the tail directly from the claim data.
A classical result in extreme value theory (see e.g., [43]) states that for a random variable X in the maximum domain of attraction of an extreme value distribution, it holds that where σ u . is a normalising constant depending on the threshold u. This limit motivates the use of a generalized Pareto distribution as a model for the distribution of X|X > u, which we will use on the level of a single municipality. For each municipality i and corresponding claim size random variable X i , we, therefore, assume that there exists a threshold u i such that the distribution of X i |X i >u i is generalized Pareto with parameter σ i and ξ i . For simplicity (and analogously to [1]) we assume that the parameter ξ i is the same for all municipalities, i.e., ξ i = ξ ∀ i. Further, we assume that the other parameter σ i is a linear function of the fraction u i of built-up area in the zone HQ100 among the complete built-up area of the municipality and also depends on the state in which the municipality is located. Finally, for the threshold we assume u i = u i u, where u is another fixed quantity, which seems reasonable as the threshold can be higher when there are more buildings at risk. Given the data and thresholds u i , the estimation of ξ and σ i is a standard procedure (cf. [44]), and we will use the function gpd.fit from the R-package ismev. For the estimation of u i there exists no standard method. We split the task of estimating u i into two subtasks. First, we will estimate the u i from the damage data and then use a Hill-type plot for different values of u. To get a complete model for the X i , we finally need the probability P(X i > u i ) of having a claim larger than the threshold, which we estimate using the data of Salzburg and Styria only (since the time series for Upper and Lower-Austria are too short for this purpose) and the distribution of X i |X i ≤ u i , for which we simply use the empirical distribution. We refer to Appendix A for further details.

Main Results
In the previous section, we have developed two (statistical) models for the assessment of flood risks in Austrian municipalities. The first (HORA) model can be applied to all municipalities in Austria. The applications of the second (EVA) model is restricted to municipalities where claim data exist. We then use Monte Carlo simulations to obtain estimates for the average yearly damages and the quantiles. On the marginal distribution level, the differences between the two approaches are discussed and analyzed in detail in Appendix A. We will then apply our dependence model, and also complement the final results by respective estimates under the assumption that flood risk of the municipalities is independent (no dependence) and that it is comonotone (complete dependence cf. [45]). This allows to illustrate best-and worst-case bounds concerning the sensitivity of the results with respect to the implemented dependence structure. The comparison with extreme dependence scenarios is also helpful to assess the presence of diversification potential of flood risk.

Results for Austria
In Table 1 we present the average yearly damage as well as the expected damage for a flood event with return period of 30, 100, 200, and 300 years for Austria as obtained by the HORA model. Note that the expected damage of a flood with a return period of 200 years is particularly interesting in terms of risk management (since the 200-year value at risk is the relevant quantity for the required safety capital for insurance in the European Union). For the average yearly damage, we get values between 80.8 mio. € and 86.7 mio. €. Theoretically, all models should provide the same average value, but the Monte Carlo estimate naturally has some variance. Since this latter variance will be larger when spatial dependence is involved, we favor here the estimate based on independent flood risk between municipalities and will, therefore, use the corresponding result of 84.2 mio. € for the yearly average damage as the benchmark for the later comparisons.
A flood with a return period of 200 years has an expected damage of about 1.9 billion €, which is approximately 22 times the average yearly damage. In other words, one has to save the average yearly damage of 22 years as solvency capital for a 200-year flood event. Comparing the results for the estimated dependence structure with the one for independence and comonotone dependence, we see that the estimated quantiles are much closer to the comonotone case, illustrating and quantifying the strong spatial dependence for flood risk in Austria. Note that for given marginal distributions for each municipality, under comonotonicity, the quantile for the aggregate risk is the sum of the individual quantiles, see e.g., [46]. In other words, since the quantile of the loss distribution is linked to the solvency capital requirement for a flood insurance activity, the respective figures in the third row of Table 1 can be interpreted as the overall amount of needed (stand-alone) solvency capital, if all municipalities would insure their flood risk on their own, whereas the figures in the first row reflect the respective capital requirement when the flood risk is pooled on the aggregate Austrian level. The ratio in the last row is hence an indication for the diversification potential in flood risk (e.g., for the relevant 200-year return level, one can reduce the amount of required solvency capital by pooling all flood risk in Austria by 38%). At the same time, given the many individual municipalities involved, this also reflects the strong spatial dependence present in flood risk, as one a priori may expect a larger reduction on the aggregate level. Table 1. Average yearly damage and estimated damage (in mio. €) for floods with return periods of 30, 100, 200, and 300 years for Austria (based on the HORA model with the estimated dependence structure, and in comparison with the case of independence and comonotone dependence). Numbers in parentheses are the 95% confidence intervals from the Monte Carlo simulation (For the calculation of a two-sided distribution-free (slightly conservative) confidence interval for a quantile q, we use the formulaF −1 q ± z 1−α q(1−q) n , where z 1−α is the (1 − α)-quantile of the standard normal distribution andF −1 is the inverse of the empirical distribution function, cf. [47] for a motivation of this formula, where we approximate the binomial distribution by a normal distribution. Note that the resulting confidence interval is not necessarily symmetric around the estimate).

Results for the States of Austria
In Table 2, we break down the results to the level of individual states in Austria. We again present the average yearly damage as well as the expected damage for a flood event with return period of (for brevity now only) 30 and 200 years for all states of Austria as obtained by the HORA model (While we did include Vienna for the results of the entire country of Austria for completeness, as mentioned before we do not provide a detailed separate analysis on its state level due to its very different characteristic in view of river regulations. Note that its average yearly damage measured akin to the other states in Table 2 would amount to 0.24, the estimated damage for a return period of 30 years is virtually 0 and for a return period of 200 an analogous estimate would be 4.35 mio. €, clearly indicating that flood risk for Vienna is negligible on the Austrian scale). We provide the results again for the three used dependence structures (independence, comonotone and the estimated dependence structure). We see from the ratios that the diversification potential within each state is now only around 20% for the 200-year event quantile levels (which can be translated into possible solvency capital reduction) compared to stand-alone management of each municipality. It was to be expected that this amount is reduced relative to pooling across all of Austria (since spatial dependence within each state is naturally stronger than across the entire country). However, the concrete numbers nicely allow to quantify how much one gains-on average-by pooling flood risk across the entire country (note that the financial management of flood risk in Austria is to quite some extent politically a competence of each state separately, so that policy recommendations towards collaboration may originate from the present considerations). Concretely, summing the 200-year return period quantiles of all the individual states (here even including the small value for Vienna), we arrive at an over-all potential solvency capital requirement of 2.32 bio. €, which can be compared with the 1.88 bio. € obtained for the entire country of Austria above. That is, pooling the risk of the states on the Austrian level leads to a solvency capital reduction of about 500 mio. €, when measured in these terms. Comparing these results further with the ones of independent risks, we may still conclude that the diversification effect for flood risk is altogether relatively small, and that flood risk in Austria is quite concentrated. When comparing to the respective damage figures of the 30-year return period, one observes that this concentration is considerably stronger in the tail (i.e., for higher quantiles).

Results for Municipalities with Claim Data
For those municipalities where claim data from the past are available, we can now also apply our second (EVA) model, which provides a sort of local refinement. In the present study, we could use data for municipalities in Lower and Upper Austria, Styria and Salzburg. Before we proceed to the presentation of the results, let us consider some empirical statistics for these data (which are normalized with respect to the building cost index) in Table 3. For purposes of comparison, we first present in Table 4 the respective results of the HORA Model for the aggregation of all the municipalities with available claim data. From the average loss amounts in that table one deduces that more than 60% of the flood risk (in terms of expected claim amounts) in Austria is in fact concentrated in those municipalities. When comparing the results of the HORA model with the empirical claim statistics in Table 3, one can sum the empirical averages of the annual claims of the four states, which in fact quite well corresponds to the model average in Table 4 (note that for the latter the historical claim sizes were only used in a very peripheral way when estimating the mean damage of individual buildings!). Also, given the almost 30 years of claim experience for Styria and Salzburg, one could interpret that the maximal claim amount over the available years is a reasonable approximation to the 30-year quantile (although with a considerable uncertainty). For Upper and Lower Austria, less than 10 years of data are available, but there we know that we have one (Upper Austria) and two (Lower Austria) historically big flood events in the data, so that the maximum there may also be an indicator towards the 30-year quantile. Despite the hand-waving nature of this interpretation, one sees that indeed the magnitude of the 30-year quantiles of the model are in a reasonable range when compared to the empirical data. Table 5 now presents the respective results for the EVA Model for the same municipalities. One sees that the estimated average annual damage is about 30% larger for the EVA model than for the HORA model. Looking at the resulting numbers for different return periods, we see that the EVA model gives considerable more weight to the tail of the distribution (which may be most easily visible when considering the numbers under the independence assumption). Yet, except for the 30-year quantile, the quantiles of the two models are relatively close (especially when one compares the ratio of the quantiles over the mean (16.5 vs. 10 times the mean for the 100 year-quantile, 24.4 vs. 21.1 times the mean for the 200-year quantile, 30.7 vs. 29.3 times the mean for the 300 year-quantile). Nevertheless, in absolute numbers the increase in expected damage as function of the return period is considerably bigger for the EVA method. Table 5. Average yearly damage and estimated damage (in mio. €) for floods with return periods of 30, 100, 200, and 300 years, for the municipalities with claim data. The damage is calculated with the Extreme Value Analysis (EVA) model and with three different dependence structures between the municipalities (the estimated dependence structure of this paper, independence and comonotone dependence). Numbers in parenthesis are the 95% confidence intervals from the Monte Carlo Simulation.

Discussion and Conclusions
In this paper, we made two main contributions to the flood risk literature. Firstly, we introduced a new approach to model the dependence structure that can be used to aggregate the information of existing risk maps to larger areas. Concretely, we suggested a spatial dependence structure that is calibrated with publicly available data of river discharges and the river network. Hence it is easily applicable to any region in Europe without the need of a hydrological model. The second contribution is to use historical claim data of municipalities for the estimation of future flood risks to obtain an alternative estimate which can either be used directly or as a benchmark for other approaches. We exemplified our approach for the country of Austria.
The choice of a dependence structure is a crucial step when aggregating individually modeled risks onto a collective level. For flood risk applications, a reliable dependence structure allows to aggregate the risk of flood damages from individual municipalities to bigger areas like states or the entire country. The resulting aggregate distribution can then be used to estimate important quantities for risk management like the required solvency capital (e.g., the expected damage of a flood with return period of 200 years in the regulatory framework of the European Union). Comparing the required solvency capital for a calibrated dependence structure with the corresponding numbers for comonotone risks subsequently quantifies diversification effects in the process of flood damage aggregation. In this paper, we developed a concrete and detailed case study for Austria. We first used the HORA risk map available in Austria and intertwined it with the building structure at a fine resolution. We then estimated loss distributions for each municipality in Austria for any flood event with a given return period. The obtained marginal flood damage distributions per municipality were then combined with a spatial dependence model based on a max-stable process to obtain the distribution of flood risks on various levels of aggregation. Here the dependence model was calibrated using the joint occurrences of river discharge levels in the Austrian river network. Monte Carlo simulations of joint appearances of return periods were then implemented to obtain flood loss distributions for each municipality and its aggregations.
The resulting model leads to a projected reduction of required solvency capital for the individual states of Austria between 29% for Lower Austria and 13% for Vorarlberg when compared to assuming the risks separately on the individual municipality levels. One also obtains that the states of Salzburg, Upper Austria, Vorarlberg and Burgenland have a higher concentrated risk than the other states. We have also seen that a further aggregation onto the country level is definitely worthwhile and leads to an overall reduction of required solvency capital of about 38% when compared to the sum of capital requirements of all individual municipalities. We additionally compare the results with the ones under an independence and a comonotone dependence assumption between municipalities, which reveals that the actual situation is much closer to comonotonic dependence than to independence, indicating that flood risk in Austria is particularly difficult to diversify and needs a substantial amount of scrutiny in its management.
As a complement to the HORA model, we also developed a model based on real historical claim data for Austria and an extreme value analysis (EVA) approach to fit the suitably normalized flood losses for each municipality. Firstly, the actual claim data allowed to challenge the findings of the HORA model as the latter is not based on historical data, suggesting consistency for three of the four states with claim data (although this comparison between the actual claims and the HORA estimates should not be over-emphasized given the limited amount of years with available claim data).
Combining the obtained marginal distributions in the EVA approach again with the same spatial dependence model as for the HORA method then led to an alternative approach to quantify aggregate flood risk on a state and country level. Compared to the earlier study [21] of flood risk in Austria, the present contribution not only uses a considerably improved dependence model, but also much more refined data for the marginal behavior, so that we were able to sharpen those previous results considerably. The new estimates suggest a reduction of the projected flood damage based on historical claim data compared to the earlier study. In comparison to the HORA model, the EVA approach leads to higher expected annual damages and more mass in the tail of the distribution. Nevertheless, in the range of flood damages for a return period of 200 years, a quantity that is relevant for risk management purposes related to Solvency 2 in the European Union, the results for the HORA and the EVA model are quite similar. This is remarkable, as the two approaches use empirical data and relations in quite a different way, indicating some robustness of the obtained magnitudes.
As mentioned in the introduction, the dependence modeling approach proposed in this paper can be applied to any other country, if marginal modeling experience and river discharge information is available, and it could be nice to see how the respective estimates for average and extreme flood losses compare to findings of other studies, often relying on hydrological models. While the present study took the non-stationarity of flood risk into account in terms of inflation correction, changes in building structure over time, and adaptations to the risk maps, it will be interesting in future research to also combine this approach with longer-term experience and projections due to climate risk. This could for instance include the consideration of time-dependent parameters in the max-stable components based on longer-term experience of flood occurrences (see e.g., [48]) or even paleo-floods (see e.g., [49]). Also, one might lift the granular dependence model approach suggested in this paper up one more level to the entire river network of Europe and in that way refine the somewhat coarse Clayton copula approach underlying the interesting findings in [4].
We would like to finish with some remarks on uncertainties. While one can quantify the error that comes from the numerical evaluation of the quantiles by Monte Carlo estimation from the specified model (and we have done so in the respective tables), the model choice itself naturally is another source of uncertainty for final estimates that cannot be quantified in a straight-forward way, and is, in fact, present in any modeling approach. Uncertainties of the HORA model are already discussed in [21], and mainly are the number of buildings affected by a flood for a given return period, the average damage of a building in a flood with a given return period and the used dependence structure for aggregating the risk of individual municipalities. Among these, the number of buildings affected by a flood of a given return period is probably the major source of uncertainty about the results, but since this error depends on the geographic characteristics of the municipality, it is hard to quantify with statistical methods. However, in this paper, the HORA model was complemented by a second, different approach using EVA techniques. The latter has considerable uncertainties itself, which mainly stem from the extrapolation nature of the method, but the comparison of the two different methods in fact assesses to some extent model risk in the calculations. Finally, adding the scenarios with comonotone and with independent risks in the analysis of this paper enabled us to establish best-and worst-case bounds for our results concerning the underlying dependence assumptions. In particular, the respective numbers can help to decide, in future studies, whether to focus more particularly on these aspects of the modeling procedure.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A.
Appendix A.1. More Comments on the Cluster Analysis The 7 clusters resulting from the clustering algorithm were shown in Figure 1. As a performance measure, let us consider the resulting probability that extreme discharge levels occur at the same time within and outside the cluster. In Table A1, the median of the probability for the joint occurrence of extreme discharge levels for arbitrary pairs of stations is provided. We can see that while in the regions north of the Alps (groups 1, 2 and 5) that median probability inside the groups is between 0.5 and 0.6 (meaning that 5, respectively 6, out of the 10 biggest discharge levels of the last 30 years happened at the same time), for the other groups the median probability is smaller, so the dependence is lower. A possible explanation is that in those areas extreme discharges are mainly forced by rather local phenomena. One can compare the resulting cluster regions to the ones of [34]. The cluster regions for the Eastern lowlands, Northern lowlands and Southern alpine regions are quite similar. On the other hand, parts of the Northern alpine regions in the middle of Austria are now in a group with the Northern lowlands. If we compare our results with the figure after Table 3.3.2.2 in [35], we see that our estimated "Group 3" looks similar to the N-Stau region, only shifted to the east. Other regions from [35] are also recovered but shifted to the east, whereas yet others are split into different sub-regions here.

Appendix A.2. Calibration of the Dependence Model to River Discharge Data
To decide about a good choice for the number m of groups and the number n of different distance kernels Γ i -s, we tried various possibilities and finally chose the best-performing among the considered options. For the number of groups, we tried with either 2 clusters (where the clusters (1,3,5) form Group 1 and clusters (2,4,6,7) form Group 2) or keep all original 7 clusters. For the number of different kernels (and hence dependence structures), we tried with either two or four (note that four is the number of different weather areas as identified in [35]). In addition, we also tried to use the coarser original process as used in [1], in order to check whether the proposed refinement is worthwhile.
For each of these choices, the results of the maximum-likelihood estimation procedure are provided in Tables A2 and A3 (note that we used the coordinate system that is also used in the Copernicus data set, while river distances are in km). Table A4 gives the corresponding negative log-likelihood and BIC values. We observe that using different groups significantly increases the likelihood and BIC when compared to the original model. Further, using 4 Γ i -s and 2 groups has a smaller likelihood than using 2 Γ i -s and 7 groups, although the latter uses fewer variables. In terms of BIC, the best model is the model with 7 groups and 4 dependence structures. If we look at the different estimated distance kernels Γ i in Table A4, we can observe that especially the values of τ differ among Γ i . A small value of τ indicates that river distance is only important for close stations, and the major weight is on the distance between catchment areas -this might indicate the importance of large events for multiple regions). In contrast, a large value of τ means that the river distance is more important, which might be the case for smaller scale events. For m = n = 2, one sees that for the stations in Group 2 the Γ i with the bigger τ have a higher weight than for Group 1, i.e., north of the Alps large scale events might be more relevant than in the south. As in [1] we want to test how well the model fits the data. For that purpose we use the extremal coefficient θ as defined earlier and its connection to the kernel Γ i . To get an estimate of θ(s, t) for two stations s and t, we can either use the madogram (see [50]) for yearly maxima, or an approximation of (1) where u is chosen such that the largest 10% of the values of each discharge station are included. Note that with the madogram it is possible to obtain values outside the interval (0,2). This estimate can then be compared with the value we obtain when substituting the fitted parameters into (2). Figure A1 plots the values of θ(s, t) across all pairs of stations using the fitted model parameters (for m = 7 and n = 4) against the respective empirical estimates as described above. Figure A2 plots the analogous points for the model with only one group and dependence structure of [1]. Comparing the two, one sees that using several Γ i -structures in the model slightly improves the correlation of the fitted and empirical θ-values. Within each figure, one observes that using Formula (1) leads to a significantly better accordance than using the madogram for the estimation of θ. Due to the better fit when using the 7 groups and 4 dependence structures, this is what we have finally chosen for the analysis.  Let us now discuss the fitting of the threshold multiplier u i to the data. The threshold depends on the fraction of buildings in the HQ100 zone. We fit the maximal yearly damage to the covariate. For that purpose, we sort the municipalities and split them into 20 groups. For each group, we estimate the median of the covariate and the median maximal damage. We then calculate a linear model for these two quantities. The adjusted squared error is 79.5, see also Figure A3.
To determine the values of u i , we first use the estimate of the linear model with the built-up area in the HQ100 zone as covariate for all municipalities and divide it by the median of the estimated values resulting in an estimateû i . In a second step, we calculate the median u s s of the fraction of the maximal yearly damage andû i for all states s. Finally, we estimate u i =û i u S s(i) u s , where s(i) is the state of municipality i and u s is the average of the values u s s . Figure A3. Median of the covariate (fraction of built-up area in the HQ100 zone against the median of the normalized damage.

Appendix A.3.2. Threshold Selection and the Estimation of Parameters
We next discuss the estimation of the threshold u and the parameters. For municipalities in the considered areas, we provide in Figures A4 and A5 a plot of the fitted generalized Pareto distribution as a function of the number of extremes used in the estimation. We apply this procedure for Styria, Salzburg, the two states combined as well as all municipalities with data. In addition to the plot, we also provide the lines when choosing the data such that the probability of extremes is 95% and 98%, respectively, the choice that we would make because of the plot. To estimate P(X i > u i ), we estimate the empirical probability P(X i > u i ) for all municipalities in Styria and Salzburg, and use the average over all empirical probabilities as the estimate for P(X i > u i ). We list the resulting u i , estimated ξ and P(X i > u i ) in Table A5. One can see that the method works well for Styria, but is also quite satisfactory for the other areas. Except for municipalities in Salzburg, the estimated ξ is between 0.85 and 0.89, while the estimation for Salzburg is 0.78, i.e., slightly lower. These results suggest that the mean of the damage exists, but not the variance, which is in line with earlier extreme value studies for flood risk in Austria (see e.g., [26]), but now with a much more refined model. The rather similar obtained values for ξ is a consistency check across the different areas. The estimation of the threshold u shows that for Styria the threshold is much smaller than for the other areas. We can see that the probability to exceed the threshold ranges from 7.2% (every 14 years) for Styria to 3.3% (every 30 years) for all municipalities.   In order to investigate how well the extreme value model describes the damage data above the threshold, we study in Figures A6-A9 QQ-and PP-plots with 95% confidence intervals under an independence assumption. We observe that the estimated models work quite well for the areas where the parameters were fitted.    HORA method. We are therefore interested to compare the implied damage distributions of the two methods on the level of individual municipalities. Figures A10 and A11 depict scatterplots on the log-scale of the expected yearly damages as well as the 30, 100, 200 and 300-year quantiles. With the exception of some outliers, there is a relatively good correspondence between the two methods for the expected yearly damages. Concerning quantiles, the two methods lead to rather different results for the 30-year level, for the 100-year return level the HORA method has a tendency to overestimate the EVA method for municipalities with higher quantiles, whereas it is remarkably close for the 200and 300-year quantiles. This is further illustrated in Figure A11, which depicts the R 2 for a regression line between the EVA and HORA quantiles across municipalities as a function of the return period. One can observe that the two methods have rather similar results for return periods between 160 and 500. At the same time, there is a strong downward peak around the 100-year quantile, which might be explained by the fact that in the HORA model, for events with a return period of 100 years or higher the damage distribution for the buildings in the HQ30 zone changes from the "small"-damage distribution to the "big"-damage distribution, whereas the EVA method (implicitly) estimates a smoother transition. Finally, we recall that in the EVA method the tail of the damage distribution for the municipalities basically depends on the fraction of buildings in the HQ100 zone, and in the HORA model the average damage of a house is fitted using the 30-year quantile of the damage distribution. It is hence not surprising that one finds a reasonable degree of linear dependence for different municipalities. However, it may be considered surprising that the absolute value of the mean (with the exception of some outliers) and the absolute value of the higher quantiles correspond so well.