Active Tectonics in Tuscany (Central Italy): Ten Years of Seismicity (2009-2019)

Strong earthquakes (moment magnitude M W ≥ 5.5) are uncommon in Tuscany and surroundings (central Italy). The last strong seismic event occurred a century ago (September 7, 1920 Garfagnana, M W = 6.53). The paucity of seismic instrumental recordings hinders the identification of the tectonic regime active in Tuscany. On the other hand, the geological and geomorphological pieces of evidence collected so far, concerning potential active and capable faults, are scarce, fragmentary and ambiguous. In this work I shed light on the active deformation of Tuscany by using two independent approaches: earthquake source mechanisms and GNSS (GPS) geodetic measurements. I have considered 41 small seismic events (M W ≤ 5.1) that occurred in the study area during the last decade. The related source mechanisms (retrieved by the Time Domain Moment Tensor method) define a relatively clear picture of the active deformation: extension along the northern Apennine watershed and strike-slip regime within inner Tuscany, up to the Tyrrhenian coast. This pattern broadly agrees with the horizontal strain field reconstructed by the geodetic velocity field. The latter has been constrained by a network of 840 GPS stations located in Italy and neighboring countries, operating in the last 20 years.


Structural Setting and Neogenic Evolution
In the study area, crustal thickness increases from about 20 km below Elba Island to 22 -26 km in central Tuscany. The Moho depth reaches 50 km beneath the northern Apennine arc, then it decreases to 26 -30 km near the Adriatic coast [16].
Geophysical survey allowed us to shed light on the architecture of this part of the Apennines [17] [18] [19] [20] [21]. It is a typical thrust-and-fold belt, formed by the progressive stacking of crustal slices. Regional, SW dipping overthrusts make the contact between the slices. They cut the sedimentary cover, dissect the metamorphic basement and perhaps prosecute in the underlying upper mantle as shear zones [20] [21] [22].
The study area belongs to the western edge of the Adriatic plate, whose deformation began in the lower Miocene with the development of eastward migrating thrust belt-foredeep systems [23] [24] After the opening of the northern Tyrrhenian basin in the late Miocene [25], the compressional regime involved an increasingly large area of the Adriatic foreland. Moreover, the related foredeeps were incorporated in the eastward migrating Apennine belt.
Two conflicting interpretations of the study area have emerged over time. The first interpretation provides that since middle Miocene, inner Tuscany was affected by an extensional stress regime. It generated a horst-and-graben (often semigraben) system, corresponding to the ridges and basins indicated in Figure  1 [8] [26]- [38]. Until the lower Pliocene, extension would have occurred by low-angle normal faults (detachments). The subsequent deformation was instead produced by high-angle normal faults, which bound the present depressions [8]. Figure 2 shows the main normal faults reported by the above works.
The alternative interpretation implies that compressional deformation dominated the study area until perhaps the middle Pleistocene [9] [40]- [52]. Indeed, the above Authors recognize the effect of normal faulting from the late Quaternary only. In this context, thrusts and reverse faults would have controlled the development of antiformal folds, corresponding to the ridges shown in Figure 1. The depressions developed as synform or thrust-top basins. Figure 3 shows the main reverse faults, thrusts and folds axes reported in the above works.
The aforementioned debate may seem irrelevant to the analysis of seismogenic processes, as it concerns past (pre-Holocene) tectonic phases. However, the gap between the two interpretations may cause serious ambiguity in mapping of tectonic features. For example, the boundary between the Mt. Cetona ridge and Chiana basin is reported as an east-dipping normal fault (VC in Figure 2 [53]), or as a west-dipping thrust (CC in Figure 3 [9]).
In both the above interpretations, most tectonic features are NW-SE directed ( Figure 2 and Figure 3). However, some NE-SW features have been suggested also. Such lineaments (e.g., Al and GP in Figure 2) would explain the interruption in the continuity of the sedimentary basins. However, these transverse fault systems are often inferred by stratigraphic and structural considerations, rather than by direct field mapping [13]. More evidence on these features is reported along the northern Apennine arc [54].
An example of the Tuscan transverse structures is the dextral strike-slip fault system located just north of Valtiberina basin (PSS in Figure 2), which however is much older than the basin itself [55].
Recently, some sinistral, transverse fault systems have been described in central Tuscany (Am, Ar, Ca, Rp, St and Vi in Figure 2). These fractures are very recent and probably active, because they are associated with growing travertine deposits [53] [56]- [61] suggest that the western sector of the Larderello geothermal field (LGF in Figure 1) is a pull-apart basin bounded by NE-SW sinistral strike-slip faults. Similar transverse faults have been described in the Pleistocene volcanoes of northern Latium (Ci, Lt and Tf in Figure 2 [62]). [63] recognize the Mt. Amiata strike-slip fault (Am in Figure 2), but they attribute it a dextral kinematics.
Although the faults and folds shown in Figure 2 and Figure 3 are often prominent morphotectonic features, they do not necessarily represent active structure or seismogenic sources. Indeed, the identification of active and capable faults [64] faces several difficulties. As regards Tuscany, it should be considered in its initial phase [63].

Geomorphological Constraints: Local and Regional Uplift
Two important uplift episodes affected the study area since the Upper Miocene. In some cases, surface uplift can be a consequence of the emplacement in the crust of plutons and magma chambers [73]. This mechanism has been adopted to explain the current elevation of the Pliocene marine deposits of southern Mts. and the adjacent Siena and Valdarno superiore basins ( Figure 1 and Figure   2) have undergone a substantial Quaternary uplift [76]. Likewise, in north-western Tuscany the Apuan Alps acquired their rugged morphology in the Quaternary only [77]. In addition, conspicuous uplift and exhumation rates are estimated both north-west and south-east of the Mugello basin, as well as in the Romagna Apennines [31] [70] [78].
Any interpretation of the active tectonics of the study area must take in due account the evidence of recent uplift, at both the local and regional level. [87]. This may be due to the fact that the most recent strong earthquakes have hit the outer, intramontane basins ( Figure 1).

Active Tectonics and Seismicity
The active structures identified (more often suggested) by the above studies usually are normal faults ( Figure 2). Active reverse faults and folds are instead reported in the Tuscan-Emilia-Romagna Apennines (e.g. CdP in Figure 3).
So far, none of the active faults proposed in the literature has revealed its seismogenic nature, because the last strong earthquake occurred in the study area dates back to a century ago (1920 Garfagnana, Figure 1). Plain earthquakes did not produce surface faulting [105].
Some complications can also occur in the presence of marked fault scarps with measurable kinematic indicators: 1) An ancient tectonic feature can be mistakenly recognized as an active fault.
For example, [106] refute the existence of a normal fault in the Val d'Agri (southern Apennines), which [107] consider as the source of the 1857 large earthquake (M W = 7.12 [6]). In this case, the detected triangular facets are not the traces of a normal active fault. Instead, they probably are flatirons, i.e. eroded remains of a limb of an ancient fold [106]. In fact, triangular facets may appear on any tectonic feature (normal fault, reverse fault and fold) and by selective erosion also [104] [108].
2) The morphological prominence of the fault scarp is not always a reliable indicator of its recent/present activity. A well-known example is the Gubbio fault (Gu in Figure 2), which bounds to the east the Gubbio intramontane basin. [104] also advises caution in assessing the recent activity of fault escarpments by morphotectonic analysis only. For example, [63] suggest the present activity of the Mt. Cetona normal fault (Ce in Figure 2), which bounds to the east the Radicofani basin. They also estimate a recent slip rate of about 0.6 mm/year. On the contrary, [53] state that this Pliocene fault is now inactive, since it is dissected by a Quaternary strike-slip fault system (St in Figure 2).
3) The so-called selective erosion, due to the different weathering of outcropping formations, can generate forms simulating normal active faults [108]. For example, the alleged active Cerbaie fault (Cb in Figure 2), located in the Valdarno inferiore basin, could actually be an erosional landscape [108].

4) A large landslide can generate a shallow deformation pattern simulating
normal, strike-slip and reverse faults in the detachment area, sides and toe of the landslide respectively [112]. [113] discuss several features misinterpreted as active normal faults in the central Apennines.   Figure A1 of the Appendix. The parameters of the nodal planes are shown in Table A1 of the Appendix. The dashed green lines delimit the three sectors discussed in the text: I) north-western Tuscany II) Tuscan-Emilian-Romagna Apennines III) inner Tuscany. The green shape identifies the territory of Regione Toscana.
In brief, the search for seismogenic structures by morphotectonic analysis should be supported by independent methodologies. In the next section I consider the use of the source mechanisms of small earthquakes that have recently occurred in the study area.

State of the Art
Global, national and local seismometric networks make possible recording many seismic events, most of which can be defined small earthquakes (M W ≤ 5). The analysis of the micro-seismicity provides basic information in order to define the seismotectonic pattern of an actively deforming zone [114]. This kind of investigation is particularly important where tectonic deformation is slow and strong earthquakes are relatively rare. Furthermore, using micro-seismicity has the obvious advantage of dealing with actual (albeit small) seismic slip events. The earthquake source mechanism (ESM), representing the geometry of the seismic fracture, is fully described by the seismic moment tensor [115]. As earthquakes are basically slip events on shear fracture, the seismic moment tensor is dominated by the double couple solution, represented by a couple of nodal planes and the related tensor eigenvectors (B, P, and T axes).
A set of ESMs allows us constraining the stress and strain tensors in the crustal volume affected by seismic activity [116] [117] [118]. Small earthquakes contribute little to the overall deformation of a given crustal volume [119]. However, the seismic strain style (i.e. the direction of the principal strain axes and the ratio between their amplitudes) is similar for large earthquakes and for small shocks. Therefore, micro-seismicity can be used to constrain the seismotectonic pattern [120]. Furthermore, small earthquakes re-distribute elastic stress on the involved faults, thus influencing the development of future seismicity [121].
Some previous works deal with the micro-seismicity of the study area, although the analysis is often limited to specific sectors. For example, a local seismometric network has provided several ESMs for the Larderello geothermal field (LGF in Figure 1 [122] [123]). In particular, [124]  [128] consider the micro-seismicity (2.5 < M < 4.8) occurred in Italy in the period 1988-1995. The EMSs were obtained by the first-arrival polarity method.
Then, [129] procedure was adopted to constrain the stress regime in various sectors of the Italian region. The two sectors named "Northern Tuscany" and "Pery-Tyrrhenian" include north-western Tuscany and the Tyrrhenian coast of Tuscany and Latium respectively. In both sectors, the estimated stress regime is extensional, with a ENE-WSW directed, sub-horizontal minimum compression axis (σ 3 ). In the stress regime map elaborated by [130] the study area is characterized by pure extension, with NE-SW σ 3 .
However, the determination of the stress regime by ESM inversion is based on some assumptions that are not easily verifiable [117]. In fact, the inversion procedures require homogeneity of rock mechanical properties and uniformity of the stress state in the crustal volume considered. Furthermore, the source mechanisms should be independent each other, which is in contrast with the fault interaction expected by the rock elastic behavior [121].
[131] considered the seismicity of the central-northern Apennines recorded in ESMs are not sufficient to clearly constrain a seismotectonic pattern.
So far, seismological evidence has not provided a widely shared interpretation of the active tectonics of Tuscany. In the next section I will discuss a new set of ESMs that could allow us to approach that objective.

New Information from Micro-Seismicity
Since 2005, the Italian Istituto Nazionale di Geofisica e Vulcanologia (INGV) has been processing seismic waveforms by the Time Domain Moment Tensor (TDMT) procedure. This concerns earthquakes having magnitude M W ≥ 3.5, whose epicenter in located in Italy and neighboring countries [133].
The TDMT procedure is based on the algorithm developed at the University of Berkeley, California (http://seismo.berkeley.edu/~dreger/mtindex.html), which inverts the three-component waveforms recorded by broadband seismic stations.
The synthetic seismograms used for seismic moment tensor inversion are calculated by means of a laterally homogeneous, 1D velocity model calibrated on a regional basis. The focal mechanism representation is obtained with the procedure developed by [134].
Using the ESMs derived by TDMT analysis has two important advantages: 1) they have been obtained from recent seismic recordings, which benefited from the modernization of the Italian seismometric network 2) the seismic moment tensor inversion methodology is common for all data considered. This eliminates the problems related to the use of ESMs having heterogeneous quality, as they were retrieved by different procedures [49] [131] [135].
Available ESMs can be downloaded from the dedicated website  Moreover, I have assigned to each datum an order number ranging from 1 to 41.
Epicenter coordinates and basic nodal plane angles are reported in Table A1 of the Appendix. For each ESM, the stereographic projection of the focal sphere, with the poles of T and P axes, is shown in Figure A1 of the Appendix. The strain style indicated by the ESMs suggests dividing data into 3 subsets, corresponding to 3 different sectors of the study area (I, II and III in Figure 4). Data numbering follows this grouping. Sector I (16 ESMs from n.1 to n.16) includes north-western Tuscany, with the Apuan Alps and the Lunigiana and Garfagnana basins (Figure 1 and Figure 2 (Figure 1). The mechanism is clearly extensional, with roughly N-S T-axis. Data from n.6 to no.15 are the aftershocks of n.4. They show very similar ESMs, analogous to those of the main shock. It is worth noting note that the above ESMs are compatible with the location, direction and kinematics of the Northern Apuan Alps fault described by [86] (NA in Figure 2). Datum n.16 is located near the Tyrrhenian coast in the Viareggio basin ( Figure 1 and Figure 2). The mechanism is still extensional, but the T axis is directed ENE-WSW.
In summary, in sector I shallow events (H < 10 km) have an extensional mechanism, while deeper shocks show strike-slip seismic strain.

Sector II
In this sector, the epicentres are located mostly along the northern Apennine arc, near the watershed between the Tyrrhenian and Adriatic Sea. Data n. 17  The mechanism of the above events is extensional, with T-axis trending from NNE-SSW to NE-SW. [127] presented similar results by considering the small earthquakes that occurred in the same area from 1997 to 2005. These ESMs contrast both the active reverse faulting suggested by [70] and the trans-pressional strain regime proposed by [49].
Datum n.30 (17/10/2014 Chiusi della Verna) falls in the Alpe di Catenaia ridge, which separates the Casentino and Valtiberina basins ( Figure 1). The strain regime is roughly strike-slip, with a NW-SE P-axis and NE-SW T-axis.
In summary, the shallow seismicity of sector II highlights an extensional style of seismic strain, with a sub-horizontal, roughly NE-SW T-axis. As noted for Sector I, shocks deeper than 10 km have a strike-slip or compressive ESM.  Figure 1). The mechanism is compressional, with a NW-SE P axis and sub-vertical T-axis. Active compression in that area disagrees with the alleged activity of the Mt. Cetona normal fault (Ce in Figure 2), suggested by [63]. Earthquake n.35 (31/05/02016 Acquapendente) is located near the Latium-Tuscany border. The mechanism is trans-tensional, with NE-SW T-axis.

Sector III
The epicenter of the shock n.36 (18/03/2013 Buonconvento) is located in the Siena basin (Figure 1). The trans-tensional mechanism shows a NE-SW T-axis. Event n.37 (04/03/2016 Colle di Val d'Elsa) is located in the Elsa basin ( Figure  1) and has a trans-pressional mechanism, as indicated by the nodal planes parameters reported in Table A1  compatible with the stress state described for that zone by [124]. Events n.40 and 41 (09/08/2014 Castelfiorentino and 25/10/2016 Certaldo) took place in the Elsa basin ( Figure 1). The mechanism is strike-slip with NW-SE P-axis.
ESMs of sector III show a remarkable similarity as regards the orientation of the nodal planes and the P and T axes (Figure 4 and Table A1). The strain style is mostly strike-slip, while no solution indicates the pure extensional regime observed in sector II.
These information cast doubt on the activity of the normal faults of southern Tuscany, suggested by [63]. Moreover, the purely extensional stress regime proposed for Tuscany [128] [130] disagrees with the above ESMs. Instead, such data suggest that the seismic deformation of Sector III is caused by shallow strike-slip faults, such as those recognized by geological survey (Am, Ar, Ca, Rp, St and Vi in Figure 2 [53] [58]). However, the epicenter of the earthquakes discussed above is quite distant from the aforementioned faults. In addition, recent seismicity affected zones as the Colline Metallifere and the Elsa and Pesa basins, where no strike-slip faults have been described so far. Thus, the link between actual seismogenic faulting and mapped active/capable faults remains weak for inner Tuscany.
In brief, the source mechanisms of recent earthquakes clearly indicate that along the northern Apennine arc (Sectors I and II in Figure 4) the strain style is extensional in the shallow crust (H < 10 km) and strike-slip or compressional at greater depth. This feature supports the conclusions drawn by [131]. Within inner Tuscany, the strain style is mainly strike-slip, with a roughly NW-SE P-axis and NE-SW T-axis.
In the next paragraph I will discuss the implications of the satellite geodesy measurements, with particular attention to the horizontal strain field.

Space Geodesy: Further Constraints on Active Deformation
Satellite geodesy networks such as the Global Navigation Satellite System (GNSS) or Global Positioning System (GPS) allow researchers to investigate plate kinematics at global, regional and local scale [137]. The geodetic velocity field has long been used to constrain the seismotectonics of the Italian region [

Analysis of the Geodetic Measurements
In this paragraph I describe the treatment of the geodetic measurements carried out in many GPS sites located in Italy and neighboring countries. Data come from 840 GPS permanent stations, whose location is shown in Figure 5. Since the treatment of GPS time series is detailed by [140] [141] [144]. I report here only basic information.  The GPS permanent stations belong to various geodetic networks, managed by research institution, professional organizations and commercial companies. However, all stations use good quality, dual frequency instruments. In addition, the adopted support always ensures tight coupling between the antenna and ground or building walls. [147] shown that the different ways of installing and managing GPS stations do not entail neither systematic noise increase in time series, nor important signal distortions. For this reason, the observations acquired by professional and commercial stations can be used for scientific studies also. Data were analyzed with the GAMIT/GLOBK software [148] [149]. The International GNSS Service (IGS) provide information about the orbit of GPS satellites, data on Earth rotation and revolution, and parameters for modeling the electronic center of antennas. Due to the large number of stations, I adopted the distributed processing proposed by [150]. The initial part of the analysis is performed by dividing the initial network (all stations), into several subnets (clusters) containing a smaller number of sites. Therefore, the network was divided into 45 clusters. I inserted in each cluster the observations of at least 6 stations belonging to the European EUREF network, which guarantees measurement quality. Then, the daily solutions of all the subnets were combined by using the position of the common stations inserted in the different clusters as a constraint. This was done on the daily solutions of the subnets, thus obtaining a general daily solution.
As the different daily solutions may not be expressed in the same system, each solution must be referred to a single reference system. For this reason, the observations of 16 stations located around the Italian region were included in the subsequent analysis. Indeed, IGS provides position and velocity of these stations, referred to the International Terrestrial Reference Frame (ITRF). Consequently, the daily solutions of the geodetic network have been included in the updated ITRF2014 version, by using the position and velocity data provided by [151].
In order to insert each station in the ITRF2014, a 7-parameter Helmert transformation (3 translations, 3 rotations and a scale factor) was carried out, by means a weighted least squares approach which compares the calculated station coordinates with those provided by IGS. Finally, I obtained the time series of the daily position of each station: north component (measured from Equator), east component (measured from the Greenwich meridian) and elevation component (height on the WGS84 reference ellipsoid).
Then, the daily time series of the above components have been elaborated by [140] [141] procedure, in order to determine the velocity of the vertices of the geodetic network. The time series of each component was analyzed independently of the others, in order to eliminate any anomalous data (outliers) by using a specific software [152]. Once the anomalous data were eliminated, I obtained a first estimate of the station velocity, along with any discontinuities eventually caused by maintenance, or by nearby earthquakes. The above parameters, estimated by means of a weighted least square approach, were used to determine the theoretical time pattern of the three components of the position of each station. This theoretical trend is subtracted from the real trend, yielding a time series of residual. This residual series is then analyzed using a Lomb-Scargle approach [153] [154] in order to estimate the series frequency spectrum. Then, the spectrum is analyzed to determine the statistically significant periods of the 5 main seasonal signals. This search is carried out in the time interval ranging from 1 month to half of the entire observation period of the station, to avoid aliasing or indetermination problems.
The time series of the 3 components (without outliers) are analyzed again to determine a new estimate of the velocity, discontinuities, amplitude and phase of the 5 seasonal signals whose periods were estimated during the spectral analysis. Then a new theoretical daily trend is estimated and subtracted from the real one, to obtain a new time series of residuals. Finally, this last series is analyzed to determine the noise content. Using this information, it was possible to realistically estimate the errors associated with the station parameters.
In summary, the above procedure provided the velocity vector of each GPS station (described by the north, east and vertical components) and the related uncertainty.

Horizontal Velocity Field
The vector sum of the north and east velocity components gives the site absolute horizontal velocity. However, it is suitable for seismotectonic analysis subtracting the motion of the Eurasian plate from absolute velocity components. Therefore, I have used the plate rotation parameters provided by [151] to remove the Eurasia plate motion at each station. The residual velocity field (with respect to the Eurasian plate) is shown in Figure 5. This kinematic pattern updates previous results, based on a shorter observation time interval [140] [144] [155]. Indeed, this work use GPS measurements gathered in the last two decades (from 01/01/2001 to 30/04/2019).
The horizontal velocity field here presented confirms previous results [140] [155]. The GPS stations located on the Adriatic side of the Apennine belt show the highest horizontal velocity, with NNE-trending vectors reaching 6 mm/year ( Figure 5). In the study area (included in the red box shown in Figure 5), the horizontal velocity clearly changes from the Adriatic to Tyrrhenian side. In the latter zone, including Tuscany, vectors are mostly NE-directed with modest amplitude (<3 mm/year). The divergence of horizontal velocity vectors between the Adriatic and Tyrrhenian sides implies active extension in the northern and central Apennines, as discussed in detail by [140] [155]. However, the study area needs an accurate reconstruction of the horizontal strain field, as described in the next paragraph.

Horizontal Strain Field
The theory of small deformations allows us calculating the horizontal strain rate from the geodetic velocity field. For this work, I used the interpolation method developed by [140] [155].
First, I superimposed a regular grid to the GPS network. By using the GPS horizontal velocity vectors ( Figure 5), a weighted least square method provided six kinematic parameters for each point of the grid: 2 velocity components and 4 velocity gradients [156]. Then, the 4 components of the horizontal strain rate tensor (and the related eigenvectors) were computed from the velocity gradients [157].
The weight factors entering in the interpolation method depend on both the where d ij is the distance between the i-th grid point and the j-th GPS site, whereas D is the so-called attenuation distance [140] [157]. To improve the reliability of the results, I imposed two further geometric constraints [158] [159]: 1) the interpolated values are acceptable only if at least 3 GPS sites are located at a distance ≤D from the grid point 2) such GPS sites must be evenly distributed in the region surrounding the grid point (at least one for each 120˚ angular sector). If both these conditions are not met, the strain rate estimate is discarded. Figure 6 shows the strain rate field obtained by using a 0.1˚ Latitude × 0.1˚ Longitude regular grid. The adopted attenuation distance is D = 50 km, about twice the average distance between the GPS stations shown in Figure 5 [140] [155].
The amplitude of the principal strain rate axes varies in the range 10 −8 -10 −7 yr −1 ( Figure 6). Such values are lower than those observed at most active plate margins (>10 −7 yr −1 ). However, the study area deforms more rapidly than the adjacent foreland areas, such as the European continent north of the Alps, where the strain rate is much lower (10 −10 -10 −9 yr −1 [160]). Values of the geodetic strain rate comparable to those shown in Figure 5 are reported by [146] for the central and southern Apennines.
In the study area, the horizontal lengthening axis always has a NE-SW direction. Consequently, the shortening axis is directed NW-SE ( Figure 6). However, one can note important changes in the relative amplitude of the two eigenvectors. Along the Tuscan-Emilian-Romagna Apennine arc, lengthening is much larger than shortening, implying a pure (up to radial) extensional strain style (see box a in Figure 6). Instead, moving towards the Tyrrhenian coast the amplitude of the two axes becomes nearly equal. This highlights a strike-slip strain style (see box b in Figure 6). Furthermore, strain rate amplitude clearly decreases from the Apennine watershed to the Tyrrhenian coast ( Figure 6).
The geodetic strain field is broadly consistent with the strain style highlighted by ESMs (Figure 4). Indeed, shallow earthquakes along the Apennine arc (sectors I and II of Figure 4) show an extensional mechanism. On the other hand, most events occurred within inner Tuscany (sector III of Figure 4) have strike-slip mechanisms. Moreover, the direction of the seismological T and P axes is compatible with the principal axis of the geodetic strain rate tensor (Figure 4 and Figure 6).
The zone represented in Figure 6 is larger than the study area. To the north, the figure shows the Po Valley zone affected by the May, 2012 seismic sequence. In this area the geodetic strain is dominated by N-S compression (see box c in Figure 6), in good agreement with the EMSs of the 2012 seismic events [3]. This agreement corroborates the close relationships between seismic and geodetic strain, even outside the study area. International Journal of Geosciences Figure 6. Horizontal strain rate field, derived from the velocity field shown in Figure 5 by the procedure described in the text (see also [140] [155]). The zone is contained in the box shown in Figure 5. The principal axes of strain rate tensor are calculated on a 0.1˚ × 0.1˚ regular grid, with attenuation distance D = 50 Km, i.e. roughly twice the average distance between the GPS stations shown in Figure 5 (see text). Converging red arrows indicate the principal shortening axis, while diverging blue arrows indicate the principal lengthening axis. The amplitude scale (expressed in 10 −8 year −1 or 3.17 × 10 −16 s −1 ) is shown at the bottom left. Boxes a, b and c are mentioned in the text.

Discussion
The similarity of the seismic and geodetic strain fields suggests what strain regime affects the study area. Such regime refers to a short time interval (the last 10 -20 years), so it cannot be extrapolated easily to past tectonic phases. However, the above results are relevant for the identification of seismogenic processes which take place in the study area.
First, seismic and geodetic strain can be used to evaluate the reliability of the active faults proposed in the literature. Along the Tuscan-Emilian-Romagna Apennine arc, the alleged active faults agree with the ESMs in two cases only. In north-western Tuscany, the Northern Apuan Alps Fault (NA in Figure 2; [86] (Figure 4). In the Mugello basin, the Ronta and Sieve normal faults (Ro and Sv in Figure 2 [49]) seem to be compatible with the mechanism of the 09/12/2019 earthquake (M W = 4.54). However, for both cases the absence of coseismic faulting precludes a definitive confirmation of this hypothesis.
On the other hand, the suggested trans-pressional [49] or compressional [70] active deformation along the Tuscan-Emilian-Romagna Apennine watershed is denied by the extensional ESMs and geodetic strain (Figure 4 and Figure 6).
For inner Tuscany, both seismic activity and geodetic data show a roughly strike-slip strain style. This casts doubt on the present activity of the normal fault shown in Figure 2, some of which are considered as seismogenic by [63].
Much less plausible is the present activity of the compressive tectonic features shown in Figure  The results provided by this work may also serve as constraints for a coherent model of the current deformation of the Apennine belt. [161] suggest that the central and northern Apennines are affected by an overall strike-slip stress regime, with sub-horizontal NW-SE σ 1 and NE-SW σ 3 axes. This stress regime would explain the main features of seismic activity recognized by [131]. [162] [163] [164] propose an evolutionary model for the northern Apennines, based on a belt-parallel compression. This model may explain several geological and geomorphological features observed in the study area, as the Quaternary uplift and intense volcanism (Section 2). Furthermore, both extensional strain along the northern Apennine arc and strike-slip strain in the inner Tuscany are compatible with the above model [162] [163] [164].

Conclusion
The source mechanism of small earthquakes (M W < 5) helps to shed light on the active deformation of a given tectonic zone. However, in a limited time window  Figure 4 (the related parameters are reported in Table A1). Each panel indicates the earthquake number, date, epicenter zone, magnitude (M) and hypocenter depth (h). The sketch shown the projection of the nodal planes and the poles of the P and T principal axes of the seismic moment tensor [115]. Note that T and P axes bisect the right dihedral angles delimited by the nodal planes. [166] first showed that the above axes coincide with the principal strain axes. Therefore the source mechanism describes the instantaneous deformation at the earthquake hypocenter, with the maximum and minimum lengthening corresponding to the T and P seismological axes respectively. Contrary to what is sometimes (improperly) reported in the literature, a single source mechanism poorly constrains the direction of the principal stress axes [124] [166] [167].