The Calabrian Arc (CA hereinafter), Southern Italy (Figure 1A), is an active structural domain where high-intensity historical and instrumental earthquakes with magnitude higher than 7 and MCS intensities up to 11 (Postpischl, 1985; Boschi et al., 1995; Barbano et al., 2005; Locati et al., 2022; Rovida et al., 2022) were recorded. In particular, some earthquakes in southern Calabria such as the seismic sequence of 1783, with the mainshock on February 5 (M 6.9–7.1; Jacques et al., 2001), the earthquake of 8 September 1905 (Mw 7.5; Presti et al., 2017), and the event of 28 December 1908 (M 7.1; Barbano et al., 2005) (Figure 1A) have been regarded as the strongest seismic events of the Italian Peninsula.

Although the seismotectonics of the CA has been explored through numerous studies (e.g., Westaway, 1993; Tortorici et al., 1995; Monaco and Tortorici, 2000), the seismogenic sources of major earthquakes occurred in the area are still debated. For instance, different interpretations were proposed regarding the location and the geometries of the causative faults responsible for the 1783 historical seismic sequence.

According to DISS (Database of Individual Seismogenic Sources, DISS 3.3.0, https://diss.ingv.it/), major extensional seismogenic sources of southern Calabria can be grouped into two categories: types A1 and A2 (Figure 1A). Type A1 consists of E-dipping normal faults bounded by the west of the Strait of Messina and the Gioia Tauro and Mesima basins (see also Valensise and D’Aaddezio, 1994; Cucci et al., 1996; Peruzza et al., 1996; Pizzino et al., 2004; Loreto et al., 2013). Type A2 structures consist of W-dipping normal faults limiting the same Quaternary basins to the east (see also Tapponnier et al., 1987; Ghisetti, 1992; Tortorici et al., 1995; Monaco and Tortorici, 2000; Jacques et al., 2001; Galli and Bosi, 2002), at the foot of the Aspromonte and Serre mountains. Unlike the Type 2 faults (i.e., Cittanova and Serre faults), which are well exposed on the surface, Type 1 (i.e., Gioia Tauro and Mesima faults) do not exhibit a clear morphological expression.

In this work, a multidisciplinary approach has been followed to reconstruct, in the MOVE Software Suite environment (granted by Petroleum Experts Limited; www.petex.com), the 3D geometry of the Cittanova fault (CF) and Serre fault (SRF), the longest tectonic structures occurring in southern Calabria (Tapponnier et al., 1987; Ghisetti, 1992; Tortorici et al., 1995; Monaco and Tortorici, 2000; Jacques et al., 2001; Galli and Bosi, 2002). These faults, belonging to the A2 category, are defined by Jacques et al. (2001) and recently, starting from macroseismic data, by Andrenacci et al. (2023) as the seismogenic sources of the 5 and 7 February 1783 seismic events, whereas other authors (e.g., Valensise and D’Aaddezio, 1994; Cucci et al., 1996) consider the Gioia Tauro and Mesima faults, belonging to the A1 category, as the most likely sources for the same seismic events. The open debate between different geometric interpretations (E-dipping vs. W-dipping) reported in DISS and different schools of thought on the sources of the 1783 seismic sequence proves that the development of a 3D model for these faults may be crucial to better characterize the seismotectonics of the area and evaluate its seismic potential. Being characterized by clear surface evidence and seismic cluster at depth, we choose the CF and SRF for our analysis. A structural survey was carried out in order to investigate the attitude of these faults and collect kinematic data. Then, a 3D model for these faults was built from their expression at surface, defined by geomorphological and structural analyses and with the available literature, and using seismological data at depth.

The geometric parameters built through the 3D modelling were then validated using the available empirical relationships (Wells and Coppersmith, 1994; Leonard, 2010). Moreover, we tested the kinematic effects of the modelled planes, simulating a uniform slip distributed on fault planes, according to the coseismic effect recognized at the base of the CF scarplet by de Dolomieu (1784) after the 5 February 1783 event. The expected total vertical displacement and its spatial distribution associated with the fault activity were also obtained to compare the current height of the fault scarps with that of the modelled faults. Finally, the potential reactivation of the modelled planes was evaluated using the slip tendency analysis method (Morris et al., 1996).

Geodetic data, measured at several GNSS stations of various permanent networks (e.g., RIGN/RDN network and TopNETlive Italy network), have been processed to obtain the velocity field of this area. Using the inversion techniques proposed by Shen et al. (1996) and Grid Strain 2D software, a strain field was also obtained for the investigated area.

The CA represents an arc-shaped terrain originated from the deformation of the European paleo-margin in the context of Africa–Eurasia plate relative convergence (Malinverno and Ryan, 1986; Jolivet and Faccenna, 2000; Rosenbaum and Lister, 2004). It is part of a structural domain characterized by a large accretionary wedge located in the Ionian Sea, the latter formed as a response to the Miocene to Quaternary subduction process of the Ionian crust (Dewey et al., 1989; Faccenna et al., 2001 Figures 1B, C). Within this geodynamic framework, CA is made up of several tectono-stratigraphic units belonging to the Kabilo-Calabride Chain, constituted by Hercynian metamorphic and magmatic rocks, Mesozoic carbonate platform, and middle–late Miocene terrigenous sequences (Ogniben, 1969; Amodio-Morelli et al., 1976). Since the Late Miocene, CA has been dominated by extensional and transcurrent tectonics, superimposed on the previous collisional context, that have given rise to the shaping of Plio-Quaternary sedimentary basins (Ghisetti and Vezzani, 1982; Tortorici et al., 1995; Monaco et al., 1996; Monaco and Tortorici, 2000; Jacques et al., 2001; Tansi et al., 2007; Scudero et al., 2017; Pirrotta et al., 2021).

Several interpretations have been proposed to explain the regional extension affecting the CA: 1) isostatic response to the progressive Ionian Slab subduction processes (Westaway, 1993; Monaco et al., 1996; Wortel and Spakman, 2000); 2) ascent of asthenosferic flow due to slab roll-back processes (Faccenna et al., 2011); 3) counterclockwise rotation of the Adria Plate (D’Agostino and Selvaggi, 2004; D’Agostino et al., 2008); and 4) pull of the Hellenic slab dragging eastward the Ionian–Calabrian domain (Goes et al., 2004). Regardless of the cause, normal faulting has given rise to cumulative escarpments that characterize the present morphology of the major reliefs of the region (i.e., Catena Costiera and Serre and Aspromonte mountains; Tortorici et al., 1995; Monaco and Tortorici, 2000). Although different source models have been proposed over time (see DISS and references therein), normal faults are considered responsible for the most destructive historical earthquakes in CA (M. 6.5–7.4, Boschi, 1997; Jacques et al., 2001; Neri et al., 2020; Figure 1A), such as the 1905 and 1908 events and the seismic sequence of 1783. Recently, the inversion of structural data collected from extensional faults and joints (see Caputo and Caputo, 1989), combined with seismological data, allowed to define the stress field acting at various depths in the Calabrian Arc (e.g., 0–15 km, 15–35 km; De Guidi et al., 2013; Scudero et al., 2017). The latter testify a predominant WNW-ESE-oriented regional extension, in agreement with previous structural studies (e.g., Tortorici et al., 1995; Jacques et al., 2001) and consistent with GPS measurements (Hollenstein et al., 2003; Caporali et al., 2003; Palano et al., 2012; Serpelloni et al., 2005, 2010; D’Agostino and Selvaggi, 2004).

The 1783 seismic sequence was characterized by the following five mainshocks (Baratta, 1901; Jacques et al., 2001): the sequence started with a catastrophic event on 5 February that destroyed towns and villages located at the western foot of the northern Aspromonte (i.e., S. Cristina D’Aspromonte and S. Giorgio Morgeto; Figure 1A). On 6 February, a second large shock struck mostly the coastal area to the southwest, bringing the cumulative damage to this area to a level almost comparable to that in the Aspromonte piedmont, causing a rockslide along the cliff to the west of Scilla which in turn triggered a tsunami wave inside the Straits of Messina. On 7 February, a third shock ruined the small towns at the western foot of the Serre Mountains, approximately 40 km northeast of the epicentral area of the catastrophic event of 5 February (Figure 1A). On 1 March, a fourth shock struck a smaller area at the western foot of the Serre Mountains, approximately 20 km north of the epicentre of the 7 February event (Figure 1A). Finally, on 28 March, the fifth and last important shock of the sequence caused destruction farther toward the northeast, closer to the Ionian Sea. According to de Dolomieu (1784), after the 5 February mainshock, a 20-km-long scarplet that developed at the western foot of the Aspromonte Mountain between the village of S. Cristina and S. Giorgio Morgeto (red dashed line in Figure 1A) was observed. Jacques et al. (2001) interpreted this as the coseismic reactivation of the west dipping CF during the 5 February mainshock. Concerning the seismic source models of 5 February, 7 February, and 1 March, some authors proposed the CF and SRF west-dipping normal faults as responsible sources (Tapponnier et al., 1987; Ghisetti, 1992; Tortorici et al., 1995; Monaco and Tortorici, 2000; Jacques et al., 2001) (Figure 1A, DSS; Debated Seismic Source). This hypothesis was supported by field geologic and morphotectonic data and observation by the distribution of the high-damage sites as exploited by the mesoseismal areas of the 1783 sequence (Baratta, 1901) and revisited macroseismic data (Andrenacci et al., 2023), the observed coseismic fracturing (de Dolomieu, 1784), and also by paleoseismological studies (Galli and Bosi, 2002; Galli and Peronace, 2015). Conversely, other authors inferred ruptures along the Gioia Tauro and the Mesima east-dipping blind low angle (∼30°) normal faults for the same events (Valensise and D’Aaddezio, 1994; Cucci et al., 1996; Peruzza et al., 1996; Loreto et al., 2013) (see also Figure 1A, ISS: Individual Seismic Source). This hypothesis was supported for the first time by Ricchetti and Ricchetti, (1991) based on stratigraphic studies and subsequently by Valensise and D’Aaddezio (1994), Pizzino et al. (2004), and Tiberti et al. (2017) based on their model with the poor evidence of recent activity along the W-dipping faults, studies regarding the sediment attitude, and also by geochemical anomalies recognized along the Gioia Tauro faults. More recently, other studies have been proposed by Cucci (2022), who, accordingly with Jacques et al. (2001), reviewed his own previous theory (Cucci et al., 1996), asserting that the distribution of the observed coseismic hydrological anomalies after the 1783 seismic sequence does not support the role of E-dipping faults as accountable sources for these events.

The CF and SRF border the western sectors of the Aspromonte and Serre mountains, respectively. These NNE-SSW-oriented, west-dipping, and 40-km-long faults (Figures 1A, C) develop cumulative scarps of up to 450 m in height (Jacques et al., 2001). Morphological and geological analyses highlighted how the faults juxtapose the Plio-Pleistocene sediments of the Gioia Tauro and Mesima basins with the crystalline rocks outcropping in the Aspromonte and Serre mountains, respectively (Monaco and Tortorici, 2000; Jacques et al., 2001). The fault activity caused the tilting of Plio-Pleistocene sediments on the hanging wall and determined the development of various sets of triangular and trapezoidal facets on the footwall (Tortorici et al., 1995; Monaco and Tortorici, 2000). The drainage network is perpendicular to the fault segments and interacted with the tectonic activity resulting in pronounced v-shaped valleys on the footwalls. Differently, the hanging-wall sectors are characterized by flat morphology and large alluvial fans.

In order to investigate the surface expression and the kinematic of the faults, a geological–structural survey was carried out. Unfortunately, due to the poorly conservative rock types cropping out in the investigated area, the kinematic indicators are difficult to observe (particularly in the SRF area). Structural data were collected in a few selected outcrops along the CF near Oppido Marmertina Village. A first outcrop (Figure 2, St. 1 in Figure 1) highlights a clear fault gauge zone along a tectonic contact that juxtaposes two facies of the metamorphic basement. Structural analysis indicates a NE-SW trending, NW-dipping fault planes. A second outcrop is located southwest of Oppido Mamertina, along the Spilinga River (Figure 3, St. 2 in Figure 1), where the fault meets the Aspromonte gneisses in the footwall and the Middle–Upper Pleistocene deposits in the hanging-wall. Pleistocene deposits are faulted and tilted toward the west (layer attitude 310/35, Figures 3A, B). Minor faults are parallel to the master fault and show steep eastward dips (see also Jacques et al., 2001). The master fault (F1) is NE-SW-oriented and shows a set of slickenlines consistent with normal-oblique sinistral kinematics (Figure 3C). Figure 4 (see also St. 3 in Figure 1) shows structural data (see also Tortorici et al., 1995) collected along a cataclastic zone near the CF plane. The master fault exhibits a mean N 30° E direction with the west-dipping plane, with sub-vertical slickensides indicating prevalent normal kinematics. Associated minor faults were also recognized and are characterized by east-dipping planes.

To investigate the recent kinematics of the faults affecting the studied area, we analysed the available focal solutions from Scarfì et al. (2021) and the focal mechanism from ISC (http://www.isc.ac.uk/iscbulletin/search/fmechanisms/) and ISIDe (http://iside.rm.ingv.it/tdmt) databases (see Figure 1C). Within the first 15–20 km depth, focal mechanisms highlight prevalent normal kinematics, with major nodal planes oriented NNE-SSW, consistent with the geometries of the faults characterizing this domain. Conversely, the seismic layer between ∼20 and ∼30 km is characterized by prevalent strike slip mechanisms. In the deepest layer (30–40 km), focal solutions indicate reverse and strike slip faulting.

Earthquakes that have been instrumentally recorded since the early 1980s were used to infer the geometry of the studied faults at depth. We selected the seismic events occurred in the Calabrian Arc from the INGV databases (https://istituto.ingv.it/it/risorse-e-servizi/archivi-e-banche-dati.html) in order to enhance the picture derived from the seismic dataset. The initial parameters (i.e., arrival time and locations) were processed using tomoDDPS software (Zhang et al., 2009), which improves the accuracy of hypocentre location through a combination of absolute and differential arrival-time readings between couples of closely -spaced earthquakes. In addition, the code allows computing of the seismic ray-tracing in a 3D velocity model; here, we used the velocity model by Scarfì et al. (2018). The final locations (approximately 8,500 seismic events with magnitudes ranging from 1 to 5.7) resulted with an average uncertainty of 0.20 ± 0.15 km in both horizontal and vertical coordinates and an average root-mean-square travel-time residual of 0.02 s.

Kernel density maps of the final dataset allowed us to identify the most seismically active zones. The highest density of crustal seismicity (0–20 km) (Figure 5A) is recognizable along the main faults, particularly at the southern tip of the SRF and at the southern tip of CF. Sub-crustal seismicity (20–45 km) (Figure 5B) is concentrated in the northern and southern sectors. Note that the numerical values of density for the latter depth layer are lower than the other values (i.e., earthquakes/km² lower than 0.37). Figure 5C shows the earthquakes between 45 and 300 km. In this case, the earthquake distribution and the corresponding high density associated with them are concentrated near the Tyrrhenian sector, which could be considered related to the Ionian Slab.

To constrain the depth geometry of Serre and Cittanova faults, a set of 10-km spaced seismological sections (five sections for each fault, Figure 5D) with a buffer projection of 5 km was created (Figure 6); the geometry at depth of the considered faults (red lines) was traced following earthquake clustering starting from the intersection of the fault on surface (black crosses in Figure 6). In particular, in the 0–20 km depth range, the instrumental seismicity shows a high concentration of events near the area of Cittanova and Serre faults; clusters are visible, especially in S1 and S6 corresponding to the southern tip areas of the CF and SRF, respectively. Despite the other sections that do not exhibit clear clusters useful to infer the attitude of the studied faults at depth, we traced their geometry considering the same trend of S1 and S6 and the earthquakes with the highest magnitude for each sections. Subsequently, using a trial-and-error approach, we geometrically tested the modelled planes with the known empirical scaling in order to find a reliable solution for fault planes capable of generating events with a magnitude of 7 (as reported by historical catalogues for this area; CPTI15 catalogue; Rovida et al., 2022). It is worth noting that other clusters may be considered to trace a possible CF plane (see clusters to the left and right of the red line shown in S2, S4, S5 of Figure 6); however, for the aforementioned reasons, we decided to consider a trace that better fits with surface data (i.e., fault dip and position). In addition, here, clusters vertically arranged in sections S1, S4, and S6 (to the right of the red line, Figure 6) are interpreted as extensional cracks forming at footwall areas in response to the dilatational strain that generates at the base of a flexurally supported upper crust during normal faulting (see Ellis and King, 1991). Nevertheless, overpressured fluids causing distributed microseismicity (as proposed by Collettini et al., 2022 in the Central Italian areas) could not be ruled out. Frequently, without fluid overpressure, faulting can generate tension fractures within the first kilometre in depth (Griffith, 1920). The development of tension fractures at greater depths may be facilitated by loading dynamics associated with slips on seismogenic faults (see also Gudmundsson, 2000; Gudmundsson, 2011).

The same seismic clusters were observed along S1–S5 and were used to infer the presence of other faults (i.e., the SEF and SF, S. Eufemia and Scilla faults, respectively, Figure 6). The S5 section is further north with respect to the SEF trace at surface (see also Figure 5D); thus, we consider this fault uncertain. Another interesting cluster is located between the tips of SRF and CF; this is also visible in the map view (Figure 5D). The cluster is roughly NNW-SSE-oriented, and it is probably generated by slips along the Nicotera–Gioiosa Fault Zone (NGFZ) (see Tripodi et al., 2018; Tripodi et al., 2022). To better investigate the origin of these seismic concentrations, we built two additional NNE-SSW trending seismological sections (S1.1 and S2.1, Figure 5D) (see Supplementary Figure S1). These sections allowed us to interpret the high concentrations of events between the SRF and CF as the result of the NGFZ activity. The fault belt is interpreted as a transfer zone that accommodates the SE migrations of the Calabrian terranes (see also Tripodi et al., 2022).

We collected GNSS data from permanent stations of the RING network (http://ring.gm.ingv.it/) and TopNETlive Italy network in order to obtain the velocity field of the area across the previously described major faults (Figure 7). GNSS data were processed using GipsyX 1.5 (Bertiger et al., 2020), with precise ephemerides and clock correction provided by the Jet Propulsion Laboratory (https://sideshow.jpl.nasa.gov). In order to correctly process the data, we used the Earth orientation parameters from the International Earth Rotation Services (https://www.iers.org), the absolute IGS antenna phase center (GPS week 1958), and the Global Mapping Function (GMF) atmospheric zenith delay models (Boehm et al., 2006). Then, we processed the RINEX (Receiver INdependent EXchanges) files using the gd2e.gy module in GipsyX to estimate the position of the GNSS stations by using the Precise Point Positioning method. The obtained coordinates for each day of the year (DOY) were used to compute time series in the ITRF2014 reference frame (Altamimi et al., 2016). Subsequently, to show the current velocity field of the area, the velocities are referred to the eastern Sicily (Carnemolla, 2021, PhD thesis, for further information see also Pirrotta et al., 2021).

In order to investigate the effects of CF and SRF on surface, we calculated the strain field of the studied area from the GNSS velocity. This analysis is used by several authors to provide useful elements for studying tectonic phenomena (e.g., Caporali et al., 2003; Anzidei et al., 2005; Serpelloni et al., 2005). The strain field was obtained using Grid Strain (2D) software (Teza et al., 2008) based on the modified least square approach (Shen et al., 1996; Shen and Jackson, 2000). Strain eigenvalues and eigenvectors were computed within a series of experimental points (EPs) near a regular grid of node constructed in function of the median distance of available GNSS stations (Figure 8). We applied a reiterative approach in order to estimate the optimal scale factor for the distance weighting and to exclude the points too far from grid nodes. The reliability of the data inversion was estimated based on a significance geometric evaluation criterion in order to verify whether the grid nodes are well distributed respect to our experimental data (see also Teza et al., 2008).

Figure 7 shows the horizontal velocities of GNSS stations resulting from our processing and error ellipses with a 95% confidence level. The velocity vectors exhibit a general south–east trend of motion with an average velocity of 2–3 mm/yr. The major change in the velocity is recognized near the Messina Strait; in this area, the velocity vectors show an increase in the module from west to east. Other minor changes in the velocity are visible across the major faults; here, a differential change in velocity is recognized near the CF and SRF traces. The WNW-ESE-oriented extensional strain is predominant, ranging from 1.2e⁻⁸ to 3.29e⁻⁸ strain (Figure 8). In the southwestern sector, near the Messina Strait, the extension rotates clockwise assuming a NW-SE direction. On the other hand, compressional strain is less evident in this domain, especially near the area surrounding major faults. The grey-shaded map shows the eigenvalue distribution that reaches the highest strain near the Messina Strait and near the footwall of CF. Maximum horizontal strain axes show an extension perpendicular to the directions of major faults. Thus, the obtained extensional dynamic is consistent with the observed tectonic framework characterizing the southern-central Calabria (Westaway, 1993; Tortorici et al., 1995; Monaco and Tortorici, 2000; Jacques et al., 2001; Ferranti et al., 2008; Palano et al., 2012).

Adopting a multidisciplinary approach based on morpho-structural field data integrated with seismological and geodetic data, we modelled the 3D geometry of the Cittanova and Serre faults in southern Calabria, performed the FRM simulation for these faults, and evaluated their propensity to slip in the MOVE environment. On the basis of our results (maximum expected magnitude, vertical slip evaluation, and slip tendency) and according to the historical report (de Dolomieu (1784)) and previous works (Jacques et al., 2001; Andrenacci et al., 2023), we found that CF and SRF (called “debated seismogenic sources” in DISS; DISS Working Group, 2021, https://diss.ingv.it/) are the best candidate sources of the 5 February, 7 February, and 1 March 1783 events. On the contrary, the Gioia Tauro and Mesima faults do not exhibit surface evidence (e.g., coseismic scarplet, cumulative morphological scarps, kinematic indicators, and tectonic control on the drainage network) that would testify a recent activity associated with them.

The presence of numerous earthquakes and geodetic differential velocity across the faults testify that the southern Calabria is an active tectonic domain. The kernel density map for crustal earthquakes (0–20 km) (Figure 5A) shows a high value of density near the southern tip of CF and SRF, highlighting an increasing stress concentration in these sectors. According to Gudmundsson et al. (2010), this feature can be related to stress concentration that can occur around the fault zone due to a sudden change in geometry or geomechanical parameters, such as a decrease in Young’s modulus within the fault zone compared with the outside zone. The earthquake concentration near these structures may also be influenced by the presence of other structures, such as S. Eufemia, Armo, Reggio Calabria, and Scilla faults in the case of the CF source, and the presence of Nicotera and Coccorino faults in the case of the SRF source (Figure 1). Some of these faults are also considered capable of generating high-magnitude earthquakes (see also Jacques et al., 2001). Earthquakes show off-fault clusters vertically arranged in the same sections (Figure 6), suggesting the occurrence of extension or well-defined hydraulic fractures, which could be related to dynamic loadings, such as those associated with slip on seismogenic faults. In our opinion, considering the dip-attitude of faults at surface (∼60°) and the known strikes (∼N35°–40°E), clusters are congruent with a west-dipping geometry.

The strain field achieved from the geodetic velocity (Figure 8), consistent with the extensional stress state obtained by structural data inversion and with the stress state previously defined by other authors (e.g., Tortorici et al., 1995; De Guidi et al., 2013; Scudero et al., 2020), testifies a prevalent long-term extensional deformation oriented toward WNW-ESE. The latter seems to be orthogonally oriented with respect to CF and SRF attitude, which probably accommodate this process. Our results also agree with previous geodetic analysis carried out in these areas (Caporali et al., 2003; D’Agostino and Selvaggi, 2004; Serpelloni et al., 2010; Palano et al., 2012; Pirrotta et al., 2021). Although the strain values are very low (max strain equal to 3.29 × 10⁻⁰⁸), the extensional strain field is clear. Low strain affecting the southern-central Calabrian is also documented by several studies (e.g., Mattia et al., 2009 ; D’Agostino et al., 2011; Palano et al., 2012 and reference therein); nevertheless, the process that determine this dynamic is still debated. In our opinion, the extensional regime is strictly related to the complex deep crustal dynamics and consequent mantle upwelling affecting the Calabrian Arc (see also Gvirtzman and Nur, 1999; Faccenna et al., 2001; Goes et al., 2004). In this scenario, the low strain in this domain vs. the high seismic moment release of the area remains an ongoing open challenge (see also Carafa et al., 2018).

Starting from the geometric parameter of the modelled planes, we investigated their reliability in relation to historical earthquakes (i.e., the 1783 events) using empirical scaling relationships (Wells and Coppersmith, 1994; Leonard, 2010; Table 1). Considering the relationship between the surface rupture length (SRL) vs. magnitude and the rupture area (RA) vs. magnitude, we found that both the CF and SRF planes are capable of generating an earthquake with M ∼7. Although the maximum length of the surface rupture observed by de Dolomieu (1784) is equal to 20 km (Figure 1C), using the empirical relationships (SRL vs. M), we calculated that the CF and SRF were activated at depth along their entire length (∼40 km) by an earthquake with magnitude ∼7. Furthermore, using the FRM module, we obtained the vertical displacement field for the modelled planes (Figure 10A). This seems to be consistent with the average slip of 3 m observed by de Dolomieu (1784) and with the cumulative escarpment associated with these faults. Note that the synthetic displacement for the FRM tests (see Section 7) was achieved from a kinematic constrain of the slip from direct observations (from de Dolomieu, 1784) and from the empirical scale by Manighetti et al. (2001) (see also Table 1).

To better understand the variability of the synthetic vertical displacement dz associated with the activity of the faults, we made a set of displacement profiles across the fault traces (Figures 10C, D). In our opinion, for the first test, the resulting synthetic dz distribution across the fault distance can depict the possible displacement field associated with an event similar to the main 5 February 1783 shock, also considering that the abrupt change in dz near the CF escarpment agrees with the average 3 m of slip on the fault plane recorded for this event. Assuming that the high-damage area corresponds to the area of maximum dz, we find a good match between the mesoseismal areas and the dz field distribution (Figure 10A). For the second test, the cumulative dz profile across the SRF and CF fits with the morphological escarpment across these faults, although this last simulation is valid only for the area enclosing the fault escarpment (i.e., 3–4 km across the fault traces, red lines in sections, Figure 10D). Moving away from the fault trace, dz can probably be overestimated or underestimated considering the presence of the other faults, such as S. Eufemia and Scilla, as well as the combination of erosion and isostatic rebound processes.

Finally, we used the slip tendency analysis (Morris et al., 1996) to evaluate a possible reactivation of the fault planes using the geometric parameters found in our 3D model and applied a remote stress obtained from structural and seismological data (see Sections 3–5). A critical review on the results of this analysis must take into account the deviation from the average geomechanical parameter values characterizing the host rocks and the fault planes. In our opinion, a density of 2,600 kg/m³, an angle of internal friction of 30°, and a Young’s module of 75 GPa represent the ideal condition for this evaluation. Furthermore, considering that these values can change with depth, the evaluation of Ts depends on change in the depth at which the resolved shear stress was computed.

A multidisciplinary study was developed to model the 3D geometry and kinematics of the Cittanova and Serre faults. Despite the scarcity of kinematic indicators in field, morpho-structural features (e.g., triangular and trapezoidal facets, drainage network, fracturing pattern, and kinematic indicators) confirm that CF and SRF are characterized by northwest-dipping planes with prevalent normal motion.

The instrumental seismicity, merged with structural field investigation and literature data, provided useful constraints to infer the geometry on surface and at depth of these faults. Earthquake concentration near the fault zone suggests that these structures are active and probably accommodate the extensional strain pattern recognized from the geodetic data. Strain pattern and geodetic velocity are also consistent with previous studies (see Palano et al., 2012; Pirrotta et al., 2021) and testify that the south-central Calabria is affected by a predominant extensional deformation.

Combining all previous data, we built, for the first time, a 3D model of the Cittanova and Serre fault planes. The geometric parameters obtained through the model of these faults are compatible with the empirical relationships (magnitude vs. rupture area and magnitude vs. fault length; Wells and Coppersmith, 1994; Leonard, 2010). Accordingly, we achieved the expected magnitude for CF and SRF and confirmed that these faults can be the probable sources of the mainshocks of the 1783 seismic sequence (5 February, 7 February, and 1 March, M = 6.9–7; Jacques et al., 2001). From the slip tendency analysis, we achieved the propensity to slip for both CF and SRF. According to this, we found that Ts values range from 0.35 to 0.85 and 0.45 to 0.80 for SRF and CF, respectively. This could suggest that both planes are astride under the stable and unstable mechanical conditions in the given regional stress state.

Considering the lack of seismic profiles onshore, the proposed model can represent an important starting point for seismotectonic modelling. The response of the modelled fault planes indicates that the simulated coseismic and cumulative vertical displacement fields agree with the historical observation of the slip along the fault planes and with actual height of the morphological scarps, respectively. Our work may lay a foundation for future data inversion to better constrain the seismogenic sources, for example, using the Okada (1992) analytical method and adopting interseismic kinematic models such as block modelling (McCaffrey, 2002; Meade and Hager, 2005) or dislocation modelling (Savage and Burford, 1973) approaches.