Macroseismic intensity attenuation models calibrated in Mw for Italy

This study aims at developing new macroseismic intensity attenuation models valid for Italy by exploiting the most updated macroseismic dataset and earthquakes catalogue, as well as the information obtained from a critical analysis of the most recent models in the literature. Several different attenuation models have been calibrated as a function of the moment magnitude (Mw) and epicentral distance from 16.260 intensity data points, that are related to 119 earthquakes occurred after 1900. According to trends and residuals analysis, the preferred calibrated intensity ( I ) attenuation function is a Log-Linear model for epicentral distance ( r in km) and a linear model for Mw as:


La Introduction
Despite the progress achieved in the last century by instrumental seismology, macroseismic intensity data continue to be fundamental for the assessment of seismic hazard and risk.The macroseismic intensity classifies, through the use of a given macroseismic scale (e.g., MCS -Sieberg, 1930;MMI -Wood and Neumann, 1931;MSK -Medvedev et al., 1964;EMS-98 -Grünthal, 1998), the severity of the effects of the earthquake shaking produced in a limited area, usually called "locality", on humans, and on the built and natural environment.For example, the European Macroseismic Scale EMS-98 (Grünthal, 1998) is compiled with theoretical descriptive frameworks, arranged in a hierarchical order.This scale assigns macroseismic intensity by considering the typology and vulnerability of damaged buildings, the relative level of damage, and comparing that data to a predefined table.
As a consequence, the macroseismic intensity is a tool in the seismological and engineering practice for classifying the degree of damage that an earthquake may cause and is a parameter that could be used to estimate expected ground shaking.Furthermore, the spatial distribution of Intensity Data Points (IDPs) is used for the characterization of the seismic source (i.e., estimates of epicentral location and magnitude) of pre-instrumental earthquakes (e.g., Bakun and Wentworth, 1997;Gasperini et al., 1999Gasperini et al., , 2010;;Provost and Scotti, 2020), that constitute the bulk of seismic catalogues in countries where the historical record is much longer than the instrumental one.In Italy, much of the information on earthquakes available has been compiled from historicalmacroseismic studies.Recent years have seen an increase in historical research, which has been incorporated into macroseismic databases (Monachesi and Stucchi, 1997;Boschi et al., 2000;Guidoboni et al., 2007;Stucchi et al., 2007;Locati et al., 2022), improving the understanding of seismic activity in the country.
Macroseismic intensity attenuation is the rate at which shaking amplitude decreases with distance from the epicentre (Musson and Cecic, 2002).Kövesligethy (1906) proposed the first intensity attenuation model, which describes the decrease in seismic energy due to absorption of the geophysical media and geometrical spreading.
The model is represented by the difference between epicentral intensity (Io) and site intensity (I) as a function of hypocentral distance (R), focal depth (h), and a free parameter.Blake (1937) simplified the model by eliminating the linear term while keeping the coefficient of the logarithm as a free parameter.However, empirical observations and macroseismic data have revealed issues with this model, such as the use of Io as a measure of earthquake size.
The accuracy of earthquake intensity predictions has been improved over the years by refining Kövesligethy's model.Refinements have included additional simplifications and considering magnitude as an earthquake size term.The electronic supplement (Table A1) shows a summary of such models grouped into families of functional types.
The models proposed in literature, mainly after 1990, are either logarithmic or non-logarithmic (linear, non linear, polynomial, etc.).These models are used to estimate the correlation between macroseismic data in isoseismal maps and seismic energy density (Howell and Schultz, 1975).Several linear, log and log-linear models have been used in different places, such as California, Northern Rhine area, Iberian Peninsula, France, Portugal, Caribbean, Andean region, Brazil, Japan and China among others (Table A2) with different spatial and source size ranges and resulting intensities.
In this study new macroseismic intensity attenuation models have been derived for Italy by exploiting the most up-to-date and publicly available datasets, including the Italian Macroseismic Database (DBMI15, Locati et al., 2022) and the Parametric Catalogue of Italian Earthquakes (CPTI15, Rovida et al., 2022a).The models were developed in the framework of a research project, leaded by the Seismic Hazard Center (Centro Pericolosità Sismica, CPS) of the Istituto Nazionale di Geofisica e Vulcanologia (INGV) and supported by the Italian Civil Protection Department (Dipartimento della Protezione Civile, DPC), aimed to produce the MPS19 seismic hazard model for Italy (Meletti et al., 2021); one task of the project was focused on the prediction models for strong motion parameters and macroseismic intensity to be used in the hazard calculations.
Leveraging the most recent research and a robust macroseismic dataset, we present a classical Log-linear attenuation model calibrated in Mw, which is the most widely used in the literature, as well as an alternative more unusual model with higher sensitivity to physical parameters such as focal depth and magnitude following Howell and Schultz (1975).To confirm the validity of our models, we also apply a validation process using an independent dataset following Mak et al. (2015) and Bakun and Wentworth (1997).

State-of-the-art in Italy
Several intensity attenuation models have been developed in Italy to provide seismic hazard evaluations in terms of macroseismic intensity or to derive earthquake parameters for pre-instrumental events from IDPs.
The published models (Table 1) have extensively focused on macroseismic intensity attenuation over whole Italy (Gasperini, 2001;Albarello and D'Amico, 2004;Gomez-Capera, 2006;Pasolini et al., 2008b) or more localised areas (e.g., Azzaro et al., 2006 for the Mt.Etna volcanic region).This multiplicity of studies is largely due to the continual improvement of Italian macroseismic databases and earthquake catalogues, as well as to a greater understanding of methodological factors that were previously overlooked.As an example, previous attenuation models proposed for Italy (Peruzza, 1996;2000), and used for seismic hazard analyses in terms of macroseismic intensity (Slejko et al., 1998;Albarello et al., 2000), did not consider aleatory uncertainty in the predicted intensity values.
The more recent models available in literature for the whole Italian territory have been proposed by Gomez-Capera (2006) and Pasolini et al. (2008b).Both models were developed using the same datasets (DBMI04, Stucchi et al, 2007 andCPTI04, CPTI Working Group, 2004) and were used to compute seismic hazard in Italy in terms of macroseismic intensity through different approaches (Gomez-Capera et al., 2010).The two models differ in methodological and computational choices, mainly concerning the functional form of intensity-distance decay and the source-size parameter (Table 1).
In particular, Gomez-Capera (2006) adopts epicentral intensity as the measure of earthquake size and the model assumed that the intensity decay (i.e., difference between epicentral and site intensity) is proportional to the cubic root of the epicentral distance (Berardi et al., 1993), independent of the earthquake's focal depth.Historical and post-1900 earthquakes were used and criteria for excluding distant IDPs were applied.Calibration was conducted using 20.873 IDPs related to 212 earthquakes that occurred from 1279 to 2002; these data were grouped according to the predominant faulting mechanism in the relevant seismic zone (Meletti et al., 2008).This model was developed for the whole Italian territory, as well as for each of the predominant style-of faulting (normal, reverse and strike-slip) and the Etna volcanic zone (Table 1).A residual analysis was reported, along with its accompanying standard deviation, as an indication of the quality of the regressions.For the Whole Italy and Normal models, the fit was good up to about 150 km; for the Strike-Slip and the Reverse model, the error increased for distances greater than 110 km and, for the Etna volcanic area, beyond 10 km.Validation process on an independent dataset was not performed and there was no calibration process with a physical parameter such as moment magnitude.
The model by Pasolini et al. (2008b) considered earthquakes after 1200 with a minimum of 10 IDPs and excluded offshore events and those that occurred in the volcanic areas of Mt.Etna and Ischia Island.To avoid any bias due to incompleteness of low macroseismic intensity values, they also disregarded all IDPs at a distance where an intensity of less than 4 MCS was expected, according to Gasperini (2001).Similarly, Gomez-Capera (2006) applied the same selection criteria, but further omitted deep events, minor earthquakes (Io <7 MCS), and those outside the seismogenic zones of the ZS9 source zone model (Meletti et al., 2008).The dependence of macroseismic intensity on source distance is modelled through a log-linear function.This choice, shared by many studies on intensity attenuation, is based on physical evidence that supports the proportionality of intensity with the logarithm of ground motion amplitude (Pasolini et al., 2008a), thereby justifying the adoption of a relationship similar to common Ground Motion Prediction Equations (GMPEs).Although the two studies adopted different functional forms to model intensity-distance decay and differed in their analysis strategies, their average attenuation curves are very similar, with the exception of the epicentre, where Gomez-Capera (2006) predicts higher values than the relevant epicentral intensity.This difference is attributed to the different source terms used by the two studies rather than to different attenuation patterns: unlike the IE parameter of Pasolini et al. (2008b), which is defined as the intensity expected at the epicentre, the "conventional" epicentral intensity Io is often not consistent with the intensity predicted by the attenuation relationship at R=0 (Table 1; Pasolini et al., 2008a).However, when examining the effect of the two predictive equations on seismic hazard, some discrepancies were observed in certain parts of Italy.D' Amico et al. (2009) suggest that these differences should be attributed more to the different parameterization of the source term (Io for Gomez-Capera, 2006;IE for Pasolini et al., 2008b; Table 1) and its consequent effect on seismicity rate calculations, rather than to the different values of standard deviation.
The models developed by Gomez-Capera et al. (2008a, b;2009) were applied to evaluate earthquake parameters of historical earthquakes from IDPs.The first model (Gomez-Capera et al., 2008a) was proposed for offshore events of the CPTI11 catalogue (Rovida et al., 2011).This model simplified the functional form by eliminating the absorption coefficient while allowing the geometric coefficient to remain a free parameter (Gomez-Capera et al., 2008b; Table 1).Gomez-Capera et al. (2009) calibrated a more recent intensity attenuation model within the context of the EC NERIES Project and the European Earthquake Catalog SHEEC (Stucchi et al., 2013).To derive epicentral parameters from IDPs, this earthquake catalogue used the Boxer method (Gasperini et al., 1999) as well as the Bakun and Wentworth method (1997), both of which had been properly calibrated.While this Italian relationship was validated on an independent dataset, it was disregarded as its accuracy could be further increased with a larger dataset.The European EPICA catalogue (Rovida et al., 2022b) retains exactly the same principles and procedures adopted for the compilation of SHEEC, of which it is an update.

Data
The input dataset considered in this work was extracted from DBMI, the Italian Macroseismic Database (Locati et al., 2016(Locati et al., , 2019(Locati et al., , 2021(Locati et al., , 2022)).The latest version of DBMI is DBMI15 v4.0 and it provides access to 123.981 IDPs, related to 3.229 earthquakes that occurred between the year 1000 and 2020.When data of instrumental origin is missing, the intensity data of DBMI is used to compile the earthquake parameters of CPTI, the Parametric Catalogue of Italian Earthquakes (Rovida et al., 2016(Rovida et al., , 2019(Rovida et al., , 2020(Rovida et al., , 2021(Rovida et al., , 2022a)).CPTI contains homogeneous macroseismic and instrumental data for earthquakes with a maximum macroseismic intensity Imax ≥ 5 MCS or a magnitude Mw ≥ 4.0 that occurred on Italian territory.The latest version of CPTI is CPTI15 v4.0 and it provides access to parameters related to 4894 earthquakes between the year 1000 and 2020.In order to ensure the accuracy of the new macroseismic intensity attenuation model, a carefully selected input calibration dataset has been compiled using data from DBMI and considering CPTI parameters.This dataset satisfies the following general criteria: • earthquakes should cover a wide range of magnitudes; • earthquakes and macroseismic intensity data should have a wide spatial distribution; • magnitude and epicentre should preferably be of instrumental origin, with a preference for Mw; • the set of IDPs associated with each considered earthquake should be as large as possible.
The compilation of the calibration dataset from CPTI15 and DBMI15 with well-defined criteria ensures the reliability of the new macroseismic intensity attenuation model in the seismic hazard frame.
To accurately calibrate a macroseismic intensity attenuation model, a careful selection of shallow tectonic earthquakes and IDPs was carried out according to various criteria (Gomez-Capera, 2006;Gomez-Capera et al., 2010): • earthquakes in the volcanic area of Mount Etna were excluded since the attenuation pattern in this zone is distinct from that of active crustal regions (Ciccotti et al., 2000); • IDPs with intensity I < 3 were discarded to avoid data incompleteness; • events with accumulated effects due to damaging aftershocks were also removed; as well as earthquakes with focal depth > 35 km were discarded; • earthquakes with an azimuthal gap in the macroseismic intensity distribution due to a lack of information from sparsely populated areas and off-shore epicentres were excluded, as were crossnational earthquakes due to the incomplete distribution of IDPs and ill-constrained epicentre location; • earthquakes characterised by low Io and/or small Mw were removed since the study will focus on strong earthquakes; • earthquakes with less than 12 macroseismic observations were removed because events with a low number of IDPs could bias the regression analysis; • macroseismic data associated with special cases identified by DBMI15 (e.g.

TE [large territory], IB [isolated building], SS [small settlement], MS [multiple settlement], DL [deserted locality], AL [absorbed locality], CQ [city quarter]
) should be evaluated carefully, since the statistical nature of macroseismic intensity may not be fulfilled.For this reason, such macroseismic data were excluded; • likewise, special care was taken when dealing with non-conventional macroseismic intensities (e.g.Felt, Damage, etc.).
The criteria are applied to capture the main physical characteristic of the macroseismic data distribution for each earthquake: the magnitude of the earthquake and its corresponding macroseismic intensity attenuation trend with distance.Following the application of the selection criteria, 16.260 IDPs related to 119 earthquakes (Figure 1) that occurred between 1908-2013 with reliable instrumentally recorded magnitude (Mw) are used to calibrate the coefficients of the macroseismic intensity attenuation model.The instrumental parameters of the earthquakes (CPTI15) used for the calibration process are provided in xlsx format in the electronic supplement (Table A4).

Method
This study analyses two attenuation functionals proposed by Howell and Schultz (1975).Such study deduces the pattern of attenuation in macroseismic intensity based on an empirical relationship between macroseismic intensity and energy, and the classical equation for energy decay: where E0 is the total energy, R is the hypocentral distance, b is a constant for geometric spreading, and c is a constant for the rate of absorption.Howell and Schultz (1975) proposed two empirical equations that relate seismic intensity, expressed as macroseismic intensity (I), to seismic energy density (E): i) The first equation is a logarithmic relationship, in which I is proportional to the logarithm of E: ii) The second equation is an alternative power relationship, in which I is proportional to a power of E: These equations were developed under the assumption that seismic energy generated during an earthquake is radiated from a point source into a space of simple geometry, such as a uniform hemisphere or layering only as a function of focal depth.Although the macroseismic intensity is conceived in terms of observable qualitative effects, these equations suggest that the steps of the macroseismic intensity scale approximate an even progression of any kind based on energy.
The state of the art of macroseismic intensity shows an empirical relationship between the intensity (I) and hypocentral distance.Most models are similar to the Logarithmic-Linear model (Log-Lin) and use the intensity at the epicentre (I0) as seen in early macroseismic models (Table A1).The following two empirical macroseismic intensity attenuation models (4) and ( 5) were derived from ( 2) and (3), respectively, using the earthquake size parameter Mw (Howell and Schultz, 1975;Gomez-Capera and Salcedo Hurtado, 2002): i) Macroseismic intensity attenuation model type 1 where an additive Normal error is implicit (Appendix A1).
ii) Macroseismic intensity attenuation model type 2 where an additive Normal error is implicit.Model ( 5) is equivalent to the multiplicative model with a2'=10 a2 and multiplicative LogNormal error (Appendix A1).
The classical functional form of macroseismic intensity attenuation, described by equation ( 4), is widely observed in literature from Kovesligethy (1906) to the most recent models calibrated in Mw (Atkinson et al., 2014;Baumont et al., 2018;Gomez-Capera et al. 2020;Mezcua et al., 2020; Table A2).This form is similar to the one used by GMPEs, taking into account geometric spreading and anelastic attenuation.
Equation ( 5) is an alternative linear functional form for the attenuation of the logarithm of macroseismic intensity, which increases with magnitude (Mw) in logarithmic form.Taking the logarithm of both sides of equation ( 6) is equivalent to equation ( 5).This functional form is less common in the literature, but it has been applied with successful regional results by Howell and Schultz (1975) for San Andreas, Cordillera, and eastern provinces in the United States and Canada; Kaila and Sarkar (1982) for the United States; Chavez and Castro (1988) for subduction zones in the south-central region of Mexico and the Trans-Mexican Volcanic Belt; Greenhalgh et al. (1989) for Australian states; and Gomez-Capera and Salcedo Hurtado (2002) for events with h<60 km in the Colombian Andes.There are a wide range of alternative equations from ( 4) and ( 6) as discussed in Howell and Schultz (1975), Cua et al. (2010), among others.
In the equations ( 4) and ( 6) for modelling macroseismic intensity, the coefficient 'a' is calibrated to reflect boundary conditions of the source; 'b' defines the rate of energy geometric spreading; 'c' is associated with anelastic attenuation; and 'd' represents magnitude dependence.The hypocentral distance 'R' is calculated as the square root of the sum of the squares of the epicentral distance 'r' and the pseudo focal depth 'h' (in km): R=(r 2 + h 2 ) 1/2 .The resulting macroseismic attenuation models (4) and ( 6) are functions of both epicentral intensity r and moment magnitude Mw, namely I=f(r,Mw), where (a1, b1, c1, d1, h1) and (a2', b2, c2, d2, h2) are coefficients to be estimated.
By assuming that all coefficients are positive except a1 which can be any number, the difference between equations ( 4) and ( 6) lies in the curvature of the variation of macroseismic intensity with epicentral distance (r) and moment magnitude Mw.The slope of the macroseismic intensity (I) with respect to r is expressed by the first order derivative, given by equation ( 7) for model ( 4) and ( 8) for model ( 6): Slope ( 7) depends solely on the epicentral distance r, whereas slope (8) depends on moment magnitude Mw as well.This implies that, unlike model ( 4), the curvature of the intensity attenuation model ( 6) changes with Mw.
Nevertheless, ( 7) and ( 8) share some properties.Both slopes are approximately zero very close to the epicentre, where the intensity I=f(r,Mw) is almost constant.As analytical functions, they asymptotically approach a constant value, precisely -c1 for model ( 4) and zero for model ( 6).This indicates that, for very large values of r, the analytical function (4) has a linear decreasing trend and assumes negative values; consequently, in our application to macroseismic data, function (4) must be truncated at the distance where the intensity reaches its minimum value (I=3 in this study).Also the model ( 6) must be analogously truncated because its analytical function tends asymptotically to zero.
As for the derivatives of models ( 4) and ( 6) with respect to Mw, they are respectively given by dI/dMw= d1 (9) dI/dMw= I*d2/Mw (10) These slopes show that the intensity of model ( 4) grows linearly at constant rate d1 as Mw increases, whereas the intensity of model ( 6) grows at a nonlinear rate if d2>1 and decreases otherwise.

Results
This section is devoted to determine which of the two models chosen in this study is the best predictor of macroseismic intensity attenuation through moment magnitude and hypocentral distance (defined by the epicentral distance r and the pseudo focal depth h).To this end, the model coefficients (a, b, c, d, h) of models ( 4) and ( 5), along with their standard errors, are estimated from the input dataset by applying a linear regression procedure based on the Levenberg-Marquardt algorithm (Synergy Software, 2021;Wolfram, 2022).The optimization procedure detects the best-fitting values of the coefficients among all real numbers.The obtained values of the coefficients, as well as their standard error, are summarised in the electronic supplement (Table A3) and the corresponding calibrated relationships in Table 2.The estimated standard deviation (σ) of the residuals provides a measure of accuracy of the model, with lower standard deviation indicating a more accurate model.

The best regressions assuming h as free parameter
The following results have been obtained for the macroseismic intensity attenuation models (4) and ( 5) considering the pseudo focal depth h as a free parameter: • Log-Lin_10 model type 1.
The regression analysis of model type 1 in Eq. ( 4) will be hereinafter denoted by Log-Lin_10 , because the pseudo focal depth h is estimated to be approximately equal to 10 Km (h= 9.87km ±5.7%).This is shown in row N=2 of Table 2 and Table A3, where the estimated coefficients of the model and their standard error are also reported.The input data and the estimated macroseismic intensity model I=f(r,Mw) are provided in Fig. 3a.Let Iobs.and Icomp.denote the observed and the predicted macroseismic intensity at sites, respectively; the histogram of the residuals (Icomp.-Iobs.) is illustrated in Fig. 3b, with estimated standard deviation σ_I=0.748derived from the root mean square of the residuals: The model is further visualised by plotting the predicted macroseismic intensity (Icomp.)and the associated uncertainty (Icomp.±σ)against epicentral distance (r) for Mw classes, in logarithmic scale (Fig. 3c).
The regression analysis of model type 2 in Eq. ( 5) will be hereinafter denoted by CRV9 , being the estimated pseudo focal depth approximately equal to 9 Km (h=8.72km±7.9%).This is shown in row N=5 of Table 2 and ESUPP5, where the estimated coefficients of the model and their standard error are also reported.These estimates also apply to the equivalent model ( 6), remembering however that a2'=10^a2.The input data and the estimated equivalent model ( 6 is given by Eq. ( 4), the macroseismic intensity I is directly proportional to Mw and decreases with the increase of R. Similarly, Fig. 4d shows the model CRV9, given by Eq. ( 6), for different Mw classes; in this case, the macroseismic intensity I is proportional to a power of Mw and decreases as R increases.
Since the estimated coefficients are all positive (Table A3), the remarks on the curvatures given in previous section 4 apply to both estimated models.Fig. 5 compares the slopes (dI/dr) versus the epicentral distance (r) of the models.As already noted, the model CRV9, unlike model Log-Lin_10, has a different slope for each Mw class.
As might be expected, the only slope of Log-Lin_10 exhibits an average behaviour with respect to the slopes of CRV9 for different Mw classes.
In both cases, near the epicentral area, the slopes rapidly and slightly decrease to a minimum (at about 10.1 Km for Log-Lin_10 and 8.1 Km for CRV9) and then just as quickly grow towards a constant value.A change in concavity of both models is also observed (at about 17.4 Km for Log-Lin_10; difficult to calculate exactly for CRV9, but likely at much less than 20 Km), from downward concavity near the epicentral area to upward concavity further away.

Simplification of the functional and the Log Model
Previous studies (Table 1; Table A1; Table A2) have demonstrated that further simplifications can be made to the Log-Linear functional to give more weight to the geometric attenuation.We conducted a one-component logarithmic regression of the hypocentral distance for both type 1 and 2 models.In Fig. 6, the type 1 Log model with only geometric attenuation (N=8 in Table 2) was seen to be comparable to the Log-Lin_10 model (N=2 in Table 2) up to a distance of 300 km, after which the model diverged and reached an epicentral distance of 700 km for macroseismic intensity I=3.In this type 1 Log model (N=8), the pseudo focal depth was modelled as h=16.6 km, with a higher standard deviation (0.751) when compared to the Log-Lin_10 model (0.748).The type 2 model with geometric coefficient only (N=9 in Table 2) is similar to the CRV9 model with both geometric and anelastic attenuation (N=5) up to an epicentral distance of 200 km (Fig. 6).Beyond this, the model diverges up to 1000 km for macroseismic intensity I=3.For this type 2 model with geometric coefficient only, the pseudo focal depth is computed as h=16.2km and its standard deviation (0.735) is slightly higher than that of the CRV9 model (0.731).
It is clear from Fig. 6 that there is an attenuation deficit in the macroseismic far field due to the lack of anelastic coefficient compared to the two best regressions, particularly for strong earthquakes.The importance of the linear attenuation component with distance (anelastic attenuation) is highlighted in the macroseismic far field (x>200 km), complementing the relevance of the geometric attenuation in the macroseismic near field.This suggests that the regressions with only a geometric spreading are not viable in the present study.However, we take h=16 km, an approximation of the estimated pseudo focal depth, as the reference threshold for the following sensitivity tests.

Macroseismic data completeness in the far field and the Log-Lincut dist I=3 model
A trial regression was conducted to reduce incomplete data following Gasperini (2001).In that paper, macroseismic data from earthquakes of the historical period and post-900 were used as input, and models that do not depend on the moment magnitude Mw were analysed (Table 1).Pasolini et al. (2008b) also applied criteria based on the distance to avoid incomplete macroseismic data while using historical period events and Mw in the regression.In this study, these criteria were also used to get a dataset as complete as possible (Gomez-Capera, (N=7 in Table 2), was Log-Linear type 1 with a slightly higher standard deviation (0.771) than Log-Lin_10 (0.748), despite the reduced data set.The pseudo focal depth was modelled as h=11.3km.The calibrated coefficients had similar values, but with greater uncertainties in comparison with the coefficients of the best Log-Linear model Log-Lin_10 (Table A3).
We recall that in this study post-900 earthquakes, macroseismic data and Mw calibrated regressions were analysed.At great epicentral distances, the macroseismic data are essentially due to strong events, as demonstrated by the IDPs between 400 km and 634 km epicentral distance (Fig. 2d) that are from 1915.01.13 Marsica (Mw=6.99), 1930.07.23 Irpinia (Mw=6.62), 1976.05.06 Friuli (Mw=6.45) and 198011.23 Irpinia (Mw=6.81)earthquakes (Table A4; Rovida et al., 2022a;Locati et al., 2022).Thus, this study found no advantage in using an epicentral distance cut-off for modelling macroseismic intensity attenuation calibrated in Mw.The macroseismic intensity attenuation model incorporating both geometric and anelastic coefficients offers the most accurate representation of macroseismic intensity attenuation for a wide range of macroseismic intensity distribution and Mw values, encompassing near and far macroseismic field, as well as moderate and strong earthquakes.
close to those obtained for the calibrated models of post-900 events by Mezcua et al., (2020), Gomez-Capera et al., (2008a), Bakun (2005) and Mezcua et al., (2004) (Table A1; Table A2).The uncertainty of the Mw coefficient was lower than those estimated for the geometric and anelastic attenuation coefficients, as well as for the independent term of all the calibrated models.The calibration coefficient of the magnitude Mw was also stable in the sensitivity analysis.Simplifications of the Log-Lin model, such as modelling only the geometric component or applying a distance cut of data in the far macroseismic field, did not provide any physical or statistical advantage for modelling attenuation (Table A3).
The results show that the Log-Lin_10 model type 1, with h as a free parameter, has the lowest standard errors of the coefficients and reasonably low standard deviation of residuals (σ) compared to the other models.Model type 1 also offers a more precise pseudo focal depth of h=9.87km±5.7%,compared to model type 2 with h=8.72km±7.9%.This suggests that model type 1 provides the most accurate representation of macroseismic intensity attenuation with epicentral distance and moment magnitude.Although the attenuation trend performance of the type 2 model ( CRV9) is not significantly different from that of type 1 (Log-Lin_10) at distances greater than 10km, CRV9 may still be better at predicting the maximum intensity of stronger earthquakes.

Validation Process
The validation process was conducted using independent observations, that is earthquakes and IDPs not considered for models' calibration.To further assess the reliability of the calibrated intensity attenuation models, we analysed the IDPs (DBMI15; Locati et al., 2022) of 15 Italian earthquakes with known instrumentally-derived earthquake parameters (location and magnitude, CPTI15, Rovida et al., 2022a;Table 3;Fig. 9).We selected the events following the criteria outlined in section 4.Among the 15 earthquakes, four events had partial IDP distributions, which we used to evaluate the performance of the calibrated attenuation models: two offshore The performance of the calibrated intensity attenuation models is validated via two prediction strategies: i) mean absolute error (MAEi) as the measure of the goodness of fit between model predictions and observed macroseismic intensities of earthquake "i" (a standard metric for model regression analysis; i.e., refer to Willmott et al., 1985, andMak et al., 2015 for some applications in geophysics); and ii) moment magnitude Mw from IDPs, using the macroseismic intensity attenuation model calibrated in Mw and the Bakun and Wentworth (1997) method.The results demonstrate the reliability of the models.

Mean Absolute Error (MAE)
Mean absolute error (MAE), as its name implies, is a measure of the difference between the forecasted and observed values.MAE was used to assess the accuracy of the macroseismic intensity attenuation model predictions and the observed macroseismic intensities for each validation earthquake.For the macroseismic field of an earthquake with a number P of macroseismic data points, MAE is defined as and resj=Iobs.j-Ipre.j where Iobs,j are the observed intensity values and Ipre,j are the predicted intensity values, for site record j of the considered earthquake.The apex N is the number of macroseismic data points for the seismic event and the wj is the weight for Iobs,j.Mak et al., (2015) did not weight the dataset (wj=1) because in the traditional macroseismic intensity database, as example DBMI15, there is no indicator for the precision of each macroseismic intensity assignment.We calculated the average of all event-based Mean Absolute Error (MAE) scores as the performance indicator for each macroseismic intensity attenuation model calibrated in Moment Magnitude (Mw) in this study.
This indicator was then used to compare with models published in the literature, such as Pasolini et al. (2008b).
The electronic supplement (Fig A1 ) of this paper provides an analysis of the correlation between residuals of the 16.260 macroseismic data points (IDPs) and their input parameters (epicentral distance, magnitude, and macroseismic intensity).The residuals for the six models appear to be random and exhibit no significant trends.
Table 3 presents the 15 earthquakes used for the validation process, along with the obtained mean event-based MAE scores for Log-Lin_10, CRV5, CRV9 and CRV16, as well as Pasolini et al. (2008b).Figures 10 and 11 display the results from Table 3, arranged chronologically and by Mw respectively, for each macroseismic intensity attenuation model.
For the 1909.04.06 Aquilano and 2016.08.24Monti di Laga earthquakes, the worst predictions were given for all attenuation models due to an incomplete distribution of far-field macroseismic data.Conversely, the best prediction for all attenuation models was given for the offshore 2002.09.06 Tirreno meridionale event.The 1919.06.20 Mugello earthquake, which had the highest Mw, maximum intensity (6.26Mw±0.17 and 10 MCS) and 565 IDPs, presented a stable performance for all calibrated models.
The four models calibrated in this study outperformed Pasolini et al. (2008b), which had a median event-based MAE score of 0.9.For comparison, the present study found a median value of all event-based MAE scores of 0.6 intensity units for the models Log-Lin_10, CRV5, CRV9, and CRV16.This result indicates that these models are better at capturing macroseismic intensity for moderate and large earthquakes than Pasolini et al. (2008b).
6.2 Moment magnitude from IDPs 6.2.1 Bakun and Wentworth (1997) method The Bakun and Wentworth (1997) method (BW97) requires a macroseismic intensity attenuation model as a function of both earthquake magnitude and distance.This is shown through equations 4 and 6, which assume a constant focal depth h=h0 km.These equations are inverted to estimate the single-site magnitude Mi (where i is an integer between 1 and P, the total number of IDPs available for the earthquake under examination) from the individual intensity values Ii observed at distances Ri: We indicate with MIi the earthquake magnitude estimated from macroseismic intensity (Ii).The location and magnitude of a given earthquake are estimated by computing the magnitude MI k i over a grid of trial source locations xk (i.e., grid nodes).For the earthquake studied, the magnitude MI k is defined as the average of the magnitudes MI k i estimated from individual macroseismic observations (IDPs) and assuming the epicentre located in the grid node xk, that is Then, considering a grid of trial source locations xk, the root mean square rms(MI k ) is computed as (16) where (17) rms0[MI k -Mi] is the minimum of all rms[MI k -Mi] over the grid of assumed trial source locations, and wi the distance weighting function (Bakun and Wentworth, 1997) The minimum of all the rms over the grid of assumed trial source locations is subtracted from each of the rms of the grid node xk.The grid node xk corresponding to the minimum rms[MI k ] is the intensity centre (IC): where the intensity magnitude MI is given by MI k evaluated at point IC.The IC corresponds to the location on the fault plane with the highest energy release (i.e., the location of the maximum fault displacement, or centroid moment; Beauval et al., 2010).Hence, this does not always match the epicentre (Bakun, 2006).In the present work, the IC is used as the macroseismic epicentre, as opposed to the classic definition of an epicentre as the point on the surface that is the vertical projection of the seismic focus where the rupture begins.RMS levels provide confidence intervals, indicating that the IC is within the area delimited by them.Typically, the 95%, 90%, 80%, 67%, and 50% are represented, and their shape is based on the number of IDPs (Bakun and Wentworth, 1999).
The magnitude MI of the IC is equivalent to the Mw of the earthquake, and the equations ( 4) and ( 6) are calibrated in Mw.However, according to Bakun and Wentworth (1997), the accuracy of the calculated Mw is dependent on the number of IDPs used.Further details are provided in Bakun et al. (2011) andGomez-Capera et al. (2022).

Mw validation results
In order to further validate the reliability of the calibration models obtained using the BW97 method, we processed 15 Italian earthquakes with known location and moment magnitude (Table 4; Fig. 9a) and a macroseismic dataset of 4.533 IDPs (Fig. 9b).The macroseismic locations were in good agreement with the instrumental ones.
Macroseismic intensity attenuation was not a critical factor for location.Reliable source location estimates depend on the quantity and quality of the IDPs, as well as the geometry of the data field relative to the earthquake source; no technique can provide accurate source locations for events outside the network of IDPs.In the calibrationvalidation process, we focused on the prediction of moment magnitude (Mw).
Moment magnitude estimates are heavily dependent upon the macroseismic intensity attenuation model, which is calibrated using events for which instrumentation is available to accurately determine earthquake parameters and IDPs.The results of this calibration are displayed in Figure 12. Figure 12a displays the instrumental Mw and the calculated Mw from macroseismic data for each event (Log-Lin_10km and CRV5,9,16km models).Figure 12b displays the instrumental Mw and calculated macroseismic Mw with Log-Linear model in comparison with Mw calculated with Pasolini et al. (2008b) model and macroseismic magnitude (bxn; boxer method; Gasperini et al., 1999;Gasperini et al., 2013) by CPT15 (Rovida et al., 2022a).
The differences between the instrumental Mw and calculated Mw values in this study, as seen in Table 4, are all within 0.0-0.3magnitude units for the Log-Lin_10 model and 0.0-0.4 for the CRV5,9,16km models, and 0.0-0.8 for Pasolini et al. (2008b).All four calibrated models exhibit similar performance, with a median absolute residual of 0.2 magnitude units (Table 4).
The Log-Linear_10km and CRV5,9,16km models computed macroseismic Mw values that align well with the instrumental ones, but less so with the magnitudes calculated with the Pasolini et al. (2008b) model and Boxer's proposal at CPTI15.This outcome was anticipated (Tab.1), as the Pasolini et al. (2008b) and Boxer method have not been calibrated on this dataset.The Log-Lin_10km and CRV5,9,16km calibration models were thus validated, confirming that the BW97 method is an effective tool for evaluating the performance of macroseismic intensity models and determining earthquake location and magnitude from macroseismic data in the present study.

Discussion and conclusions
In this study we propose new macroseismic intensity attenuation models calibrated in Mw for the Italian territory, using the most recent releases of the CPTI15 earthquake catalogue (Rovida et al., 2022a) and the associated DBMI15 macroseismic database (Locati et al., 2022).Two attenuation functionals, originally proposed by Howell and Schultz (1975), are analysed.These models are given by equations ( 2 A carefully selected set of Italian shallow tectonic earthquakes and IDPs were used (119 events from CPTI15 and 16.260 IDPs from DBMI15) to calibrate the macroseismic intensity attenuation models.This set of data was selected based on criteria such as magnitude range, spatial distribution, quality of instrumental moment magnitude, epicentral locations, and intensity range.The two new attenuation models calibrated in Mw are analysed, along with a set of alternative models.We found that the best models are Log-Lin_10 and CRV9, the former having a slightly higher standard deviation, but lower error in the calibrated coefficients than the latter.The estimated pseudo-focal depths are (h=9.87km±5.7%)for Log-Lin_10 and and (h=8.7km±7.9%)for CRV9, showing that both results of the calibration of coefficient h are in good agreement with the mean observed in the input dataset, which is h=10km.
These findings lead us to conclude that the two best models, Log-Lin_10 and CRV9, are both essentially valid to describe the macroseismic intensity attenuation in Italy with hypocentral distance and moment magnitude.
Following a principle of parsimony, a preference can be expressed in favour of the mathematically and computationally simplest model, which in our case is the regression Log-Lin_10 model for macroseismic intensity I.The CRV9 model type 2 still remains an alternative attenuation model calibrated in Mw: it may offer an alternative solution to the more frequently used type 1 model, particularly for strong (Mw>6) and shallow (h<10 km) events, as it can more accurately predict the maximum intensity in these cases.
As a result of the sensitivity analysis, the CRV model is more affected than the Log-Lin model to changes in h and Mw, especially in terms of maximum intensity for events of Mw>6.0.It is however noteworthy that the standard deviation of both type 1 and type 2 calibrated models (Table 2; Table A3) remains relatively unchanged under variations in pseudo focal depth h.agreement with the instrumental ones, with differences of çResMwç≤0.3and of çResMwç≤0.4respectively.These results indicate that the intensity attenuation models obtained in this study are reliable and can be used to accurately estimate macroseismic intensity at a site from Mw and epicentral distance, or to assess Mw from IDPs for a given earthquake.
In conclusion, this study provides updated macroseismic intensity attenuation models that can accurately predict the macroseismic intensity and/or Mw, and can be applied to evaluate seismic hazard in terms of macroseismic intensity in Italy.Gasperini et al. (bxn;1999) Table 1 The more recent intensity attenuation models referenced in literature for Italy

Figures Fig 1a Location
Table 3 Instrumental parameters of the earthquakes used for validation process and mean absolute error (MAE) for each event and macroseismic intensity attenuation model Table 4 Instrumental parameters of the earthquakes used for validation process and Mw estimations among the macroseismic intensity attenuation models and BW97 method for each event Figure 2c and 2d respectively.The intensity ranges between 3 and 11 MCS and shows a median value of 5 MCS and standard deviation of 1.6.The average epicentral distance is 63 km, the median is 43 km and the standard deviation is 63 km.The distribution of the epicentral distances is asymmetric and the majority (3 standard deviation) of IDPs are located within 252 km from the epicentre, reflecting the decay of the number of intensity estimates at long distances.The following table summarises the range of the calibrated model: ) are represented in Fig. 4a.The histogram of the residuals (Icomp./Iobs.)for model (6) are shown in Fig. 4b, with a standard deviation of 0.149 computed from the entire dataset and 0.15 obtained from the theoretical formula (A1.4) in Appendix 1.For easy comparison with Fig. 3b, the histogram of the residuals (Icomp.-Iobs.) is also shown in Fig. 4c, with a standard deviation of 0.731.Fig 4d illustrates the intensity I defined by model (6) versus epicentral distance r, in logarithmic scale, for magnitude Mw classes.The values of a, b, c, d, and pseudo focal depth h are different in both calibrated models (Log-Lin_10 and CRV9), which indicates that each equation has a different way of relating the variables R and Mw to the macroseismic intensity (I).Fig. 3c illustrates the model (Log-Lin_10) for different Mw classes; according to this model, which 2006): for each macroseismic intensity at the epicentre (I0), the local distance for I=3, called Dist_I3 was determined using the relationship by Pasolini et al. (2008b); IDPs with D≥Dist_I3 were disregarded.This restriction significantly reduced the input data to 12,587 IDPs.The resulting model, named Log-Lincut dist I=3
) and (3) which are combined with the classical equation for energy decay (1) and the earthquake size parameter (Mw) to propose two empirical macroseismic intensity attenuation models (4) and (5), respectively named Log-Lin type 1 and CRV type 2. These two model types have different ways of relating the variables R, Mw to the macroseismic intensity I; e.g., a direct proportionality is assumed between I and Mw in the Log-Lin model and between their logarithms in the CRV model.This also implies different expressions of their curvatures, which are related to the pseudo focal depth in the Log-Lin model and to both pseudo focal depth and magnitude of the earthquake in the CRV model.

Further
analysis has been performed to test the models on a separate set of data, called the validation dataset, including 15 earthquakes from CPTI15 and 4.533 IDPs from DBMI15.The results of the validation process show that the four calibrated intensity attenuation models Log-Lin_10, CRV5, CRV9, and CRV16 performed better than the Pasolini et al. (2008b) model, with a median value of all event-based MAE scores equal to 0.6 (intensity unit).The calculated Mw values from IDPs obtained with the Log-Lin_10 and CRV9 models are in good

Fig. 7 Fig. 8
Fig. 7 Sensitivity analysis of the Log-Linear Model (Type 1) to pseudo focal depth h=5, 9.87, and 16 km, for different values of Mw