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

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
(

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,

To estimate the probability

In order to evaluate the optimal distance threshold to characterize

Number of neighboring localities within 10

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

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

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

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.

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

The conditional probability

Two different analyses were performed in order to estimate the frequency
distribution

Values of

Relative frequency of

In both analyses, we then verified the impact on the results of the distance
among localities. Figure 4 shows the effect of distance on

Relative frequency of

To test the effectiveness of this procedure, the probability distributions in
Eq. (

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

List of the 28 test localities with their maximum intensity (

The probability

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

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

using as

Observed (

Table 3 shows that, for analysis c, the differences between observed and
predicted values are less than 10

Observed (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

Differences 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

Differences between the observed intensity values, as
reported in DBMI15 (Guidoboni et al., 2007), for the 1980 Irpinia earthquake
and the intensity values computed with

For the 2009 L'Aquila earthquake, Fig. 6a shows that for 32 out of 315 sites
(10

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

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

Modal values of the probability distribution

Intensity 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

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

With the aim of integrating the intensity distribution of this earthquake, 968
localities of DBMI15 within a radius of 20

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

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

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

CPTI15 is available at

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.