Reappraisal and Analysis of Macroseismic Data for Seismotectonic Purposes: The Strong Earthquakes of Southern Calabria, Italy

: In tectonically active areas, such as the Italian peninsula, studying the faults responsible for strong earthquakes is often challenging, especially when the earthquakes occurred in historical times. In such cases, geoscientists need to integrate all the available information from historical reports, surface geology, and geophysics to constrain the faults responsible for the earthquakes from a seismotectonic point of view. In this paper, we update and review, according to the EMS-98 scale, the macroseismic ﬁelds of the ﬁve main events of the 1783 Calabria sequence (5, 6, and 7 February, 1 and 28 March, M w 5.9 to 7.1), two other destructive events within the same epicentral area of the 1783 sequence (1791, M w 6.1 and 1894, M w 6.1), plus the Messina Strait 1908 earthquake (M w 7.1). For the 1783 seismic sequence, we also elaborate an updated and new catalog of coseismic effects. The new macroseismic ﬁelds were analyzed using a series of MATLAB algorithms to identify (1) the unitarity of the ﬁeld or its partitioning in sub-sources and (2) the ﬁeld and sub-ﬁelds’ main elongation. A collection of earthquake scale laws from literature was used to compute the average source parameters (length, width, and area) with their range of variability, and an elliptical map-view representation of the source geometry was calculated and made available. The analyses of such data allow us to speculate on the earthquakes/faults association, as well as propose new interpretations and reconstruct the space–time evolution of the signiﬁcant southern Calabria seismic sequences in the last ﬁve centuries.


Introduction
Since historical times, the Italian peninsula has been the scene of destructive (M w up to 7.1) extensional earthquakes [1]. In the instrumental age, these major events and other minor-magnitude seismic sequences are recorded with extreme precision, allowing the seismogenic sources to be delineated in detail (e.g., [2][3][4][5][6][7][8][9][10][11][12]). However, when it comes to pre-instrumental seismicity, the challenge is quite different. Historical earthquakes in areas with a long history of destructive events are essential for seismotectonic and seismic hazard purposes, especially when instrumental seismicity is moderate or almost missing. This is the case of the Southern Calabrian Arc, in the area between the Catanzaro Basin and the Messina Strait, which, after a paroxysmal activity, with six M w 6.5-7 events in 125 years from 1783 to 1908 and other two M w 6 events in 1791 and 1894 [1,[13][14][15], has only experienced small events with a few minor isolated sequences (M w < 4.5) [16]. Therefore,

Seismic Sources in the Literature
The attribution to the respective seismogenic sources of the earthquakes studied herein has long been debated, mainly due to the often absent or scarce coseismic effects of primary faulting. Most historical sources describe the 5 February 1783 earthquake (M w 7.1) as long with multiple events, possibly rupturing several fault segments [72]. In 1784, Dolomieu [87] described a~18 km long and 1-m wide surface fracture that occurred during the 5 February 1783 event along the Cittanova fault. Jacques et al. [72] and Galli and Bosi [26] considered it as clear evidence of surface faulting, attributing the earthquake to the W-dipping Cittanova fault [22,26] or the Cittanova and Sant'Eufemia faults [72,78,88] (Figure 2a). Some authors identify the source of this earthquake in a blind E-dipping low-angle normal fault on the west side of the Gioia Tauro basin [83,89,90] (Figure 2a). Similarly, the 7 February event (M w 6.7) was associated with either an E-dipping large normal fault (Mileto Fault in Figure 1) bounding the Mesima graben [91,92] (Figure 2b) or the W-dipping Serre fault [93,94]. The 6 February and 1 March events are not localized in the Italian parametric catalogs [95] because earthquakes with less than ten intensity data points or incomplete intensity distributions were discarded from the dataset [34]. An S-dipping transversal strike-slip fault located, between the uplands of the Serre Mts. and the Gioia Tauro Plain and the Cittanova and Serre fault has been proposed for the 6 February 1783 earthquake [83] ( Figure 2c).
The source of the 28 March event is questioned even more. This event is the second most energetic of the overall 1783-SS (M w~7 .0), although it may have been overestimated because it was the last shock of a long-lasting destructive sequence. The event generated destruction at the hanging wall of the Serre fault, but also largely expanded toward the Ionian coast. Some authors [85,96] propose the activation of the S-dipping LCFS along the northern border of the Catanzaro basin. This fault would represent the surface expression of an E-W deep transform [97], allowing the SE-ward advancement of the Calabrian arc in the overriding plate. Other authors hypothesize NS-striking [98] or NNW-SSE-striking [99] sources on the Ionian slope of the Catanzaro basin ( Figure 2b). Pirrotta et al. [81] proposed the NNW-dipping Caraffa fault system as the source responsible for the event.  [1,15]) from 1000 to 2020 are noted in the yellow squares and indicate the historical earthquakes studied in this paper. (c) Location map; red, blue, and green lines represent the crustal (≤40 km) deformation field from the P and T axes of focal mechanisms present in the TDMT (Time Domain Moment Tensor) catalog [86] (1985-2023).  (9) Vallefiorita faults (described in the next section). Each minor circle within the rose diagram represents a step count of 5 items. (b) Major historical earthquakes from the Italian Parametric Catalogue (CPTI15 v4.0, [1,15]) from 1000 to 2020 are noted in the yellow squares and indicate the historical earthquakes studied in this paper. (c) Location map; red, blue, and green lines represent the crustal (≤40 km) deformation field from the P and T axes of focal mechanisms present in the TDMT (Time Domain Moment Tensor) catalog [86] (1985-2023).
The 1791 (M w 6.1) and 1894 (M w 6.1) earthquakes are less represented in the literature of pre-instrumental earthquakes. The former struck the Serre chain in the same epicentral area as the event of 7 February 1783; in fact, it is associated with a possible triggering phenomenon related to the 1783-SS. The epicenter of the 1894 earthquake was located on the northern slope of Aspromonte, in an area just south of the 5 February 1783 mesoseismic zone. Two N-dipping transverse strike-slip faults between the uplands of the Serre Mts. and the Gioia Tauro plain and the Cittanova and Serre faults, one onshore and one offshore of Scilla in the Tyrrhenian Sea, have been suggested for this earthquake [83] (Figure 2c).
Regarding the seismogenic fault of the 1908 earthquake (M w 7.1), the scientific community is divided into two groups: one supporting a high-angle NW-dipping onshore Calabrian fault [84,100,101], the other an offshore low-angle SE-dipping major fault in northern Sicily [30,[90][91][92][102][103][104]. A large amount of the literature discusses this problem and an updated overview may be found in Barreca et al. [74]. Mulargia et al. [105] considered the possibility of simultaneous activation of both faults, particularly of the southern portion of the W-dipping Reggio Calabria fault and the northern portion of the Messina E-dipping fault [84] (Figure 2d). In this paper, we consider it useful to provide a shapefile containing all the seismic sources of the various earthquakes in the literature previously discussed (Supplementary Data S1).

Data and Methods
In this study, we present an innovative methodology that utilizes EMS-98 macroseismic field intensity points as the basis. By incorporating a statistical analysis of the intensity data points, examining relevant scale-laws from existing literature, and utilizing information on the coseismic and Quaternary fault pattern, our approach enables the identification of potential associations between earthquakes and faults. We have compiled a workflow diagram in the Supplementary Materials ( Figure S1) to summarize the methodological steps described below.

Data and Methods
In this study, we present an innovative methodology that utilizes EMS-98 macroseismic field intensity points as the basis. By incorporating a statistical analysis of the intensity data points, examining relevant scale-laws from existing literature, and utilizing information on the coseismic and Quaternary fault pattern, our approach enables the identification of potential associations between earthquakes and faults. We have compiled a workflow diagram in the Supplementary Materials ( Figure S1) to summarize the methodological steps described below.

Tectonic Trends Determination
The determination of tectonic trends of the potentially seismogenic faults of the study area is the first step of the procedure, as it serves to validate the trends obtained from the analysis of macroseismic fields. This step entails compiling fault traces from the literature (i.e., the faults described in the "Tectonic setting" section above) on a GIS platform and establishing the average trends, taking into account variations along the strike of any given structure. To generate directional rose diagrams, depicted in Figure 1, we first used ArcGIS functions to resample the fault traces with a sampling step of 1 km. Then, we selected the fault traces, taking into account the segmentation of the structures. Finally, we plotted the trends in rose diagrams to highlight the mean directions and, by distinguishing the plots from each other, the bends in the main fault alignments (Figure 1a). Table 1 shows the average trends and dip directions of the fault groups represented in the rose diagram plots. Table 1. Average direction and dip-direction of the most relevant fault systems in the study area. Id refers to rose diagram plots in Figure 1a.

Id
Fault Name Mean Direction Mean Dip-Direction Serre main fault and internal splay N029 Cittanova and outer splay N032 Armo and Reggio Calabria Scilla and Sant'Eufemia dell'Aspromonte Messina and West fault N033

EMS-98 Intensity Fields
We reanalyzed the historical sources reported in the CFTI5Med catalog [106] by using the original historical accounts reported in the Supplementary Data S2. We reassessed the intensity of the selected strong earthquakes using EMS-98 [38], which estimates the intensities by considering the buildings' vulnerability class and damage level. This scale has recently been used by many geoscientists (e.g., [37,107,108]) because the Mercalli-Cancani-Sieberg (MCS) scale usually used in the historical catalogs [1,106] does not take into account the building typology. We considered that almost all the buildings were of vulnerability class A and few were of Class B, as described in most of the contemporaneous sources. The use of poor-quality construction materials, often rubble stones, and the widely adopted construction technique known as "a sacco", which used bare stones, poor-quality mortar, and slight stone facades (buildings of vulnerability class A) (e.g., [109]), were usually the cause of the widespread collapse of many buildings. Only a few masonry buildings (simple or massive stone, i.e. buildings of vulnerability class B) were constructed in the towns, especially as noble palaces and religious buildings, such as monasteries and churches (e.g., [109,110]). When it was not possible to estimate the percentage of the damaged buildings and/or the level of damage, the intensity was assigned using uncertain values (i.e., 5-6 to 10-11 EMS-98).
Regarding the 1783-SS, many coeval reports are available. However, the overlap of numerous shocks in a very limited time and space makes interpreting the historical data difficult, since historical accounts describing damage in several localities do not distinguish between the shocks.
To separate the effects of the different shocks, we analyzed the descriptions of the event from single localities and compared the corresponding times that the shocks were felt with the times the shocks were felt at the other sites. After having unraveled a tangle of multiple dates, hours, and spurious connections between several earthquakes, we identified collections of historical data dealing with single seismic events and assessed the intensities considering the cumulated effects of different shocks. In particular, we reanalyzed the events with M w >~6 recorded along the southern Calabria ridge, between the Catanzaro and Messina straits, i.e., the five main shocks of the 1783-SS, the 1791 event located close to the third 1783 event, the 1894 event located close to the first 1783 event, and the 1908 earthquake (see Figures S2-S9 in the Supplementary Materials).

Coseismic Features
In this paper, for the 1783-SS, we compiled an updated database of the coseismic effects (e.g., CEDID [111] and the EEE catalog (http://eeecatalogue.isprambiente.it/viewer.php, accessed on 26 May 2023)), starting from the references reported in the CFTIMed catalog [106] and integrating it with Cucci [112]. We reanalyzed the primary contemporaneous sources and relocated the coseismic effects distinguishing ground fracturing and openings, landslides, liquefaction phenomena, and hydrogeological anomalies (map reported in the Results section). For completeness, we provide the dataset as a shapefile for further use (Supplementary Data S3). The coseismic dataset provides information on the maximum ground shaking and possible coseismic faulting. Most of the coseismic effects, such as ruptures or openings, align along the main faults of the area or at their hanging wall. In the Supplementary Data S2 (see Tables S1-S4), we compiled detailed tables for each seismic event, reporting a description of the earthquake, the coseismic effects reported by historical sources in the original language and translated, all historical sources, and the revisited macroseismic field. There are also other papers that assemble and describe the coseismic effects of the Calabria earthquakes (e.g., [13,20,[113][114][115][116][117]) that used the same historical sources for different purposes, such as paleoseismological research and application of the ESI intensity scale based on earthquake environmental effects (i.e., EEE) [118,119].

Macroseismic Field Partitioning
Based on the area of the maximum macroseismic intensity of an earthquake, it is possible to assume that the hypocenter of the earthquake and the seismic source are located near the most affected locations [36,38]. However, this approximation must consider uncertainties arising from the analysis of the data and historical sources used to reconstruct the macroseismic field associated with the earthquake (lack or uncertainty of historical sources, absence of inhabited locations, local amplification effects, cumulative damage, vulnerability of buildings, etc.) and geological structural complexities, such as the involvement of multiple seismogenic sources.
The recent central Italy 2016 extensional seismic sequence has taught us that individual earthquakes may have a complex source, with synchronous activation of distinct fault segments or different faults. This is the case of the 24 August event (M w 6.2), consisting of two neighboring twin sources [120][121][122], and the 30 October (M w 6.5) event, characterized by a major off-fault deformation that also affected a different structure from the main seismogenic one [6].
To identify potential macroseismic sub-fields representing individual sub-sources during major historical events, a comprehensive analysis of the spatial distribution of the intensity field is essential. Therefore, we have developed a MATLAB code structured as follows:

1.
Intermediate macroseismic intensities within the given macroseismic field (i.e., uncertain classes) are reassigned to the lower class (e.g., intensity points IX-X are assigned to class IX).

2.
Points corresponding to the highest intensity class (I max ) are selected, as indicated by the points within the dashed line in Figure 3a. If there are only a few points, additional classes can be selected.

3.
Using the selected intensity points, the KMEANS function is applied to identify two sub-clusters and their centroids (yellow dots in Figure 3a).

4.
If the number of selected points is ≥10 and the two centroids are more than 10 km apart, the intensity points within the overall macroseismic field are partitioned into sub-fields according to their proximity to one of the centroids (Figure 3b).

5.
If the conditions specified in step 4 are not met, no partitioning of the macroseismic field occurs, and it remains unmodified.
The pursuit of these steps clearly cannot disregard the importance of conducting a firstorder geological assessment to validate the meaningfulness of any sub-clusters that may be statistically obtained. This is an additional and preliminary step that serves as a crucial test to ensure the reliability and significance of the identified sub-clusters. As in other cases in the literature [36,123], our methodological workflow cannot take into account such very detailed complexities (i.e., at the scale of the seismic microzonation [107,124]); it follows a statistical path that does not consider the local effects (well highlighted, for example, by Teves-Costa and Batlló [125]) or other problems arising from the intensity assessment for historical earthquakes; for these reasons, the control exerted by the geoscientist is vital. Specifically, as regards our study area, in which we obtained a partition of the macroseismic field, we verified that there was no correspondence with clear variations of the geological substrate typologies (rock types) on the most up-to-date existing geological maps.
3. Using the selected intensity points, the KMEANS function is applied to identify two sub-clusters and their centroids (yellow dots in Figure 3a). 4. If the number of selected points is ≥10 and the two centroids are more than 10 km apart, the intensity points within the overall macroseismic field are partitioned into sub-fields according to their proximity to one of the centroids (Figure 3b). 5. If the conditions specified in step 4 are not met, no partitioning of the macroseismic field occurs, and it remains unmodified. The pursuit of these steps clearly cannot disregard the importance of conducting a first-order geological assessment to validate the meaningfulness of any sub-clusters that may be statistically obtained. This is an additional and preliminary step that serves as a crucial test to ensure the reliability and significance of the identified sub-clusters. As in other cases in the literature [36,123], our methodological workflow cannot take into account such very detailed complexities (i.e., at the scale of the seismic microzonation [107,124]); it follows a statistical path that does not consider the local effects (well (a) Clustering of I max points of the whole macroseismic field: points with maximum intensity within the dashed line (intermediate classes are downgraded to the lower class) are divided into two clusters and the position of the corresponding centroid is calculated (yellow dots); (b) partitioning of the macroseismic field with respect to cluster centroids: points and squares represent the macroseismic subfields separated with respect to the proximity of the two centroids (indicated by the red arrows), respectively; (c) magnitude and epicenter calculation of the whole macroseismic field (yellow star) and representation of the axes of elongation. Red dashed ellipses are the 80% probability contour calculated using the selected points and squares of the macroseismic subfields; red solid ellipses are the 80% probability contours calculated excluding the outlier points (marked with black crosses); (d) representation of the elliptical seismic sources. Light yellow ellipses represent the elliptical seismic sources of the macroseismic subfield.

Macroseismic Epicenter and Equivalent Magnitude
We obtained new EMS-98 macroseismic intensity data for the largest earthquakes of southern Calabria and, following the Boxer software (Version 4.2.1, 2022) approach [36,123], we estimated the centroid of the largest intensities, assumed to correspond to the earthquake epicenter position (yellow star in Figure 3c). With the same code, we evaluated the equivalent magnitudes of the overall macroseismic field and of the macroseismic subfields where we partitioned it.

Macroseismic Field Elongation Axis
The elongation of a macroseismic field may be calculated by applying different functions (see e.g., [123]). In this study, the macroseismic field's elongation trend was identified by calculating the best-fitting longitudinal and transverse axes delineating an ellipse. To compute the ellipse, we used a MATLAB function (modified from "Statistical Methods in the Atmospheric Sciences"; https://waterprogramming.wordpress.com/2016/11/07 /plotting-probability-ellipses-for-bivariate-normal-distributions/ accessed on 1 February 2023). This function calculates and plots the probability ellipses for bivariate normal distributions. Given a set of points, the code plots the constant probability elliptical contour (i.e., 50%, 80%, 90%, and 99%) within which a given percentage of the distribution lies. The ellipses are constructed with respect to a central point representing the distribution's centroid relative to the averages of the macroseismic point locations used in the probability calculation.
In particular, to compute the direction of the macroseismic elongation axes, we developed a MATLAB code that (1) selects the macroseismic field points that have intensity ≥ I max − 1 (intensity points from 9.5 to 10.5 in Figure 3c); (2) calculates the probability ellipses for bivariate normal distributions, i.e., the constant probability contour (50%, 80%, 90%, and 99%) within which a given percentage of the distribution lies and its centroid; (3) represents the probability contour of 80% (red dashed ellipses in Figure 3c) and selects the points that lie within the ellipse; (4) recalculates the probability ellipses for bivariate normal distributions, as described in step (2), using only reselected points (black crosses in Figure 3c indicate the outlier points excluded from the calculation); and (5) selects the 80% probability contour ellipse and calculates the direction and length of the minimum and maximum axes of the ellipse (red ellipses in Figure 3c).
The orientation of the maximum axes of the ellipses provides valuable constraints into the orientation of the seismic sources responsible for the earthquakes. Such information and the possibility of subdividing a total macroseismic field into statistically constrained sub-clusters must be evaluated in light of the geological and instrumental constraints on the orientation of the potentially seismogenic sources of the area. To test the reliability and validity of the method, we ran it using macroseismic data from the 2009 L'Aquila earthquake of central Italy (M w 6.2), obtaining results comparable to what was observed with instrumental data. As shown in Figure S10 in the Supplementary Materials, the macroseismic axes of the ellipse obtained with our method deviate less than 10 • compared to the trend of the seismogenic structure constrained with instrumental seismology and source models from various authors [84,[126][127][128][129].

Scale Law in the Literature
In the case of historical earthquakes, at best, the available parameters to define the geometry of the fault rupture are the earthquake equivalent magnitude together with possible information on (i) the fault surface rupture length (SRL) from historical reports, (ii) the long-term outcropping extent of possibly seismogenic faults and their average dipangle, and (iii) the thickness of the seismogenic layer and the fault kinematics from local instrumental data. Since the pioneering research by Wells and Coppersmith [53], many authors have tackled the scaling of earthquake faults in different tectonic contexts, relating magnitude (M w ) and seismic moment (Mo) to different parameters, such as SRL, subsurface rupture length and width (L and W), rupture area (A), and coseismic displacement (D). As a result, many empirical laws and corresponding regressions are available in the literature, developed for specific tectonic contexts and with different types of data input. For example, most authors use seismological data, such as the spatial pattern of early aftershocks of seismic sequences and the thickness of the seismogenic layer, or geological data, such as the extension of the outcropping fault and the coseismic rupture at the surface (e.g., [21,[44][45][46]53]). Others use seismic and geodetic data to reconstruct the seismogenic source parameters using finite-fault earthquake source inversion methods ( [43,47,50,51]).
In this paper, we consider a selection of scale laws computed for extensional settings around the world [43,44,[47][48][49][50][51][52][53][54], specifically in Italy [21,50] (Table 2). In Figure 4, we separately represent, in logarithmic scale, the regression curves computed in the literature for scaling the relationships of M w with SRL, L, W, and A. Table 2. List of regression laws from the literature for extensional tectonic regimes. SRL: surface rupture length; L: sub-surface rupture length; W: rupture width; A: rupture area.

SRL (km) L (km) W (km)
A (km 2 ) W&C94 [53] Log ( Table 2. Yellow stars represent an earthquake with Mw 6.5; the vertical black line is the average value of seismic source geometrical parameters (SRL, L, W, A); vertical dashed black lines represent the variability of seismic source geometrical parameters (based on standard deviation values).

Determination of Seismic Source Size
The analyzed scaling laws show the high range of variability of the source size that may be calculated for the same earthquake with a given magnitude. To qualitatively associate a reference value with its range of variability to any magnitude, we prepared a function in MATLAB that returns the average values (SRL, L, W, and A) with their  Table 2. Yellow stars represent an earthquake with M w 6.5; the vertical black line is the average value of seismic source geometrical parameters (SRL, L, W, A); vertical dashed black lines represent the variability of seismic source geometrical parameters (based on standard deviation values).
We have made available a "user-friendly" Excel spreadsheet that allows for the calculation and comparison of all the different parameters (See Supplementary Data S4).
We observed systematic variations in the scaling law values due to the typology and quality of the data input used to perform the linear regression, e.g., geological or seismological data. Generally speaking, scaling laws computed with earthquake data, from both a hypocentral database and source models [45][46][47]51], correspond to the predicted size of the geometric parameters (A, L, W, SRL) for any given M, with respect to the geological data input. Additionally, the extent of the earthquake catalog significantly influences the results.
Most of the M-SRL regression lines compiled in this paper ( Figure 4a) have a similar slope and close positions, except for Pavlides and Caputo [49] and Wesnousky [54], which are significantly different and more conservative because they assign smaller SRL than others for the same M.
The M-L regression lines (Figure 4b) show consistent trends, except for the Martin Mai and Beroza [47] law, which tends to underestimate L. Greater differences are observed at M < 6.0, where the corresponding L shows greater variability for the same magnitude.
The M-W regression lines (Figure 4c) are less common in the literature, and the variability of the results is large for both high and low M.
The M-A regression lines (Figure 4d) all have similar slopes, but a wide range of variability between the lowest and highest values for all the magnitude ranges (i.e., more than a half point in M).

Determination of Seismic Source Size
The analyzed scaling laws show the high range of variability of the source size that may be calculated for the same earthquake with a given magnitude. To qualitatively associate a reference value with its range of variability to any magnitude, we prepared a function in MATLAB that returns the average values (SRL, L, W, and A) with their standard deviation between the highest and lowest values. As an example, in Figure 4, we represent an M w 6.5 earthquake (yellow star) projected on the regression lines within the space delimited by the vertical dotted lines representing the computed standard deviation.
Our specific aim in this paper is to determine the size of the herein investigated Calabria historical earthquakes starting from the newly computed equivalent magnitudes, considering both the total and the partitioned sources.
Once the extensional kinematics from independent long-term geology and instrumental seismology were determined, we derived the average A, L, and W parameters and their variability (represented by the standard deviation) from the regression laws previously discussed.

Map-View Earthquake Source Representation
Most commonly, individual earthquake sources are represented in map-view as rectangular boxes (e.g., [81]), dimensioned by L and width surface projection (Ws). In order to calculate the Ws, we assume a dip-angle of 60 • is suitable with extensional settings and fault/slip data available for the southern Calabria region [20,26,73,75,93,94,130,131].
An alternative graphical representation we propose in this paper is that of ellipticalshaped sources, which better approximate the physical model of the seismic rupture propagation on the fault plane [132], also highlighted by earthquake interferometry (e.g., [133][134][135]). The ellipse must have the same areal extent as the rectangular source and an ellipticity controlled by a long-to-short axes ratio equal to the L/W ratio, directly derived from the scale laws. The following equations summarize the steps described above: Surface semi-minor axis = b cos 60 • Based on these formulas, we developed a MATLAB code that calculates the average elliptical seismic sources oriented and centered with respect to the ellipse representing the elongation of the macroseismic field (Figure 4d). The devastating 1783-SS, which started on 5 February 1783 with a M w 6.8 event, was followed by a long series of shocks that migrated NNE-ward from the Messina Strait, along the axis of the Calabrian Arc. More than 180 localities were nearly destroyed and there were more than 30,000 victims. The aftershocks lasted for about ten years. There were five earthquakes in 2 months with I max X-XI EMS and M w between 6.0 and 6.8. The sequence caused a geomorphological crisis that changed the Aspromonte and Serre Mts. landscape ( Figure 5), triggering ground-fracturing, landslides, and liquefactions, and inducing a significant change in the fluvial network, with deviation and damming of the fluvial channels, as well as the formation of numerous lakes [13] (see Table S1 in Supplementary Data S2 and Supplementary Data S3).  [106] integrated with [112,113]). Openings and fractures, soil liquefaction, hydrological anomalies, and landslides are reported on the same structural map as in Figure 1a. The shapefile of coseismic effects is made available in Supplementary Data S3. Figure 5. Updated map of the coseismic effects of the 1783-SS (data from the CFTIMed catalog [106] integrated with [112,113]). Openings and fractures, soil liquefaction, hydrological anomalies, and landslides are reported on the same structural map as in Figure 1a. The shapefile of coseismic effects is made available in Supplementary Data S3.

Results
The area with the strongest effects of the 5 February shock is located on the Tyrrhenian side of the Calabrian Apennines (Figures 6a and S2). The effects of this earthquake brought a few changes. The maximum intensity assigned is X-XI I EMS for two reasons: (i) to take into account the buildings' vulnerability and (ii) to account for the cumulated effects of more shocks in a short time. Most of the buildings were of vulnerability class A and a few were class B, and some buildings were left standing, indicating that the intensity was certainly X (most buildings of vulnerability class A sustained damage of grade 5 (collapse), but it could have been intensity XI (destruction: very heavy structural damage of most buildings of vulnerability class B) [42]. Without knowledge of the precise number of buildings, we decided to use this uncertainty. Mainly for historical earthquakes, the macroseismic scales, including EMS-98 from degrees X to XII, demonstrated saturated levels of damage, losing their diagnostic value (see [118]). Furthermore, most of the accounts report that the data were gathered after that all the sequences had occurred, making it possible that the destruction was the result of the cumulated effect of more shocks. Moreover, in the area of the Strait of Messina on both the Sicilian and Calabrian sides, the intensities of some localities with higher values than the neighboring ones were decreased (Table S5 in Supplementary Data S2). These localities have been added to the earthquake of 6 February (Figures 6b and S3). In fact, this earthquake is known above all for the tsunami that struck Scilla and Torre Faro [136]. Nevertheless, several historical sources (see Table S1 in Supplementary Data S2) describe this shock as being as violent as the previous one, but of shorter duration than the first. This earthquake mainly caused damage in Messina and some nearby towns and villages, as well as in Reggio (I EMS = VIII-IX). Since this shock occurred about twelve hours after the initial shock, information was reported for only a few localities. By using the new intensity values (Table S6 in Supplementary Data S2), the epicenter is moved southernmost with respect to the previous one [103].
The next day, on 7 February, the earthquake moved NE, towards the Mesima valley (Figures 6c and S4) on the Tyrrhenian side of the Serre chain (I EMS = X-XI). Most of the localities (Table S7 in Supplementary Data S2) were already damaged by the previous shocks.
On 1 March, at 08.30 Italian time (01.40 GMT), an earthquake was felt across a wide region between Messina and the northern border of the Calabria region, with the area of the maximum effects (I EMS = VIII-IX), which moved further north (Figures 6d and S5) along the Apennines. The most affected villages were Polia (now known as Trecroci), Poliolo, Castelmonardo, and Francavilla Angitola (Table S8 in Supplementary Data S2). The 28 March earthquake (X-XI I EMS ) was felt throughout southern Italy. The most affected area was in the Catanzaro basin, with damage both on the Ionian and the Tyrrhenian sides (Figures 6e and S6). This is the northernmost greatly damaged area of the entire seismic sequence (Table S9 in Supplementary Data S2). The extended area of damage was likely due to the cumulated effects of more shocks; indeed, on the Tyrrhenian side, some localities were already badly damaged by another minor shock that occurred on 28 February, as well as the 1 March earthquake.

1791 Central Calabria, 1894 Southern Calabria
The 13 October 1791 (M w 5.9) earthquake damaged the area hit by the 7 February 1783 earthquake (Figure 7a and Figure S7). The seismic sequence lasted until 20 October. The I max was VIII-IX I EMS and the estimated magnitude was equal to 5.89. There were 15 victims. The most damaged towns were Spadola, Simbario, Serra San Bruno, and Soriano (Table S10 in Supplementary Data S2). The buildings that were impacted in the last village were strongly conditioned by the previous earthquakes. Geosciences 2023, 13, x FOR PEER REVIEW 16 of 35  (Figures 7a and S7). The seismic sequence lasted until 20 October. The Imax was VIII-IX IEMS and the estimated magnitude was equal to 5.89. There were 15 victims. The most damaged towns were Spadola, Simbario, Serra San Bruno, and Soriano (Table S10 in Supplementary Data S2). The buildings that were impacted in the last village were strongly conditioned by the previous earthquakes.  The 16 November 1894 (Mw 6.1) earthquake affected southern Calabria and eastern Sicily and had destructive effects in an area of 80 km 2 between the Aspromonte plains and the Tyrrhenian coast (Figures 7b and S8). San Procopio and Sant'Eufemia d'Aspromonte were the most damaged villages (IEMS = IX, Table S11 in Supplementary Data S2), with many buildings collapsed totally or in large parts. In 17 other localities in the province of Reggio Calabria (including Bagnara Calabra, Palmi, and Seminara), the earthquake caused a few collapses (see Table S3 in Supplementary Data S2), but most of the buildings suffered serious damage, structural failures, and partial collapses. In Messina and Reggio Calabria, many buildings were damaged, some heavily. There were two foreshocks before the mainshock and several aftershocks until 16 June 1895. The Imax was IX IEMS and the Mw was 6.1. There were 96 victims.

1908 Messina Straits Earthquakes
The 28 December 1908 earthquake (Mw 6.8) (Figures 7c and S9) and subsequent tsunami devastated southern Calabria and northeastern Sicily. The double catastrophe destroyed Messina, Reggio di Calabria, and dozens of nearby coastal towns (Table S12 in The 16 November 1894 (M w 6.1) earthquake affected southern Calabria and eastern Sicily and had destructive effects in an area of 80 km 2 between the Aspromonte plains and the Tyrrhenian coast (Figures 7b and S8). San Procopio and Sant'Eufemia d'Aspromonte were the most damaged villages (I EMS = IX, Table S11 in Supplementary Data S2), with many buildings collapsed totally or in large parts. In 17 other localities in the province of Reggio Calabria (including Bagnara Calabra, Palmi, and Seminara), the earthquake caused a few collapses (see Table S3 in Supplementary Data S2), but most of the buildings suffered serious damage, structural failures, and partial collapses. In Messina and Reggio Calabria, many buildings were damaged, some heavily. There were two foreshocks before the mainshock and several aftershocks until 16 June 1895. The I max was IX I EMS and the M w was 6.1. There were 96 victims.

1908 Messina Straits Earthquakes
The 28 December 1908 earthquake (M w 6.8) (Figures 7c and S9) and subsequent tsunami devastated southern Calabria and northeastern Sicily. The double catastrophe destroyed Messina, Reggio di Calabria, and dozens of nearby coastal towns (Table S12 in Supplementary Data S2). About 90% of the buildings in Messina were destroyed [131], with the worst damage in the central and northern parts of the city, which were built on soft soil and, therefore, were possibly susceptible to liquefaction. Damage was also reported as less severe in the western part of the city, particularly for structures built on more compact terrain, such as around the Gonzago Castle (Table S4 in Supplementary Data S2). The large number of damaged buildings highlighted the vulnerable nature of the building stock in Messina at that time [109]. Poor-quality construction materials were blamed for the widespread collapse of many buildings. Buildings constructed with better quality materials or practices were less prone to collapse [137]. The main event lasted 30-40 s with I max = X-XI I EMS and an M w = 6.83 calculated with the Boxer code. The aftershocks lasted until 1913, and several shocks probably reached an intensity of VIII I EMS , but it is not easy to estimate the intensity when most of the buildings were nearly destroyed. Considering the historical reports, Baratta [109] identified two areas defined as "centers of shaking", primary and secondary, respectively, obtained through observations on the direction of the shock. Observing the seismogram, Oddone [138] distinguished a second shock, which occurred five minutes after the first one. Rizzo [139] observed that the earthquake was characterized by a first pulse NNE-SSW orientation, followed by a stronger shock, ESE-WNW directed, and a third, even stronger series of vibrations along the NNE-SSW direction. This description would suggest the occurrence of multiple events related to the activation of different faults or fault segments. The number of casualties due to the Messina earthquake remains uncertain (between 60,000 and 80,000) [109].

Earthquake Source Parameters
Starting from the EMS-98 revised macroseismic fields and following the previously described methodological steps, for each earthquake, we analyzed and/or computed the following:

•
Tectonic trends of neighboring Quaternary tectonic structures; • Macroseismic field partitioning; • Equivalent magnitude and macroseismic epicenter; • Direction and length of the elongation axes of the macroseismic area; • Scale law-derived source size according to regression laws; • Elliptical seismic sources parameters.
Using the Boxer program (version 4.2.1, 2022), we first computed the equivalent magnitude and location of the macroseismic epicenters for all the newly estimated EMS-98 macroseismic fields [123]. Subsequently, we used the KMEANS algorithm in MATLAB to identify the presence or lack thereof of intensity point clusters (see Section 3.4.1), and in the case of a partitioned field, we recomputed the equivalent magnitude from Boxer. We specify that the equivalent magnitudes of the partitioned fields are not to be intended as sub-seismic events, but are necessary for calculating the respective seismic source's size.
For all the identified fields, following the previously explained procedural steps, we calculated the directions and sizes (length and area) of the seismogenic sources and represented them in a GIS platform with both an elliptical and a rectangular map-view shape (Supplementary Data S5). In Table 3, the earthquake source parameters are also compared with the corresponding values in the CPTI15 v4.0 parametric catalog.
Subsequently, we proceeded to validate the seismic source parameters inferred from the macroseismic data by comparing them with the tectonic trends of the potentially associated faults. This comparison enabled us to assess the consistency and compatibility between the observed macroseismic patterns and the known tectonic characteristics of the fault systems in the region.
A single elliptical source was evident for most events, except the first event of the 1783 sequence and the 1908 earthquake. We refer the reader to Table 3 for details on the earthquakes with the single sources computed. Table 3. Earthquake parameters. Key: I 0 : epicentral intensity (calculated as in [36]); M w : equivalent magnitude with error; MDP: macroseismic data points; 1 MDP used on MDP selected; 2 axis parameters obtained using selected points and unselected points (in the brackets); 3 value ± the standard deviation; 4 surface width calculated considering a seismic source with a dip angle of 60 • . For the 5 February earthquake (M w 6.8), the code identified two clusters with centers 17 km apart within the Gioia Tauro plain, both consisting of intensity points X and X-XI. This allowed us to divide the entire macroseismic field into two sub-fields: a northern one with 215 points and a southern one with 152 points ( Figure S11). For the northern macroseismic field, the algorithm selected 22 points with intensities between X-XI and IX-X; for these points, a maximum elongation axis of~23 km in the direction of N25 • was determined (Table 3 and Figure S13a). On the other hand, the southern macroseismic field has an axis of maximum elongation 13 km long and oriented N59 • , which was also calculated with 22 points selected (Table 3 and Figure S13a). Based on these data, we calculated, for the northern macroseismic field, a magnitude of 6.56 and an elliptical source with an average size of 392 km 2 and an average length along the major axis of 27 km, oriented NNE-SSW. For the southern macroseismic field, the equivalent magnitude is~6.6, and the resulting seismic source has an average size of 400 km 2 and an average length of 27 km, rotating along the Aspromonte chain to NE-SW (Figure 8a). For the 5 February earthquake (Mw 6.8), the code identified two clusters with centers ~17 km apart within the Gioia Tauro plain, both consisting of intensity points X and X-XI. This allowed us to divide the entire macroseismic field into two sub-fields: a northern one with 215 points and a southern one with 152 points ( Figure S11). For the northern macroseismic field, the algorithm selected 22 points with intensities between X-XI and IX-X; for these points, a maximum elongation axis of ~23 km in the direction of N25° was determined (Table 3 and Figure S13a). On the other hand, the southern macroseismic field has an axis of maximum elongation 13 km long and oriented N59°, which was also calculated with 22 points selected (Table 3 and Figure S13a). Based on these data, we calculated, for the northern macroseismic field, a magnitude of 6.56 and an elliptical source with an average size of 392 km 2 and an average length along the major axis of 27 km, oriented NNE-SSW. For the southern macroseismic field, the equivalent magnitude is ~6.6, and the resulting seismic source has an average size of 400 km 2 and an average length of 27 km, rotating along the Aspromonte chain to NE-SW (Figure 8a). For the 28 December 1908 earthquake (M w 6.8), the code identified two clusters located north and south of the macroseismic epicenter with centers 14 km apart, both consisting of intensity points X and X-XI. This allowed us to divide the macroseismic field into northern and southern subfields of 497 and 317 points, respectively ( Figure S12). The total equivalent magnitude associated with this earthquake is~6.8, while, considering the separate partial macroseismic fields, we obtained magnitudes of 6.6 and 6.5 for the northernmost and southernmost fields, respectively. The estimated northern elliptical seismic source is~ENE-WSW oriented, while the southern one is~NNE-SSW, with an average area of 427 km 2 and 345 km 2 and a length of 28 km and 25 km, respectively (Figure 9c). For the northern source, the algorithm selected 51 points from the macroseismic subfield with intensities between X-XI and IX-X, resulting in a maximum elongation length of 23 km and a N80 • trend. For the southern source, the macroseismic subfield can be described as having a length of 17 km, oriented N20 • , using 35 selected points (Table 3 and Figure S14c). The obtained values can be underestimated because the areas are separated by the sea; therefore, the macroseismic field is not complete. Among the analyzed fields, especially interesting is the 28 March 1783 earthquake, with a recalculated magnitude of ~6.8 and located in the middle of the Catanzaro basin Among the analyzed fields, especially interesting is the 28 March 1783 earthquake, with a recalculated magnitude of~6.8 and located in the middle of the Catanzaro basin (Figure 6e). Based on 23 macroseismic field points with intensities between X-XI and IX-X, the algorithm calculated an axis of maximum elongation of 25 km, N26 • -oriented (Table 3 and Figure S13e). We obtained an elliptical seismic source with an average area of 715 km 2 and an average length of 37 km, which crosses the Catanzaro basin and aligns approximately NNE-SSW with the Serre chain (Figure 8e). For a comparison of the seismic sources described here and trends of the tectonic structures shown in Figure 1, see Section 4.3.

Earthquake Sources Quality Parameters
For quality control, we assigned the quality parameters to the seismic sources obtained herein. The first quality parameter (defined as Qm in Table 4) indicates the degree of variability of the possible direction of the seismic source. This is calculated as the difference between the directions of the preliminary and final elongation axes (Figure 3c). Quality A corresponds to an angular difference ≤20 • , B if between 20 • and 40 • , C if between 40 • and 60 • , and D if >60 • . The second quality parameter (Qε in Table 4) considers the degree of ellipticity (ε) in terms of the shape ratio (the higher the ellipticity, the higher the quality); Quality A corresponds to ε ≥ 2, B if between 1.5 ≤ ε ≤ 2, and C if ε < 1.5. The third parameter describes the comparability with the seismogenic structures present in the area (see Figure 1 and Table 1). This parameter, defined as Qf in Table 4, corresponds to the difference between the direction of the seismic sources obtained in this paper and the faults' trends of Table 1. In addition, in this case, quality A corresponds to an angular difference ≤20 • , B is between 20 • and 40 • , C is between 40 • and 60 • , and D is >60 • . As shown in Table 4, in most cases, we obtained high qualities (A and B), while we obtained medium-low qualities for the parameter Qm related to the 6 February 1783 event and Qf related to the northern sub-event of 28 December 1908. In the first case, this can be attributed to the presence of the Messina strait. In the second case, the low quality is due to the strike variability of the Messina fault and West fault (from NNE-SSW to ENE-WSW). Table 4. Earthquake seismic sources' quality parameters. (1) Direction and (2) shape ratio of major and minor axes lengths of the preliminary macroseismic field elongation calculated in step 3 of Section 3.4.3 (see also Figure 3c); (3) direction and (4) shape ratio of major and minor axis lengths of the final (after outlier points exclusion) macroseismic field elongation calculated in step 5 of Section 3.4.3; ID as in Table 1; (5) mean direction of fault systems, as in Table 1; (Qm) quality of the macroseismic field elongation computed as the difference between (1) and (3); (Qe) quality of the ellipse computed as the ratio between the major and minor axes of the ellipse; (Qs) quality of the seismic source obtained comparing its direction (3) with the directions of the fault system on the area (5); see main text for description of qualities.

Comparison of EMS-98 Intensities and MCS with Implications for Source Size
When studying historical earthquakes, the use of different seismic intensity scales can result in substantially different assessments of the damage that structures sustained as a result of shaking. Unlike the MCS scale, the EMS-98 scale classifies residential buildings into six vulnerability classes and assesses their damage distribution in five classes. This leads to a more specific attribution of earthquake damage among different buildings in a single location, and also takes into account the possible differences in construction techniques among different residential centers in a given earthquake-affected area [107]. The MCS intensity assignments are generally higher than EMS for the same initial data [140].
In this study, we observed an apparent reduction in the value of macroseismic intensities between evaluations with the EMS-98 scale and the macroseismic fields present in CPTI15, where the MCS scale is used. For the earthquakes analyzed, the epicentral intensity (I 0 ) shows a reduction of half a degree, except for the events of 7 February 1783 and 16 November 1894 (Table 3). Such a difference, which is due to the use of an uncertain intensity value (e.g., X-XI), together with the review of historical reports and the re-attribution of the recorded damage to the respective earthquake events and sub-events, leads to the calculation of equivalent magnitudes lower than those generally attributed in the literature. In all eight earthquakes analyzed, we obtained a reduction in the equivalent magnitudes of even a few tenths, especially for the higher energy events (5 February 1783 and 28 December 1908), as shown in Table 3. Although this reduction may seem insignificant, an increase of 0.2 additional magnitude points corresponds to a doubling of the energy released by the earthquake [141]. This then affects the size of the seismic sources calculated by the scaling laws from the relative magnitude values (as shown in Section 3.5.2). In fact, a 0.2 increase in the equivalent magnitude results in a seismic source with larger length and width values of about 1.3 and 1.2, respectively, as well as an increase in the source area by a factor of about 1.5. We are aware that, especially for pre-instrumental earthquakes, representing the magnitudes with two decimal places is excessive, but in this way, we maintained the correct relationship with the sizes of the sources. Similar results were obtained by Blumetti et al. [113] for the 5 February 1783 earthquake and by Comerci et al. [115] for the 1908 earthquake using the ESI scale [118]. For the 5 February 1783 earthquake, from the extent of normal surface faulting (rupture length about 35 km; maximum displacements at least in the order of 0.8 m) and the total area affected by secondary effects (about 3500 km 2 ), Blumetti et al. [113] estimated an ESI epicentral intensity equal to X. Possibly due to local characteristics of the territory and to dense secondary effects, such as landslides and liquefactions, in four localities (Santa Cristina di Aspromonte, San Giorgio Morgeto, Molochio, and Oppido Mamertina), the ESI local intensity values reached an intensity degree = XI, higher than the epicentral intensity [113]. For the 1908 earthquake, the MCS local intensities were often larger than the corresponding ESI values obtained by Comerci et al. [115], especially in the near-field area, at the intensity interval IX-XI. Indeed, the earthquake severity estimated by the MCS intensities in the highest degrees might be influenced by both the poor quality of the buildings and the weakness of the buildings already damaged by the earthquakes of 1894, 1905, and 1907 that occurred in the same area [137,[142][143][144][145], as suggested by the authors.

Seismotectonic Implications and Earthquake/Fault Association
Although most of the southern Calabria outcropping active faults dip westward, for all the destructive Calabrian earthquakes, a discussion is open in the literature about the E-dip or W-dip direction of the extensional seismogenic major fault or, rather, the activation of a strike-slip transfer fault. The E-dip model for the 1783 earthquakes fits well with the geochemical evidence from local springs and well observations [82], but it contrasts with the long-term and coseismic features. In fact, the 5 February earthquake caused coseismic surface ruptures along at least 23 km of the Cittanova fault at the base of the mountain slope, between the localities of Santa Cristina and Cinquefrondi [22,87,110]. The same fault shows evidence of Holocene faulting in paleoseismic analyses [22,26].
In this paper, we have provided more detailed and complete information on the coseismic deformation that, according to our reconstruction, involves both the west-dipping Cittanova and the Serre faults.
The sources that we reconstructed for the 5 February event fit with the attribution of the Cittanova fault for the northern portion, as previously proposed by [22,26], and of the Cittanova and Sant'Eufemia d'Aspromonte faults for the southern portion [72,78,88] ( Figure 8a). Our results fit with the models that associated the source of this quake with W-dipping faults. Similarly, for the 7 February event, we propose the attribution to the W-dipping Serre fault [93,94] or, alternatively, to an easternmost alignment, which also fits with our seismogenic source (see Figure 8c).
For the 6 February event, our results highlight the possibility of the involvement of a W-dipping source located on-land in Calabria, possibly the Reggio Calabria or the Armo sources (Figure 8b), while for the 1 March event, we advance the hypothesis of the activation of the central segment of the Serre fault system (Figure 8d).
The source of the 1791 event would fit with the central-easternmost portion of the Serre fault (Figure 9a). Based on our results, the 1894 earthquake could be related to the NNW-dipping portion of the Sant'Eufemia fault (Figure 9b).
In the case of the 28 March 1783 earthquake, we also highlight the difficulty in identifying a well-constrained source. Nevertheless, considering the best-fitted NNE-SSW direction of the revisited macroseismic field and the macroseismic ellipse's position at the Serre fault's footwall (Figure 8e), we advance two alternative hypotheses. One sees the possible involvement of the NW-dipping Caraffa fault system as a NNE-ward continuation of the Serre fault [81]. The second is a new suggestion that speculates on the possibility of its association with a basal portion of an ESE-dipping extensional detachment fault, developed along with the down-dip prosecution of the Gulf of Sant'Eufemia SE-dipping major fault responsible for the 8 September 1905 earthquake (M~7) [143] (Figures 8e and 10). This source recalls, although with a different magnitude range, the L'Aquila 2009 earthquake of central Italy, which, two days after the release of the main event (M w 6.2) at upper crustal depths on the SW-dipping Paganica fault, released a third major event (M 5) along the east-dipping basal detachment at the footwall of the main fault and about 5 km deeper [127].
Our analysis of the macroseismic field of the 1908 earthquake shows that the seismogenic source derived from the complete macroseismic field has an overall N-S direction and a vast extended area that does not fit with the outcropping extensional fault pattern in the area (Figures 1 and 7c). The difficulty in identifying a likely source is evident when looking at the wide range of proposed solutions that largely variates in trend, size, and position (see Figure 2d and Supplementary Data S1). This could be also related to the fact that the epicenter area falls offshore in the Strait of Messina, leading to the absence of macroseismic data inside the most affected area. developed along with the down-dip prosecution of the Gulf of Sant'Eufemia SE-dipping major fault responsible for the 8 September 1905 earthquake (M~7) [143] (Figures 8e and  10). This source recalls, although with a different magnitude range, the L'Aquila 2009 earthquake of central Italy, which, two days after the release of the main event (Mw 6.2) at upper crustal depths on the SW-dipping Paganica fault, released a third major event (M 5) along the east-dipping basal detachment at the footwall of the main fault and about 5 km deeper [127].  The possibility of interpreting the event as a double rupture arises when considering the herein analyzed configuration of the macroseismic field that may be statistically partitioned into two confining clusters, one on the Sicilian side and the other on the Calabrian side (Figures 7c and S12). In the light of this observation, and in consideration of the events that occurred two five minutes apart, as recognized by Oddone [138] in the seismogram and by Baratta [109] in the two "centers of shaking" (see areas a and b in Figure S15), we advance the hypothesis that the 1908 earthquake either activated in sequence two fault segments or nucleated from a central hypocenter and simultaneously propagated in different directions, activating two distinct faults. Taking into consideration both the macroseismic field geometry reconstructed in this paper and the on-shore and off-shore fault pattern recognized in the literature (Figure 1 with references), two new alternative scenarios may be added to the previously proposed ones: (1) he earthquake activated a SSE-dipping fault segment offshore of Messina and an antithetic WNW-dipping segment possibly corresponding to the Armo fault and its northward continuation; (2) It activated the SSE-dipping fault segment offshore and a synthetic ESE-SSE-dipping offshore fault corresponding to the West fault of Barreca et al. [74] (Figure 9c). These two faults may be activated by the normal slip on the foreland-dipping discontinuity (the West fault [74]), inducing additional stress and promoting failure in the overlying brittle faults.
We recall that the possible activation of two opposite dipping faults has already been hypothesized by Mulargia et al. [105] and by Barreca et al. [74], even if the sources in these papers partially differ in direction, position, and size. Furthermore, the involvement of two synthetic west-dipping faults (Armo and Sant'Eufemia di Aspromonte) has been suggested by Aloisi et al. [146]. Earthquakes generated on off-faults, which branch from the main fault or form in the surrounding medium due to stress concentrations around the rupture tip, have been recognized in instrumental recording times, as in the case of the Norcia 2016 earthquake [6,121]. Off-faults may play a key role in accommodating the large slip and enhancing the seismic radiation of this event, nearly simultaneously activating neighboring structures. Evidently, the new earthquake source associations proposed herein for the 1908 earthquake are preliminary hypotheses that need to be tested and validated, especially with respect to the available geodetic data.

Time-Space Evolution of Southern Calabria Earthquakes
In historical times, but not in instrumental times, Calabria has been affected by numerous seismic events (Figure 1b), in particular, by at least 40 earthquakes ( Figure S16) with M > 5.5, according to the information available from the CPTI15 catalog since 1172 (Table S13). For these events, we studied the time-completeness using a magnitudefrequency method [147] and compared the cumulative frequency with time (in years). With the scatter diagram obtained, we estimated the time-completeness for events with M > 5.5 for data fitted by a straight line [148], i.e., starting from 1609 ( Figure S17). To analyze the time-space distribution of seismicity in southern Calabria, we projected the elliptical-shaped source of all the events with M w > 5.5 released since 1609 along the trace of a vertical ENE-WSW transect ( Figure 10).
The size of the earthquake sources, which are not analyzed in this paper, was calculated considering the magnitude values in the CPTI15 catalog ( Figure S18; Table S14) and scaling them according to the dataset of regression laws for seismic source parameters, as shown in Paragraph 3.5, while for the eight earthquakes studied in this work, we used our results ( Table 3). As a simplification, we assume that the considered earthquakes are associated with a single seismic source whose preferred direction is~N30 along the section profile (Figures 10 and S18).
From the space-time diagram, at least three earthquake cycles that inter-correlate everỹ 125 years are observed ( Figure 10). These are characterized by large seismic sequences that have released considerable energy across the overall region, involving most of the main fault systems (i.e., the Serre-Mileto fault system, the Armo-S. Stefano-Delianuova-Cittanova fault system, and the S. Eufemia dell'Aspromonte-Reggio Calabria-Messina fault system from NE to SW). The diagram also shows a seismic gap in the southernmost part of Calabria and the Strait of Messina from 1600 to 1749, with a lack of earthquakes of significant magnitude (M w ≥ 5.5).
Considering the earthquakes studied herein, they all fall into the two most recent cycles (cycles 2 and 3 in Figure 10). Cycle 2 consists almost entirely of the events of the 1783-SS, on which the area activated by the 1791 earthquake is superimposed. These earthquakes activated a total of~125 km of faults along the strike, not considering the overlapping portions. The events mentioned took place about 1.5 centuries after cycle 1, which activated the northern~75 km of the study area in~50 years, leaving the southerñ 100 km substantially inactive. Finally, cycle 3 developed between the end of the 1800s and the beginning of the 1900s, activating~85 km of faults in a period of~50 years, both in the central-northern portion of the study area and in the southern one, divided by a gap of 25 km and with the 1894 and 1908 earthquakes both in the southern portion.
A similar consideration can be made for the cumulative slip rate distribution, derived from the value of the seismic moment using Hanks and Kanamori's [141] formulas (see lower panel of Figure 10). The displacement values obtained are shown in Supplementary  Table S14. These values have been used to plot the cumulative displacement for each cycle (in red) and the total cumulative displacement (in blue) in the lower panel of Figure 10. In the latter, Cycle 1 shows a maximum cumulative displacement of~2.3 m, which affected the northern area, while Cycles 2 and 3 have a more homogeneous distribution of slip along the transect with a maximum displacement of~1.7 m. The total cumulative displacement curve shows two different patterns of vertical movement related to the northern and southern areas, where the former has a maximum cumulative slip of 5.5 m, and the latter has a slip of 2.5 m with a slip deficit area (~1.5 m) corresponding to the Gioia Tauro basin. The variability between the two sectors can be attributed to several factors: when working with historical information, we oftentimes have incomplete data and uncertainty of the available data; therefore, the estimates of magnitude can be affected by errors. It must also be taken into account that during cycle 1 (1600-1720), in the area between the Strait of Messina and the Gioia Tauro basin, there were at least 12 earthquakes with magnitudes between 4.1 and 5.1, which we excluded from the analysis because the completeness magnitude we used is 5.5.

Conclusions
Southern Calabria has experienced two of Italy's most powerful and destructive earthquakes: the 1783 sequence and the 1908 single event. However, the seismotectonic interpretation of these events remains a subject of extensive debate in the scientific community due to the limited availability of energetic earthquakes in instrumental times.
Our paper aims to provide valuable insight and perspectives from a modern seismotectonic approach to the analysis of macroseismic data. It comprises the following key contributions: (1) We have revised and documented the macroseismic fields of the previous events along with those of neighboring earthquakes (Supplementary Data S2) using the EMS-98 scale. This revised information is made accessible and accompanied by an updated catalog of the 1783-related coseismic effects (Supplementary Data S3). (2) We have developed a methodological approach to extract geometric information on seismogenic sources from the macroseismic field. This approach combines the conventional Boxer method with new statistical tools and with scale-law algorithms.
Our findings indicate that the average trends of the main seismogenic faults in the area (see Table 1) are consistent with the orientations of the elliptical sources obtained in our study (Table 3). The 1908 earthquake appears to be a single event with the simultaneous activation of a major WSW-dipping fault on the Calabrian side and a SSE-dipping fault on the Sicilian side. Alternatively, the SSE-dipping source may have acted in conjunction with the synthetic ESE-SSE-dipping West fault described by Barreca et al. [74]. As for the 28 March final event of the 1783 sequence, we propose its nucleation on a SE-dipping detachment located at the footwall of the NW-dipping complex multi-segmented network along the Serre-Cittanova-Delianuova west-dipping normal fault alignment responsible for the previous events of the 1783-SS.
In conclusions, it is important to note that these hypotheses are preliminary and require further investigation and validation, particularly in light of the available geodetic data and earthquake fault modeling. This ongoing scientific investigation is crucial for a better understanding of seismic hazard assessment and developing effective strategies for mitigating earthquake risks in southern Calabria and similar tectonic regions.  (Tables S13 and S14).