the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

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

### Andrea Antonucci

### Andrea Rovida

### Vera D'Amico

### 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.

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).

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 *l*th event,
the discrete probability density distribution *p*_{l}(*I*_{s}|*I*_{v}) is
computed to associate with each possible intensity value *I*_{s} at the site
*s* a probability value conditioned by the occurrence of effects of intensity
*I*_{v} at any other site *v*:

Here *p*_{l}(*I*_{s}) represents the “prior” probability density which is
deduced from an IPE by using the epicentral parameters (location, epicentral
intensity, magnitude, etc.) of the *l*th event. In general, the most
common IPEs (in their probabilistic formulation) have the form

(Albarello and D'Amico, 2004), where the average *μ* is a function of the
epicentral distance *D* and the intensity at the epicenter *I*_{e} (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

where *g* is the probability distribution which expresses the uncertainty
affecting the epicentral intensity and *I*_{max} is the upper bound of the adopted macroseismic scale (e.g., 12 for the MCS scale). The probability density *p*_{l}(*I*_{s}) can be computed in the form

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*(*I*_{v}|*I*_{s}), 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*(*I*_{v}|*I*_{s}) expresses the constraining
power of *I*_{s} on *I*_{v}. According to Albarello et al. (2007), *q*(*I*_{v}|*I*_{s}) 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.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 (https://emidius.mi.ingv.it/ASMI/index_en.php, 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, https://emidius.mi.ingv.it/CPTI15-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 (https://emidius.mi.ingv.it/CPTI15-DBMI15/query_eq/, 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, $\mathrm{6}\phantom{\rule{0.25em}{0ex}}\mathrm{1}/\mathrm{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
*I*^{′}–*I*^{′′}, e.g., 6–7). In the case intensity *I*_{v} is uncertain between two
contiguous values ${I}_{\mathrm{v}}^{\prime}$ and ${I}_{\mathrm{v}}^{\prime \prime}$ (e.g., *I*_{v}=6–7),
Eq. (1) becomes

where an equal probability is assigned to the hypotheses ${I}_{\mathrm{v}}={I}^{\prime}$ and ${I}_{\mathrm{v}}={I}^{\prime \prime}$. This is the way uncertain intensity values contained in DBMI15 are treated in the following analysis.

## 3.2 Results

To estimate the probability *q*(*I*_{v}|*I*_{s}) 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\left({I}_{\mathrm{v}}\right|{I}_{\mathrm{s}})\approx q({I}_{\mathrm{v}})$; i.e., *I*_{s} becomes
non-informative about *I*_{v}. The closer the sites considered are, the higher
the informative power of *I*_{s} on *I*_{v} 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*(*I*_{v}|*I*_{s}) 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*(*I*_{v}|*I*_{s}), 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*(*I*_{v}|*I*_{s}).

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.

The conditional probability *q*(*I*_{v}|*I*_{s}) in Eq. (1) was estimated from the relative frequencies of the differences between *I*_{v} and *I*_{s}
computed for each of the 546 selected earthquakes. In this analysis *I*_{v}
and *I*_{s} represent any pair of intensity values observed at neighboring
localities (i.e., within a distance of 20 km). If the intensity values
*I*_{s} and *I*_{v} 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 *I*_{s} and *I*_{v}
are certain values (e.g., 7), the difference is counted four times
(Albarello et al., 2007). In general, one has

where $\mathrm{\Delta}I={I}_{\mathrm{v}}-{I}_{\mathrm{s}}$ and the dependence of both *I*_{v} and
*I*_{s} is due to the lack of defined metrics for intensity degrees. As a
preliminary step, we assume that

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

In both analyses, we then verified the impact on the results of the distance
among localities. Figure 4 shows the effect of distance on *q*(Δ*I*). 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.

## 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
*p*_{l}(*I*_{s}|*I*_{v}), derived at each *j*th site for each
*l*th earthquake, were used to calculate the predicted number of
occurrences for each intensity degree *I*_{s} (*N*_{pred}) over the total of the *M* sites and *N* earthquakes considered:

The predicted values (*N*_{pred}) can be compared with the observed number of occurrences (*N*_{obs}) for the same intensity degree *I*_{s}. If the observed intensity value is certain (e.g., 6), it can be expressed as

In the case the *I*_{s} 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:

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

Equation (10) can be used to evaluate the statistical significance of
the discrepancy between predicted (*N*_{pred}) and observed
(*N*_{obs}) values. According to Albarello and D'Amico (2005), when
$\left|Z\right|$ 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 *p*_{l}(*I*_{s}|*I*_{v}) with Eq. (1) using the probability
distribution *q*(Δ*I*) 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 *I*_{s} for all sites and earthquakes through
Eq. (8) and compared these estimations with the observed occurrences (Eq. 10).

The probability *p*_{l}(*I*_{s}|*I*_{v}) has been computed with three different
analyses for the prior distribution *p*_{l}(*I*_{s}) and for *q*(Δ*I*)

- a.
using a uniform distribution over the intensity range 2–11 for

*p*_{l}(*I*_{s}) and the intensity observed at the nearest locality within 20 km (i.e., probability*q*(Δ*I*)_{near}in Table 1); - b.
using a uniform distribution over the intensity range 2–11 for

*p*_{l}(*I*_{s}) and iteratively considering in Eq. (1) the intensities observed at all the localities within 20 km (i.e., probability*q*(Δ*I*)_{all}in Table 1); - c.
using as

*p*_{l}(*I*_{s}) the probability computed through an IPE and the intensities observed at all the localities within 20 km (i.e., probability*q*(Δ*I*)_{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 (*N*_{obs})
for a given intensity value *I*_{s} was compared with the predicted number
(*N*_{pred}) derived for the three possible choices (a, b, c) of the
prior distribution *p*_{l}(*I*_{s}) and *q*(Δ*I*) for the 28 selected
localities (Table 3). The differences in percentage between predicted and
observed values were computed and expressed as
$\left[\right(\mathrm{1}-{N}_{\text{pred}}/{N}_{\text{obs}})\times \mathrm{100}]$. *Z* represents the
standardized Gaussian statistics: if $\left|Z\right|>\mathrm{2}$, the resulting discrepancy can
be considered statistically significant (probability 0.05).

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 (*N*_{pred}) is
consistent with the number of observed occurrences (*N*_{obs}) 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).

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 *M*_{w} 6.3
event that occurred on 6 April 2009 in the L'Aquila area (central Italy) and the
*M*_{w} 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.

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.

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.

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 *M*_{w} 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 *p*_{l}(*I*_{s}|*I*_{v}) 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).

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.

DBMI15 is available at https://doi.org/10.13127/DBMI/DBMI15.2 (Locati et al., 2019).

CPTI15 is available at https://doi.org/10.13127/CPTI/CPTI15.2 (Rovida et al., 2019).

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.

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.

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, https://doi.org/10.1785/0120040062, 2005.

Albarello, D. and Mucciarelli, M.: Seismic hazard estimates from ill-defined macroseismic data at a site, Pure Appl. Geophys., 159, 1289–1304, https://doi.org/10.1007/s00024-002-8682-2, 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, https://doi.org/10.1785/0120000058, 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: http://esse1.mi.ingv.it/d10.html (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, https://doi.org/10.1785/0220200056, 2020a.

Albini, P.: Nine Major Earthquakes in the United States of the Ionian Islands, 1815–1864, Seismol. Res. Lett., 91, 3595–3611, https://doi.org/10.1785/0220200205, 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, https://doi.org/10.1007/s10950-018-9730-4, 2018.

Allen, T. I., Wald, D. J., and Worden, C. B.: Intensity attenuation for active crustal regions, J. Seismol., 16, 409–433, https://doi.org/10.1007/s10950-012-9278-7, 2012.

Ambraseys, N. N. and Douglas, J.: Magnitude calibration of north Indian earthquakes. Geophys. J. Int., 159, 165–206, https://doi.org/10.1111/j.1365-246X.2004.02323.x, 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, https://doi.org/10.1111/j.1365-246X.2003.02073.x, 2003.

D'Amico, V. and Albarello, D.: SASHA: A computer program to assess seismic hazard from intensity data, Seismol. Res. Lett., 79, 663–671, https://doi.org/10.1785/gssrl.79.5.663, 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., https://doi.org/10.13127/QUEST/20090406, 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, https://doi.org/10.1785/0120090330, 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, https://doi.org/10.4401/ag-4264, 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: http://storing.ingv.it/cfti4med/ (last access: 12 April 2021), 2007.

Hough, S. E. and Martin, S. M.: Which Earthquake Accounts Matter?, Seismol. Res. Lett., 92, 1069–1084, https://doi.org/10.1785/0220200366, 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], https://doi.org/10.13127/DBMI/DBMI15.2, 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, https://doi.org/10.1007/BF00879501, 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, https://doi.org/10.1007/s10518-017-0236-1, 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, https://doi.org/10.4401/ag-8579, 2021.

Mucciarelli, M., Albarello, D., and D'Amico, V.: Comparison of probabilistic seismic hazard estimates in Italy, B. Seismol. Soc. Am., 98, 2652–2664, https://doi.org/10.1785/0120080077, 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, https://doi.org/10.2312/GFZ.NMSOP-2_ch12, 2012.

Olea, R. A.: Geostatistics for engineers and earth scientists, Kluwer Academic Publishers, Dordrecht, https://doi.org/10.1007/978-1-4615-5001-3, 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, https://doi.org/10.1785/0120070021, 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, https://doi.org/10.1785/0220200064, 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, https://doi.org/10.1007/s10950-017-9724-7, 2018.

Rotondi, R., Varini, E., and Brambilla, C.: Probabilistic modelling of macroseismic attenuation and forecast of damage scenarios, B. Earthq. Eng., 14, 1777–1796, https://doi.org/10.1007/s10518-015-9781-7, 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], https://doi.org/10.13127/asmi, 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], https://doi.org/10.13127/CPTI/CPTI15.2, 2019.

Rovida, A., Locati, M., Camassi, R., Lolli, B., and Gasperini, P.: The Italian earthquake catalogue CPTI15, B. Earthq. Eng, 18, 2953–2984, https://doi.org/10.1007/s10518-020-00818-y, 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, https://doi.org/10.1785/0120080299, 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, https://doi.org/10.1785/0120050176, 2006.

Stucchi, M., Albini, P., Mirto, M., and Rebez, A.: Assessing the completeness of Italian historical earthquake data, Ann. Geophys.-Italy, 47, 2–3, https://doi.org/10.4401/ag-3330, 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, https://doi.org/10.1785/0120100130, 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, https://doi.org/10.1007/s10950-012-9335-2, 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, https://doi.org/10.1007/s10518-015-9795-1, 2015.