Articles | Volume 21, issue 8
Research article
03 Aug 2021
Research article |  | 03 Aug 2021

Integrating macroseismic intensity distributions with a probabilistic approach: an application in Italy

Andrea Antonucci, Andrea Rovida, Vera D'Amico, and Dario Albarello

The geographic distribution of earthquake effects quantified in terms of macroseismic intensities, the so-called macroseismic field, provides basic information for several applications including source characterization of pre-instrumental earthquakes and risk analysis. Macroseismic fields of past earthquakes as inferred from historical documentation may present spatial gaps, due to the incompleteness of the available information. We present a probabilistic approach aimed at integrating incomplete intensity distributions by considering the Bayesian combination of estimates provided by intensity prediction equations (IPEs) and data documented at nearby localities, accounting for the relevant uncertainties and the discrete and ordinal nature of intensity values. The performance of the proposed methodology is tested at 28 Italian localities with long and rich seismic histories and for two well-known strong earthquakes (i.e., 1980 southern Italy and 2009 central Italy events). A possible application of the approach is also illustrated relative to a 16th-century earthquake in the northern Apennines.

1 Introduction

Characterizing earthquake effects on the anthropic environment is of paramount importance for estimating seismic risks and planning prevention policies. This characterization is performed by classifying earthquake effects according to macroseismic scales. Each macroseismic scale considers a set of scenarios, 12 in the most used scales in Europe (i.e., MCS – Sieberg, 1932; MSK – Medvedev et al., 1964; EMS-98 – Grünthal, 1998), ordered in terms of increasing severity of the effects. Through macroseismic scales, observed seismic effects concerning human behavior, damage to buildings and geomorphological phenomena at a site are compared with the scenarios proposed in the scale to assess an intensity value. An intensity value, referring to a specific earthquake and a specific place identified through its geographic coordinates, defines the intensity data point (IDP in the following). The spatial distribution of IDPs is considered for the characterization of the earthquake source (i.e., estimates of epicentral location and magnitude) in the absence of instrumental data (e.g., Bakun and Wentworth, 1997; Gasperini et al., 1999, 2010; Provost and Scotti, 2020). Collecting these parameters in homogeneous seismic catalogues (e.g., Fäh et al., 2011; Stucchi et al., 2013; Manchuel et al., 2018; Rovida et al., 2019, 2020) is a key element to providing a seismic characterization of a region, and this information represents a basic tool for seismic hazard estimates (e.g., Stucchi et al., 2011; Woessner et al., 2015; Meletti et al., 2021). In particular, in countries with rich macroseismic data (e.g., Italy and France) the histories of documented earthquake effects at a site can be consistently used for local seismic hazard assessment (e.g., Albarello and Mucciarelli, 2002; D'Amico and Albarello, 2008). Moreover, macroseismic intensity can be useful to check the outcomes of probabilistic seismic hazard assessments (Stirling and Petersen, 2006; Mucciarelli et al., 2008; Rey et al., 2018), especially in countries where the historical record is much longer than the instrumental one.

Retrieving seismological data from documentary information requires specific methodologies and expertise (e.g., Guidoboni and Stucchi, 1993) and presents several criticalities mainly due to the incompleteness of the documentation (e.g., Albarello et al., 2001; Swiss Seismological Service, 2002; Stucchi et al., 2004; Hough and Martin, 2021). The probability of retrieving such documentation depends on the period, size and location of the event and is hampered by the survival of sources and the capability of retrieving and analyzing them (e.g., Albini and Rovida, 2018; Albini, 2020a, b). This implies that intensity distributions of historical events may present important gaps, which depend also on the density and importance of the settlements affected by the earthquakes.

To fill these gaps, documented seismic effects may be integrated with “synthetic” intensities, which can be estimated in different ways. Until the second half of the 20th century, qualitative contouring procedures were used to draw isoseismal maps (e.g., Shebalin, 1974; Barbano et al., 1980; Postpischl et al., 1985; Ferrari and Postpischl, 1985; Patané and Imposa, 1985), in which hand-drawn isoseismals bounded areas enclosing sites with intensity overcoming any given threshold (Musson and Cecić, 2012). This form of regularization aims at reconstructing a general radiation pattern for historical earthquakes but is affected by biases induced by the conceptual background of the “tracer”. To overcome this drawback, some authors (Ambraseys and Douglas, 2004; Rey et al., 2018) proposed geostatistical approaches (e.g., Olea, 1999) to identify areas affected by similar seismic effects. They applied the kriging spatial interpolation technique to compute the expected values of macroseismic intensity through a semivariogram that describes the correlation between neighboring IDPs. This kind of approach, however, disregards the inherent ordinal and discrete nature of intensity data, which requires specific formalizations to account for uncertainty affecting intensity estimates (see, e.g., Magri et al., 1994; Albarello and Mucciarelli, 2002).

An alternative approach to obtaining synthetic intensities makes use of intensity prediction equations (IPEs), which provide the possible intensity values at any site as a function of the epicentral distance and maximum or epicentral intensity or magnitude (e.g., Pasolini et al., 2008; Sørensen et al., 2009; Allen et al., 2012; Rotondi et al., 2016). The limitation of this approach is the hypothesis that the radiation pattern of seismic waves from the source is the only one responsible for the intensity at a site, disregarding lateral heterogeneities induced by the fracture process and geological and/or geomorphological features.

To account for these features, we present an alternative probabilistic approach, which improves the one proposed by Albarello et al. (2007) and D'Amico and Albarello (2008). The key element is a combination, through a Bayesian approach, of probabilistic estimates provided by an IPE constrained by observed intensities that are spatially close to the site of interest. The proposed procedure is described first; then it is tested on a set of localities and macroseismic fields included in the Italian Macroseismic Database DBMI15 (Locati et al., 2019).

2 Methodology

In the frame of a coherent Bayesian formalization, the proposed procedure combines intensities estimated at a site with an IPE with observed intensities at neighboring localities for the same earthquake, taking into account the inherent uncertainty. To this purpose, considering any lth event, the discrete probability density distribution pl(Is|Iv) is computed to associate with each possible intensity value Is at the site s a probability value conditioned by the occurrence of effects of intensity Iv at any other site v:

(1) p l ( I s | I v ) = p l ( I s ) q ( I v | I s ) I = 1 12 p l I q ( I v | I ) .

Here pl(Is) represents the “prior” probability density which is deduced from an IPE by using the epicentral parameters (location, epicentral intensity, magnitude, etc.) of the lth event. In general, the most common IPEs (in their probabilistic formulation) have the form

(2) S ( I s | I e D ) = prob [ I I s | I e , D ] = 1 σ 2 π I s - 0.5 e - J - μ I e , D 2 σ 2 d J

(Albarello and D'Amico, 2004), where the average μ is a function of the epicentral distance D and the intensity at the epicenter Ie (or, eventually, the estimated magnitude). Both the average μ and the standard deviation σ are determined from the statistical analysis of the available information (e.g., Pasolini et al., 2008, for Italy). To account for the uncertainty affecting epicentral intensity, the marginal probability can be computed as

(3) P l I s = I e = 1 I max g ( I e ) S ( I s | I e , D ) ,

where g is the probability distribution which expresses the uncertainty affecting the epicentral intensity and Imax is the upper bound of the adopted macroseismic scale (e.g., 12 for the MCS scale). The probability density pl(Is) can be computed in the form

(4) p l I s < I max = P l I s - P l I s + 1 p l I s = I max = P l I max .

It is worth noting that Eq. (1) can be iteratively applied when an increasing number of neighboring sites are considered to constrain intensity. This can be simply performed by substituting the prior distribution with the output of the preceding estimate. The key element of Eq. (1) is the conditional probability density q(Iv|Is), which expresses the correlation between intensity values at neighboring localities. In other terms, such a probability density represents the “smoothness” of the macroseismic field and plays the role of covariance in classical geostatistics. More specifically, q(Iv|Is) expresses the constraining power of Is on Iv. According to Albarello et al. (2007), q(Iv|Is) can be estimated empirically by considering observed intensity distributions. To this purpose, data provided by the Italian Macroseismic Database DBMI15 (Locati et al., 2019) were considered for a case study.

3 Assessing the spatial variability of macroseismic data in Italy

3.1 The Italian macroseismic database DBMI15

The long tradition of historical macroseismic investigation in Italy has produced a wealth of studies and data on the seismic history of the country and neighboring areas. All such studies are collected and organized in the Italian Archive of Historical Earthquake Data – ASMI (, last access: 12 April 2021; Rovida et al., 2017) – which grants access to information on more than 6200 earthquakes that have occurred in the Italian area from 461 BCE to 2019 CE. The data gathered in ASMI are of several typologies and formats and provide a large number of intensity data from different sources, such as macroseismic bulletins, online databases, and many scientific papers and reports. The different information collected in ASMI for each earthquake requires a careful comparison in order to identify the reference study among those available for the compilation of the Italian Macroseismic Database DBMI. The current release of the latter (DBMI15,, last access: 12 April 2021; here considered in its version 2.0; Locati et al., 2019) is the result of the specific selection of these data according to the content and quality of each study and to the number and spatial distribution of intensity data. DBMI15 makes available 123 756 IDPs related to 3219 Italian earthquakes in the time window 1000–2017 and referring to 20 000 localities, of which 15 332 are in Italy. DBMI15 results from 189 different studies and the intensity data they provide are not homogenous as regards the geographic coordinates and the standards used for assessing macroseismic intensities. For this purpose, a series of operations were performed in order to obtain a homogenous set of intensity data: (i) a unique gazetteer, covering the whole national territory was adopted in order to match the position of a locality with the macroseismic observation, and (ii) a standard to express the macroseismic intensity (e.g., 6, 6–7, 7) was defined and non-conventional descriptive codes (e.g., D for damage or F for felt) were adopted when the available information was not sufficient for assessing an intensity value. DBMI15 allows direct access to seismic histories of Italian localities and provides data upon which the macroseismic parameters of the Parametric Catalogue of Italian Earthquakes – CPTI15 (, last access: 12 April 2021; Rovida et al., 2019, 2020) – are built.

Although the guidelines of the European Macroseismic Scale EMS-98 (Grünthal, 1998) recommend the users preserve the integer character of the intensity scale and avoid forms such as 6.5, 61/2 or 6+, in many studies intensity data are listed as intermediate values in order to express uncertainty affecting the intensity estimate. This is also the solution adopted in DBMI15. Following D'Amico and Albarello (2008), intensity data can be classified into two categories: “certain data” (one single intensity value I, e.g., 6) and “uncertain data” (pair of values II′′, e.g., 6–7). In the case intensity Iv is uncertain between two contiguous values Iv and Iv′′ (e.g., Iv=6–7), Eq. (1) becomes

(5) p l ( I s | I v ) = 0.5 p l ( I s | I v ) + 0.5 p l ( I s | I v ′′ ) ,

where an equal probability is assigned to the hypotheses Iv=I and Iv=I′′. This is the way uncertain intensity values contained in DBMI15 are treated in the following analysis.

3.2 Results

To estimate the probability q(Iv|Is) in Eq. (1), one can consider the relative frequencies of the differences between intensity values at pairs of sites affected by the same event. Such probability is expected to monotonically decrease with the distance between the sites, and above any distance threshold q(Iv|Is)q(Iv); i.e., Is becomes non-informative about Iv. The closer the sites considered are, the higher the informative power of Is on Iv is expected to be because closer sites possibly also share similar seismostratigraphical and geomorphological conditions. On the other hand, since the selected sites correspond to urbanized areas (settlements, villages, towns, etc.), each conventionally represented with the geographical coordinates of a single point, there is a lower limit to the distances between considered sites, which depends on the size and density of urbanized areas. Moreover, distances in the estimate of q(Iv|Is) should be large enough to include at least two sites. In other words, the search radius is selected by balancing the needs for maximizing the number of intensity data within the radius (related to the average density of settlements in Italy) and minimizing the possible geological heterogeneities present in the same area.

In order to evaluate the optimal distance threshold to characterize q(Iv|Is), the geographic distribution of the 15 332 Italian localities in DBMI15 has been investigated. In particular, for each locality, the number of localities within a set of possible distance thresholds has been computed. Figure 1 shows that there are significant parts of the Italian territory (mainly in southern areas) where the mutual distance of localities is larger than 10 km. On the other hand, for all the sites (except for nine localities on small islands), there is at least another locality within 20 km. For this reason, the latter was selected as the reference distance threshold for the characterization of q(Iv|Is).

Figure 1Number of neighboring localities within 10 km (a) and 20 km (b) for all the Italian localities in DBMI15. The 28 sites selected for testing (Sect. 4.2) are shown.

To this purpose, a dataset derived from DBMI15 has been defined by selecting 546 earthquakes with at least 10 IDPs with intensity greater than or equal to 5. These earthquakes occurred in the period 1117–2017 CE and are well distributed over the whole Italian territory. From the intensity distributions of these earthquakes, we discarded

  • i.

    non-numerical macroseismic observations (e.g., “Damage” or “Felt”);

  • ii.

    data related to unidentified localities or large areas (see Locati et al., 2019, for details);

  • iii.

    macroseismic observations related to earthquakes with epicenters inside the active volcanic areas (i.e., Mt. Etna and Campanian volcanoes), due to the faster attenuation observed in these zones with respect to the rest of Italy (e.g., Carletti and Gasperini, 2003).

We obtained 58 062 IDPs with intensity ranging from 1–2 to 11 (Fig. 2), referring to more than 12 500 Italian localities.

Figure 2Frequencies of selected intensity values related to the 546 earthquakes used in the analysis.


The conditional probability q(Iv|Is) in Eq. (1) was estimated from the relative frequencies of the differences between Iv and Is computed for each of the 546 selected earthquakes. In this analysis Iv and Is represent any pair of intensity values observed at neighboring localities (i.e., within a distance of 20 km). If the intensity values Is and Iv are both uncertain (e.g., 6–7), the two adjacent integer degrees (i.e., 6 and 7) are considered equiprobable and the differences between the four intensity values are computed; if both Is and Iv are certain values (e.g., 7), the difference is counted four times (Albarello et al., 2007). In general, one has

(6) q I v | I s = q ( Δ I | I v , I s ) ,

where ΔI=Iv-Is and the dependence of both Iv and Is is due to the lack of defined metrics for intensity degrees. As a preliminary step, we assume that

(7) q Δ I | I v , I s = q Δ I ,

which corresponds to the assumption of linear metrics for intensity values. This hypothesis will be tested below.

Two different analyses were performed in order to estimate the frequency distribution qI) from the residuals (IvIs) considering (i) only the nearest IDP within 20 km and (ii) all the IDPs within 20 km. Table 1 reports the values of qI) expressed as the relative frequency of cases where, for each earthquake, site intensity Is differs by ΔI from the intensity Iv observed at the nearest locality (qI)near) and at all localities (qI)all) within 20 km. The evident difference between the two analyses is that for qI)all the probability distribution results peaked less at ΔI=0 than qI)near. Figure 3 shows the effects of releasing the assumption in Eq. (7). To evaluate this aspect, we tested the dependence of qI) on Is. The results show that the probability of having the same intensity in the nearest locality is slightly higher for Is being equal to 7 and 8, while this probability tends to decrease for Is being equal to 6, 4, 9 and 5. The outcomes of this analysis seem to be almost independent from the intensity Is. As a consequence, we consider the approximation in Eq. (7) reliable.

Table 1Values of q(Iv|Is) expressed as the probability that site intensity Is differs by ΔI from the intensity Iv observed at the nearest locality within 20 km (qI)near) and at all the neighboring localities within 20 km (qI)all).

Download Print Version | Download XLSX

Figure 3Relative frequency of qI) as a function of intensity Is for the nearest locality within 20 km.


In both analyses, we then verified the impact on the results of the distance among localities. Figure 4 shows the effect of distance on qI). This effect is quite weak and only concerns the probability of observing the same intensity at the two sites (ΔI=0). Considering only the nearest locality within 20 km (Fig. 4a), the frequency of ΔI=0 decreases from around 52 % for the distance range 0–5 km to 42 % for the range 15–20 km. When all localities within 20 km are considered (Fig. 4b), this relative frequency decreases from 48 % to 36 % for the ranges 0–5 and 15–20 km, respectively.

Figure 4Relative frequency of qI) as a function of different distance ranges for (a) the nearest locality within 20 km and (b) all the localities within 20 km.


4 Testing

4.1 Testing procedure

To test the effectiveness of this procedure, the probability distributions in Eq. (1) were computed for a set of Italian localities and then compared with available observations. The estimated probabilities pl(Is|Iv), derived at each jth site for each lth earthquake, were used to calculate the predicted number of occurrences for each intensity degree Is (Npred) over the total of the M sites and N earthquakes considered:

(8) N pred = j = 1 M l = 1 N p l j I s | I v .

The predicted values (Npred) can be compared with the observed number of occurrences (Nobs) for the same intensity degree Is. If the observed intensity value is certain (e.g., 6), it can be expressed as

(9) N obs = j = 1 M l = 1 N I s l j .

In the case the Is value is uncertain (e.g., 6–7), an equal probability (0.5) is assigned to the two adjacent integer degrees. The central limit theorem was used to check the consistency between predicted and observed values. The statistical Z test follows the standardized Gaussian distribution:

(10) Z = N obs - N pred σ pred .

The standard deviation (σpred) associated with the predicted values was estimated as described in Albarello and D'Amico (2005):

(11) σ pred = j = 1 M l = 1 N { p l j I s | I v [ 1 - p l j I s | I v ] } .

Equation (10) can be used to evaluate the statistical significance of the discrepancy between predicted (Npred) and observed (Nobs) values. According to Albarello and D'Amico (2005), when |Z| is greater than 2, the resulting discrepancy can be considered statistically significant at the 5 % confidence level.

4.2 Application

The above approach was applied to 28 localities with at least 40 intensity data in DBMI15 homogeneously distributed over the Italian territory (Fig. 1; Table 2). As reported in Table 2, the number of intensity data, the maximum intensity and the time coverage of the seismic histories of each considered site are different. For each locality and for each earthquake, we computed the probability pl(Is|Iv) with Eq. (1) using the probability distribution qI) in Eq. (6) derived from all localities within 20 km (Table 1), by excluding the intensity observed at the site of concern. We then estimated the predicted number of occurrences for each intensity degree Is for all sites and earthquakes through Eq. (8) and compared these estimations with the observed occurrences (Eq. 10).

Table 2List of the 28 test localities with their maximum intensity (Ix), time coverage of seismic histories in terms of starting (Start) and ending (End) years, and number of intensity data (No.) as reported in DBMI15.

Download Print Version | Download XLSX

The probability pl(Is|Iv) has been computed with three different analyses for the prior distribution pl(Is) and for qI)

  • a.

    using a uniform distribution over the intensity range 2–11 for pl(Is) and the intensity observed at the nearest locality within 20 km (i.e., probability qI)near in Table 1);

  • b.

    using a uniform distribution over the intensity range 2–11 for pl(Is) and iteratively considering in Eq. (1) the intensities observed at all the localities within 20 km (i.e., probability qI)all in Table 1);

  • c.

    using as pl(Is) the probability computed through an IPE and the intensities observed at all the localities within 20 km (i.e., probability qI)all in Table 1), with the IPE defined for Italy by Pasolini et al. (2008) and recalibrated by Lolli et al. (2019) with IDPs from DBMI15 and earthquake parameters provided by CPTI15 (Rovida et al., 2019, 2020).

Using Eq. (10) the number of observed occurrences (Nobs) for a given intensity value Is was compared with the predicted number (Npred) derived for the three possible choices (a, b, c) of the prior distribution pl(Is) and qI) for the 28 selected localities (Table 3). The differences in percentage between predicted and observed values were computed and expressed as [(1-Npred/Nobs)×100]. Z represents the standardized Gaussian statistics: if |Z|>2, the resulting discrepancy can be considered statistically significant (probability 0.05).

Table 3Observed (Nobs) and predicted (Npred) number of intensity values for each analysis (a, b, c) with their differences as a percentage (Diff) and the results of the Z test (Z).

Download Print Version | Download XLSX

Table 3 shows that, for analysis c, the differences between observed and predicted values are less than 10 % for intensity 3, 6, 7 and 8. For analysis a and b, this is verified only for intensity 3 and 7 and for intensity 4 and 6, respectively. Results for intensity 2 and 11 cannot be considered significant due to the strong data incompleteness for the former (Fig. 5) and to the lack of data (only one observation available) for the latter. Furthermore, for intensity 5 there is a significant underestimation of the observed intensities in all the analyses, whereas for intensity 4 analysis b and c tend to overestimate the observed values (see Fig. 5 and Z values in Table 3). These outcomes indicate that the number of predicted values (Npred) is consistent with the number of observed occurrences (Nobs) at the 28 test localities. Among the three analyses, analysis c is more effective than the others, although the discrepancies expressed with the Z test are statistically significant for intensity 4 and 5. However, this may depend on the selected dataset because DBMI15 contains only earthquakes with a maximum intensity greater than or equal to 5 (Locati et al., 2019).

Figure 5Observed (blue bar) and predicted number of intensity values for analysis a (orange bar), b (yellow bar) and c (violet bar).


To verify the impact of using this procedure rather than the IPE alone to predict intensity values, a comparison test was carried out for two well-documented recent Italian earthquakes, i.e., the Mw 6.3 event that occurred on 6 April 2009 in the L'Aquila area (central Italy) and the Mw 6.8 Irpinia (southern Italy) earthquake of 23 November 1980. For both earthquakes we computed the differences between the observed intensity values as reported in DBMI15 (Figs. 6 and 7, respectively) and (i) the intensity values computed with the IPE by Pasolini et al. (2008) recalibrated by Lolli et al. (2019) and (ii) the intensity values estimated with the proposed procedure following analysis c described above. In both cases, the modal value of each probability distribution computed by Eq. (1) for all the sites was considered.

Figure 6Differences between the observed intensity values, as reported in DBMI15 (Galli and Camassi, 2009), for the 2009 L'Aquila earthquake and the intensity values computed with (a) the IPE alone (Pasolini et al., 2008; Lolli et al., 2019) and (b) the proposed procedure.


Figure 7Differences between the observed intensity values, as reported in DBMI15 (Guidoboni et al., 2007), for the 1980 Irpinia earthquake and the intensity values computed with (a) the IPE alone (Pasolini et al., 2008; Lolli et al., 2019) and (b) the proposed procedure.


For the 2009 L'Aquila earthquake, Fig. 6a shows that for 32 out of 315 sites (10 %) the values predicted with the IPE alone are equal to the observed ones, and for 183 sites (58 %) the predicted intensities differ by more than 1 intensity degree from the observations. The results obtained with our procedure (Fig. 6b) show a higher predictive performance because 218 sites (69 %) present the same predicted intensity as the observed value and, for 288 sites (91 %), the differences are within 1 intensity degree. For the 1980 Irpinia earthquake, Fig. 7 shows that the intensity values predicted with the IPE alone are equal to the observed ones for 652 out of the 1202 considered sites (54 %), whereas using the proposed methodology these sites become 822 (68 %). A difference of 1 intensity degree between the predicted values and the observed ones is shown at 478 sites (40 %) with the IPE alone (Fig. 7a), whereas it is shown at 350 sites (29 %) with our procedure (Fig. 7b).

This test demonstrates that the intensity values obtained by means of the proposed procedure better reproduce the observed intensities than using the IPE alone. In fact, more than 90 % of differences between predicted and observed intensity values are within 1 intensity degree, which is the uncertainty associated with any macroseismic intensity assessment.

Figure 8Intensity distribution of the 1542 Mugello earthquake assessed by Guidoboni et al. (2007) and reported in DBMI15.

Figure 9Modal values of the probability distribution pl(Is|Iv) computed at the 968 considered localities (small dots) for the 1542 Mugello earthquake; colored circles bound areas of different intensity values predicted by the IPE.

Figure 10Intensity distribution of the 1542 Mugello earthquake at the 968 considered localities represented as the probability of being greater than or equal to intensity 6 with (a) the IPE alone and (b) this procedure and of being greater than or equal to intensity 8 with (c) the IPE alone and (d) this procedure.

The results shown in Figs. 9 and 10 demonstrate the impact of the proposed methodology in reproducing the pattern of observed intensities with respect to the simple isotropic decay of intensity with distance predicted by IPEs.

5 Case study

How this procedure may serve the purpose of reconstructing the macroseismic fields of past earthquakes, especially those with scattered IDPs, is shown by means of the case study of an earthquake that occurred on 13 June 1542 in the Mugello area (northern Apennines), with Mw 6 and an epicentral intensity equal to 9 according to CPTI15. In DBMI15 there are 45 IDPs with a maximum intensity equal to 9, as assessed by Guidoboni et al. (2007). As reported in Fig. 8, the effects of this earthquake are primarily known in the epicentral area with 31 localities with an intensity greater than or equal to 8, whereas the macroseismic information at the localities far from the epicenter is extremely scattered.

With the aim of integrating the intensity distribution of this earthquake, 968 localities of DBMI15 within a radius of 20 km from each of the 45 IDPs were considered. Figure 9 shows the modal values of the probability distribution pl(Is|Iv) computed at each of the 968 localities assuming as prior distribution the probability derived through the IPE (Pasolini et al., 2008; Lolli et al., 2019) and using the intensities observed at all the localities within 20 km (analysis c). Such values are compared with the intensities (expressed as modal values) predicted by the IPE alone, represented with colored circles that bound areas of different intensity values. Figure 9 shows that the intensity values estimated by the two approaches are quite different, particularly in the epicentral area. For example, focusing on the area where the IPE alone predicts intensity 7, the intensities computed by the proposed procedure are equal to 6, 7 and 8. On the contrary, moving away from the epicentral area, the two approaches provide similar results for intensity 5.

Figure 10 displays the geographical distribution of the predicted intensities at the 968 localities represented as the probability of being greater than or equal to intensity 6 and 8, computed through (i) the IPE alone and (ii) the proposed procedure. As shown in Fig. 10a and b, the probability of an intensity greater than or equal to 6 is more than 90 % for 278 localities using the IPE alone, while the same probability extends to 488 localities using the second procedure. The differences between the two approaches become more evident in the case of localities where the probability of an intensity greater than or equal to 8 is higher than 50 %, that is 96 localities using the IPE alone (Fig. 10c) and 212 using our procedure (Fig. 10d).

6 Conclusions

The procedure proposed in this article estimates the probability distribution for a given intensity value at the considered site through a Bayesian approach. The procedure takes into account (i) region-dependent empirical relations to model macroseismic intensity attenuation with source distance (i.e., IPEs), (ii) probability distributions resulting from the in-depth analysis of the spatial variability in intensity data collected in the Italian Macroseismic Database DBMI15 (Sect. 3.2), and (iii) the discrete and ordinal nature of macroseismic intensity and its uncertainties. This procedure allows improvement of the macroseismic intensity distributions of historical earthquakes constraining the intensity values calculated at a site through an IPE with intensities observed at neighboring localities for the same earthquake.

The results obtained in the application part (see Sect. 4.2) emphasize that this method reproduces well the observed values for intensities greater than or equal to 6 and equal to 3. On the other hand, the outcomes for intensity 4 and 5 show an overestimation and underestimation, respectively, that could be linked to both (i) the incompleteness of the analyzed dataset due to the input threshold of DBMI15 (intensity ≥5) and (ii) the incompleteness of historical documentation for lower intensity degrees. These outcomes demonstrate that the intensities predicted with the proposed procedure match the observed values better than those obtained using the IPE alone.

This procedure is thought to integrate incomplete and scattered intensity distributions while avoiding the isotropic decay of intensity with distance resulting from existing IPEs. Through a more realistic modeling of the pattern of predicted intensities, this procedure takes into account the spatial distribution and variability of observed intensity data to constrain the results. Not unexpectedly, the obtained results are dependent on the spatial distribution of the data observed for the selected earthquake and on the number of intensity values available in nearby localities.

The proposed procedure aims at the integration and enrichment of both the intensity distributions of individual earthquakes and the seismic history of single localities. Together with suggestions to further document the spatial distribution and severity of effects in the framework of historical seismological investigation, the outcomes provided by this procedure can be used for local seismic hazard assessment, as well as for planning activities aimed at risk mitigation.

Data availability

DBMI15 is available at (Locati et al., 2019).

CPTI15 is available at (Rovida et al., 2019).

Author contributions

AA edited most parts of the paper, performed the statistical analyses and tested the results. AR, VD'A and DA contributed to the manuscript and supervised the research.

Competing interests

The authors declare that they have no conflict of interest.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors wish to thank Vasiliki Kouskouna and Paolo Gasperini for their thorough reviews and to thank the editor Filippos Vallianatos. The authors are grateful to Paola Albini for her valuable suggestions and comments on the text and to Rodolfo Puglia for his support in writing MATLAB codes.

Review statement

This paper was edited by Filippos Vallianatos and reviewed by Paolo Gasperini and Vasiliki Kouskouna.


Albarello, D. and D'Amico, V.: Attenuation relationship of macroseismic intensity in Italy for probabilistic seismic hazard assessment, B. Geofis. Teor. Appl., 45, 271–284, 2004. 

Albarello, D. and D'Amico, V.: Validation of intensity attenuation relationships, B. Seismol. Soc. Am., 95, 719–724,, 2005. 

Albarello, D. and Mucciarelli, M.: Seismic hazard estimates from ill-defined macroseismic data at a site, Pure Appl. Geophys., 159, 1289–1304,, 2002. 

Albarello, D., Camassi, R., and Rebez, A.: Detection of space and time heterogeneity in the completeness of a seismic catalog by a statistical approach: An Application to the Italian Area, B. Seismol. Soc. Am., 91, 1694–1703,, 2001. 

Albarello, D., D'Amico, V., Gasperini, P., Pettenati, F., Rotondi, R., and Zonno G.: Nuova formulazione delle procedure per la stima dell'intensità macrosismica da dati epicentrali o da risentimenti in zone vicine, Progetto DPC-INGV S1, available at: (last access: 12 April 2021), 2007. 

Albini, P.: Documenting Earthquakes in the United States of the Ionian Islands, 1815–1864, Seismol. Res. Lett., 91, 2554–2562,, 2020a. 

Albini, P.: Nine Major Earthquakes in the United States of the Ionian Islands, 1815–1864, Seismol. Res. Lett., 91, 3595–3611,, 2020b. 

Albini, P. and Rovida, A.: Earthquakes in southern Dalmatia and coastal Montenegro before the large 6 April 1667 event, J. Seismol., 22, 721–754,, 2018. 

Allen, T. I., Wald, D. J., and Worden, C. B.: Intensity attenuation for active crustal regions, J. Seismol., 16, 409–433,, 2012. 

Ambraseys, N. N. and Douglas, J.: Magnitude calibration of north Indian earthquakes. Geophys. J. Int., 159, 165–206,, 2004. 

Bakun, W. H. and Wentworth, C. M.: Estimating earthquake location and magnitude from seismic intensity data, B. Seismol. Soc. Am., 87, 1502–1521, 1997. 

Barbano, M. S., Cosentino, M., Lombardo, G., and Patané, G.: Isoseismal maps of Calabria and Sicily earthquakes (Southern Italy), CNR-PFG, Catania, 116 pp., 1980. 

Carletti, F. and Gasperini, P.: Lateral variations of seismic intensity attenuation in Italy, Geophys. J. Int., 155, 839–856,, 2003. 

D'Amico, V. and Albarello, D.: SASHA: A computer program to assess seismic hazard from intensity data, Seismol. Res. Lett., 79, 663–671,, 2008. 

Fäh, D., Giardini, D., Kästli, P., Deichmann, N., Gisler, M., Schwarz-Zanetti, G., Alvarez-Rubio, S., Sellami, S., Edwards, B., Allmann, B., Bethmann, F., Wössner, J., Gassner-Stamm, G., Fritsche, S., and Eberhard, D.: ECOS-09 Earthquake Catalogue of Switzerland Release 2011 report and database, Swiss Seismological Service, ETH Zurich, Report SED/RISK/R/001/20110417, 2011. 

Ferrari, G. and Postpischl, D.: The Mugello earthquake of June 29, 1919, in: Atlas of isoseismal maps of Italian earthquakes, edited by: Postpischl, D., Quaderni della Ricerca Scientifica, 114, 124–125, 1985. 

Galli, P. and Camassi, R.: Rapporto sugli effetti del terremoto aquilano del 6 aprile 2009, Rapporto tecnico QUEST, DPC-INGV, Roma, 12 pp.,, 2009. 

Gasperini, P, Bernardini, F, Valensise, G., and Boschi, E.: Defining seismogenic sources from historical earthquake felt reports, B. Seismol. Soc. Am., 89, 94–110, 1999. 

Gasperini, P., Vannucci, G., Tripone, D., and Boschi, E.: The Location and Sizing of Historical Earthquakes Using the Attenuation of Macroseismic Intensity with Distance, B. Seismol. Soc. Am., 100, 2035–2066,, 2010. 

Grünthal, G.: European Macroseismic Scale 1998 (EMS-98), Cahiers du Centre Européen de Géodynamique et de Séismologie, 13, 1–99, 1998. 

Guidoboni, E. and Stucchi, M.: The contribution of historical records of earthquakes to the evaluation of seismic hazard, Ann. Geophys.-Italy, 36, 201–215,, 1993. 

Guidoboni, E., Ferrari, G., Mariotti, D., Comastri, A., Tarabusi, G., and Valensise, G.: CFTI4Med, Catalogue of Strong Earthquakes in Italy (461 B. C.–1997) and Mediterranean Area (760 B.C.–1500), INGV-SGA, available at: (last access: 12 April 2021), 2007. 

Hough, S. E. and Martin, S. M.: Which Earthquake Accounts Matter?, Seismol. Res. Lett., 92, 1069–1084,, 2021. 

Locati, M., Camassi, R., Rovida, A., Ercolani, E., Bernardini, F., Castelli, V., Caracciolo, C. H., Tertulliani, A., Rossi, A., Azzaro, R., D'Amico, S., Conte, S., Rocchetti, E., and Antonucci, A.: Database Macrosismico Italiano (DBMI15), versione 2.0, Istituto Nazionale di Geofisica e Vulcanologia (INGV) [data set],, 2019. 

Lolli, B., Pasolini, C., Gasperini, P., and Vannucci, G.: Prodotto 4.8: Ricalibrazione dell'equazione di previsione di Pasolini et al. (2008), in: Il modello di pericolosità sismica MPS19, rapporto finale, edited by: Meletti, C. and Marzocchi, W., Centro Pericolosità Sismica, Istituto Nazionale di Geofisica e Vulcanologia, maggio 2019, Roma, 168 pp., 2 App., 2019. 

Magri, L., Mucciarelli, M., and Albarello, D.: Estimates of site seismicity rates using ill-defined macroseismic data, Pure Appl. Geophys., 143, 618–632,, 1994. 

Manchuel, K., Traversa, P., Baumont, D., Cara, M., Nayman, E., and Durouchoux, C.: The French seismic CATalogue (FCAT-17), B. Earthq. Eng., 16, 2227–2251,, 2018. 

Medvedev, S., Sponheuer, W., and Karník, V.: Neue seismische Skala Intensity scale of earthquakes, in: 7. Tagung der Europäischen Seismologischen Kommission vom 24.9. bis 30.9.1962, in Jena, Veröff. Institut für Bodendynamik und Erdbebenforschung in Jena, Deutsche Akademie der Wissenschaften zu Berlin, 77, 69–76, 1964. 

Meletti, C., Marzocchi, W., D'Amico, V., Lanzano, G., Luzi, L., Martinelli, F., Pace, B., Rovida, A., Taroni, M., Visini, F., and the MPS19 Working Group.: The new Italian seismic hazard model (MPS19), Ann. Geophys.-Italy, 64, SE112,, 2021. 

Mucciarelli, M., Albarello, D., and D'Amico, V.: Comparison of probabilistic seismic hazard estimates in Italy, B. Seismol. Soc. Am., 98, 2652–2664,, 2008. 

Musson, R. M. and Cecić, I.: Intensity and Intensity Scales, in: New Manual of Seismological Observatory Practice 2 (NMSOP-2), edited by: Bormann, P., Deutsches GeoForschungsZentrum GFZ, Potsdam, 1–41,, 2012. 

Olea, R. A.: Geostatistics for engineers and earth scientists, Kluwer Academic Publishers, Dordrecht,, 1999. 

Pasolini, C., Albarello, D., Gasperini, P., D'Amico, V., and Lolli, B.: The attenuation of seismic intensity in Italy, Part II: Modeling and validation, B. Seismol. Soc. Am., 98, 692–708,, 2008. 

Patané, G. and Imposa, S.: The Linera earthquake of May 8, 1914, in: Atlas of isoseismal maps of Italian earthquakes, edited by: Postpischl, D., Quaderni della Ricerca Scientifica, 114, 120–121, 1985. 

Postpischl, D., Branno, A., Esposito, E., Ferrari, G., Marturano, A., Porfido, S., Rinaldis, V., and Stucchi M.: The Irpinia earthquake of November 23, 1980, in: Atlas of isoseismal maps of Italian earthquakes, edited by: Postpischl, D., Quaderni della Ricerca Scientifica, 114, 52–59, 1985. 

Provost, L. and Scotti, O.: QUake-MD: Open-Source Code to Quantify Uncertainties in Magnitude–Depth Estimates of Earthquakes from Macroseismic Intensities, Seismol. Res. Lett., 91, 2520–2530,, 2020. 

Rey, J., Beauval, C., and Douglas, J.: Do French macroseismic intensity observations agree with expectations from the European Seismic Hazard Model 2013?, J. Seismol., 22, 589–604,, 2018. 

Rotondi, R., Varini, E., and Brambilla, C.: Probabilistic modelling of macroseismic attenuation and forecast of damage scenarios, B. Earthq. Eng., 14, 1777–1796,, 2016. 

Rovida, A., Locati, M., Antonucci, A., and Camassi, R.: Italian Archive of Historical Earthquake Data (ASMI), Istituto Nazionale di Geofisica e Vulcanologia (INGV) [data set],, 2017. 

Rovida, A., Locati, M., Camassi, R., Lolli, B., and Gasperini, P.: Catalogo Parametrico dei Terremoti Italiani (CPTI15), versione 2.0, Istituto Nazionale di Geofisica e Vulcanologia (INGV) [data set],, 2019. 

Rovida, A., Locati, M., Camassi, R., Lolli, B., and Gasperini, P.: The Italian earthquake catalogue CPTI15, B. Earthq. Eng, 18, 2953–2984,, 2020. 

Shebalin, N. V.: Atlas of isoseismal maps, III, UNDP-UNESCO Survey of the seismicity of the Balkan region, Skopje, 275 pp., 1974.  

Sieberg, A.: Geologie der Erdbeben. Handbuch der Geophysik, Gebrüder Bornträger, Berlin, 2, 550–555, 1932. 

Sørensen, M. B., Stromeyer, D., and Grünthal G.: Attenuation of macroseismic intensity: a new relation for the Marmara Sea region, northwest Turkey, B. Seismol. Soc. Am., 99, 538–553,, 2009. 

Stirling, M. and Petersen, M.: Comparison of the historical record of earthquake hazard with seismic hazard models for New Zealand and continental Unites States, B. Seismol. Soc. Am., 96, 1978–1994,, 2006. 

Stucchi, M., Albini, P., Mirto, M., and Rebez, A.: Assessing the completeness of Italian historical earthquake data, Ann. Geophys.-Italy, 47, 2–3,, 2004. 

Stucchi, M., Meletti, C., Montaldo, V., Crowley, H., Calvi, G. M., and Boschi, E.: Seismic hazard assessment (2003–2009) for the Italian building code, B. Seismol. Soc. Am., 101, 1885–1911,, 2011. 

Stucchi, M., Rovida, A., Gomez Capera, A.A, Alexandre, P., Camelbeeck, T., Demircioglu, M. B., Gasperini, P., Kouskouna, V., Musson, R. M. W., Radulian, M., Sesetyan, K., Vilanova, S., Baumont, D., Bungum, H., Fäh, D., Lenhardt, W., Makropoulos, K., Martinez Solares, J. M., Scotti, O., Živcic, M., Albini, P., Batllo, J., Papaioannou, C., Tatevossian, R., Locati, M., Meletti, C., Viganò D., and Giardini, D.: The SHARE European Earthquake Catalogue (SHEEC) 1000–1899, J. Seismol., 17, 523–544,, 2013. 

Swiss Seismological Service: ECOS – Earthquake Catalog of Switzerland ECOS report to PEGASOS, version 31.03.2002 SED, Zürich, 2002. 

Woessner, J., Danciu, L., Giardini, D., Crowley, H., Cotton, F., Grünthal, G., Valensise, G., Arvidsson, R., Basili, R., Demircioglu, M. B., Hiemer, S., Meletti, C., Musson, R., Rovida, A., Sesetyan, K., and Stucchi, M.: The 2013 European Seismic hazard model: key components and results, B. Earthq. Eng., 13, 3553–3596,, 2015. 

Short summary
We present a probabilistic approach for integrating incomplete intensity distributions by means of the Bayesian combination of estimates provided by intensity prediction equations (IPEs) and data documented at nearby localities, accounting for the relevant uncertainties. The performance of the proposed methodology is tested at 28 Italian localities with long and rich seismic histories and for the strong 1980 and 2009 earthquakes in Italy. An application of this approach is also illustrated.
Final-revised paper