1 Introduction

Climate sciences foresee a future where extreme weather events could happen with increased frequency and strength as a consequence of accumulating greenhouse gases in the atmosphere. Specifically, climate change would favour extreme precipitations to occur with higher frequency and intensity, causing more riverine, coastal and flash floods, which are already the main source of losses in the world (UNISDR 2015). The Mediterranean basin is considered a hot spot in terms of climate vulnerability and risk (Blöschl et al. 2019). The higher probability of these events to happen, coupled with the effects of land-use change and, in particular, of urban growth that increases soil sealing and, hence, water runoff, in turn determines increased potential impacts (Kundzewicz et al. 2014; Slater et al. 2015; Viero et al. 2019; Winsemius et al. 2016). The expected consequences are an increase of fatalities and injuries, but also of economic losses in urban areas, commercial and productive sites, infrastructures and agriculture, and they depend on the specific vulnerability of the area hit by the event.

The definition of vulnerability can vary significantly depending on the community addressing it. Two major conceptualizations of vulnerability can be found in the literature (USAID 2014): contextual vulnerability, focusing on the factors that define the ability to withstand and recover from a shock, and outcome vulnerability combining “information on potential climate impacts and on the socio-economic capacity to cope and adapt.” (Füssel 2010; O’Brien 2007).

The IPCC with its report on managing the risk of extreme events (IPCC-SREX 2012) made an effort for harmonizing the definition of vulnerability among the communities of Disaster Risk Reduction (DRR) and Climate Change Adaptation (CCA) and defined it as “the propensity of predisposition of exposed receptors to be negatively affected by hazard events”.

Having defined the concept and relevance of vulnerability, the need for methods for its quantification emerges. Vulnerability assessment (VA) has been covered extensively in the literature, however without a consensus on what constitutes a best practice (USAID 2014). Four major conceptual frameworks for VA can be identified (USAID 2014): (i) the Intergovernmental Panel on Climate Change (IPCC) framework with vulnerability as a function (i.e. outcome) of exposure, sensitivity and adaptive capacity (IPCC 2001; Parry 2007); (ii) the extended vulnerability framework (Turner et al. 2003; Birkmann 2006), which develops upon the IPCC framework including a broader array of place-based contextual factors, conceptualizing also the main feedbacks among elements; (iii) the livelihood framework (Carney 1998) developed by the United Kingdom’s Department for International Development, which describes five capitals (natural, social, financial, human, physical); (iv) the SREX framework mentioned above (IPCC 2012), building on contextual vulnerability, which becomes one of the three variables determining risk, together with exposure and hazard.

In this work we adopt the IPCC-SREX definition and assessment approach to vulnerability, following the KULTURisk (KR) framework, which considers vulnerability as one of the three independent variables of risk (Giupponi et al. 2013; Mojtahed et al. 2013). KR considers both the physical and human dimensions of vulnerability (see Figure S1). The physical dimension is given by the susceptibility of exposed man-made structures, namely their predisposition of being negatively affected by hazards. The social dimension is made of two components: adaptive capacity (AC) (ex-ante) that is “the ability to anticipate and transform structure, functioning, or organization to better survive hazards” and coping capacity (CC) (ex-post) that is “the ability to react to and reduce the adverse effects of experienced hazards” (Gain et al. 2015).

To assess social vulnerability through its two components of adaptive and coping capacities, several socio-economic characteristics need to be taken into account and measured by means of indicators, such as age, gender, race, overcrowding, ethnicity, social class, unemployment rate, immigrant status, density and quality of the built environment, land use, housing tenancy and the presence of informal support networks (Borden et al. 2007; Burton and Cutter 2008; Cutter et al. 2000, 2003; Fekete 2009; Finch et al. 2010; Masozera et al. 2007; Rygel et al. 2006).

The results of VA are quite often presented as maps showing the spatial variability of social vulnerability as a non-dimensional index, highlighting community needs at different level of disaster response and mitigation, useful for emergency management or risk reduction planning (Boruff and Cutter 2010; Cardona 2005; Chakraborty et al. 2005; Wood et al. 2010).

More rarely vulnerability is assessed as a dynamic variable, which vary not only across space, but also in time. In general, being vulnerability the outcome of the combination of variables, which evolve over time, it should be expected to vary over time too. Moreover, the dynamics of vulnerability is intertwined with the occurrence of extreme events, with the direct experience of flooding influencing the perception of future risk and the responsibility to act, increasing the propensity to believe to be affected again in the future, the sense of personal duty and willingness to adapt (Adger et al. 2013; Di Baldassarre et al. 2015). Especially when there is a frequent occurring of flooding, this leads to a decrease of social vulnerability, i.e. the so-called adaptation effect (Mechler and Bouwer 2015; Wind et al. 1999), with the society gaining AC and CC through experience or by means of flood risk management policies put in place by local authorities, such as early warning systems, flood risk awareness programmes or land-use planning regulations (Johnson et al. 2005; Penning-Rowsell et al. 2006). On the contrary, in the case of non-occurrence of frequent flooding thanks to the presence of flood protection structures, an increase of social vulnerability can be observed, i.e. the so-called levee effect (Montz and Tobin 2008). In fact, measures taken in order to prevent flooding may lead to a shift from frequent but small flooding, to rare but catastrophic ones (Burton and Cutter 2008; Di Baldassarre et al. 2013; Kates et al. 2006; Ludy and Kondolf 2012).

Notwithstanding the dynamic nature of vulnerability, methodologies for vulnerability and flood risk assessment available in the literature are usually static. In some cases, changes in risk are assessed by comparing scenarios driven by climate change and socio-economic development (Apel et al. 2009; Winsemius et al. 2013) or by flood policies (Di Baldassarre et al. 2015).

One of the main limitations to dynamic assessment of flood risk and vulnerability lays in the difficulties in obtaining time series of the spatial indicators used in the assessments. In recent years, earth observation technologies (EO) have contributed to overcome such limitations by providing accurate and prompt data for different purposes, including rapid assessment of impacts, support emergency management, DRR and CCA (De Sherbinin 2014; Giupponi and Biscaro 2015; Marconcini et al. 2013; Wolters and Kuenzer 2015).

The EO data revolution started recently, with the launch of the Copernicus programme by the European Commission in 2014, providing a systematic weekly observation of the Earth in multi-spectral optical and radar mode. Together with the USGS’s Landsat data, available from 1984 at 30 m resolution, these are constituting the so-called EO big data (Guo et al. 2014, 2015), an unprecedented tool for the observation of global environmental changes. While several spatial and non-spatial parameters required for detecting and quantifying flood risk can be extracted directly from RS imagery, many have to be elicited indirectly. This is especially true for VA, where the inherent characteristics of several socio-economic urban elements cannot be extracted directly from the images (Ghaffarian et al. 2018). Therefore, using proxies have become central and the predominant approach in RS-based risk assessment studies (Taubenböck et al. 2009).

Ghaffarian et al. (2018) provide an exhaustive review of EO-based proxies for disaster risk management, including vulnerability. Most of the proxies they found in the literature refer to the physical component of vulnerability, the one we call susceptibility. For instance, detecting buildings material allows to classify their level of vulnerability (Geiß et al. 2016; Yuan and Wang 2004) in case of a shock affecting their structure, or detecting the elevation of the buildings compared to the street level allows to understand their vulnerability to floods (Müller 2013). The analysis of road networks and road characteristics allows for VA in relation to flows of people, both in a frame of risk reduction and impacts assessment (Hu et al. 2017; Kumagai 2012). They found also several studies focusing on the socio-economic components of vulnerability, where for instance the income of population is inferred by analysing the slope where buildings are located or by the percentage of green in the urban environment (Ebert et al. 2009).

Further contributions to VA and its dynamics are coming from machine learning techniques, helping to extract information from different sources of data. For instance, Schwarz et al. (2018) used machine learning, remote sensing (RS) and census data to predict the socio-physical vulnerability to floods and dynamically deliver that information to decision makers through a web platform in the case of Senegal.

Given the current status of the literature briefly outlined above, this work proposes a novel method for calculating a spatial and multi-temporal index of social vulnerability to floods (flood vulnerability index, FVI), integrating EO, census data and spatial statistics about flood events. The method is flexible enough to easily allow decision makers and stakeholder to customize the way of integrating the different indicators by changing their overall weight within the framework, based on local knowledge and experience. Vulnerability is assessed through three components—AC, CC and susceptibility—by combining physical properties of urban areas (extent, growth, imperviousness, etc.) with economic and societal features (age of population, income classes, commercial activities, flood memory, etc.) by means of an Ordered Weighted Averaging method.

2 Methodology

Figure 1 shows the workflow followed to compute the FVI. From the input data, several indicators are extracted, which are later normalized by means of value functions and weighted using a simple weighting method (SW). Finally, an ordered weighted averaging (OWA) is applied to obtain the final FVI. The next two sections explain in detail all the steps.

Fig. 1
figure 1

Workflow with hierarchical aggregation of indicators and simple weighting for calculation of the flood vulnerability index. Indicators highlighted in blue are computed from census data, in green from earth observation data and in red from flood record

2.1 Flood vulnerability indicators

Fourteen indicators have been defined and derived from the input data. Table 1 shows their definition, their description and to which component of vulnerability they belong to.

Table 1 List of indicators used in the assessment of vulnerability

They have been selected from the generalized list of indicators proposed by Giupponi et al. (2013) and Mojtahed et al. (2013) for the three components of vulnerability. The selection has been based on data availability and on characteristic and spatial extent of the area of study.

The indicators have been stored in GIS layers with highest possible resolution: at census cell level for census data, municipality for flood records, 30 × 30 m grid for data deriving from RS. All vector data were later rasterized at 30 m resolution to perform a pixel-based analysis.

Indicators derived from census data have been computed for the years 1991, 2001 and 2011 corresponding to the years of the census surveys with the exception of the indicator “Conservation Status” not computed for 1991 and 2001 for the lack of data in the two censuses, and the indicator “Number of Multi-Storey building” not computed for 1991 for the lack of data in the census.

The three indicators derived from EO data have been computed from 1991, 1996, 2001, 2006, 2011 and 2016, with exception of the indicator “Newcomers” not computed for 1991 given that data for urban areas prior to 1991 were not available.

The “Memory Effect” indicator, which is based on the experience of the society to cope with a flood based on the occurrence of past events, was computed for each of the years of the analysis adding a dynamic component to the FVI. The indicator has been computed only for the Veneto region since the data on past flood were available only for this region. For a detailed description of the employed model, see Sect. 1 of Supplementary Materials (SM).

2.2 Indicators normalization, weighting and aggregation into the FVI

The extracted indicators are aggregated in two steps to define the social and physical dimensions of vulnerability, first weighted through a SW method, and later aggregated by means of OWA, a spatial multi-criteria analysis method (Yager 1988), giving as a result the final FVI. To allow the aggregation of the indicators, values expressed in different units were normalized by means of value functions expressing the contribution of each indicator to the three dimensions of vulnerability, thus obtaining values between 0 and 1 as indicated in Table 2.

Table 2 Value functions adopted for indicators requiring normalization

Figure 1 shows in detail how the different components of vulnerability are considered and weighted in the SW step. At the stage of methodological proposal, we considered that weights could be equally distributed within each group of indicators to be aggregated, while in operational applications they should be defined with competent decision makers and stakeholders and subject to specific sensitivity analysis (see Bojovic et al. 2018 or Ceccato et al. 2011, for the techniques of participatory weighting with stakeholders). The physical dimension, which weights one third, is the susceptibility and it is given by the properties of the buildings and the value of imperviousness (percent impervious surface) of the area. The social dimension is made of AC and CC (each weighted one third of the FVI). AC takes into consideration the skills of the society, such as employment and income level. CC takes into consideration demography, the urban growth and the characteristics of the urban environment.

After the SW, OWA is applied where the values obtained per each pixel are ordered and a second weight vector is applied, in which weights are not attributed to each indicator, but are instead applied to the ordered sequence of weighted values to be aggregated. This second weighting step allows to overcome the full compensation of simple weighted averages. Variations in the level of skew in the order weights results in solutions with different levels of risk. Therefore, the balancing between ANDness (risk averse) and ORness (risk taker), and the values of the weight vectors allow for much better consideration of the risk attitude and thus improve significantly the opportunity to meet the preferences of the decision/policy makers. Table S1 and Figure S2 show the weight order and value for the three decision makers’ profiles considered in this study: balanced, optimistic (risk taker) and pessimistic (risk averse).

Let \( \bar{\varvec{x}} \) be the vector containing all the n indicators; let \( \bar{\varvec{w}} \) be the vector containing all the n weights for each corresponding indicator, and let \( \bar{\varvec{w}} \) be the vector containing the 3 weights corresponding to the three components of vulnerability, AC, CC and S.

Then, the vector containing all the simple additive weighted indicators will be:

$$\overline{{\varvec{S}}{\varvec{W}}{\varvec{I}}}= \left(\begin{array}{c}\begin{array}{c}{SWI}_{1}\\ \vdots \\ {SWI}_{k+1}\end{array}\\ \vdots \\ {SWI}_{h+1}\\ \vdots \\ {SWI}_{n}\end{array}\right)= \left(\begin{array}{c}\begin{array}{c}{x}_{1}*{w}_{1}*{W}_{1}\\ \vdots \\ {x}_{k+1}*{w}_{k+1}*{W}_{2}\end{array}\\ \vdots \\ {x}_{h+1}*{w}_{h+1}*{W}_{3}\\ \vdots \\ {x}_{n}*{w}_{n}*{W}_{3}\end{array}\right)$$
(1)

with k being the number of indicators for the component CC, h the number of indicators for the component AC, and n the total number of indicators.

Let \( \bar{\varvec{h}} \) be the vector containing all the n OWA weights with \( h_{n} \in \left[ {{\text{0}},{\text{1}}} \right]{\text{and}}\sum _{{i = 1}}^{n} h_{i} = 1 \) and \(\bar{\varvec{b}} \) the ordered vector of the SWI indicators, where bi is the ith largest and bn is the smallest of the SWI indicators.

Then, the flood vulnerability index (FVI) is:

$$FVI\left({SWI}_{1},\dots , {SWI}_{n}\right)= \sum_{i=1}^{n}{h}_{i}{b}_{i}$$
(2)

Finally, the FVI values are masked using the built-up area layer of the year of interest. FVI is obtained with a score between 0 and 1, where 1 represents high vulnerability and 0 no vulnerability.

A further step is performed for the sake of presenting the results at a census level in a suitable form for a potential policy maker using this analysis. Specifically, FVI values have been weighted over a census cell (i) based on the percentage of built-up area on the total area of the cell:

$${FVI}_{wgt}=\mathrm{FVI }* \frac{{Builtup}_{i}}{{Area}_{i}}.$$
(3)

3 Data and case study

3.1 Census data

Census data for Italy come every ten years. We decided to consider the census data of 1991, 2001 and 2011 to have an overlap with the EO dataset, in particular the Landsat data available from 1984 until today (as presented in the next subsection). The data have been retrieved from the Italian National Institute of Statistics (ISTAT) (ISTAT, Dati Censimento).

These data provide a wide range of information about demography, education, employment, characteristic of buildings. Data are aggregated at census cell level, which are very detailed subdivisions of municipalities. Their dimensions can go from few hundreds square metres in case of densely populated city centres, to several tens of square kilometres for mountainous area with little population.

It has to be noted that census surveys changed slightly in the years, with new indicators added in the most recent census (as explained in Sect. 2.1) and with minor changes in the division of municipalities and census cells. This latter problem is overcome by the fact that our analysis is pixel-based and not area-based; therefore, the different boundaries of the cells in different years are not affecting our analysis.

Income of each family is also collected by the census and it would be a useful indicator of the capacity of each family to take adaptation actions. Unfortunately, these data are not freely accessible due to privacy reasons and therefore they could not be used. Nevertheless, given the level of wealth in the area analysed, this is not a limitation and it has been hypothesized that a similar discriminant in terms of AC can be the percentage of rented houses.

3.2 EO data

EO data have been used to derive three vulnerability indicators, as specified in detail in the Methodology section. The EO product (produced by the Smart Cities and Spatial Development team of the German Aerospace Centre (DLR)) on which these indicators are based are: (i) the settlements extent maps for year 1991–2016 at 5 years interval at 30 m resolution, (ii) the corresponding maps estimating the percent impervious surface (PIS) per pixel and (iii) spatial networks whose nodes are the centroids of the detected settlements and edges connect each node with all neighbour nodes located within a given distance buffer (used to ultimately estimate the relevance of difference settlements at regional level for each target year).

The settlement extent maps have been extracted from the World Settlement Footprint (WSF) evolution dataset (Marconcini et al. 2019), a novel global layer outlining—on a yearly basis—the settlement growth from 1985 and 2015, derived by exploiting temporal statistics of specific spectral indices computed from Landsat archived imagery. On top of these, the corresponding PIS has been calculated according to the method presented by Marconcini et al. (2015). Specifically, each pixel is ultimately associated with the estimated percentage of the corresponding surface at the ground covered by buildings (intended as structures having a roof supported by columns or walls and intended for the shelter, housing or enclosure of any individual, animal, process, equipment, goods or materials of any kind) or paved surfaces (intended as any level horizontal surface covered with paving material, i.e. asphalt, concrete, concrete pavers or bricks but excluding gravel, crushed rock and similar materials).

Finally, spatial networks have been produced starting from the abovementioned settlement extent maps following the approach described in Esch et al. (2014); in particular, these have been generated for multiple distance buffers, namely 5, 4, 3, 2 and 1 km. Different relevance measures are computed to characterize the importance of different nodes in the network; among these, in this study we considered the mean Euclidean distance among the given node and its 3 nearest neighbours as a proxy for the settlement compactness.

3.3 Memory of past floods

A database of historical information about floods and landslides occurred in Italy has been maintained by the National Research Council, known as AVI archive (Aree storicamente Vulnerate da calamità Idrogeologiche (areas affected by hydrogeologic vulnerability)). The archive contains information starting from 1917 until 2002 (Guzzetti and Tonelli 2004) reporting past flood events, the municipalities affected and the duration of the events. Unfortunately, the archive is not maintained anymore, and data are not up to date.

For a VA study, Roder et al. (2017) and Sofia et al. (2017) updated the archive for the Veneto region up to 2015 by mean of several sources of information, such as newspapers, civil protection reports, etc. From this updated archive, we used the list of flood events occurred in the Veneto region from 1950 until today, in particular the municipalities affected, the dates of the events and their duration, to model the memory effect, which affects the social component of vulnerability. In fact, once a community experiences an adverse event (which affected people and caused damages), it builds a memory of the event that contributes to increase its coping capacity for a future adverse event. The memory effect will have a higher contribution to the coping capacity right after an adverse event, which will decrease with time until the community has no more memory of any past event.

The model employed for deriving this indicator is explained in detail in the supplementary materials–“Sect. 1. Vulnerability Dynamics: Memory Effect”.

3.4 Case study

For our analysis, we selected the three administrative regions that belong to the macro-region of Northeast Italy, Veneto, Trentino-Alto Adige and Friuli-Venezia Giulia, an area of about 18,400 km2. The area is densely populated especially in the southern and south-eastern part, characterized by a flat and highly built-up land. The northern area is covered by the Prealps and Alps, making it barely populated with the exception of few towns and cities, such as Trento, Bolzano and Belluno.

The whole area has been affected by numerous flood events in the past. Many of these occurred in the last decade causing enormous impacts, both in terms of people affected and in terms of economic losses, making this area of high interest in terms of flood risk assessment and risk reduction, for which understanding vulnerability is fundamental (Cian et al. 2018a,b).

3.5 Analysed scenarios

The time steps analysed are six and correspond to the years 1991, 1996, 2001, 2006, 2011 and 2016. They can be divided in two main groups based on their different nature: the “census” and “inter-census” scenarios, depending if they fall on the year of a census survey or not.

The choice to analyse inter-census scenarios is based on the fact that RS data are available on a yearly basis and allows for an update on certain characteristics of the vulnerability (mainly in relation to the built environment) without waiting for a new census. To exploit this advantage, particularly important in urban areas affected by rapid and significant changes, it has been chosen to create scenarios in between two census surveys, therefore at year 1996, 2006 and 2016.

The census data have been kept constant (to the previous time step) in the “inter-census” scenarios.

4 Results

4.1 FVI values across space and time

The FVI has been computed for six different years (from 1991 to 2016 at 5-year time step) for three different decision marker’s profiles (balanced, risk taker, risk averse). The resulting maps are reported in Figures S3–S38. For each profile, we present the FVI and FVIwgt results for a total of 12 maps each.

Figure 3 shows the FVIwgt map for 2016 balanced profile on a colour scale from 0 to 0.7 (the max FVI value reached in the region), with 0 meaning no vulnerability and 1 full vulnerability, and the evolution of FVIwgt from 1991 to 2016 for the subarea comprising the cities of Venice, Treviso, Padova and Vicenza, one of the most populated zone of Triveneto. The white areas, indicating no vulnerability, correspond mainly to mountainous areas and water bodies, where there is no presence of population. We observe higher values of vulnerability in correspondence to the urban areas, especially in bigger cities in their historical and productive centres. We can also notice a sensitive decrease of FVI values in time, indicating that some of the changes occurred to society or to the physical environment had a positive effect on vulnerability. The introduction of the indicator of “Number of multi-storey buildings” from the 2001 census contributes to lower the vulnerability especially in the city centres where most of the buildings have more than one floor.

Nevertheless, it has been included in the assessment to improve the quality of the most recent maps. For operational implementations, sensitivity analysis could be applied to decide about its inclusion, or, for example, the use of 2001 data also for earlier dates.

Analysis of results obtained in specific areas helps understanding the significance of the results obtained (see the zooms in Fig. 3). For example, in certain areas of the historical centre of Padova, a much higher number of empty buildings and the increase in the memory of population, results in a consistent decrease of vulnerability. This should be a combined result of previous flooding (occurred in 2010 and 2012) and the effect of the economic crisis, pushing people to other areas, cheaper and safer. A second example can be taken from the San Giobbe neighbourhood of Venice (north-west of the city), and the improved capacity of society is driving the decrease of vulnerability, specifically by a better dependency ratio, an increase in employment and education level and the increase of owned property. The last factor could be taken as an example where the knowledge of the local stakeholder can be important in the tailoring of the FVI, an idea that drove our design of the index. In fact, it is known that in cities like Venice, with a high number of students and tourists, there are cases of illegal renting in properties where people instead declare to be resident. If this is known to be a significant phenomenon, the indicator could be dropped or assigned to a smaller weight. A last example can be found nearby the historical centre of Vicenza, to the southeast, where a significant decrease of FVI is due to the increase of memory (the city in 2010 has been severely hit by a flood (Viero et al. 2013)), an increase of empty building, a smaller percentage of children and/or elderly people, better education and more employment. The latter may also due to smaller percentage of illegal worker, known to be higher in the 90s in the region.

Table 3 shows the mean values of FVI for the whole Triveneto area from 1991 to 2016 for the three risk profiles. We can observe a continuous decrease of the FVI values for the three profiles with an almost identical magnitude, about 17% relative decrease from 1991 to 2016.

Table 3 Average values of the flood vulnerability index in the Triveneto region for the three different risk profiles analysed from 1991 to 2016

The inter-census values (i.e. 1996, 2006 and 2016) are slightly higher (2–4%) than their corresponding census values in case of 1996, practically identical in 2006 and slightly lower (−0.5 to −1.3%) in 2016. This tells us that the main driver of change is not due to the information derived from EO data, i.e. urban growth as proxy of newcomers, imperviousness and compactness.

Compared to 1991 values, in 1996 the risk taker profile sees an increase of 2% of mean FVI value compared to the 4% increase of the risk averse profile. This tells us that the indicators derived from EO performed relatively bad and therefore were taken into account with a greater magnitude by the risk averse profile, suggesting a bigger impact of the physical development of the areas (i.e. construction of new buildings).

With each new census (i.e. 2001 and 2011), we see a decrease of FVI of about 7–8%, suggesting an average trend of improvement in terms of flood vulnerability.

What the EO data is telling us is how the land-use change is affecting the overall FVI values. In fact, all the EO-derived indicators are connected to built-up areas in the region. The increase of FVI from 1991 to 1996 is suggesting that the growth in built-up areas (as can be noticed from the built-up growth data itself) are increasing the value of vulnerability, assuming a constant capacity of the society (given by census data that are static). Vice versa, the decrease that we observe from 2011 to 2016 is suggesting that the built-up growth has slowed down (contributing to a smaller extent to the value of vulnerability), with census-based indicators considered with a greater weight in the OWA.

The improvement coming from urban growth in the last years (i.e. a slowing down in urban growth and therefore the relative indicators contributing less to the final value of FVI) is one of the reasons of a higher decrease of the FVI values in the last census update (i.e. 2011). Social factors are also playing an important role, by improving the capacity of society in coping and adapting to floods.

4.2 FVI dynamics and territorial clusters

In order to derive synthetic information about the combinations of FVI values in the study area and their temporal variability, we used the ISODATA cluster analysis technique (iterative self-organising data analysis technique), which is a consolidated k-clustering method used for identifying land classes from stacks of multiple images in RS studies (Johnson and Wichern 2007; Richards 2012). The map of ISODATA clusters provides a synthesis of multivariate spatial variability of the evolution of FVI over time (six dates of FVI) and can be considered as a preliminary support for the identification of a series of different zones characterized by relative internal homogeneity, thus requiring specific approaches in terms of policies and measures for managing flood risk in the area. Cluster analysis has been applied to the multi-temporal stack of the six final FVI maps reported in Fig. 2.

Fig. 2
figure 2

Multi-temporal maps of flood vulnerability index

Cluster 1 identifies areas with no or very limited presence of settlements, i.e. the vast majority of the study area, where the FVI is very low because of limited presence of people, such as mountain areas (Table S2 provides the average values of the FVI in each cluster over time). Cluster 5 is similar to 1, and the FVI shows average values around 0.02 in 1991, slowly increasing up to 0.03 in 2015. In practice, this cluster identifies rural areas at the periphery of cities. Cluster 3 and 4 are similar, showing increasing FVI over time and the identify areas of development of new settlements in and around the historical urban settlements. Cluster 3 has values increasing over time from 0.06 to 0.19, while cluster 4 increases from 0.14 to 0.29. Cluster 6 identifies spots where the FVI has decreased over time. Figure 3 presents a zoom of the ISODATA map in the central area of the Veneto region, including the cities of Venezia, Treviso and Padova.

Fig. 3
figure 3

Clusters of flood vulnerability index, resulting from ISODATA analysis of the six multi-temporal maps

4.3 Analysis of sensitivity of FVI under the effect of aggregation

In order to explore the sensitivity of the results of the FVI to the OWA algorithm and in particular to the ordered weight vector, three parallel calculations of the AC, CC, susceptibility and FVIs were performed per each date. The results discussed so far were obtained for the balanced profile, giving relative higher weight in the final aggregation to both indicators with good performances (higher values obtained with SAW weights) and bad ones (lower values of weighted indicators). That produced a U-shaped set of weights with highest weights being approximately 40% higher than the lowest ones, with gradual changes between the two extremes, applied on a pixel per pixel basis (see Figure S2). Two other cases were considered: the case of a risk taker (or optimistic) decision maker who is satisfied with those results in which at least a subset of indicators have good performances; the case instead of a decision maker who is risk averse (or pessimistic) and thus weights more the indicators with bad performances. In the first case, the OWA weights were defined with a linear increase from bad indicators to those with best performances, thus weighting the best performing indicator approximately as much as twice the worst. The opposite was applied in the case of risk averse attitude.

For each date, “difference maps” were calculated (risk averse minus risk taker map values) to explore the effects of changing attitudes to risk on vulnerability mapping. The case of 2011 FVI is selected here to exemplify the sensitivity of results: the maximum difference value of FVI is 0.17, with an average difference of 0.0028 (the first map having an average FVI value of 0.034, while the second having 0.031), and a sample standard deviation of 0.01). In operational implementation of the proposed approach, these sensitiveness will be of great interest for decision makers, because they allow for the identification of the areas in which the inherent subjectivity of risk attitudes may generate higher variability of results. Risk attitude is only one of the sources of uncertainty and such possibility of efficient sensitivity analyses would become very important for decision makers when the proposed methodology would be adopted for flood risk management planning and zoning.

5 Conclusions

A method for assessing multi-temporal flood vulnerability to people has been presented and demonstrated in Northeast Italy. One of the novelties of the method is the proposed multi-temporal combination of both census and EO data, which contributes to the understanding of the dynamic evolution of vulnerability over time and makes it possible to go beyond the limited time resolution of census surveys, which for Italy comes every ten years. EO data allows to create intermediate time steps and thus contribute to explore and understand, the evolution of the vulnerability.

The proposed method and the obtained results can be an important instrument for local stakeholders and policy makers. In particular, the proposed vulnerability index could be used to analyse the effect of societal evolution and past urban planning and/or DRR, measured through the multi-temporal sequence of vulnerability maps. The vulnerability maps should then be combined with coherent spatial descriptions of hazard and exposure to allow for comprehensive assessments of flood risk over space and time. Very importantly, the method could be easily adapted for analysing not only the past, but also future scenarios and thus contribute to assess the expected effectiveness of future policies or measures under consideration by policy makers. Moreover, desired targets in the social and/or physical domain could be defined, allowing the stakeholder to understand where it is more effective to invest in order to reduce the vulnerability of people and subsequently the risks to which they are exposed.

The current implementation of the proposed approach has an academic feature. In order to move towards a consolidated FVI mapping for planning purposes, a new substantial component should be added to the procedure, i.e. the involvement of relevant stakeholders and decision makers. Only the implementation of a sound approach for the identification and involvement of stakeholders would allow for the consideration of their preferences in terms of indicator weighting and risk attitude. As a result of the participatory process, one should expect a multitude of parallel runs, producing subjective declinations of the same algorithm. The proposed approach, being fully implemented in a transparent and reproducible GIS macro, would thus allow for efficient treatment of different sources of uncertainty through the parallel execution of multitude of runs. Sensitivity analysis supported by statistical and data mining techniques would eventually allow to understand the possibility of adopting average results and their robustness, identifying also crucial elements of sensitivity, in terms of critical indicators and weights, and those areas in which the subjectivity of judgements may affect the orientation of planning instruments.