## Appendix A  
The KF inversion method allows us to retrieve fault information of historical earthquakes based on macroseismic data. The model should be as simple as possible, assuming that the seismic source is a dislocation propagating horizontally on a rupture plane of unit width. The kinematic function, \(KF\), in [23] is  
\[KF(P,I)=\frac{R(P,I)}{D(P,I)[1-\left(\frac{V_{s}}{V_{s}}\right)\cos\theta(P,I)]} \tag{1}\]  
where \(KF\) is the contribution of a source dislocation point propagating along on a horizontal line on a rupture plane of unit width at depth H and distance \(l\) from the nucleation point to the displacement-related ground motion at the receiver point \(P\) on the surface. \(KF\) considers body wave and \(R(P,I)\) is the radiation pattern of S-waves [31]. \(D(P,I)\) is the distance between the source and the receiver \(P\) and \(\theta(P,I)\) is the angle between the ray reaching \(P\) and the direction of the rupture propagation. To summarize, on the denominator, \(D(P,I)\) is the geometrical spreading of the waves and the other term is the directivity. \(KF\) is consistent with the asymptotic approximation [29]. In the asymptotic assumption, only the far-field term of the impulse response ([PERSON]'s function) in the representation theorem is used (see [31] for the representation theorem and [29] for [PERSON]'s function). The \(KF\) model uses 11 source parameters adopting the convention of the positive direction of the strike ranging from \(0^{\circ}\) to \(360^{\circ}\), with the fault plane dipping to the right. According to this, rake anglesranging between \(1^{\circ}\) and \(179^{\circ}\) indicate faults with reverse (compressive) components. In our convention, the total length of the rupture is the sum of the absolute values of the parts along the strike (considered positive L+) and antistrike (negative L\(-\)). The same goes for the Mach number (Mach+, and Mach\(-\)). In the asymptotic approximation, only the high-frequency radiation of the source is considered, where high frequency means wavelengths shorter than the shortest observer \(P\) source. This means that \(KF\) uses as a propagation medium an elastic half-space in the distance range of about 5 to 100 km from the source. In the version in [53] with the niching genetic algorithm (NGA) technique [54], it was reduced to 10 parameters because L is derived from \(M_{0}\) by the empirical relationship M versus RLD, the subsurface rupture length, [33] for all types of faults. So, the length parameter is the percentage of L in the strike direction. \[\log(\text{RLD})=0.59\text{ M}-2.44 \tag{10}\]  
The scale parameter (11 th in the new version) is the seismic moment \(M_{0}\), which is used in the calculation to convert the values of \(KF\) into pseudo-intensities \(I_{calc}\) by the semi-empirical data-fitting function A3, [17]:  
\[ \begin{array}{rl}I_{calc}(x,y)=\int[118.488-logKF\times 0.9456+(logKF)2\times 0.8101-logM_{0}\times 9.686\\ \qquad\qquad+0.2575\times logKF\times logM_{0}+0.2152\times(logM_{0})2]\end{array}  \tag{11}\]  
In Equation (11), \(I_{calc}(x,y)\) is pseudo-intensity at location (\(x\),\(y\)), \(KF(x,y)=\max KF\) (\(P\),\(I\)) is referred to a Cartesian plane, and \(M_{0}\) is the seismic moment in dyne cm. It is worth stressing that the ANOVA test shows that the overall quality of the model (11) is rather good; in fact, the probability of the null hypothesis is \(10^{-5}\) (see [17] for details). The objective function for the fitness of the inversion is \(\sum R_{s}^{2}\), where \(R_{s}\) is the difference between pseudo-intensity \(I_{calc}\) and the intensity \(I\) observed in the same site. It must be remembered that \(\sum R_{s}^{2}\) is now computed taking into account the partition coefficients of non-numerical intensities, according to [2], and those of uncertain intensities, according to Appendix B.