Natural Hazards and Earth System Sciences Local scale multiple quantitative risk assessment and uncertainty evaluation in a densely urbanised area ( Brescia , Italy )

The study of the interactions between natural and anthropogenic risks is necessary for quantitative risk assessment in areas affected by active natural processes, high population density and strong economic activities. We present a multiple quantitative risk assessment on a 420 km2 high risk area (Brescia and surroundings, Lombardy, Northern Italy), for flood, seismic and industrial accident scenarios. Expected economic annual losses are quantified for each scenario and annual exceedance probabilityloss curves are calculated. Uncertainty on the input variables is propagated by means of three different methodologies: Monte-Carlo-Simulation, First Order Second Moment, and point estimate. Expected losses calculated by means of the three approaches show similar values for the whole study area, about 64 000 000C for earthquakes, about 10 000 000 C for floods, and about 3000C for industrial accidents. Locally, expected losses assume quite different values if calculated with the three different approaches, with differences up to 19 %. The uncertainties on the expected losses and their propagation, performed with the three methods, are compared and discussed in the paper. In some cases, uncertainty reaches significant values (up to almost 50 % of the expected loss). This underlines the necessity of including uncertainty in quantitative risk assessment, especially when it is used as a support for territorial planning and decision making. The method is developed thinking at a possible application at a regional-national scale, on the basis of data available in Italy over the national territory.


Introduction
The management and the reduction of natural and technological risks are strategic problems for disaster prone communities in the last century.The frequency and the consequences have increased worldwide (Munich Re Group, 2004).For this reason, the reduction of risk is recognized as an integral component of both emergency management and sustainable development, involving social, economic, political, and legal issues (Durham, 2003).
A necessary condition for risk prevention, mitigation and reduction is its analysis, quantification, and assessment.Quantitative risk assessment (QRA) is demonstrated to be effective in the resolution and mitigation of critical situations, and it poses some interesting challenges to the scientific communities.As a consequence, different QRA methodologies have been proposed in the literature for different natural and technological risks, based on different approaches for the calculation of societal, individual and economic risk (e.g.Kaplan and Garrick, 1981;CCPS, 1989;Jonkman et al., 2003;Bell and Glade, 2004;Calvo and Savi, 2008;Agliardi et al., 2009;Fuchs and Brundl, 2005).
In some contexts, e.g.where the presence of sources of natural hazards is combined with a strong industrial and urban development, the integration of different risks can help in the design of mitigation measures (Lari et al., 2009).In this context, multiple QRA allows the comparison and integration of risks deriving from different threats.
Many approaches have been proposed to assess specific natural and technological risk, but only few studies combine multiple sources of hazard to obtain an overall analysis.Some transnational studies have been performed to compare risk levels in different countries (e.g.Cardona et al., Published by Copernicus Publications on behalf of the European Geosciences Union. 2004; UNDP, 2004;ESPON, 2005) or to identify key "hotspots" where the risks of natural disasters are particularly high (Dilley et al., 2005).However, these studies generally define risk in a qualitative way, based on relatively simplified approaches that make use of national-level indicators, without a detailed spatial analysis of hazard and element at risk patterns.Only a few local scale multi-risk analyses have been proposed including multiple sources of natural (Granger and Hayne, 2001;Middelmann and Granger, 2001;Van Westen et al., 2002;Kappes et al., 2012b) and natural/technological hazards (Barbat and Cardona, 2003;Ferrier and Haque, 2003;Lari, 2009).
The analysis and the management of technological and natural risks are affected by large uncertainties (Patè-Cornell, 1995), reflecting incomplete knowledge (epistemic uncertainty) or intrinsic randomness of the processes (aleatory uncertainty).QRA is also used to model rare events, where experimental verification of its validity is missing, and uncertainties can be significant.For these reasons, it is important that uncertainties in the results of the QRA are correctly characterized and interpreted (Parry, 1996), in order to support meaningful decision making.
In this paper, we present a local-scale heuristic quantitative risk assessment for a set of flood, earthquake and industrial accident scenarios in a high risk area (Brescia and lower Val Trompia, Lombardy, Northern Italy).We analysed risk scenarios derived from existing regulatory maps.These scenarios are easily available for the entire country, but they can be affected by large uncertainties because they are produced with a simplified approach at small scale.
In order to assess uncertainty, we applied and compared Monte-Carlo-Simulation, (MC), First Order Second Moment, (FOSM), and point estimate, (PE) techniques.
The aim of the analysis is to compare the values and distributions of different risks coexisting in the area, and to quantify the uncertainty introduced in the analysis by testing different methodologies.The identification of the most relevant risk in different zones is the basis for optimal allocation of resources in risk mitigation strategies, and for a more detailed investigation at the local scale.At the same time the proposed method is developed thinking at a possible application at a regional or national scale, by direct use of available datasets.The analysis does not consider domino effects, where some of the components of risk assessment (i.e.probability of occurrence and vulnerability) are not resulting from the simple sum of their values for single threats (Kappes et al., 2012a).

Study area
The study area includes the city of Brescia, its neighbourhoods, and lower Val Trompia (Lombardy, Northern Italy).It covers an area of 420 km 2 (Fig. 1), including a plain zone (90 m a.s.l.), and a pre-alpine zone, with a maximum elevation of 1360 m a.s.l.The population (722 000 people), is Fig. 1.Risk hot spot area (in black) as obtained from multi-risk analysis at the regional scale (Lari et al., 2009).Flooding zones as defined by the basin-scale regulatory map zonation (PAI -Piano stralcio di Assetto Idrogeologico, 2007).Three flooding areas are identified: zone A is defined as the area which contains 80 % of the discharge for a 200-yr flood; zone B is defined as the flooding area for a 200 yr flood; zone C is defined as the flooding area for a 500-yr flood.
distributed in 36 municipalities, with a maximum density of 2068 inhabitants km −2 in the urban area of Brescia and along Val Trompia.
The area is economically strongly developed and highly industrialised since the beginning of the 20th century: iron and steel, mechanic, chemical and foundry industries are diffused.Zootechnical and agricultural activities are relevant in the southern part of the area.
The study area was selected since classified as a risk hot spot in a regional risk assessment study performed for Lombardy Region (Lari et al., 2009), by integrating natural (rockfalls, shallow landslides, debris flows, floods on alluvial fans, deep-seated landslides, floods, earthquakes, wildfires) and technological sources of hazard (car crashes, industrial accidents, work injuries).Hot spot risk areas were detected according to the number of coexisting natural and technologic threats, and to their level.

Context for risk analysis, data availability and scenarios description
The major threats menacing the area, as resulting from the regional scale analysis (Lari et al., 2009), are floods, earthquakes and industrial accidents.Floods occur along the floodplain of the Mella river (basin, 311 km 2 in size), which flows through Val Trompia and the Brescia urban center.Six flooding events in the last century (1928,1959,1960,1966,1993; AVI-Aree Vulnerate Italiane da frane ed inondazioni, 2007), caused the destruction of bridges, roads and buildings in the area.The basin-scale regulatory map zonation (PAI -Piano stralcio di Assetto Idrogeologico, 2007) identifies three flooding areas (A to C). Zone A is defined as the area that contains 80 % of the discharge for a 200-yr flood.This area is flooded very frequently, and corresponds to the river bed.In the study area, zone A is poorly relevant for risk analysis since only few structures are located in it.Zone B is defined as the flooding area for a 200-yr flood.In the study area, this zone is missing due to hydraulic works which efficiently contain the 200-yr river discharge within the A zone.Zone C, defined as the flooding area for a 500-yr flood, embraces large parts of the plain area, with a potential impact on a large number of residential and industrial buildings, and crops (Fig. 1).In the area, only flooding zones A and C are included, with a return period of respectively 200 and 500 yr.
Seismic risk is quite relevant in the area of Brescia, being located 50 km from a system of active faults along Garda lake.According to the Italian Seismic hazard map (INGV -National Institute of Geophysics and Vulcanology -MPS, 2004), the peak ground acceleration (PGA) with exceedance probability of 10 % in 50 yr ranges from 0.125 to 0.15 g (Fig. 2).This value refers to reference site conditions with stiff soils characterised by shear wave velocities higher than 800 m s −1 .According to the PGA values, the area is classified in Italian law by ministerial decree OPCM 3274 20/3/03 as one of moderate risk (class 3d).Six events with magnitude larger than 5 were registered in a radius of 50 km in the last century (1901,1918,1932,1951,1979,2004), and 27 events since 1065, with a maximum historical magnitude of 5.7 Richter in 1901 AD.For the assessment of seismic risk, we considered three scenarios with different return periods (75, 475 and 2500 yr) corresponding respectively to an exceeding probability of 50 % in 50 yr, of 10 % in 50 yr and of 2 % in 50 yr (Meletti and Montaldo, 2006), as defined in the National Institute of Geophysics and Vulcanology (INGV).
Industrial risk in the study area is due to the presence of hundreds of industrial plants, mainly related to manufacturing.We focused our analysis on eight of them, classified as major risk plants by legislative decree D.Lgs.238/05, according to Council Directive 96/82/EC on the control of major-accident hazards involving dangerous substances, also known as the Seveso II Directive.For these plants, safety plans are available for different accident scenarios.The plans include the probability of occurrence, the area of impact, the intensity of the event and the effects on structures and on human life.Due to incomplete documentation about the productive processes and the accident scenarios, only explosionrelated accidents were considered in this work, neglecting the release of toxic gases and pollutants.Overall, 8 accident scenarios were taken into account, related to 3 industrial plants.For each of the scenarios, we considered only damages related to impacts exceeding the plant perimeter.Actually, damages within plants, which can also reach high values, are covered by the companies themselves, and do not have any impact on public funds.
The choice of risk scenarios to be considered for the analysis was based mainly on the availability of data, having the perspective of the applicability of the approach on different study areas.The use of scenarios defined by national law and zonations makes the approach potentially applicable on the whole national territory, providing an homogeneous methodology for multiple QRA.The number of scenarios for each threat, and the related return time are different.This problem is overcome by the use of curves which express the expected losses as a function of exceedance probability, deriving from F − N curves (Ale et al., 1996;Vrijling and Van Gelder, 1997;Jonkman et al., 2003).These curves allow to integrate losses related to all the considered scenarios for each threat, and to provide a total value of risk as the subtended area (Ale et al., 1996;Jonkman et al., 2002).

Methodology for risk assessment
Natural and technological risks are generally agreed to be defined as the measure of the probability and severity of an adverse effect to life, health, property, environment, and can be expressed as the probability of an adverse event times the consequences if the event occurs (ISSMGE Glossary of Risk Assessment Terms, http://www.engmath.dal.ca/tc32/) or the probability that a given loss will occur (e.g.Kaplan and Garrick, 1981;Baecher and Christian, 2003;CEDIM, 2005).Kaplan and Garrick (1981) define risk by a combination of the expected consequences of a set of scenarios, each with a probability and a consequence.Therefore, risk encompasses three aspects: hazard, vulnerability of the affected element and the asset of exposed elements at risk (Kleist et al., 2006).
We perform a local-scale multiple quantitative probabilistic risk assessment, aimed at calculating the consequences of adverse events, in terms of expected economic annual loss for a set of scenarios (2 for floods, 3 for earthquakes and 8 for industrial accidents) with different probabilities and intensities.
The aggregation of annual expected losses, E(L), obtained for each of the considered scenarios, is performed with an approach analogous to F −N curves, by representing the annual probability of exceedance of a certain level of loss, x 0 , as a function of the economic loss.The expected value of E(L) can be derived from the probability density function (pdf) of the economic loss, f L (x) (Ale et al., 1996;Vrijling and Van Gelder, 1997;Jonkman et al., 2003): (1) Ale et al. (1996) propose the area under the F − N curve as a measure for societal risk, analogous to the F − N curve and the expected number of fatalities.It can be shown that the area below the curves representing the annual probability of exceedance of a certain level of loss as a function of the economic loss, equals the expected value E(L) (Jonkman et al., 2002).The curves obtained combining the scenarios for each threat allow the quantification and the comparison of the economic impact of the threats in an area.
The census parcel (2460 parcels) with an area ranging from 0.0001 to 19 km 2 has been adopted as terrain unit.The choice was driven by the fact that detailed socio-economic data, provided by the National Institute of Statistics (ISTAT, 2001), are based on this territorial unit.
In this paper, we apply a methodology based on Fell et al. (2005), in which the frequency of the event, its probability of reaching the element at risk, the temporal spatial probability of the element at risk (i.e.exposure), its vulnerability and value are considered.
For each scenario we evaluate the expected annual loss E(L) as where P (event) is the annual exceedance probability of the event of a given intensity I , calculated with a Poissonian model of the type: where t is the period on which the probability is calculated (here, 1 yr); µ is recurrence interval of the event.P (event) is also referred as Hazard.
The potential loss given the event, E(L|event), is calculated as where P (I |event) is the probability that an event, once occurred, could impact a parcel; P (L|I ) is the vulnerability of a parcel, here intended as the physical vulnerability, i.e. the degree of loss to a given element or set of elements within the area affected by a hazard, as defined by ISSMGE Glossary of Risk Assessment Terms; W is the total economic value of a parcel, dependent on the number and typology of buildings.
The probability of impact given the event, P (I |event) is calculated in GIS environment by means of a geometric analysis as the portion of built area of each parcel potentially impacted by an event (Fig. 3).This probability expresses, for each parcel, the exposure to the specific threat.use map (DUSAF, 2007) was used for mapping the built areas.P (I |event) is always equal to1 for seismic risk.
In the literature, a few risk assessment studies provide some examples or the quantitative estimation of values of the exposed elements (e.g.MURL, 2000;IKSR, 2001;Dutta et al., 2003;Grunthal et al., 2006).In some cases, the value of residential buildings was estimated considering the mean insurance value for the buildings, which in general represent the replacement costs of the buildings (MURL, 2000;Grunthal et al., 2006).Dutta (2003) used economic replacement values for structures and other land uses, distinguishing buildings according to their use.This approach, based also on the use of aerial photographs, is adequate only for limited areas.Blong (2003) used construction costs published by Australian authorities.
In this paper, the value W of each census parcel has been calculated by using minimum and maximum market values (C m −2 ) distinguishing among residential and nonresidential buildings, according to data provided by Italian Observatory of Real Estate Market (Osservatorio Mercato Immobiliare, 2010) (Fig. 4).

Methodology for uncertainty evaluation
Monte Carlo analysis (MC), or its variants such as the Latin hypercube sampling, is the most common technique for uncertainty evaluation, and it has been applied to many different fields (e.g.Fishman, 1997;Breuer et al., 2006;Hall, 2006;Calvo and Savi, 2008).In this approach, the uncertainties on parameters are generally described as probability distributions (e.g.normal, log-normal, uniform, triangular, discrete), which could be continuous or discrete.Probability distributions provide the range of values that the variable could take, and the likelihood of occurrence of each value within the range.Monte Carlo methods are a class of algorithms based on repeated random sampling, and include a wide range of engineering and scientific simulation methods based on randomized input.Monte Carlo methods are used when it is infeasible or impossible to compute an exact result with a deterministic algorithm.Probability distributions are a suitable way of describing uncertainty in variables of a concerning risk analysis.
The First Order Second Moment (FOSM) method consists in the truncation of Taylor's series expansion, considering only the first order terms.It provides an estimate of the mean and the variance (first and second moment) of the output through computation of its derivative to the input at a single point (e.g. in Yen et al., 1986;Baecher and Christian, 2003;Uzielli et al., 2006).The discarded terms are functions of the second and higher order derivatives, the variances and shapes of the probability density functions of the input variables, and the correlations among input variables (El-Ramly et al., 2003).The method requires the evaluation of partial derivatives.In complex problems, this could not be possible, or too complicated and time costly.
In point estimate (PE) method (Rosenblueth, 1981), random variables are replaced with point estimates (e.g.Harr, 1989;He and Sallfors, 1994;Christian and Baecher, 1999;Ellingwood, 2001).Each variable is replaced with a central value (expectation, median, or mode), or with one which is consciously biased, and the estimates are deterministic.The aim of PE is to calculate the first two or three moments of a random variable, which is a function of one or more random variables.In decision making, a calculation of the first moment (expected value) of functions of the random variables would often suffice (Rosenblueth, 1981).The method overcomes the deficiencies of a deterministic treatment, sacrificing the accuracy of a rigorous probabilistic analysis (Rosenblueth, 1981).PE reduces the computational effort of propagating uncertainty through a function, by eliminating the calculation of derivatives or use of Monte Carlo sampling.The pdf of each random variable is just represented by discrete points, located according to the first, second and eventually third moments.
In general, the uncertainty evaluation performed by means of MC is computationally more demanding, and time consuming.The advantage of the method consists in the fact that the output function (i.e. the risk) is defined in its distribution, and all the moments can be calculated.The propagation of uncertainty is more accurate in Monte-Carlo-Simulation, as a variety of distributions can be chosen for the input variables, whereas this is not possible in the deterministic approaches.On the other side, FOSM and PE are simpler and faster.FOSM, however, requires the calculation of many derivatives, which in complex problems with many input variables can become problematic.In that case, a PE analysis should be preferred (Harr, 1989).
Many authors (e.g.Hoffman and Hammonds, 1994;Paté-Cornell, 1996;Parry, 1996;Faber, 2005;Merz andThieken, 2005, 2009;DerKiureghian and Ditlevsen, 2009) underline the necessity of keeping the distinction between aleatory and epistemic uncertainty, since the two types assume different meanings in risk assessment.The aleatory uncertainty is a fundamental and integral part of the structure of the QRA, and it is related to the intrinsic uncertainty in the occurrence and in the effects of phenomena.It is an uncertainty related to the variables on which risk is dependent, and it is not eliminable.Epistemic uncertainty, instead, is related to the level of knowledge and understanding we have of the scenarios we are modelling, and it can be reduced.
Other analysts (Hora, 1996;Hofer, 2001) have found that sometimes it is difficult to distinguish between the two types of uncertainty, especially when modelling the occurrence or the impacts of extreme physical phenomena, which can be outside our direct experience.In this work, epistemic and aleatory uncertainties are analysed jointly, being strongly connected at this scale of analysis.The use of regulatory maps and studies for risk assessment hardly allows to distinguish among the uncertainty deriving from lack of data related to processes or territory assets due to the scale (epistemic uncertainty), and the one related to intrinsic variability of the phenomena on such wide scenarios (aleatory uncertainty).Hence, the uncertainties expressed for the single variables and for the obtained risk values represent an overall unpredictability, lack of data and knowledge about the variables.Only in some cases, where possible, a distinction is maintained just in the phase of assessing uncertainty to the independent variables, as described below.
A relevant level of uncertainty was perceived in some of the steps of risk assessment.For the hazard sources, we assigned uncertainties to both the annual probability and the intensity of the event (e.g.depth of water for floods, peak ground acceleration for earthquakes).
For the elements at risk, we assigned uncertainties to the number of buildings of different typologies within each parcel, the value for unit area (C m −2 ), and the vulnerability.Data related to number of buildings of different typologies refer to 2001 (ISTAT, 2001).In more than 10 yr new areas have been urbanised, and some buildings can have changed morphology and structural characteristics.This introduces an uncertainty that is not negligible.The economic value of buildings with different use destinations (Osservatorio Mercato Immobiliare, OMI, 2010) is susceptible to market and local fluctuations.The use of a fixed mean value C m −2 , for each census parcel has been associated to a measure of uncertainty calculated on the basis of the range of values provided for each census parcel by OMI (2010).
The uncertainty for these variables is represented by coefficients of variation (COV, the ratio between the standard deviation and the expected value).A high coefficient of variation means a higher dispersion around the expected value and, thus, a larger uncertainty.
The assignment of COV values was performed according to expert knowledge.No values found in the literature (e.g.; Harr, 1996;Beck et al., 2002;Merz et al., 2004;Kaynia et al., 2008) were adaptable to the specific case study, since uncertainty is strictly related to the characteristic of each phenomenon and to the quality of the available data for the present case.In this work, both these aspects were taken into account for assignment by expert knowledge of COV values.
In absence of more specific information, variables used for uncertainty analysis are supposed to be mutually independent.This assumption, is considered realistic given the Nat.Hazards Earth Syst.Sci., 12, 3387-3406, 2012 www.nat-hazards-earth-syst-sci.net/12/3387/2012/nature of the variables (e.g.number of buildings, their economic value, probability of impact of a census parcel), except for P (E) and the intensity of the scenario.However, the choice of independence considerably simplifies the analysis, especially for FOSM and PE.On the other hand, the fact that P (E) and the intensity are inversely correlated, leads to a conservative approach, where a small risk attenuation introduced by the correlation is neglected.
In MC simulation, in absence of evidences leading to other assumptions, Gaussian distributions have been assigned to all the input uncertain variables, assuming that this type of distribution is generally the most probable.The only exception was for the number of buildings, for which a Pareto distribution was observed.Symmetric distributions have been always used in FOSM and PE.
MC simulation was performed with Latin hypercube sampling, and 5000 iterations.In Latin hypercube sampling, the input probability distributions are stratified, which means that the cumulative curve is divided into equal intervals on the cumulative probability scale (0 to 1).From each interval, a sample is then randomly taken.In this way, sampling is forced to represent values homogeneously in each interval over the entire range of the distribution.
For the application of FOSM, assuming that the variables are mutually independent, the calculation of output mean and variance is significantly reduced to a simpler form, since the correlation is neglected.Given a function of n variables: if θ is a sum of mutually independent variables x i , the first order approximation becomes exact, leading to and For the product of random variables, the approximation leads to where µ θ , and For the product of random variables, the approximation leads to: ) where µ θ , Ϭ θ and COV θ are the mean, the variance and the coefficient of variation of the output function, µ Xi, Ϭ Xn and COV i are the mean, the variance and the coefficient of variation of the input variable x i .
In the propagation of uncertainty by means of PE, the Roesenblueth's procedure (1981)  to the expected value of each variable x , increased or reduced of one standard deviation 13 Xn and COV i are the mean, the variance and the coefficient of variation of the input variable x i .
In the propagation of uncertainty by means of PE, the Roesenblueth's procedure (1981) for non-skewed independent variables was used.The procedure calculates the first and second moments of each function θ = f (x 1 , x 2 , ...x n ) using 2 n points with coordinates corresponding to the expected value of each variable x i , increased or reduced by one standard deviation (Christian and Baecher, 1999).
In this work, the number of buildings of different types (n. of floors, building materials) for each census parcel (IS-TAT, 2001), dates back to 2001.Considering that in the last decade a portion of about 12 % of agricultural surface was built up (Pileri and Maggi, 2010), the number of buildings is significantly underestimated.In order to characterise the frequency distribution of the urbanisation rate in the census parcels, 100 parcels were randomly sampled and analysed into detail through comparison of 1998 and 2007 aerial photographs.This allowed us to recognise a power-law frequency distribution of the urbanisation rate, with an average value of 11.3 %.This increase of urbanised areas can occur only for those parcels which are not completely built yet.For these parcels, a Pareto pdf is used in Monte-Carlo-Simulation in order to describe the power-law behaviour.For the other methods (FOSM and PE), we used a symmetrical distribution, with a COV of 0.1.For the parcels that are completely built, the uncertainty was assumed to be null.
In some cases, the assignment of COV values was based on the available distribution of data.In other cases, it was performed based on expert knowledge.This introduces a subjective judgement, which could be reduced in presence of a larger amount of data, but never eliminated.However, this does not hamper the repeatability of the analysis once COV values are assigned and explicitly declared.Moreover, the methodology here proposed makes use of regulatory maps and data, which are available for the whole national territory: this suggests that homogeneous COV values could be used for the whole nation.Obviously, they could be discussed by the collaboration of experts of different fields.

Flood risk assessment
To calculate the intensity (i.e.water depth) for flood scenarios, we reconstructed the flood surface by using a 20 m × 20 m digital terrain model and assuming the flood as a plane intersecting the outer border of the hazard zone (A, B, C).For each parcel, the mean water depth value was used in the computations.Data related to flow velocity, bed shear stress, dynamic forces, rate of flood rise were not available.
The damage to buildings and to their content was calculated by using depth damage vulnerability models for different typologies of building.Among the available data related to building structure, just the presence of one or more floors was considered as significant in order to assess flooding  expected loss.The building material (e.g.masonry or reinforced concrete) was not considered to influence the vulnerability, because a collapse scenario is rather unlikely to happen.
Among the models proposed in the literature (CH2M Hill, 1974; Sangrey et al., 1975;Smith, 1994;Torterotot et al., 1992;Hubert et al., 1996;Blong, 1998;ICPR, 2001;Dutta et al., 2003;USACE, 2003;Penning-Rowsell et al., 2005;Büchele et al., 2006;Kreibich and Thieken, 2008), we selected three models (ICPR, 2001;Dutta et al., 2003;USACE, 2003) which account for both damage to buildings and to their content.USACE (2003) also distinguishes, among the others, curves for one-storey and two or more storey residential buildings with basements, which is the typology of building diffused in the study area.In order to quantify the uncertainty related to the choice of the vulnerability model, we performed the risk analysis for each set of vulnerability curves.The damages related to buildings content are modeled in USACE (2003) as a percentage of structure value.In the use of the other vulnerability models, the value of content was considered equal to the one of the structure.This is considered a balanced choice, introducing an approximation by excess in case of residential buildings, where content value is in most cases lower, and an approximation by defect in case of non-residential buildings (which include productive plants, sanitary structures, and technology facilities), where content value is generally high.An example of the compared use of vulnerability models for a 1-storey residential building can be found in Table 1.
For the assessment of flood risk, the uncertainty assigned to each random variable is reported in Table 2.For economic value of residential and non-residential buildings, COV were calculated from OMI (2010) data, which provides a range of market values expressed as C m −2 for each census parcel.For the other variables, in absence of specific data, expert knowledge criteria was used.To the water depth, a COV equal to 0.3 was assigned.Actually, water depth was calculated by water surface elevation obtained starting from regulatory mapping of the scenarios, realised at the regional scale.The precision of data related to number of buildings of different type was higher (COV = 0.1), but not free from changing for new urbanisation, as explained before.The probability of event for each scenario was assumed to have an uncertainty of 0.1, being determined by National Basin Authority on the basis of a wide observation of past events.
The expected losses obtained by means of the three uncertainty methods are quite similar, although they reach their highest values in MC, and lowest in PE (Table 3).200-yr scenario, expected losses are limited because only a few buildings are potentially impacted.For both scenarios the highest level of expected losses is observed in small parcels close to the river, where the expected water depth is higher (Fig. 5).On the study area, the overall expected loss calculated by the integration of the scenarios is almost 10 000 000 C. Regarding the uncertainty, COV associated with expected loss ranges between 10 and 27 %, with small differences between the two scenarios.For all the three uncertainty methods, small parcels with high flood intensity show lower COV values (Figs. 6 and 7).

Seismic risk assessment
Seismic hazard has been defined in terms of peak ground acceleration, PGA (MPS, 2004), corrected in PGA c by a soil amplification factor, S s , for the horizontal acceleration spectrum based on soil categories (Table 4) defined according to velocity of shear waves in the first 30 m of depth, V s , 30 (DM 14 gennaio 2008).
To account for the effects of earthquakes on different building types, data related to the number of buildings for each parcel, built with different materials (masonry, 52 % of the total, and reinforced concrete, 48 % of the total), were considered (ISTAT, 2001).Vulnerability of buildings was calculated using four sets of fragility curves as a function of PGA (Kostov et al., 2004;Hancilar et al., 2006;Rota et al., 2008;Kappos et al., 2006).As for floods, risk analysis was performed for each set of curves in order to compare the results (Table 5).Kappos et al. (2006) propose curves for concrete frames with unreinforced masonry walls, mid rise, and low level of seismic design, here used to model fragility of masonry buildings, and reinforced concrete systems, low level of seismic design, here used for reinforced-concrete.Among the large number of curves proposed by Kostov et al. (2004) for Bulgary, the most similar typologies are masonry buildings constructed before 1945 and masonry buildings with reinforced concrete floors constructed after 1945.Hancilar et al. (2006), provide fragility curves only for buildings with reinforced concrete frame and shear walls, in Turkey.In this case, the damage related to masonry buildings has been underestimated (Table 3).Rota et al. (2008) propose for Italy susceptibility curves for different building typologies.Among them, 1-3-storey reinforced concrete buildings without a seismic design, and multiple storey masonry buildings have been selected.The percentages of damage for the structure associated to each damage state were taken from Kappos et al. (2006).
Uncertainties for the variables used in seismic risk assessment are reported in    3 for description of the categories.
-an epistemic uncertainty due to the estimation of amplification, and to approximations introduced in the definition of the soil categories, in the order of 0.2 (Table 2).
An epistemic uncertainty arises also from the averaging of peak ground acceleration on the whole census parcel area, which resulted to be negligible at the scale of analysis.COV values for the number of buildings of different types, and for the economic value of buildings, were calculated as explained for flood risk assessment.The probability of event presents a low COV value (0.001), being derived from a very detailed hazard analysis performed on the Italian territory (MPS, 2007).
Seismic risk affects the whole study area (Fig. 9), showing a predictable correlation with the level of urbanisation and productive activities.For all the scenarios, the highest losses occur in parcels with soil categories C (Medium density coarse soil or medium stiff fine soil, depth > 30 m,    4), in presence of a high number of masonry buildings.The highest level of expected annual loss is obtained for the 75-yr scenario (Table 6).For 475 and 2500 yr scenarios, the larger potential damages, due to the higher earthquake intensity, are compensated by a longer return period, resulting in an overall lower annual loss.On the whole study area, the three uncertainty methods provide similar results (Table 6).
Figure 10 represents the exceedance probability as a function of expected losses for seismic risk obtained by means of MC simulation, both for single parcels, (grey lines are some random examples, all ranging between minimum and maximum, i.e. in the grey area) and the whole study area (black line).Total expected loss, considering all seismic scenarios on the study area, was derived as the area under the curve.
On average, the uncertainty propagated by the different approaches is about 40 % of the expected value, and does not show any dependence on risk value (see upper quadrants in Fig. 11, where COV values are shown with respect to the level of risk, since census parcels are ordered).COV values resulting from FOSM are higher and more scattered, while MC provides lower and more constant values (Fig. 11).

Industrial risk assessment
A high value of COV (0.4) was assigned to the loss to buildings provided in safety plans of the industrial plants.Actually, the possible effects of the scenarios are described in a qualitative way.The problem of converting descriptions in quantitative damage percentages can be affected by a high uncertainty.For the probability of occurrence of the events, a COV of 0.2 was assigned, considering that safety plans, financed by the same companies and compiled by private consulting groups, can be affected by an underestimation of hazard.
Industrial risk is localised only in few census parcels (Fig. 12).The probability of occurrence and the impacted area for each scenario were derived, together with vulnerability, from the available safety plans, provided by the companies (Table 7).
FOSM, MC and PE show quite similar values of expected losses (Table 8), orders of magnitudes smaller than flood and seismic risk.This is due to the fact that the impacted area is small, and the probability of occurrence from the safety plans is extremely low.As discussed in the previous chapter, industrial risk is underestimated because we analysed only explosions for major risk plants, thus neglecting other type of scenarios (toxic releases, etc.) and small industries.Due to the typology of source of hazard (e.g.tanks and parts of the plants), damages produced in case of accident affect a small area around the plants.In this case, the use of census parcels as territorial units spreads the effects over a larger area, which is not always realistic.COV values have similar patterns with the three methods, reaching about 69 % of the expected value.This result is justified by the assignment of high COV values to the vulnerability of exposed elements, because this vulnerability is described qualitatively in the safety plans and the position of elements within the parcels is largely unknown.

Risk assessment
The presence of multiple sources of hazard in areas with a strong industrial and urban development makes multiple QRA important to detect coexistence of threats and to establish priorities for mitigation strategies.The comparison of risks that are different in nature (natural vs. technological, spatially distributed vs. spot, frequent vs. extremely rare) demands a quantitative risk assessment (QRA) in order to deal with a common and comparable measure (e.g. the expected annual loss).
For our case study, we observe that seismic risk affects the whole study area with higher expected losses, E(L) equal to 64 000 000 C (Fig. 13).On the contrary, industrial accidents produce limited expected consequences, mostly due to a low occurrence probability and a very limited impact area outside the plant, which was the only type of impact considered in this study.This result reflects the characteristics of the specific study area and of its industrial context.expected losses for single parcels, seismic risk results to be dominant over 93.8 % of the parcels (corresponding to about 96 % of the area).Flood dominates over the rest 6.2 % of the parcels, which are located in proximity to the Mella river, mainly in downtown Brescia.The total expected losses obtained using MC, FOSM and PE are similar when considering the whole study area (Tables 3, 6, and 8).However, the comparison performed for losses at each single parcel shows significant differences among the methods.In order to highlight these differences, pair-wise comparisons of the results provided by the three methods were performed for each parcel, and the frequency distributions of the expected loss differences was analysed (Fig. 14).The differences observed between MC and the other two methods always show a bimodal distribution, which is due to the way the uncertainty related to the number of buildings of different type for each parcel is treated.For completely urbanised parcels (49 % of the total), we assumed no uncertainty on the number of buildings, because a further urbanisation would be impossible or extremely limited, or eventually just connected to changes in the building function.For the other parcels, the uncertainty was modelled by means of a Pareto pdf in MC, or a COV of 0.1 in PE and FOSM.Being the Pareto distribution strongly skewed toward positive values, expected losses for these parcels in MC (i.e. the mean value of the loss distribution) are systematically higher than in the two other methods, for which the underneath distribution is always assumed to be symmetrical.The use of skewed distributions in a probabilistic risk analysis provides strong differences in the results obtained using different methods.In our analysis, this effect is undetectable on the overall result because a skewed distribution was used only for one variable and 51 % of the parcels.In other circumstances, this effect could be much greater.This would make the MC preferable, because it is more suitable for accounting for skewness.Implementing skewness in PE is complicated, and impossible in FOSM.On the other hand, we should consider that dealing with non-normal distributions in risk analysis at regional scale is not so common, due to a general lack of information about the specific shape of the distribution.In conclusion, the choice of the most appropriate method of uncertainty propagation should be based on data availability, also dependent on the scale of the analysis, and on sustainability of time-computational efforts.
The calculated expected losses are probably underestimated because they include only buildings (i.e.lifelines and agricultural/natural resources are not included), without a specific discrimination according to use (e.g.hospital vs. private house).
In this work, no domino effects were considered.Nevertheless, due to the presence of a high density of industrial plants, also belonging to those defined by Council Directive 96/82/EC as Major Risk Plants, the possible consequences of domino effects are not beforehand negligible.
Multiple risk assessment could represent an important tool for decision making and territorial planning, providing the information related to the typology of threats present in each place, to the related level of risk, and to their simultaneous presence, which could generate domino effects and increase the effects of each single scenario.An effective territorial planning should account for data provided by this type of study, compatible with economic and environmental needs.For instance, new buildings could be planned with technical features according to the threats they will be exposed to, at a certain place.The planning of new residential or productive areas, and most of sensitive structures (e.g.schools, hospitals) could be performed choosing less exposed areas.Protective structures or non-structural mitigation measures could be implemented where their benefits could be maximised.

Uncertainty
The propagation of uncertainty in multiple QRA has the purpose of showing to the decision maker the effect of a lack of knowledge or precision in the calculation of risk.This is fundamental, also considering that often decisions need to be made on the basis of the available data, even if they are not detailed or complete at best.Uncertainty propagation, furthermore, provides a heuristic quantitative estimation of the level of reliability of the obtained results.This is useful in decision making, to decide whether a decision can be taken in presence of a certain degree of uncertainty (which is in this case acceptable), or if uncertainty must be reduced before coming to a decision.In this case, efforts to obtain a more precise result can be addressed to the right target.The three methods for uncertainty assessment provide a similar coefficient of variation of the output, ranging in mean from 0.1 to 1.2 (Figs. 6 and 7, Tables 6 and 8).
For flood and seismic scenarios, uncertainty results to be lower for small, densely urbanised parcels.In these parcels, the uncertainty related to the number of buildings is null, and the uncertainty related to the value of residential and nonresidential buildings ( C m −2 ) is low, since values are quite uniform in small areas.On the other side, in large less urbanised parcels, the value of buildings can significantly vary, introducing a higher uncertainty.The size of census parcels is a variable that must be considered with great care when QRA is performed.In general, smallest census parcels are located in highly urbanized areas where the number and value of exposed elements can be high.On the contrary, the largest census parcels are located in the suburban areas where the number of exposed elements can be relatively low but also where there is the maximum potential for future developments or for uncontrolled or rapid urbanization.This causes large uncertainty relative to number and type of exposed elements.Moreover, large census parcels are affected by uncertainty for the variables describing the event intensity (e.g.water depth, PGA, landslide velocity or size).Census parcels are the most commonly adopted land units for collecting and representing statistical information/data available at a regional/national scale (e.g.number of buildings of different types per parcel).Nevertheless, for many processes, census parcel has not always the optimal spatial resolution to account for the shape of the source areas, the propagation of the events, the exact area of impact.
The effects of different vulnerability models on flood and seismic risk is shown in Fig. 15.For floods, the use of different vulnerability models leads to large differences in total risk value, with respect to USACE (2003): 15 % for Dutta et al. (2003), and38 % for ICPR (2001).For earthquakes, the  Kostov et al. (2004), Kappos et al. (2006), Hancilar et al. (2006), andRota et al. (2008) vulnerability models, respectively.differences with respect to Kappos et al. (2006) vulnerability model amounts to 19 % for Kostov et al. (2004), 38 % for Rota et al. (2008), and 72 % for Hancilar et al. (2006).In some parcels, local differences can reach values up to 86 %.The observed results suggest that the choice of vulnerability models has a strong impact in the whole risk assessment, introducing significant uncertainty.In general the choice of a model should be based on (1) the geographic context in which the model was developed, which should be as similar as possible to the study area; (2) the detail of the model: if available, a model considering highly detailed data (e.g.different typologies of buildings, different number of floors, number of damage curves etc.) could provide a better estimate of the level of damage; (3) the conservativeness of the model.

Conclusions
Multiple quantitative risk assessment performed with three different methods for uncertainty propagation (MC, FOSM and PE) leads to similar overall results for the whole study area.However, significant differences of expected value of total loss have been observed for single parcels (up to 19 %).This difference could be relevant when using the results of risk assessment for decision making, or territorial planning, at local scale.In particular, large differences are observed when skewed distributions are used in the risk analysis, because of different capability of the three methods to account for skewness.In this case MC analysis is preferable.
Although the use of census parcels works well with societal and statistical data, we show that this type of terrain unit introduces large uncertainty related to number and type of exposed elements, representative values for the variables describing the event intensity, degree of exposure of the elements to the threat.
The proposed methodology for quantitative risk and uncertainty assessment can be applied all over the national territory, being based on scenarios defined by national law and zonations and data available at the national scale.However, since these scenarios can be produced with simplified approach at small scales, the analysis can be affected by large uncertainties, which has to be considered.
Multiple QRA associated with uncertainty evaluation can represent an important tool for -territorial planning and development according to the type of threat, the coexistence of different risks, and their level; -priority assessment in fund allocation and mitigation measures implementation; -detection of risk hot spots where a more detailed risk assessment at a higher scale could be performed.
In all these situations, a quantification of risk must be supported by uncertainty: all the assumptions and lack of knowledge introduce an error which has to be accounted for.

Fig. 2 .
Fig. 2. Values of peak ground acceleration PGA, in g, with exceedance probability of 10 % in 50 yr, and location of historical earthquakes.Magnitude is expressed as Maw, momentum magnitude.

Fig. 3 .
Fig. 3. Example of calculation of flood P (I |event) for census parcels.

Fig. 4 .
Fig. 4. Economic value in euro m −2 , W , for residential buildings (a, c), non-residential buildings (b, e).Frequency of surface of census parcels in the study area (d).
for non-skewed independent variables was used.The procedure calculates the first and second moments of each function points with coordinates corresponding to the expected value of each variable x i , increased or reduced of one standard deviation θ and COV θ are the mean, the variance and the coefficient of variation of the output function; µ Xi , θ and COV θ are the mean, the variance and the coefficient of variation of the 7 output function, µ Xi, Ϭ Xn and COV i are the mean, the variance and the coefficient of variation 8 of the input variable x i .9In the propagation of uncertainty by means of PE, the Roesenblueth's procedure (1981) for 10 non-skewed independent variables was used.The procedure calculates the first and second

Fig. 10 .
Fig.10.Expected losses as a function of exceedance probability for seismic risk(Kappos et al., 2006 vulnerability model).Grey curves for a random selection of 100 census parcels, ranging in the grey area between the minimum and the maximum level of risk.Curve for the whole study area is shown in black.

Fig. 11 .
Fig. 11.Values of seismic risk for the 2500-yr scenario, and of the related uncertainty, for the three methods, calculated by means of Kappos et al. (2006) vulnerability model.Parcels are ordered according to the level of risk.Light lines show the risk level increased and diminished of the standard deviation, dark line indicates the value of risk.All the 2460 parcels of the study area are impacted.

Fig. 13 .
Fig. 13.Expected losses as a function of exceedance probability for the analysed scenarios.Expected losses are those provided by MC.Value of flood risk under the curve corresponds to 9 643 805 C (US-ACE, 2003 vulnerability model).Seismic risk equals 63 897 207 C (Kappos et al., 2006, vulnerability model), industrial risk is 2956 C.

Fig. 14 .
Fig. 14.Frequency of the differences in expected loss values between couples of methods for (a) flood 500-yr scenario, and (b) seismic 475-yr scenario.Arrows indicate the mean value of the differences for the whole study area.

Table 1 .
Example of the application of the three flood vulnerability models to a census parcel containing one 1-floor non-residential building completely impacted by the flooding corresponding to a return period of 500 yr.

Table 2 .
Uncertain variables, and values for coefficients of variation (COV) adopted in the analysis.P (event) is the annual exceeding probability of the scenario, W the economic value.

Table 5 .
Example of calculation of seismic vulnerability of a census parcel containing one 1-floor non-residential building shaken by a T r 475 yr earthquake.

Table 6 .
Kappos et al., 2006ses and COV values for seismic scenarios with return time (T r ) 75, 475 and 2500 yr, vulnerability model fromKappos et al., 2006.Range of values for all the census parcels, and total on the study area.Comparison of the three methods: Monte-Carlo-Simulation, First Order Second Moment, and point estimate

Table 7 .
Vulnerability, P (L|I ) to industrial risk as from safety plans.

Table 8 .
Expected annual losses E(L) (C yr −1 ) and COV values for industrial scenarios.Comparison of the three methods: Monte-Carlo-Simulation, First Order Second Moment, and point estimate.