Uncertainty analysis for seismic hazard in Northern and Central Italy

In this study we examine uncertainty and parametric sensitivity of Peak Ground Acceleration (PGA) and 1-Hz Spectral Acceleration (1-Hz SA) in probabilistic seismic hazard maps (10% probability of exceedance in 50 years) of Northern and Central Italy. The uncertainty in hazard is estimated using a Monte Carlo approach to randomly sample a logic tree that has three input-variables branch points representing alternative values for bvalue, maximum magnitude (Mmax) and attenuation relationships. Uncertainty is expressed in terms of 95% confidence band and Coefficient Of Variation (COV). The overall variability of ground motions and their sensitivity to each parameter of the logic tree are investigated. The largest values of the overall 95% confidence band are around 0.15 g for PGA in the Friuli and Northern Apennines regions and around 0.35 g for 1-Hz SA in the Central Apennines. The sensitivity analysis shows that the largest contributor to seismic hazard variability is uncertainty in the choice of ground-motion attenuation relationships, especially in the Friuli Region (∼0.10 g) for PGA and in the Friuli and Central Apennines regions (∼0.15 g) for 1-Hz SA. This is followed by the variability of the b-value: its main contribution is evident in the Friuli and Central Apennines regions for both 1-Hz SA (∼0.15 g) and PGA (∼0.10 g). We observe that the contribution of Mmax to seismic hazard variability is negligible, at least for 10% exceedance in 50-years hazard. The overall COV map for PGA shows that the uncertainty in the hazard is larger in the Friuli and Northern Apennine regions, around 20-30%, than the Central Apennines and Northwestern Italy, around 10-20%. The overall uncertainty is larger for the 1-Hz SA map and reaches 5060% in the Central Apennines and Western Alps. Mailing address: Dr. Anna Maria Lombardi, Istituto Nazionale di Geofisica e Vulcanologia, Via di Vigna Murata 605, 00143 Roma, Italy; e-mail: lombardi@ingv.it


Introduction
It has become common practice to apply the Probabilistic Seismic Hazard Analysis (PSHA) method to develop seismic hazard studies for input to various aspects of public and financial policy.PSHA calculates ground motions with various probabilities of exceedance during a certain time interval.This information is used in the formulation of seismic design maps, as well as to develop national recommendations for building codes (e.g., Youngs et al., 1987;Petersen et al., 1996;Building Seismic Safety Commission, 1997;Wong and Olig, 1998).Since they influence policy decisions on issues ranging from building codes to science funding, an appreciation for the uncertainties and limiting assumptions underlying such maps is valuable for the user and decision-making communities.
One controversial issue is the nature of uncertainty and how to express it in hazard maps.It is useful to differentiate between two sources of uncertainty (or variability).The first, aleatory variability, is due to the inherent randomness of natural processes and is not reducible with better data or models (e.g., locations of major earth-quakes, rupture directions, occurrence times).The second, epistemic variability, results from lack of data or uncertainties in knowledge (e.g., long-term rates of seismic activity, maximum event magnitudes, ground-motion response), and is reflected in differences between models.The random variability is explicitly taken into account in Probabilistic Seismic Hazard Analysis (Reiter, 1990), but the epistemic uncertainty must be incorporated by a suitable method.In this analysis we consider a logic tree of model-variable alternatives.In order to evaluate the ground motion level with a fixed probability of being exceeded during a certain interval time, exceedance rates for each alternative model are added with weights assigned in the logic tree.When the number of alternative models is high or some parameters have a continuous distribution, then a series of Monte Carlo runs can be used.
The problem of assessing the effects of variability in ground-motion predictions has been discussed for some decades.The first quantitative results on this subject were presented by McGuire (1977) and McGuire and Shedlock (1981).These authors pointed out the ambiguities concerning the definition of some input variables for hazard calculations, and quantified their effects on ground-motion predictions for the Eastern United States.They provided a computational method for determining variability in seismic hazard calculations which result from statistical uncertainty in the assumed model or parameters.Moreover they demonstrated that, for the best estimates of seismic hazard by logic tree, only mean values of uncorrelated and symmetrical parameters need to be considered, and that explicit inclusion of their uncertainty is unnecessary.This consideration was confirmed later by Bender and Perkins (1993), who provided important suggestions about the use of different attenuation relationships in hazard calculations, choice of weight for each of them, and procedures to evaluate their sensitivity on ground-motion predictions.
The first complete description of a methodology to estimate the overall uncertainty and parameters sensitivities in hazard assessment was provided by Cramer et al. (1996) and Cramer (2001).This approach allows one to account for parameters with continuous distributions by the introduction of a Monte Carlo series of calcula-tions.An advantage of this approach is the assessment of individual branch-point sensitivity on the overall uncertainty.Their results show that the overall variability of ground motions predictions is up to 50% in California and the New Madrid seismic zone, due to many parameters related to seismic sources, activity rate and attenuation relations.Moreover, attenuation relation is identified as the variable with the largest effect on ground-motion predictions, together with parameters related to geological information (rupture models, slip rate and Mmax of faults, magnitude-frequency distribution).The same procedure has been extended by Cramer et al. (2002) to the smoothed seismicity based hazard methodology, developed by Frankel (1995).Considering that this method does not need to take into account geological information about single faults, the uncertainty in hazard calculations is due only to variables related to seismicity and ground-motion (activity rate, attenuation relationship, minimum magnitude, Mmax and completeness period).In the southern Illinois Basin, they found that the overall variability of PGA values with 2% probability of exceedance in 50 years was up to 70%.In their calculations, the most sensitive branch points were catalogue lower magnitude, activity rate and ground motion attenuation relation.
Mean hazard maps of PGA with 10% probability of exceedance in 50 years were generated by Akinci et al. (2004) following the smoothed seismicity approach described by Frankel (1995) for three zones of Italy (Eastern Alps, Western Alps and Apennines).New regional ground-motion relationships were considered along with established models developed by Ambraseys et al. (1996) and Sabetta and Pugliese (1996).They were incorporated into hazard calculations by a logic tree of model-variable alternatives which include three parameters: attenuation relationship, b-value and M max.Since it is the first time for Italy that we introduce new inputs (e.g., regionalized attenuation relationships) into the hazard calculations and use a smoothed seismicity approach, we believe that an uncertainty analysis is necessary to quantify the variability of both peak ground acceleration and 1-Hz spectral acceleration for each studied macro zone.To do so, we use a Monte Carlo approach and fol-low the procedure used by Cramer et al. (1996Cramer et al. ( , 2002) ) and Cramer (2001).

Seismic hazard analysis
The procedure described by Frankel (1995) was used by Akinci et al. (2004) to generate hazard maps for Northern and Central Italy.The approach of using spatially-smoothed historical seismicity for the hazard calculations is different from the one used previously by Slejko et al. (1998), Romeo et al. (2000), and by Working Group-2004 for the National Seismic Hazard Maps for Italy (Working Group, 2004), in which source zones are drawn around the seismicity and the tectonic provinces.In its purest form the smoothed-seismicity method simply assumes that patterns of historical earthquakes predict fu-ture activity, but it can easily be supplemented by tectonic or geodetic-based zones or other model elements if seismicity catalogues are insufficient.
To apply Frankel's (1995) procedure we have used cells equal to 0.1°in latitude and longitude, and chosen a threshold magnitude equal to 4.0.The Gaussian correlation distance considered is equal to 25 km, as estimated by Console and Murru (2001).This parameter permits the best smoothing of seismic activity: larger values spread out seismicity too much, whereas smaller values emphasize effects due to clustering.
The historical catalogue used is CPTI, Catalogo Parametrico dei Terremoti Italiani (Gruppo di Lavoro CPTI, 1999).Considering the different crustal and tectonic characteristics of the region (Meletti et al., 2000;Montone et al., 2003;Akinci et al., 2004) as well as the completeness properties of the seismic catalogue Fig. 1.Regionalization of the study area.The maximum magnitude, Mmax is assumed 7.5 for the Central Apennines, 5.5 for volcanic regions and 7.0 elsewhere.The b-value is assumed constant and is estimated separately in each macro zone: values obtained using a Maximum Likelihood Method are 0.75 (Northwestern Italy), 0.71 (Northeastern Italy) and 0.72 (Central Italy).Circles show locations of historical events of CPTI catalogue (Gruppo di Lavoro CPTI, 1999) with magnitude larger than or equal to 5.0.
Figure 2a,b shows comparisons between the three regional attenuation relations and the empirical regressions developed by AMB96 and SP96 for a rock site; each is plotted as a func-Fig.2a,b.Comparison between a) Peak Ground Acceleration (PGA) predicted by regional attenuation relationships (Malagnini et al., 2000(Malagnini et al., , 2002;;Morasca et al., 2004) for M=5.0 and M=7.0 (solid line) with the results of the empirical regressions by Ambraseys et al. (1996) (dotted curves) and Sabetta and Pugliese (1996) (dashed curves); b) 1-Hz Spectral Acceleration (SA).a b (Camassi et al., 2002), the study region has been divided into three macro zones: Northwestern Italy, Northeastern Italy, and Central Italy (see fig. 1).
Each zone has been characterized by three variables: b-value, maximum magnitude and attenuation relationship.The magnitude distribution is assumed to follow the exponential Gutenberg-Richter relationship (Gutenberg and Richter, 1949) in which b indicates the relative frequency of small and large magnitudes.Here it has been assumed constant and estimated using the maximum-likelihood method (Bender, 1983), for each zone.The values obtained are 0.75 for Northwestern Italy, 0.71 for Northeastern Italy and 0.72 for Central Italy.Even though these three values are very similar, hazard calculations are quite sensitive to the b-value, so that we prefer to use these three values separately at each region, instead of using one single parameter for three regions.
In the CPTI catalogue different types of magnitudes are available, but the only magnitude reported in the database for all earthquakes is Ma.It is similar to moment magnitude for magnitude greater than 5.5 (see for details Gruppo di Lavoro CPTI, 1999;and Akinci et al., 2004).We set M max equal to 7.0 for most of Italy.Since the Apennines'axis region has a history of large earthquakes of magnitude around 7.0, and has the capacity of producing larger earthquakes (events with magnitude M>6.0 mostly occurred in this zone), the maximum magnitude is taken as 7.5 in this zone (fig. 1, dashed zone).Following Model PS4, used for the previous National Hazard Map for Italy (Albarello et al., 2000), we set M max equal to 5.5 for volcanic regions (Alban Hills, Amiata, Roccamonfina) (see fig. 1).
New regionalized predictive ground-motion relationships had been introduced into the maps, developed by Malagnini et al. (2000Malagnini et al. ( , 2002) ) and by Morasca et al. (2004) for the Apennines, Northeastern Italy and the Western Alps, respectively.In each macro-zone the appropriate regional relationship has been integrated with ones derived from strong motion data by Ambraseys et al. (1996, hereafter AMB96) and Sabetta and Pugliese (1996, hereafter SP96) using a logic tree model.Assigned weights are 0.5 for the regional models and 0.25 for the other two relations.tion of distance for earthquakes of magnitude M=5.0 and M=7.0, for PGA (fig.2a) and 1-Hz SA (fig.2b).The standard deviation reported by authors for each attenuation relationship has been considered in the hazard calculation.
The PGA and 1-Hz SA hazard maps, for 10% probability of exceedance in 50 years are shown in fig.3a,b, respectively.We used three parameters for hazard assessment: b-value (calculated for each macro-zone separately), Mmax (taken for each single region separately) and the three attenuation relationships weighted as 0.25 (SP96), 0.25 (AMB96) and 0.5 (regional relations).Details on methodology, historical catalogue, tectonic and geological setting are reported in Akinci et al. (2004).

Methodology for uncertainty analysis and alternative parameters
We investigate the overall uncertainty and parametric sensitivity in the ground-motion for Northern and Central Italy.Since the seismic hazard maps are based on historical seismicity (Akinci et al., 2004) and not on the geological information coming from single faults, we study only three branches of the logic tree: ground motion attenuation relationship, b-value and maximum magnitude.
Even though the correlation distance, the minimum magnitude and the catalogue completeness are three of the input parameters for the smoothed seismicity approach, we do not consider uncertainty due to them.The optimal correlation distance of 25 km was already obtained by Console and Murru (2001) for Italy using a trial-and-error procedure.Moreover good information about the historical seismicity of Italy disclosed the minimum magnitude and the catalogue completeness periods with good reliability (Camassi et al., 2002).Assuming no epistemic uncertainty for these three parameters and for the underlying Gutenberg-Richter Law, the computed activity rate depends on the b only and we do not include uncertainty in the other smoothed-seismicity parameters.
The results of our analysis are synthesized in two maps: the overall uncertainty maps provide a confidence interval for the PGA and 1-Hz SA values and the parameter uncertainty maps determine the sensitivity of hazard to variability of each logic tree branch.
In order to generate parametric uncertainty maps for b and M max, we use the well known Monte Carlo approach (Cramer et al., 1996).First we vary the examined parameter in a series of Monte Carlo calculations while the other branches of the logic tree are fixed.Then, the ground-motion predictions of PGA and 1-Hz SA (with 10% probability of exceedance in 50 years) at each grid point are sorted into ascending order, and the 2nd, 50th and 98th percentile values are chosen.In fig. 4 we show the distribution of 100 simulated PGA values at grid point 46.5N, 13.0E used in order to obtain the b-value uncertainty.The uncertainty maps are generated by selecting the larger difference between the 50th and 2nd percentiles and the 98th and 50th percentiles.This represents a 95% confidence band that is a value that can be added/subtracted to the median to compute 95% confidence limits of a distribution; it indicates variability of the hazard predictions by changing values of the examined input parameter.The 95% confidence interval in a normal distribution is given by the range between the 97.5 and 2.5 percentiles.If the simulated ground motions have a normal distribution, the 95% confidence band is an estimator of two standard deviations (Cramer et al., 1996).
Since often the 95% confidence band increases with the average ground motion, uncertainty results are often better indicated by the Coefficient Of Variation (COV).We compute it as the ratio of half of the 95% confidence band divided by the mean.This value, multiplied by 100, indicates a percent variability of the ground motion.The smaller the COV, the more reliable is the estimate.
For generating uncertainty maps, we performed 100 calculations, which is the minimum number of simulations needed to obtain stable estimates of uncertainty (Cramer et al., 1996).
Some trials with more simulations (up to 1000) do not show significantly different results.
For each macro zone (Northwestern Italy, Northeastern Italy, and Central Italy) the b-value is represented by a normal random variable with mean equal to the estimated value and standard deviation of 0.1.The histograms in fig.5b show 100 simulations of b-values for each zone.The Mmax is assumed as a discrete uniform random variable for the three values that are considered here: (5.5, 7.0 and 7.5) with ranges of variability of [5.3, 5.7], [6.8, 7.2] and [7.3, 7.7], respectively, based on the observed seismic history in each region.Five values are considered for every range (one decimal digit approximation) with an equal probability of 0.2. Figure 5a shows the histograms of 100 simulated values considered for each Mmax value.We use a uniform distribution because we want assign the same probability to five values in each range.
In order to perform the uncertainty analysis of the attenuation relationships, we have chosen an alternative method, which considerably reduces the computation time.If only the variation of this branch is taken into account, predicted PGA at a point with coordinates (i, j) can assume one of three values , provided by the three attenuation relationships.The probability for each of these three values is the weight (indicated by pk, k=1, 2, 3) assigned to the corresponding attenuation relationship in the logic tree.Then PGA is a discrete random variable, with standard deviation given by ( , , x x x ( , ( , ( , ) ) ) (3.1) where (mean while).We consider the attenuation relationship 95% confidence band equal to 2⋅s.
To generate overall uncertainty maps we perform 100 Monte Carlo simulations, and all three branches are simultaneously varied (using the distributions given in fig.5a for Mmax and 5b for b-value).Since the branches of logic tree are statistically independent, each parameter map represents the contribution to the overall uncertainty.

Results
Applying the above-described methodologies, we perform the overall and parametric uncertainty analysis for hazard assessment in 95% Confidence band maps -The map of mean PGA with 10% probability of exceedance 0.08-0.1 g (fig.6a), which is slightly higher than the b-value uncertainty, 0.06-0.08g (fig.7a).
In Northern Apennines contributions to overall uncertainty by attenuation relationships and b-value are very close (0.04-0.06 g).
A larger difference between b-value and attenuation relationship 95% confidence band is found in the Central Apennines region.The bvalue variability is similar to the Friuli Region (0.08 g) and is higher than attenuation relation-   ship variability (0.05 g).This is because, as other investigations have shown, large earthquakes have much leverage on hazard in both zones, and the regional attenuation relationship used in Central Italy (Malagnini et al., 2000) is closer to AMB96 and SP96 than the one for Northeastern Italy (Malagnini et al., 2002) for large magnitude (see fig. 2a).Moreover, the PGA 95% confidence band of attenuation relationship (fig.6a) is slightly lower in the Northern  Apennines than in Friuli, although the same regional relation, b-value and maximum magnitude are used for both zones.This is because the Northern Apennines has more modest seismic activity (no historical event with magnitude larger than 6.0) than Friuli (see fig. 1) and the differences between the three attenuation relationships decreases with magnitude, especially for distances shorter than 50 km (fig.2a).
The overall PGA 95% confidence bands (fig.9a) are about 0.02-0.04g in Northwestern Italy, where the lowest PGA predictions are observed.Comparing the ground-motion relationship and the b-value 95% confidence band maps (figs.6a and 7a) in this zone, it is clear that the sensitivity on hazard predictions of the first parameter variability is generally larger; only in the Western Alps do these two input parameters have comparable effects.
The overall 95% confidence band map for 1-Hz SA (fig.9b) shows larger variability than PGA in the Central Apennines, up to 0.35 g (figs. 2b and 6a,b).
In all three regions, the contribution of Mmax to the seismic hazard variability is negligible (fig.8a,b).Similar results were obtained also by Cramer et al. (2002) in the Southern Illinois region (U.S.A.).
Coefficient of variation maps -COV is often a clearer estimator of uncertainty than 95% confidence band.According to the overall COV map for PGA (fig.12a) the uncertainty in hazard is largest in the Northern Apennines, in Friuli and in Northwestern Italy (COV=20--30%).In the first two regions the overall COV is the sum of about equal contributions (10--15%) provided by the ground-motion relationships (fig.10a) and b-value (fig.11a).
In Northwestern Italy the overall COV is essentially due to attenuation relationship uncertainty, which is larger here than in all the other zones.This confirms the results shown in the 95% confidence band maps.Unlike the other zones, in this region the predicted ground-motion values are strongly controlled (more than 70%) by earthquakes with magnitude less than 5.5, occurring at middle to long distance.Comparing trends of PGA attenuation relationships in the logic tree for Northwestern Italy and for M=5.0 (fig.2a), the differences between predicted PGA values for epicentral distance greater than 20 km are clear, and they increase  with increasing distance.The remaining 30% of PGA predictions are essentially due to stronger earthquakes occurring at larger distance: the three attenuation relationships predict rather different PGA values for these events (fig.2a).In all other zones (Western Alps, Tuscany and Central Apennines), the percent variation of PGA is about 10-20% and the impacts of attenuation relationship (fig.10a) and b-value (fig.11a) are essentially comparable.
The largest values of overall COV for 1-Hz SA occur in the Central Apennines and Western Alps, up to 50-60% (fig.12b).Figures 10b and  11b show roughly equal contributions from attenuation and b-value uncertainties in the Central Apennines, while attenuation uncertainty dominates in the Western Alps.The Mmax COV maps of 1-Hz and PGA confirm that this parameter has a negligible effect on hazard assessment in the studied region.They are not reported because the percent variation is less than 10%.

Conclusions and discussion
A Monte Carlo logic tree approach has been applied to the uncertainty analysis of a PSHA model for Northern and Central Italy, for 10% probability of exceedance in 50 years only.PGA uncertainty, as represented by overall COV maps, is equal to 20-30% in the Friuli region, the Northern Apennines, and in Northwestern Italy near Milan, and it is about 10% elsewhere.Overall COV in 1-Hz SA hazard ranges from 50-60% in the Central Apennines and the Western Alps to 20-30% elsewhere.The largest contributions to variability in both PGA and 1-Hz SA arise from uncertainties in groundmotion relationships and b-values, while the effect of M max is negligible at 10% probability of exceedance in 50 years.
Knowledge of the most critical branch points can be a guide for future research to reduce uncertainty in these calculations.Our results suggest that the greatest reduction in uncertainty in seismic hazard calculations for the studied regions could be achieved by continued research on strong-ground motion.
Since we do not introduce geological information for large earthquakes, our results con-cern variability due only to factors related to seismicity and ground-motion.Obviously fault parameters are very important sources of uncertainty in hazard assessment.Further studies will concern introduction of faults parameters in hazard model and estimation of sensitivity of ground motion predictions by them.Introduction of geological information will require modeling of further sources of epistemic uncertainty such as the magnitude-frequency distribution (Gutenberg-Richter, characteristic magnitude model) or the recurrence model (Poisson model, time-dependent models, etc.).
Seismic hazard maps are useful as long as their limitations are recognized.In order to make the maps useful for engineers, insurance analysts and emergency planners as well as for public policy, input data, assumptions, and model limitations should be forthrightly discussed.We believe that it will be useful to consider multiple hazard estimates developed by various governmental, academic and commercial groups under different assumptions.This is one of the reasons why we used a different and new approach, with new inputs for the hazard calculations in Italy, and discussed the variability of predictions and ground motions.
Northern and Central Italy.The mapped values are the 95% confidence bands (figs.6a,b to 9a,b), which indicate a confidence interval about the expected value of PGA and 1-Hz SA, and the coefficients of variation (figs.10a,b to 12a,b), which quantify a percent variation of ground-motion predictions.

Fig
Fig. 5a,b.a) Logic tree alternatives for maximum magnitude.It is assumed a discrete random variable uniformly distributed on the five values reported.Histograms show the 100 simulated values used in uncertainty analysis.b) Histograms of 100 simulated values for b-value uncertainty estimate in each of three zones.In every macro-zone this parameter is represented by a normal random variable with standard deviation 0.1 around the estimated value.

Fig
Fig. 8a,b.Parameter 95% confidence band map in Northern and Central Italy for maximum magnitude, a) PGA and b) 1-Hz spectral acceleration.

Fig
Fig. 9a,b.a) Overall PGA 95% confidence band map in Northern and Central Italy.b) Overall 1-Hz spectral acceleration 95% confidence band map in Northern and Central Italy.

Fig
Fig. 10a,b.a) PGA and b) 1-Hz SA Coefficient Of Variation (COV) maps in the Northern and Central Italy for ground-motion relationships.

Fig
Fig. 11a,b.a) PGA and b) 1-Hz SA Coefficient of Variation (COV) maps in the Northern and Central Italy for b-value.

Fig
Fig. 12a,b.a) The overall PGA Coefficient of Variation (COV) map in the Northern and Central Italy.b) The overall 1-Hz SA Coefficient Of Variation (COV) map in the Northern and Central Italy.