Modern earthquakes as a key to understanding those of the past: the intensity attenuation curve speaks about earthquake depth and magnitude

The Italian historical earthquake record is among the richest worldwide; as such it allows the development of advanced techniques for retrieving quantitative information by calibration with recent earthquakes. Building on a pilot elaboration of northern Italy earthquakes, we developed a procedure for determining the focal depth of all Italian earthquakes from intensity data alone. In a second step the procedure calculates their magnitude, taking into account the inferred depth. Focal depth exhibits a substantial variability countrywide, but so far received little attention: pre-instrumental earthquakes 5 were routinely “flattened” at upper crustal level (∼ 10 km), on the grounds that the calculation of focal depth is heavily dependent on the largely unknown local propagation properties. We gathered a learning set of 42 earthquakes documented by reliable instrumental data and by numerous intensity observations. Rather unexpectedly we observe (1) that within 50 km from the epicenter the ground motion attenuation rate is primarily controlled by focal depth and largely independent of magnitude; (2) that within this distance the fluctuations of crustal attenu10 ation properties are negligible countrywide; and (3) that knowing both the depth and the expected epicentral intensity makes it possible to estimate a reliable magnitude. 1 https://doi.org/10.5194/nhess-2022-30 Preprint. Discussion started: 17 February 2022 c © Author(s) 2022. CC BY 4.0 License.

across the full earthquake magnitude range. The new procedure is summarized in Figure 2 with reference to a specific recent earthquake (20 May 2012; ID 13 of Table 1). After drawing the first 10 km-wide, circular search area centered in the epicenter, and the following nine 10 km-wide ring-shaped search areas, each one shifted by 5 km from the previous one (0-10, 5-15, 150 10-20, 15-25, 20-30, 25-35, 30-40, 35-45, 40-50 and 45-55 km: red and blue lines in Figure 2b, up to the distance of 55 km from the epicenter, the resulting 10 intensity estimates of averaged MDPs are used to build the attenuation diagram ( Figure   2c).  Table 1) from the HSIT database (available from: https://e.hsit.it/772691/index.html); b) map showing the first 50 km from the epicenter and the ten ringshaped search areas centered in the epicenter (shown by red and blue lines), each one shifted by 5 km with respect to the previous one; c) plot of the 10 intensity values obtained averaging the MDPs falling in each of the rings: #1 reports the average intensity calculated for the 0-10 km search area, #2 the average intensity for the 5-15 km search area, and so on. 7 https://doi.org/10.5194/nhess-2022-30 Preprint. Discussion started: 17 February 2022 c Author(s) 2022. CC BY 4.0 License.

Data selection criteria
To compose our learning set (Table 1) we searched the whole Italian territory, selecting all instrumentally well documented earthquakes that also feature good quality macroseismic information. As in Sbarra et al. (2019a) we built a learning-set comprising macroseismic data obtained either from a direct survey or collected through the web, so as to gather information to be used as a sort of "Rosetta stone" for obtaining the parameters of historical earthquakes. For pre-2007 events we used the data stored in the DBMI15 v2.0 catalogue , a compilation of macroseismic intensities for Italian earthquakes that occurred in the time window 1000 CE-2017, whereas for more recent events of M w ≤ 5.9 we either used intensity 160 data from the web-based HSIT catalogue (Tosi et al., 2007;De Rubeis et al., 2019;Sbarra et al., 2019b), or used the MDPs collected by targeted post-earthquake surveys conducted by experienced INGV-QUEST personnel (QUick Earthquake Survey Team; http://www.ingv.it/quest/index.php/rilievi-macrosismici; last accessed on 20 December 2021). Only for the M w 5.5, 18 January 2017 event we had to make an exception to this rule, due to the incompleteness of HSIT data caused by the evacuations following the M w 6.0 mainshock of 24 August 2016. 165 The events comprising the learning set were further selected based on the following criteria:  4. only for pre-2012 earthquakes, the event must not be an aftershock occurring within a week of the mainshock, or a foreshock that occurred less than 24 hours before the mainshock; 5. all earthquakes having M w ≤ 5.8 must not be aftershocks of the central Italy sequence of 2016; 6. the earthquake must be documented by at least 100 MDPs (as in Sbarra et al., 2019a), at least 60 of which must fall within the first 55 km from the epicenter; 175 7. the MDPs falling at a 10-55 km distance from the epicenter must be distributed in an azimuth range > 180 • ; 8. the attenuation slope must be calculated based on six or more averaged points; 9. the standard error of the estimated attenuation slope must be ≤ 0.01.
All 42 earthquakes listed in Table 1 fulfill these rather strict criteria with the only exception of #6 and #17 (MDP < 60), two deeper events that were already included in the learning set of Sbarra et al. (2019a) as they are crucial for characterizing lower 180 crustal and subcrustal seismicity.
Notice that the selection criteria 1-4 had already been adopted by Sbarra et al. (2019a). The additional criteria (5 through 9) were added in consideration that the present work deals with the entire Italian territory, and hence with a much larger diversity of the potentially concerned earthquakes. More specifically: central Italy sequence, due to the widespread evacuations following the M w 6.0 mainshock of 24 August 2016 and to the superposition of the effects of subsequent shocks; criterion #6 was added after various experimental tests, in order to achieve more reliable and stable estimates of the attenuation; criteria #7 and #8 were introduced to discard earthquakes located offshore or near the coastline, whose epicentral location 190 generally exhibits greater uncertainty; criterion #9 was adopted to retain only earthquakes for which we could calculate a reliable attenuation slope.

Analysis of the learning set
We analyzed separately two data subsets, respectively comprising only earthquakes located in Northern Italy, and earthquakes located in the rest of the Italian peninsula. We made this choice because the dataset used in Sbarra et al. (2019a) included only 195 earthquakes from a region whose lithospheric structure and wave propagation properties are rather homogeneous; conversely, in this work we wanted to evaluate the possible influence of variable attenuation properties resulting from the full range of tectonic and geodynamic diversity characterizing the Italian peninsula.
Due to the intervening minor updates in our methodology -and specifically in the calculation of the starting point of our moving window, which implies a slightly different slope for the first 50 km of the attenuation curve -we first recalculated 200 the attenuation slope for all the 20 earthquakes comprising the learning set used by Sbarra et al. (2019a) (Table 1, ID from 1   to 20; see Figures 1 and 3). We added to this dataset the 1996, Emilia earthquake (ID 23), originally rejected because its depth from the ISIDe catalogue (ISIDe Working Group, 2007) was fixed at 10 km; for this event we now use the depth evaluated by Selvaggi et al. (2001). We then analyzed the earthquakes we selected for the rest of Italy and calculated their intensity attenuation slope (Figures 1, 4; Table 1, from #21 to #42, except for #23).

205
As discussed earlier, in both datasets, which together form our new learning set, we observed a distinct break in slope at an epicentral distance of about 50 km. In describing this feature of the Italian attenuation curves, Gasperini (2001) contended that within a ∼ 50 km epicentral distance the ground shaking is dominated by direct seismic phases, whose propagation is highly sensitive to earthquake depth, whereas Moho-reflected phases dominate at larger distances. According to this hypothesis, the exact distance of the transition would be controlled by the average Moho depth along the source-receiver path.

210
For all the earthquakes of the learning set we then plotted the steepness (S), i.e. the absolute value of the slope, versus focal depth (D), and found two separate but very similar best-fitting logarithmic functions ( Figure 5). For northern Italy we found S = (−0.020 ± 0.003) ln D + 0.093 ± 0.009 (1) whereas for central and southern Italy we found The coefficients of both these functions fall within their combined confidence limits, suggesting that our method does not detect any statistically significant change in the attenuation properties between the two domains, at least over the first 50 km of epicentral distance. This finding also suggests that an approach based on averaging the intensity values distributed over circular search areas has the ability to smooth out most of the inevitable azimuthal differences in crustal propagation properties.
We decided to calculate a new logarithmic function using all 42 earthquakes of the learning set, so as to obtain a law that 220 may be used over the whole Italian region (green line in Figure 5): The Pearson coefficients of the three regressions are quite good, corresponding to a significance level of correlation better than 10 −3 . As we expected, Eq. 3 is similar to the previous two equations, and exhibits a narrower 95% confidence interval resulting from the larger number of available data points. Eq. 3 can be applied for a steepness interval 0.058 ≤ S ≤ 0.012, 225 which corresponds to a depth interval 5 ≤ D ≤ 73 km. Notice that the function is not constrained beyond these limits, and hence should not be used for shallower or deeper events. Notice also that for epicentral distances > 50 km the curves shown in Figures  (Table 1, from #21 to #42, except for #23; red symbols in Figures 1 and 5). Individual intensity data points were obtained by averaging the intensity values as shown in Figure 2.
We obtained the linear fit for the first 50 km of each curve and calculated the resulting slope.
3 and 4 exhibit different attenuation between the northern and the central-southern Italy datasets, but are remarkably similar in the first 50 km, in agreement with the observations made by Gasperini (2001). In particular, the attenuation of earthquakes occurring in northern Italy, where the crust is generally thicker than in the rest of the country, shows a characteristic, very 230 gently sloping plateau that has been interpreted as due to Moho reflected seismic waves by Bragato et al. (2011).
A further element to be taken into account is the difference in seismic wave propagation between the Tyrrhenian and Adriatic sides of Italy, most likely resulting from a rather different efficiency of the crust-upper mantle system (Mele et al., 1997;Lolli et al., 2015;De Rubeis et al., 2016;Di Bona, 2016). This difference, however, can only be appreciated in the far-field, i.e.
beyond a 50 km distance from the epicenter, implying that most likely it does not show in our analyses, nor does it affect their 235 results.
Finally, using the new learning set we obtained a new Intensity Prediction Equation (IPE) for Italy (4): where I is the intensity and r is the hypocentral distance in km. This equation rests on the assumption that the macroseismic fields used to build it contain fairly well-distributed data, both in the near-field and in the far-field.  the central and southern Italy datasets, respectively: the corresponding best fitting logarithmic functions are shown in blue (Eq. 1) and red (Eq. 2), respectively, along with their 95% confidence intervals. The best fitting function obtained for the whole dataset (Eq. 3) is shown in green. Each earthquake is labelled with a unique identifier (see Table 1) and is plotted along with its standard error (shown by a vertical bar).

Independence of earthquake depth from magnitude
The slope of the first 50 km of the attenuation curves calculated for the earthquakes of our learning set (Figures 3, 4) is independent from magnitude, as already experimentally observed by Sbarra et al. (2019a) for a smaller sample of events. To prove this statement we plotted in Figure 6a the attenuation curves for four earthquakes falling in the M w range 4.8-6.5, but having a similar instrumental depth, in the range 7.3-8.7 km (#4, 16, 24, 39; see Table 1). Figure 6a shows that all the calculated 245 steepness values fall in a rather narrow range (0.045-0.051), regardless of magnitude. Figure 6b shows that the same behavior is observed also for four deeper earthquakes (#8, 10, 11, 26; see Table 1, which share a similar instrumental depth (24.5-29.2 km) but exhibit a different M w (4.0-5.6).  Table 1 for further details). The slope of the bestfitting line in the first 50 km is similar among the four events reported for each group, providing experimental evidence for the independence of depth from magnitude.
The invariance of the attenuation slope with magnitude is a key point as it makes our approach suitable for analyzing historical earthquakes even if their size is not well constrained. On the contrary, nearly all the methodologies developed in 250 the past to calculate earthquake depth use magnitude as an essential input parameter (Jánosi, 1906;Kövesligethy, 1907;Blake, 1941;Sponheuer, 1960;Ambraseys, 1985;Burton et al., 1985;Levret et al., 1996;Musson, 1996;Traversa et al., 2018;Provost and Scotti, 2020).

A logarithmic law to determine earthquake depth
We analyzed the possibility of reproducing the experimental trend of Figure  both. To explore the variation of the attenuation slope with depth we therefore used three of the models that explicitly include this parameter: our new IPE (Eq. 4), the IPE by Musson (2005), and the GMPE by Cauzzi and Faccioli (2008). We chose these functions because they feature a simple functional form which determines a magnitude-independent attenuation slope, as suggested by real earthquake data. Conversely, a functional form containing a magnitude-distance mixed term would lead to a change in the shape of the attenuation curve with distance and to a variation of the slope for a variable magnitude. 265 We then used two different conversion equations (Faenza and Michelini, 2010;Masi et al., 2020) to turn the PGA values obtained with the adopted GMPEs into macroseismic intensities, so as to also test the influence of the conversion process.
We used these equations to compute the macroseismic intensities caused at several epicentral distances by a hypothetical earthquake located at variable depth. We then applied the same 10 km-moving window average method used for the analyzed earthquakes and calculated the regression line within a distance of 50 km.
270 Figure 7a shows some of the average intensity values obtained using the IPE proposed by Musson (2005), along with their regression lines. We remark that, although the macroseismic intensity is proportional to the logarithm of the hypocentral distance, the linear regression of intensity versus epicentral distance gives statistically significant results in the adopted distance range. Figure 7b shows the steepness of the regression lines thus calculated as a function of the earthquake focal depth, and the values, calculated with the same method, derived from PGA (using the GMPE proposed by Cauzzi and Faccioli (2008), 275 converted into macroseismic intensity. It is worth noting that the differences caused by the use of the IPE in place of the GMPE are comparable to the differences caused by the use of two different conversion equations. At any rate, in all three cases the trend of the values as a function of the depth is similar to that found experimentally. Also for these synthetic data, the Pearson coefficients of the four linear regressions between the steepness and the logarithm of focal depth (R = 0.995) are significant, with a probability of less than 10 −3 to find by chance the same regression in uncorrelated data. Having been obtained with a 280 completely different kind of data, this result suggests that the approach followed for deriving Eq. 3 is adequate and reliable.

Reliability and validation of the depth estimation method
The reliability of the slope of the first 50 km of the attenuation curve depends on the quality and spatial distribution of the available MDPs and on the accuracy of the epicentral locations. Italian macroseismic data are systematically stored in the DBMI database v2.0 ; as a rule of thumb, the older is the earthquake, the less complete and reliable are the 285 historical sources from which macroseismic intensities were derived (e.g., Guidoboni and Ebel, 2009).
To test our procedure we investigated the minimum number of MDPs of the macroseismic field that are needed to obtain a reliable estimate of the attenuation slope. To this end we intentionally depleted the macroseismic field of the 20 May 2012, M w 5.8, Emilia earthquake (#13), a well-recorded event for which over two hundred spatially well-distributed MDPs are available, using data from the HSIT database (De Rubeis et al., 2019); Figure 2a. For each of the ten ring-shaped search areas (see Figure   290 2b, 2c) we performed a gradually increasing reduction of the number of MDPs (from 1% to 99%), the same for each step and for all areas. The regression of the attenuation trend was calculated 1,000 times for each depletion step, so as to evaluate the steepness variability through its standard deviation.  Our depletion test demonstrates that we may obtain a reliable attenuation slope -and hence a reliable depth estimationeven for historical earthquakes for which only few MDPs are available, provided that they are homogeneously distributed for 300 each distance window.
We remark that the reliability of our estimates improves significantly when the value of each of the variables indicated by our selection criteria #6 to #9 (number of MDPs within 55 km, azimuthal range, number of used circular windows, and standard error: see Sect. 3.1) rises above the average value calculated for all 206 selected events. Therefore, the steepness values obtained through our methodology are more reliable if the macroseismic field includes more than 99 MDPs within the 305 first 55 km from the epicenter, distributed over an azimuthal range > 270 • , its attenuation curve is described by 10 averaged points and the standard error is < 0.0058, implying that the data points are well aligned on a straight line (see Table S1). https://emidius.mi.ingv.it/ASMI/event/20120520_0203_000), which occurred nearby. These circumstances misled our method, causing a drastic overestimation of the earthquake depth (36.8 km). Conversely, thanks to the rapidity in the response given by citizens and to the ensuing lack of contamination, the HSIT dataset method returned a depth ≤ 5 km, much closer to instrumental estimate (8.1 km).

A two-step method for estimating magnitude based on intensity and depth
While the focal depth is independent of magnitude and can be obtained simply based on the steepness of the line that best fits the first 50 km of the attenuation curve, the estimation of the magnitude itself affects the y-intercept (the expected intensity at the epicenter, I E ) of the linear fit ( Figure 6): for a constant depth the y-intercept increases for an increasing magnitude ( Figure 6), and decreases if depth increases. Therefore, a reliable magnitude determination based on macroseismic data must 320 necessarily take into account earthquake depth.
The method used in this work to estimate magnitude differs from that described by Sbarra et al. (2019a), where magnitude was calculated simply by reversing the IPE proposed by Tosi et al. (2015). To determine magnitude more reliably we devised a two-step procedure where depth is estimated first (Step 1), then M w is estimated using our learning set data to derive a standard least squares regression equation among D (depth), I E and M w (Step 2). Figure 9 shows the learning set data and the contour Magnitudes obtained through this procedure are referred to as y-intercept M w . In this perspective the attenuation curve becomes a sort of 'earthquake identity card', as it contains all the elements needed to retrieve magnitude and depth from the observed intensities, provided that a reliable calibration scheme is available. Such calibration can be regarded as an application 330 to Seismology of the principle of actualism popularized by British naturalists in the late 18th century: "Observing modern earthquakes to understand those of the past".
The method summarized by Eqs. 3 and 5 is simpler and more intuitive than the methods based on a joint inversion of magnitude and depth (Musson, 1996;Bakun and Wentworth, 1997;Sirovich et al., 2013;Traversa et al., 2018;Provost and Scotti, 2020), plus it may allow a geological verification of the depth before estimating the magnitude. Our Step 2 uses a method 335 similar to that proposed by Gutdeutsch et al. (2002), who applied it to carefully selected datasets only, so as to minimize the bias caused by a poorly constrained depth or by an incomplete macroseismic field.
In conclusion, starting from our experimental observations of the independence of the attenuation slope from magnitude, we were able to mitigate the trade-off between magnitude and depth when estimating both these parameters from macroseismic data.

Dealing with larger magnitude earthquakes
Our approach works well if the size of the seismic source is negligible relative to the epicentral distances, but it may not be immediately applicable to estimate the attenuation of the macroseismic intensity for a large magnitude earthquake (Gasperini, 2001;Albarello and D'Amico, 2004;Pasolini et al., 2008). To test the validity and possible limitations of this assumption we evaluated the maximum magnitude for which the use of a point-source approximation is granted, using both our learning set 345 and our analyzed set.
We used the empirical relationships proposed by Wells and Coppersmith (1994) to calculate the rupture area and the expected length of the seismogenic source based on the y-intercept M w (Eq. 5). Assuming a dip angle of 45 • for every fault, irrespective of its kinematics and tectonic setting, we calculated the surface projection of each rectangular source and the radius R e of the equivalent circle (i.e. a circle having the same projected area as the fault: thus R e is a function of M w ). We found R e greater 350 than 10 km only for earthquakes of M w ≥ 6.75; but not having any such earthquakes in the learning set (Table 1), we used a geometric correction only to infer the depth of 21 earthquakes of the analyzed set (see discussion at the end of this section and in Table S1).
Then we applied to this group of larger magnitude earthquakes a procedure that we call "variable moving windows". More specifically, we used as the first search area a circular window of radius R e , inside which we averaged the MDPs intensities, 355 while for the subsequent windows -each one shifted by 5 km, as usual -we adopted the standard 10 km radius increase.
For the 13 January 1915, M w 7.1, Marsica (central Apennines) earthquake, one of the largest in Italian history, we made a test using the R JB distance (Joyner and Boore, 1981) and calculating the average of the MDPs inside a rectangular rather than a ring-shaped area. We singled out this earthquake because its macroseismic field should not have been contaminated by the effects of previous significant events, which make it difficult to separate the individual contribution of a specific shock to the 360 cumulative damage (Grünthal, 1998;Grimaz and Malisan, 2017;Graziani et al., 2019); a circumstance that would ultimately affect the attenuation slope and hence contaminate the inferred earthquake depth. This is a recurring problem in historical earthquake catalogues; a condition that is hard to overcome even for modern earthquakes, and even if a very rapid damage survey is carried out, because the first large shock inevitably causes an increase in the vulnerability whose effects on later shocks are virtually impossible to identify.

365
The comparison of the attenuation slope calculated using the R JB or using the moving window/variable moving window methodology proposed here shows only modest fluctuations; in fact, they are comparable to the error arising from the uncertainties in the epicentral location. For the 1915 earthquake we found a steepness of 0.044, 0.044 and 0.047, respectively using the R JB approach, the moving window approach and the variable moving window approach. These steepness fluctuations imply a difference of about 1.5 km in the expected depth. Given the uncertainty about our knowledge of the source geometry 370 for historical earthquakes and the limited impact of using the R JB distance in the window approach, we decided to recalculate all distances for earthquakes having M w ≥ 6.75 using the variable moving windows method only, which analyses the MDPs over circular search areas.
In conclusion, using an extended source approach for the largest earthquakes has a minimal influence on the steepness.
Conversely, the effect on the y-intercept (I E ) is not negligible. The correct way of calculating their magnitude would be using 375 R JB distances, but due to the lack of information on source geometry we use the variable moving windows method and apply a geometric correction to the intercept value. As a result, for 21 earthquakes having M w ≥ 6.75 (see Table S1) we assumed an extended source with a radius of R e . Consequently, the distances of the relevant MDPs were systematically reduced by R e , leading to a geometric correction of the regression line and of its intercept: Finally, we recalculated the magnitude of these 21 earthquakes using Eq. 5.

380
3.8 A method for inferring the depth and magnitude for earthquakes of the CPTI15/DBMI15 catalogue We applied our methodology to the pre-1984 earthquakes of the DBMI15 v2.0 catalogue . We analyzed only pre-1984 events because their parameters were computed from intensity data as their instrumental location is generally unreliable, although there are notable exceptions (see discussion in Rovida et al. (2021)).
We first selected the earthquakes to be analyzed: they must meet all criteria listed in Sect. 3.1 "Data Selection criteria" except 385 for #6, which we relaxed by reducing the minimum number of MDPs from 60 to 30, based on the conclusions drawn in Sect.
3.5 "Reliability and validation of the depth estimation method". These criteria were passed by 206 out of 2,679 earthquakes ( Figure 10 and Table S1), which comprise the analyzed set of this work. Unfortunately, most of the events listed in DBMI15 (80%) exhibit less than 30 MDPs within the first 55 km from the epicenter, and therefore had to be discarded. To estimate the depth of the 206 events that were retained we first calculated the steepness of the line that best fits the first 50 km of the 390 attenuation curve of each event, then we used Eq. 6, which is simply the reverse of Eq. 3: In Sect. 3.2 we clarified that we can calculate a reliable depth only for events whose S falls in the interval 0.058 to 0.012 ( Figure 5), which corresponds to the depth interval 5.0-73.0 km, respectively (Table S1). This implies that an inferred 5.0 km depth must be intended as ≤ 5.0 km, and similarly, a 73.0 km depth stands for ≥ 73.0 km. Notice that the uncertainty 395 associated with the inferred depths is determined by the confidence interval shown in Figure 5, and is hence larger for deeper earthquakes. In addition, Eq. 3 and Eq. 6 are affected by the accuracy of the instrumental location of the learning set earth- quakes, on the basis of which the logarithmic curve data are fitted. Furthermore, we must keep in mind that the reliability of the depth determination ultimately depends on the quality of the available macroseismic data.
Once the depth of the 206 selected earthquakes is known, we can estimate their magnitude using Eq 5. All estimated depths 400 and magnitudes are shown in Figure 10 (see Table S1).

Reliability of the magnitude estimation method
As a countercheck we used our method to calculate the depth first (Eq. 6) and then the magnitudes (Eq. 5) of the 42 events of our learning set ( Figure 11, Table S2), so as to analyze the departure from the instrumental magnitudes listed in Table 1.
We obtained differences in the range 0.68 to -0.41 magnitude units, respectively, with an average of -0.03 and a mean squared 405 deviation of 0.28. Figure 11. Correlation between the Mw calculated with the y-intercept approach proposed in this work and the instrumental Mw reported in We then compared the macroseismic magnitudes calculated through our method with those calculated through the Boxer method, using the very same intensity dataset from DBMI15. Notice, however, that the parameters of the earthquakes comprising our learning set were computed using also data from other sources, such as HSIT, CFTI5Med etc. (Table 1). Table   2 lists the magnitude of all learning set earthquakes for which the comparison was possible. Notwithstanding the significant 410 differences between the two methods, the mean squared deviations between instrumental magnitudes and those estimated with Boxer and our method are comparable, as they are 0.38 and 0.35 respectively.
We wish to stress once again that the reliability of the magnitude and depth determinations shown in Figure 10 and Table   S1 depends on both the quality of the macroseismic data and the accuracy of the epicentral locations. For completeness of information, Table S1 reports also the full details of the processing for each of the selected events. Some of the results may 415 appear unrealistic, due to inherent uncertainties that are reflected in the determination of the steepness and of the y-intercept (as defined in Sect. 3.7): these may include the cumulation of damage from subsequent shocks, unpredictable anomalies in wave propagation, strong source directivity and site amplification effects, all of which may also cause a sizable shift in the epicentral location.
Finally, the reliability of y-intercept M w is a function of the accuracy of the estimated depth. For instance, we examined the 420 13 January 1909, Northern Italy earthquake (https://emidius.mi.ingv.it/ASMI/event/19090113_0045_000), whose macroseismic field is suggestive of a rather deep source. We obtained a depth of 44 km, yielding a M w 5.5; should this earthquake be much shallower (e.g. 5 km), Eq. 5 would return a substantially smaller M w (5.1).

Comparison between y-intercept M w and Boxer M w
Before comparing the M w estimates obtained with our approach and those listed in the CPTI15 catalogue we must recall 425 that all M w estimates supplied by this catalogue are inherently hybrid, because the decision to adopt the macroseismic or the instrumental value as 'preferred' is made on a case-by-case basis. To minimize the ambiguities that may arise from these circumstances our analyzed set includes only 206 pre-1984 events that satisfy all our selection criteria; for the vast majority of these events CPTI15 adopted as preferred the intensity-based magnitude (Rovida et al., 2021). Figure 12 shows a comparison of the magnitude obtained with our methodology with the corresponding magnitude listed 430 in CPTI15 (Table S1). The two sets of estimates are generally consistent, yet on average the magnitudes calculated in the present work show a difference of +0.25 magnitude units. Moreover, Vannucci et al. (2021) stated that the magnitudes of all pre-instrumental earthquakes in CPTI15 might be overestimated by 0.1-0.2 units due to differences in the response of pre-1960 seismographs relative to the response of more recent and better calibrated electromagnetic sensors. If this is the case, the difference between our estimates and the CPTI15 estimates summarized in Figure 12 would be even larger.

435
The calculated M w may also vary if we consider macroseismic intensities assigned using the MCS or the EMS scale; according to Vannucci et al. (2021), using one or the other may cause differences in the macroseismic location. It is important to be aware that the calculation of M w from macroseismic data is controlled by a number of variables whose relative weight is critical: assigning proper weights, however, is not an easy task, regardless of the quality of the data and of the reliability of the adopted algorithm.  The +0.25 magnitude unit difference we found implies that on average our seismic moments are 2.3 times larger than those obtained with conventional methods; a conclusion that may have strong implications for the assessment of seismic hazard. In this study we proposed a two-step procedure for deriving the depth and magnitude of Italian earthquakes from official, publicly accessible macroseismic intensity data, i.e. from the "classic" DBMI15 database and from the innovative, web-based 445 HSIT dataset. The main merit of the proposed methodology is its objectivity and ease of application.
Web-based macroseismic systems allow a large amount of data to be collected through crowd-sourcing, and are often the only available source of information concerning the effects of low-to-medium magnitude earthquakes. In fact, HSIT data were critical to performing this work because -especially for deeper earthquakes (> 30 km) -they were almost always the only available observations.

450
We proved that the initial 50 km of the attenuation curve contain all the elements needed to retrieve not only the depth, but also the magnitude of any given earthquake.
The first step of our procedure involves the calculation of the earthquake depth (Eq. 6). Based on our experimental observations we have shown that the steepness of the attenuation curve in the first 50 km from the epicenter is independent of magnitude and is solely a function of the source depth. We also found that the steepness does not vary much as a result of any 455 regional differences in the propagation properties. This finding implies that these properties do not change much countrywide, despite the well-known complexity of Italian geodynamics and despite the ensuing geological heterogeneity, and that our new relations are valid for the whole Italian territory ( Figure 5).
The second step involves estimating the magnitude using an empirical law obtained from a regression function that relates the expected epicentral intensity to the depth and magnitude of the 42 earthquakes comprising our learning set (Eq. 5).

460
We applied this methodology to 206 earthquakes from the CPTI15 catalogue, after removing all events whose macroseismic field is too sparse or too inhomogeneous to return reliable results.
Our approach allowed us to verify that the inferred depth is consistent with the presumed earthquake-causative tectonic structures, and is essential to obtain a well-calibrated magnitude value. We contend that the new methodology may be crucial for mitigating the trade-off between earthquake depth and magnitude; this is a pre-condition for calculating reliable depth 465 estimates -and hence reliable magnitudes -for earthquakes of the pre-instrumental era.
The historical record is the main pillar of any seismic hazard analysis, conducted at any scale and using any approach.
We maintain that the revised framework discussed in this work may ultimately serve for exploiting more systematically the enormous potential of historical earthquake data, and ultimately for providing inherently more reliable input data for seismic hazard assessment.