Reﬁning Rates of Active Crustal Deformation in the Upper Plate of Subduction Zones, Implied by Geological and Geodetic Data: The E-Dipping West Crati Fault, Southern Italy

: We investigate crustal deformation within the upper plate of the Ionian Subduction Zone (ISZ) at different time scales by (i) reﬁning geodetic rates of crustal extension from continuous Global Navigation Satellite System (GNSS) measurements and (ii) mapping sequence of Late Quaternary raised marine terraces tectonically deformed by the West Crati normal fault, in northern Calabria. This region experienced damaging earthquakes in 1184 (M 6.75) and 1854 (M 6.3), possibly on the E-dipping West Crati fault (WCF) which, however, is not unanimously considered to be a seismogenic source. We report geodetic measurements of extension and strain rates across the strike of the E-dipping WCF and throughout the northern Calabria obtained by using velocities from 18 permanent GNSS stations with a series length longer than 4.5 years. These results suggest that crustal extension may be seismically accommodated in this region by a few normal faults. Furthermore, by applying a synchronous correlation approach, we reﬁne the chronology of understudied tectonically deformed palaeoshorelines mapped on the footwall and along the strike of the WCF, facilitating calculation of the associated fault-controlled uplift rates. Raised Late Quaternary palaeoshorelines are preserved on the footwall of the WCF indicating that “regional” uplift, likely related to the deformation associated either with the subduction or mantle upwelling processes, is affected by local footwall uplift. We show that GIS-based elevations of Late Quaternary palaeoshorelines, as well as temporally constant uplift rates, vary along the strike of the WCF, implying normal faulting activity through time. This suggests that (i) the fault slip rate governing seismic hazard has also been constant over the Late Quaternary, over multiple earthquake cycles, and (ii) our geodetically derived fault throw rate for the WCF is likely a more than reasonable value to be used over longer time scales for an improved seismic hazard assessment. Overall, we emphasize the importance of mapping crustal deformation within the upper plate above subduction zones to avoid unreliable interpretations relating to the mechanism controlling regional uplift. moreover, the location of this historical seismic event lies where we estimate the maximum co-seismic subsidence [29]. kinematic analysis) and E.S.; S.C.; data


Introduction
Destructive earthquakes with a magnitude higher than nine, such as the 1960 Chile earthquake (M 9.5) [1][2][3], 1964 Alaska earthquake (M 9.2) [1,4], 2004 Sumatra earthquake (M 9.1) [5,6], and 2011 Japan earthquake (M 9.1) [7], have occurred on tectonic plate boundaries associated with subduction of oceanic crust. Processes of subduction produce uplift of the overriding plate that can be recorded by flights of marine terraces that are cut into the land by the sea [8,9]. These terraces record the history of tectonics and sea-level changes over the Quaternary [10]. Data of uplift from overriding plates above active subduction  (Table 1). Green-coloured stars show that historical earthquakes occurred within the investigated area. Pink-coloured dots show the network of continuous GNSS stations. Inset A shows Quaternary active normal faults within the Calabrian-Peloritani Arc and the trace of the Ionian Subduction Zone (ISZ).
It is still unclear whether these seismic events have ruptured (or not) the West Crati Fault or other active faults across the strike. Empirical scaling correlations between fault size and expected earthquake magnitude may yet be used to suggest that the ~40 km mapped fault length of the E-dipping West Crati Fault could be capable of producing seismic events as large as Mw 6.7-7 [68,69]. However, this potential seismogenic source is still considered a "Debated Seismogenic Source" within the Database of Individual Seismogenic Sources (DISS) [30,31], even if it shows very similar fault geometry and size to other mapped "Individual Seismogenic Source", deforming recent sedimentary deposits [9,63]. Indeed, we bear in mind the concept that any evidence of fault movement since the Middle/Late Pleistocene has been accepted by some as the definition of an "active fault" [70,71]. We report that the West Crati Fault could be capable of producing earthquake ruptures affecting the topographic surface because previous authors show kinematic measurements taken from striated fault scarp with dip-slip normal movement, indicating slip at the surface [36].

Raised Middle-Late Quaternary Palaeoshorelines and Existing Age Controls in
Northern Calabria, in the Tyrrhenian Side Regional uplift is affecting the Calabrian-Peloritani Arc above the ISZ, likely due either to the collision/subduction processes or mantle upwelling beneath the crust of Calabria or a combination of these processes [38,39,72]. The main geological signature of in-  (Table 1). Green-coloured stars show that historical earthquakes occurred within the investigated area. Pink-coloured dots show the network of continuous GNSS stations. Inset A shows Quaternary active normal faults within the Calabrian-Peloritani Arc and the trace of the Ionian Subduction Zone (ISZ) [29].

Raised Middle-Late Quaternary Palaeoshorelines and Existing Age Controls in Northern Calabria, in the Tyrrhenian Side
Regional uplift is affecting the Calabrian-Peloritani Arc above the ISZ, likely due either to the collision/subduction processes or mantle upwelling beneath the crust of Calabria or a combination of these processes [38,39,72]. The main geological signature of interaction between regional uplift and active faulting within the Calabrian-Peloritani Arc is the presence of dramatic flight of marine terraces, formed since the Middle Pleistocene [9,63,[73][74][75][76][77][78][79][80][81]. We here focus on evidence that tectonically deformed sequences of palaeoshorelines have been recording the effects of Middle-Late Pleistocene normal faulting activity within the Calabrian-Peloritani Arc [9,63,80,82,83]. Moreover, as previous authors proposed in the past, we recognize on DEM, using 10 m high resolution DEMs [84], seaward-sloping marine surfaces as terraced surfaces cut by wave action, updip bounded by palaeo sea cliff, with their shoreline angle (or inner edge) identifying the palaeo mean sea level at the formation of the terrace [85][86][87][88]. In particular, within the footwall of the West Crati Fault, these palaeoshorelines can be traced parallel to the present-day coastline between the towns of Coreca and Cetraro (Figures 2 and 3). These marine terraces are carved either into calcarenites and clays deposits of Upper Miocene age or in places into the Palaeozoic crystalline basement [41]. ward-sloping marine surfaces as terraced surfaces cut by wave action, updip bounded by palaeo sea cliff, with their shoreline angle (or inner edge) identifying the palaeo mean sea level at the formation of the terrace [85][86][87][88]. In particular, within the footwall of the West Crati Fault, these palaeoshorelines can be traced parallel to the present-day coastline between the towns of Coreca and Cetraro (Figures 2 and 3). These marine terraces are carved either into calcarenites and clays deposits of Upper Miocene age or in places into the Palaeozoic crystalline basement [41].  Age controls associated with this sequence of marine terraces are limited, and in this section, we review existing constraints as shown in Table 1. Furthermore, these ages are slightly located north compared to the investigated region where local tectonic effect from active faults is less prominent, in Scalea and Grotta del Prete towns ( Figure 1). Indeed, these age controls lie between the southern fault tip and the hanging wall of the "Maratea Fault" [89], likely suggesting a quasi-negligible "hangingwall subsidence" local effect. For instance, in Grotta del Prete town, an age of 329 (+196 and −64) ka for the associated marine terrace at 14-22 m.a.s.l has been firstly proposed, implied by U/Th dating on corals (LB 96 sample- Table 1 from Carobene et al., 1991) [90]. However, we discard this age because new investigations have demonstrated that those corals showed re-crystallization, suggesting a younger age of 83.8 ± 3.6 ka (GRP2 sample from [56]) for the same marine terrace [91]. Furthermore, close to Scalea town, a marine terrace mapped at 16 m is assigned to the Marine Isotope Stage 5c (100 ka) implied by U/Th dating measurements on corals (SLC sample- Table 2 from [56]). The latter samples of U/Th absolute dating (GRP2 and SLC) permit suggestion of an averaged uplift rate of~0.45 mm/yr over the Late Pleistocene [91] for a region just a few kilometres away, north of our investigated area ( Figure 1). Furthermore, these latter samples show robust activity ratios for U/Th dating on corals (δ234 Ui) with values of 190 (GRP2) and 335 (SLC), in agreement with present-day oceanic δ234 Ui value of~145 [19,73].
Remote Sens. 2022, 14, x FOR PEER REVIEW 6 of 34 Figure 3. Evidence of raised marine terraces and their associated palaeoshorelines: in (a) three prominent sea-level highstands such as the 125 ka (MIS 5e), the 240 ka (MIS 7e), and the 340 ka (MIS 9e), which are mapped from Google Earth just south of Paola town, close to our Profile 11; in (b) a modified photo taken from Westaway (1993), showing a sequence of marine terraces from field close to Profile 1, close to Coreca town. Inset A shows a sketch of conceptual evolution when marine terraces form.
Age controls associated with this sequence of marine terraces are limited, and in this section, we review existing constraints as shown in Table 1. Furthermore, these ages are slightly located north compared to the investigated region where local tectonic effect from active faults is less prominent, in Scalea and Grotta del Prete towns ( Figure 1). Indeed, these age controls lie between the southern fault tip and the hanging wall of the "Maratea Fault" [89], likely suggesting a quasi-negligible "hangingwall subsidence" local effect. For instance, in Grotta del Prete town, an age of 329 (+196 and −64) ka for the associated marine terrace at 14-22 m.a.s.l has been firstly proposed, implied by U/Th dating on corals (LB 96 sample- Table 1 from Carobene et al., 1991) [90]. However, we discard this age because new investigations have demonstrated that those corals showed re-crystallization, suggesting a younger age of 83.8 ± 3.6 ka (GRP2 sample from [56]) for the same marine terrace [91]. Furthermore, close to Scalea town, a marine terrace mapped at 16 m is assigned to Figure 3. Evidence of raised marine terraces and their associated palaeoshorelines: in (a) three prominent sea-level highstands such as the 125 ka (MIS 5e), the 240 ka (MIS 7e), and the 340 ka (MIS 9e), which are mapped from Google Earth just south of Paola town, close to our Profile 11; in (b) a modified photo taken from Westaway (1993), showing a sequence of marine terraces from field close to Profile 1, close to Coreca town. Inset A shows a sketch of conceptual evolution when marine terraces form.
This lack of absolute age constraints for outcropping yet mapped marine terraces within the footwall of the West Crati normal fault is a common problem for regions affected by relatively low uplift rates [9,19,63,88]. For this reason, it is important to use the knowledge of dated palaeoshorelines to estimate ages of other undated palaeoshorelines. In the past, several authors have tried to estimate uplift rates by applying a "sequential correlation" approach [85,[92][93][94][95][96]. This approach is based on the concept that, given a dated palaeoshoreline, the next higher and older palaeoshoreline should belong to the next Remote Sens. 2022, 14, 5303 7 of 31 older glacio-eustatic sea-level highstand; this is shown in previous investigations using raised marine terraces with relatively low uplift rates [95][96][97]. Yet, this approach is prone to erroneous uplift rate derivations because it does not take into account the well-known "re-occupation or overprinting problem" [38] where older marine terraces are re-occupied by subsequent younger and higher palaeo-sea-level highstands, identifying a new marine terrace [9,19,20,63,73,88]. If this is not taken into account, it may lead to incorrect age assignments for undated marine terraces, with likely erroneous uplift rate derivations through time [9]. However, previous investigations have shown a different approach, named "synchronous correlation" approach, whereby given an uplift rate driven by age controls, it allows calculations of all expected/predicted elevations of glacio-eustatic sea-level highstands identified by marine terraces [19,98]. These predicted elevations are correlated with all GIS-based mapped and measured palaeoshoreline elevations synchronously rather than sequentially; this permits overcoming of the "overprinting problem", allowing more robust scenarios for uplift history to be provided, as widely shown by some within the Mediterranean realm [9,19,20,63,73].

Methods
In this section, we present how we investigate (i) refining changes in space and time of uplift and the associated rates along and across the strike of the West Crati Fault and (ii) deriving strain and extension rates from GNSS measurements for investigated region in northern Calabria. This was aimed at better understanding if the regional uplift is affected by local active faulting and if different timescale deformation rates may lead to a refined seismic hazard assessment for the West Crati normal fault and across the strike of the northern Calabria. Below, we describe methodologies used in this study.

Elevation Data from DEM-Based Topographic Analysis of Mapped Marine Terraces
We have carried out detailed GIS-based topographic and geomorphological analysis by using 10 m high resolution digital elevation model (DEM) provided by INGV TINITALY [84]. This approach has been widely proposed where raised marine terraces are used for tectonic and crustal deformation investigations in the Mediterranean realm [9,20,63,73,87]. In particular, we have produced several topographic profiles perpendicular both to the actual coastline and the strike of the fault over the~40 km extent of the fault (Figures 2 and S1-S4), where each profile intercepts a sequence of marine terraces outcropping on the fault footwall. We have mapped inner edges (or shoreline angles) identifying palaeoshorelines ( Figure 3 and Inset A), representing glacio-eustatic palaeo-sea level associated with a highstand. These are confirmed in many locations by the existence of prominent breaks of slope identifying palaeoshoreline inner edges on the DEM, where scarp-like palaeo-sea cliffs bound a seaward sloping flat marine updip. Palaeoshoreline indicators and the topography are used to depict the location of each topographic profile; for instance, areas affected by fluvial incision caused by the coupled action of rivers and uplift are avoided to ascertain that we investigated geomorphological features associated with marine and not fluvial geomorphological features [63]. Indeed, our topographic profiles are produced along uplifting interfluvial areas where marine terraces are better preserved and mapped.

Synchronous Correlation Approach to Model Multiple Glacio-Eustatic Sea-Level Highstands and Multiple Raised Marine Terraces
This method has been widely described by previous studies where uplifted marine terraces are used for tectonic and seismic hazard investigations in the Mediterranean realm and elsewhere, outlined in [9,19,20,63,73,86,88,98]. This approach is based on the concept that glacio-eustatic sea-level highstands, which combined with tectonic uplift are thought to produce marine terraces [10], are unevenly spaced in time, indicating that for a constant uplift rate over time, uplifted marine terraces will be preserved in the landscape unevenly spaced in elevation [19,98]. It is important to highlight that this approach also works if uplift rates change through time as demonstrated by previous studies in Greece and southern Italy [19,99]. Furthermore, re-occupation of marine terraces within regions showing relatively low uplift rates may occur for both constant and fluctuating uplift rate scenarios. Indeed, this occurs where palaeo-sea level at a highstand lies beneath the level of present-day sea level because of reduced volume of global ocean at that particular time, combined with a low uplift rate [73]. These coupled processes imply that younger sealevel highstands with larger volume of global ocean and higher palaeo-sea-level elevation overcome and re-occupy older palaeoshorelines during sea-level rise [38,73]. We report that the synchronous correlation technique overcomes this problem due to the fact that all palaeoshoreline elevations are calculated on topographic profiles, yet sequential correlation tends to mislead for uplift rate calculations in this scenario if the "re-occupation" problem is not taken into account [19,73]. Furthermore, age constraints drive the uplift rate iteration, and it is desirable to have at least one dated palaeoshoreline across the investigated area. This allows us to test the first and simplest assumption of having a constant uplift rate through time and evaluate if the iteratively calculated glacio-eustatic sea-level highstand elevations match the measured and mapped palaeoshoreline elevations [73]. If a robust correlation is not found for the simplest scenario, a changing uplift rate through time is investigated, with iteration of uplift rates driven by age constraints aimed to find the best match with DEM-mapped and measured palaeoshoreline elevations [73]. For the calculation of glacio-eustatic sea-level elevations, well-known sea-level curves for our palaeoshoreline modelling are used (Table 2) [100,101], with negligible differences in terms of uplift rate derivation if different sea-level curves are used [20].
Our topographic analysis of marine terraces is carried out by constructing topographic profiles intercepting the palaeoshorelines from DEM data ( Figure 2). In order to perform the synchronous correlation, iteration of uplift rates is enabled by the so-called "Terrace Calculator" built in Excel already available in the literature where data from sea level highstands derived by well-known sea level curves are input [9,63,73]. As previously shown, this "Terrace Calculator" allows using uplift rate (UR) as input data, which is iterated to compute the predicted elevations (PE) of all sea-level highstands using the age of the highstands (T) and the sea-level elevations of each highstand (SLE) using the formula PE = (UR × T) + SLE [9,20,63,73]. In more detail, rates of uplift are iterated for the entire topographic profile, constrained where available by dated palaeoshorelines, in order to find the best match between multiple measured palaeoshoreline elevations and multiple predicted sea-level highstand elevations (within a range of ±10 m). This allows for assigning predicted ages of glacio-eustatic sea-level highstands to undated palaeoshorelines. We then perform linear regression analysis, forcing the R 2 value to maximize, to quantify if a correlation exists between the "measured" (or mapped) elevations of all palaeoshorelines from DEM data and the "predicted" (or expected) elevations identifying sea-level highstand elevations implied by iteration of uplift rates driven by age controls and by means of wellaccepted sea-level curves for the last~1 Ma (Table 2) [100,101]. Finally, we replicate this iteration for each topographic profile, constructed across the strike of the marine terraces and along the entire strike of the fault in order to determine the extent of the deformation caused by the fault.

GNSS Data Analysis and Strain Rate Computation
We have studied the crustal deformation for the investigated area by using GNSS measurements. In northern Calabria, there are 32 regionally well-distributed permanent GNSS stations. However, only 19 of them can be considered for our study because they must show a time series of at least 4.5 yr, which is the minimum time-series length compatible with reliable results, as suggested by [49]. Moreover, one station, KROT, is placed on an unstable slope; thus, it cannot be used to estimate strain rates and crustal extension, as already outlined by previous studies [47,52]. For these reasons, we have used 18 GNSS stations for our study, whose main details are summarized in Table 3, and, including the WGS84 coordinates (EPSG 4326), in Table S1 (Supplementary Material). We also report that for the investigated region in northern Calabria, there is no reliable time series from non-permanent GNSS stations to our knowledge.
The time series used for this study were the 24 h final solutions available at Nevada Geodetic Laboratory (NGL) database [102]. These time series have been processed by NGL by means of JPL GipsyX v1.0 software on the basis of several parameters such as nominal troposphere, troposphere delay models from VMF data server, elevation-weighted observations, higher order ionospheric calibrations, improved JPL Repro 3 orbits, and the global reference frame IGS14. Since northern Calabria is close to the southern margin of the Eurasian plate, two plate-fixed reference frames were considered: (i) Eurasian Plate (EU) and African Plate (AF). In both cases, the plate rigid motion has been removed using the Euler poles defined previously by [103], allowing the recognition of plate rotation effects.
We have processed the EU and AF time series for each considered station after having obtained daily positions from NGL, in order to obtain horizontal velocity components in east and north directions and the associated standard deviations (SDs). In particular, for each station and each direction, we have estimated parameters such as (i) the seasonal and half-seasonal components and (ii) models for chosen noises such as white noise and power law noise by means of the maximum likelihood estimation (MLE) method implemented in HECTOR software package [50]. It is important to note that before estimating the abovementioned parameters, we carried out offset recognition and outlier removal. Specifically, we defined measurements outside three times the interquartile range (IQR) as outliers and therefore removed them. Table 3. Some data about GNSS stations in the studied area. More data are shown in Table S1 which can be downloaded as Supplementary Material.

GNSS Station
Time Series UTM33N Coordinates (m) Once the velocities of the GNSS stations and the corresponding standard deviations were obtained, we then computed the horizontal strain rate field on a regular grid by means of the modified least square (MLS) approach implemented in grid_strain MATLAB toolbox [104]. The contribution of a GNSS station to the strain-rate estimation in a grid point is weighted by exp(−d/d 0 )/σ 2 v , where σ v is the standard deviation (SD) of the station velocity, d is the station-node distance, and d 0 is a scale factor introduced to define the level of locality of the estimation. This allows for carrying out different analyses focusing, for instance, on the area across either one or more of the faults or on the entire region. We provide, for each scale factor and for each node of the chosen grid: (i) the strain rate, i.e., maximum and minimum principal strain rates ∆ neglects nature and geometry of the strain tensor and therefore allows an evaluation of the deformation pattern inside an area having a specific behaviour (extensional or compressional). Instead, the normalized engineering shear rate, . γ, also provides information on the sign of both the principal strains and therefore is particularly suitable to the recognition of the boundaries between areas characterized by different kinematics. Our analysis also provides the uncertainties of the computed parameters and a map of the geometric significance of the results. In particular, the results obtained in a grid node are considered to have high geometric significance if the GNSS stations within the scale factor are distributed on all four 90 • quadrants around this node, whereas they are considered to have mid geometric significance if no more than three quadrants are filled. In all the other cases, the results are considered to have low significance. Since the strain rate does not depend on the choice of a specific reference frame, it was computed with the EU velocities only.

Vertical Co-Seismic Deformation Calculation to Estimate Earthquake Scenarios for the West Crati Fault
In order to (i) investigate possible earthquake scenarios associated with the co-seismic rupture of the entire length (~40 km) of the West Crati Fault and (ii) test if a correlation may exist between co-seismic and long-term crustal deformation, we carried out elastic half-space modelling using Coulomb 3.4 software [105]. We used a modified trace of the fault from [37], according to the geomorphology mapped using Google Earth and high resolution DEMs. The mapped fault trace was used as input information into the modelling using a new MATLAB code [106] into Coulomb 3.4 software proposed in the past for other investigations where curvilinear faults are used to derive co-seismic vertical deformations and Coulomb stresschange calculations, producing corrugated 3D fault surface [73,87,106,107]. Indeed, this code allows for calculations of vertical and horizontal displacements arising from seismic events rupturing active faults with variable-strike geometry [107]. It is important to note that it is well-recognized that geometry of active faults affects Coulomb stress transfer [106,108,109], consequently influencing strain and displacements nearby the causative fault. We then used an averaged value of fault dip angle of 45 • for the West Crati Fault, based on several field data provided in the literature [28,36,37,110]. We also used a concentrical slip distribution according to well-known scaling relationships between earthquake magnitude, slip at depth, and fault length [68], as already proposed by some in southern Italy [73,107], to obtain a modelled earthquake magnitude of~Mw 6.7-6.8, which are values that match magnitudes of the historical earthquakes, within a margin of error, shown in Figure 2. Finally, we correlate co-seismic and long-term crustal deformation along the N-S-oriented West Crati normal fault.

Results
In the following, we firstly show that our synchronous correlation technique advocated in the previous section allows us to replicate measured palaeoshoreline elevations by finding the best match with predicted sea-level highstand elevations. We attempt to correlate elevations of palaeoshorelines mapped on the uplifting footwall and along the strike of the West Crati Fault. We therefore use the derived temporal and spatial constraints related to the geometries of palaeoshorelines to derive long-term rates of uplift over the Late Quaternary and how these are associated with displacements on the West Crati normal fault. We finally show refined short-term rates of extension and strain from new GNSS measurements and how these rates relate to the long-term fault-controlled crustal deformation implied by investigating raised marine terraces.

Synchronous Correlation Modelling to Refine Ages for Undated Marine Terraces
We present our synchronous correlation between DEMs-based palaeoshoreline elevations and those predicted by iterating rates of uplift (Figures 4 and S5-S10). It is important to note that this iteration is driven by age constraints ( Table 1) identified north of the northern fault tip of the West Crati Fault ( Figure 1); they lie between the hanging wall yet very close to the southern fault tip of the Maratea Fault and over its southern tip, suggesting minimum to negligible deformation from local upper-plate faulting. Indeed, we have iterated rates of uplift to find the best match between mapped and predicted elevations of marine terraces following a few criteria already proposed by previous studies [9,63], such as (i) prominent mapped marine terraces were iteratively assigned to the most common sea-level highstands mapped worldwide (125 ka-MIS 5e, 240 ka-MIS 7e, and 340 ka-MIS 9e, Figure 3) [9,63,73,111] and (ii) maximizing the coefficient of determination, R 2 value, showing how well other less clearly mapped marine terraces fit the predicted sea-level elevations. For the first time, refined ages are proposed (Table 4) for mapped yet undated raised marine terraces outcropping on the footwall of the West Crati Fault by applying a synchronous correlation approach, which takes into account the "re-occupation" problem for relatively low uplifting regions such the Calabrian-Peloritani Arc [9,63]. This permits us to investigate marine terrace elevations and associated uplift rates along the strike of the West Crati Fault. Note that we have used 10 m high-resolution DEMs as proposed previously for other studies within the Calabrian-Peloritani Arc [9,63,112] where excellent correlations have been shown between marine terraces mapped on DEMs and in the field, suggesting a reliable robustness of our GIS-based marine terraces mapping for this study. undated raised marine terraces outcropping on the footwall of the West Crati Fault by applying a synchronous correlation approach, which takes into account the "re-occupation" problem for relatively low uplifting regions such the Calabrian-Peloritani Arc [9,63]. This permits us to investigate marine terrace elevations and associated uplift rates along the strike of the West Crati Fault. Note that we have used 10 m high-resolution DEMs as proposed previously for other studies within the Calabrian-Peloritani Arc [9,63,112] where excellent correlations have been shown between marine terraces mapped on DEMs and in the field, suggesting a reliable robustness of our GIS-based marine terraces mapping for this study.  Table 4. Table 4. Mapped inner edges from DEM analysis with refined ages for marine terraces assigned via synchronous correlation.
Profile Number (Sea-Level Highstand Referred to Figure 2, Figure S1- Figure 2, Figure S1- Figure 5 shows a very robust synchronous correlation between our multiple mapped palaeoshoreline elevations and the predicted ones ( Figure 4, Table 2, and Spreadsheet S1) within a margin of error, suggesting we have obtained consistent uplift rate estimates. This is a crucial key observation because it also suggests that rates of uplift mapped on the footwall of the West Crati Fault have been constant through time. Furthermore, we have mapped 17 order of palaeoshorelines, in agreement with previous studies, for the Calabrian-Peloritani Arc [37,38,42,113] associated with the sea-level highstands as old as 100-980 ka (not all mapped within a single topographic profile), as previously recognized elsewhere within the Mediterranean realm [9,19,20,63,73,88]. Evidence in places of marine terraces as old as 980 ka confirm that the investigated region has been uplifting since then, in agreement with previous investigations [37,38]; however, we cannot rule out that the onset of the uplifting process may likely be older. This is in contrast with the previous proposed onset at 700 ka for the "regional" uplift of the Calabrian-Peloritani Arc [38].  Figure 5 shows a very robust synchronous correlation between our multiple mapped palaeoshoreline elevations and the predicted ones ( Figure 4, Table 2, and Spreadsheet S1) within a margin of error, suggesting we have obtained consistent uplift rate estimates. This is a crucial key observation because it also suggests that rates of uplift mapped on the footwall of the West Crati Fault have been constant through time. Furthermore, we have mapped 17 order of palaeoshorelines, in agreement with previous studies, for the Calabrian-Peloritani Arc [37,38,42,113] associated with the sea-level highstands as old as 100-980 ka (not all mapped within a single topographic profile), as previously recognized elsewhere within the Mediterranean realm [9,19,20,63,73,88]. Evidence in places of marine terraces as old as 980 ka confirm that the investigated region has been uplifting since then, in agreement with previous investigations [37,38]; however, we cannot rule out that the onset of the uplifting process may likely be older. This is in contrast with the previous proposed onset at 700 ka for the "regional" uplift of the Calabrian-Peloritani Arc [38].  Figure 6a shows spatial changes of uplift along the strike of the West Crati Fault, which have been recorded by the geometry of the marine terraces. In particular, our synchronous correlation shows that rates of uplift are constant over the Middle-Late Pleistocene yet change spatially, with values of uplift rates between 0.87 mm/yr in the centre of the fault footwall, diminishing towards the fault tips and beyond where the uplift rates are 0.4 (southern tip) and 0.52 (northern tip) mm/yr (Figure 6b). The latter value is not surprising if we consider that [56] show an averaged value of ~0.45 mm/yr where a quasi- Figure 5. Graph showing a linear regression analysis between iteratively calculated sea-level highstand elevations and our mapped palaeoshoreline elevations on DEMs. Note that the predicted elevations have been derived by defining a constant uplift rate through time and iterating this value to find the best match to the measured and mapped palaeoshorelines. Figure 6a shows spatial changes of uplift along the strike of the West Crati Fault, which have been recorded by the geometry of the marine terraces. In particular, our synchronous correlation shows that rates of uplift are constant over the Middle-Late Pleistocene yet change spatially, with values of uplift rates between 0.87 mm/yr in the centre of the fault footwall, diminishing towards the fault tips and beyond where the uplift rates are 0.4 (southern tip) and 0.52 (northern tip) mm/yr (Figure 6b). The latter value is not surprising if we consider that [56] show an averaged value of~0.45 mm/yr where a quasi-negligible effect may be considered from the hanging-wall tectonic subsidence due to the Maratea Fault activity in the northern area. Furthermore, our refined uplift rates are in agreement with previous studies where values of~0.86 mm/yr had been proposed for this region since the end of the Early Pleistocene (700 ka) [114]. It is important to note that beyond the fault tips, no active faults are mapped, yet a "regional" signal of uplift is recorded (Figure 6b), suggesting that mapped and uplifted marine terraces on the footwall of the West Crati Fault represent a combination of the "regional" uplift signal and the "local" fault-controlled uplift. Moreover, we have tested whether the tilting process of palaeoshorelines (Figure 6a) has occurred over the Middle-Late Pleistocene or anachronously to the formation of all marine terraces. Following an approach from previous investigations [9,20,63], we assess tilting values along the strike of the fault and if these change for different ages of marine terraces; if faulting activity by the West Crati Fault has occurred with time, we would assume older marine terraces to be more tilted compared to younger ones, along the strike of the fault. Indeed, Figure 6c shows that higher and older marine terraces present higher tilt angle values, implying that they have experienced a longer history of faulting activity with a constant rate, a key observation shown by others when marine terraces are used for deriving faulting activity through time [9,20,63].  showing that older palaeoshorelines have higher tilt angles; this suggests that they have experienced a longer history of differential uplift and that differential uplift was ongoing progressively during the Late Quaternary. Note that values of tilt angle for each investigated marine terrace have been calculated as a tan −1 of a gradient m of a straight line (equation y = mx), as proposed by previous studies [20,63,73].

Uplift and Uplift Rates along the Strike of the West Crati Normal Fault
We interpret the spatially changing (i) palaeoshoreline elevations, (ii) uplift rates and the temporally varying tilting angle values along the strike as evidence to indicate that the West Crati Fault has been active, with a constant rate of deformation, over time during the progressive formation of the investigated sequence of marine terraces. Moreover, it is not surprising to show that higher uplift rates allow the preservation of a higher number of marine terraces (Figure 6a,b), which is a diagnostic stratigraphic evidence to show variations along the strike in fault activity rate, in agreement with other regions deformed by active normal faults [9,19,20,63,85]. Figure 7a shows the MLE-based horizontal velocity vectors for both EU and AF plate fixed-reference frames, including the corresponding error ellipses, whereas Figure 7b shows a contour plot of the MLE-based vertical velocity. We interpret the spatially changing (i) palaeoshoreline elevations, (ii) uplift rates and the temporally varying tilting angle values along the strike as evidence to indicate that the West Crati Fault has been active, with a constant rate of deformation, over time during the progressive formation of the investigated sequence of marine terraces. Moreover, it is not surprising to show that higher uplift rates allow the preservation of a higher number of marine terraces (Figure 6a,b), which is a diagnostic stratigraphic evidence to show variations along the strike in fault activity rate, in agreement with other regions deformed by active normal faults [9,19,20,63,85]. Figure 7a shows the MLE-based horizontal velocity vectors for both EU and AF plate fixed-reference frames, including the corresponding error ellipses, whereas Figure 7b shows a contour plot of the MLE-based vertical velocity.  As already stated, the station KROT (point 32 in Table 3 and Figure 7a) is affected by a motion mainly due to local slope instability; therefore, it cannot be used for the strainrate computation at a regional scale even if the corresponding time series is 9 yr long. In all figures, the WGS84/UTM33N (EPSG 32633) coordinate system is used. Moreover, in Figures 7 and 8, only those stations whose time series have a duration of at least 4.5 years are shown.  The distance between the stations varies from 5.6 km to 131 km, and the mean distance between the Crati normal faults varies from 12 km to 16 km. Moreover, the stations across the area between these faults are~20-35 km away from it. In order to have a spatial resolution of the strain-rate field compatible with the minimum distance, a 4 km side grid was chosen. Moreover, the chosen scale factor was 40 km in order to highlight the effect of the stations close to the fault and to reduce that of very distant stations. In Figure 8a, we show the principal directions of the strain rate in the area characterized by high geometric significance, whereas in the large-sized Figure S13 (Supplementary Material), besides the strain-rate tensors in the high geometric grid nodes (solid vectors), the ones in the mid geometric significance nodes are also shown (dotted vectors). The strain rate is nearly everywhere characterized by extension, prominently E-W-oriented. In particular, we show that in the northern part of the Crati Valley the maximum strain rate is~19 ± 3 nstrain/yr with an angle with respect to the north direction in the range between −88 and −90 • ; whereas in the southern part of this valley the maximum strain is~18 ± 2 nstrain/yr, with an angle with respect to the north direction from −81 to −83 • . The extension rate in the area between the two faults is 0.2-0.3 mm/yr.

Strain Rate, Extension Rate, and Vertical Rates from GNSS Measurements
The map of the field of the rate of change in area, also called dilatation (i.e., the sum of principal strain rates, which is the trace of the strain rate tensor and is the first invariant of strain), and of the field of engineering shear rate normalized to the change in area are shown in Figure 8b,c, respectively, where the faults are also shown. Both the plots show the results in the high geometric significance areas. This map highlights that, in the region, there are no areas characterized by peculiar kinematics and that the behaviour is generally homogenous within the investigated region, with subtle variations. In particular, the kinematics of the Crati Valley seems to be similar if compared to neighbouring areas, even if in a N-W portion of the valley the results are characterized by no more than mid geometric significance (i.e., there are stations within the scale factor in three 90 • quadrants).
Finally, since the maximum width of the studied area, throughout northern Calabria, is~90 km, we choose a 100 km scale factor, which allows us the estimation of the corresponding regional strain rate. Indeed, we estimate a mean strain rate of~21 nstrain/yr, which results in an extension rate at a regional scale of~2.3 mm/yr in the direction 81-83 • .
Finally, we carried out the analysis of the distribution of the vertical velocity components of the 18 considered stations. The corresponding map, computed by linear interpolation of these velocities, is shown in Figure 7b. We highlight that even though the expected behaviour over the Late Quaternary in the area is a widespread uplift [38,110], short-term subsidence is observed in the western coastline of northern Calabria and within the Crati valley. However, it should be noted that, as shown in Table 1, the uncertainties on the vertical velocity components are relatively high. In particular, the uncertainties related to the stations GRI1, AQSA, LUZZ, CELI, CCRI, STSV, and CUTR (affected by uplift) and ROS4, CAR8, and AMA1 (affected by subsidence) have the same order of magnitude of the corresponding values.
This does not happen for the horizontal components, for which the relative error is on average 7%, ranging from 1.5% to 15%.

Discussion
In this study, we have refined rates of crustal deformation, at different scales, within the uplifting Calabria-Peloritani Arc, throughout the northern Calabria region, by (i) synchronously correlating multiple raised Middle-Late Pleistocene marine terraces and (ii) computing new strain, extension, and vertical rates from GNSS measurements (Figures 4-8). This has permitted us to study the relationship between the "regional" signal of uplift and the "local" effect of faulting activity along and across the strike of the West Crati Fault. We stress that the obtained results agree with the ones obtained by D'Agostino et al. (2011). In particular, there is an agreement between the strain rates computed in the Crati valley both in intensity (18-19 nstrain/yr) and direction (almost E-W). However, our study has used a higher number of stations (18 instead of 12) and longer time series. In particular, as already stated (Table 3), all the time series had a duration of at least 4.5 yr, and 10 out of 18 had a duration of at least 8 yr. The mean duration of all the used time series was 10.8 yr. The percentage of gaps of the used time series was generally in the range 5-12%. Previous authors have demonstrated that a time series should have at least 4.5 yr of duration to provide meaningful time series for geodetic investigations, and 8 yr of minimum duration should be recommended [49]; this is an advance from the viewpoint of GNSS data analysis.
In this section, we now discuss our results in terms of local and regional tectonic implications associated with refined throw rates, T mean , correlation between potential megathrust earthquakes investigations and intracrustal upper-plate seismic deformation, and possible co-seismic deformation scenario for the West Crati normal fault.

Local Tectonic Implications: Estimating Geodetic Fault Throw Rate and ERI for the West Crati Fault
The investigated area shows clear evidence that rates of uplift have changed along the strike of the West Crati Fault since the Middle Pleistocene, indicating that this fault has to be considered active, if the definition of "active fault" by some [70,71] is followed. Furthermore, this is not surprising given that this region has been historically hit by damaging seismic events such as the 1184 "Valle del Crati" earthquake (M 6.75) and the 1854 "Cosentino" earthquake (M 6.3) (Figure 1), which may have likely ruptured the West Crati normal fault. However, it is important to note that we cannot rule out that the abovementioned seismic events may have occurred on the antithetic, W-dipping, East Crati Fault. This would suggest that more investigation is needed, even though for this study this is beyond the scope of the paper.
In this study, it is important to highlight that we are not able to derive direct measurements of long-term fault throw rates from faulted marine terraces, and for this reason we estimate geodetically derived fault throw rates from our geodetic refined strain and extension rates throughout northern Calabria (Figures 7 and 8). Indeed, we know that a seismically active crustal extension accommodated by normal fault systems, between tectonic plate boundaries, needs to accommodate regional deformation rates [115]. This implies that summed deformation rates (throw rates on faults) across all the involved fault segments (i) will equal the regional deformation rates associated with the plate motions and (ii) will be temporally constant if all the same number of faults remains active over time [115]. However, it has also been demonstrated that throw rates on faults can change through time within normal fault systems due to fault interaction processes; this would imply that some faults would accelerate and others would decelerate, changing their throw rates through time to equal the regional deformation [19,32,115,116]. It is important to note that in places within the Calabrian-Peloritani Arc, rates of crustal deformation have been studied at three different timescales, with negligible discrepancies: those measured (i) over a few tens of years (from GNSS measurements), (ii) over a few thousand years (Holocene palaeoseismology), and (iii) over tens to hundreds of thousands of years (deformed Quaternary marine terraces) [9,46,117,118]. This would suggest long-term constant rates of crustal deformation confirmed by some [9,63]. Our results show a constant fault throw rate through time for the West Crati Fault from refined constant fault-controlled uplift rate evidence by modelling uplifted and deformed marine terraces, a crucial key observation for detailing the geography of the seismic hazard in northern Calabria. This suggests that fault interaction across the strike of active faults within the investigated region may be likely ruled out. This may be in contrast with the view of some proposing fault interaction in northern Calabria from fault-controlled uplift rates between 100 and 300 ka, where the timing of changes may be relatively imprecise, as results came from an inversion of river profile data and not absolute dating of deformed horizons [119]. If we then assume that strain and extension are homogenously accommodated across the strike of active faults such as the West and East Crati Faults and the Lakes Fault (Figure 9) with fault throw rates constant through time, deforming northern Calabria, estimates of geodetically derived fault throw rates may likely be more than reasonable values to be used over longer time scales. Note that estimates of geodetic fault throw rates are derived by using a qualitative trigonometric geometrical approach where the heave (H, horizontal component) and the throw (T, vertical component) are equal if the fault dip angle (α) is assumed to be 45 • (Figure 9). This would suggest that our regional, geodetically derived extension rate of 2.3 mm/yr accommodated by all three before-mentioned active normal faults with an averaged fault dip angle of 45 • would imply throw rates of~0.8 mm/yr (Figure 9, inset A). This value is not surprising, taking into account that throw-rate values measured for other active normal faults that have been accommodating the Plio-Pleistocene crustal extension along the Italian Apennines and the Calabrian Arc are within the range 0.3 to 2.0 mm/year [9,116,120,121].  [29]. Inset A shows a no-scaled cartoon where homogenous heave rates (HR) from computed extension rates are assigned among the active faults such as the West Crati Fault, East Crati Fault, and the Lakes Fault. Inset B shows a scheme of a normal fault with the three components where the deformation is measured and the fault dip angles: heave, throw, and displacement. Note that for a 45° fault dip, angle, throw, and heave will be equal.

Regional Implications for Upper-Plate Crustal Deformation above Subduction Zones
Our results from marine terraces modelling to derive refined long-term fault-controlled uplift rates show that the "regional" uplift from either ongoing subduction or deep-seated (mantle upwelling) processes has been spatially affected by "local" faulting activity; this has been observed throughout the Calabrian-Peloritani Arc [9,63] and should be considered when attempts for seismic hazard assessment for subduction-related megathrust earthquakes are studied. This is crucial if we consider that some studies show that the Ionian subduction may be active [124,125] and then capable of producing megathrust earthquakes, even though it is important to highlight there is an ongoing debate regarding if this subduction zone is either seismically coupled or aseismic. For instance, a previous study has shown an apparent "decoupling" about the crustal deformation occurring on the Calabrian upper plate and deformation on this subduction zone; this is also suggested by the fact that there is no record of strong subduction earthquakes, with a relatively aseismic "zone" at least in the first 30-60 km [112,124,126]. This would agree with the idea proposed by [39,127] of having (i) rising asthenosphere just beneath the Calabrian crust producing uplift and (ii) particular thermal conditions in the crust suggesting a poor coupling between the Ionian subducting plate and the upper plate of the Calabrian arc. From this, some speculation may be carried out on the maximum compressive earthquake (Mw  [29]. Inset A shows a no-scaled cartoon where homogenous heave rates (HR) from computed extension rates are assigned among the active faults such as the West Crati Fault, East Crati Fault, and the Lakes Fault. Inset B shows a scheme of a normal fault with the three components where the deformation is measured and the fault dip angles: heave, throw, and displacement. Note that for a 45 • fault dip, angle, throw, and heave will be equal. Furthermore, we have also iterated geodetically derived fault throw rates for the East Crati Fault and Lakes Fault, hypothetically sharing 66% of our total refined computed regional extension (~0.8 × 3 faults~2.3 mm/yr), iterating their fault dip angles (30 • to 60 • ) (Spreadsheet S2). For instance, given a homogeneously partitioned geodetic heave rate of~0.8 mm/yr for each fault, we show that for fault dip angles of 60 • and 58 • , respectively, for the East Crati Fault and Lakes Fault, values in agreement with previous field studies [35,122,123], we obtain geodetic throw rates of~1.3 for the East Crati Fault and 1.2 mm/yr for the Lakes Fault. These latter geodetic-derived values are in agreement with geological long-term throw rates obtained by previous studies for the same active faults [37,117]. For this reason, we stress that our geodetically derived fault throw rate of 0.8 mm/yr for the West Crati Fault (with an averaged dip angle of 45 • from literature data) may be the most likely scenario for seismic hazard modelling investigations, even for a longer timescale. Furthermore, it is important to stress that debated geological fault throw rates had been proposed within a poorly constrained range between 0.4 and 1.4 mm/yr [28,[35][36][37], therefore misleading estimations of T mean ; our proposed geodetically derived throw rate of 0.8 mm/yr agrees with one of these previously proposed geological throw rates [37]. This would also suggest a refined value for the T mean if well-known empirical correlations between maximum expected magnitude, fault length, and maximum expected displacement are applied [68,69]. Indeed, we are able to estimate a T mean for the West Crati Fault, considering the length of the fault (~40 km) capable of producing earthquakes with a maximum magnitude (M) of 6.8. Given a maximum co-seismic slip event of~1500 mm in M 6.8, in the centre of the fault and assuming a rupture of the entire fault length, to produce a fault throw rate of~0.8 mm/yr would imply a T mean~1 875 yr; this value is not unusual if compared to those describing active normal faults along the Italian territory [9,116].

Regional Implications for Upper-Plate Crustal Deformation above Subduction Zones
Our results from marine terraces modelling to derive refined long-term fault-controlled uplift rates show that the "regional" uplift from either ongoing subduction or deep-seated (mantle upwelling) processes has been spatially affected by "local" faulting activity; this has been observed throughout the Calabrian-Peloritani Arc [9,63] and should be considered when attempts for seismic hazard assessment for subduction-related megathrust earthquakes are studied. This is crucial if we consider that some studies show that the Ionian subduction may be active [124,125] and then capable of producing megathrust earthquakes, even though it is important to highlight there is an ongoing debate regarding if this subduction zone is either seismically coupled or aseismic. For instance, a previous study has shown an apparent "decoupling" about the crustal deformation occurring on the Calabrian upper plate and deformation on this subduction zone; this is also suggested by the fact that there is no record of strong subduction earthquakes, with a relatively aseismic "zone" at least in the first 30-60 km [112,124,126]. This would agree with the idea proposed by [39,127] of having (i) rising asthenosphere just beneath the Calabrian crust producing uplift and (ii) particular thermal conditions in the crust suggesting a poor coupling between the Ionian subducting plate and the upper plate of the Calabrian arc. From this, some speculation may be carried out on the maximum compressive earthquake (Mw ca. 6.0), where strong shallow crustal earthquakes are limited to the upper plate, while aseismic creep dominates along the subduction interface [112]. Yet, if we assume the Ionian subduction zone as an active one potentially producing megathrust earthquakes as suggested previously, permanent and long-term uplift may be ascribed to several seismic and aseismic processes [128], one of which is the slip along the subduction interface with a locked (or coupled) section [128,129]. If we then assume coupling between the subducting slab and the uplifting and extending upper plate within the Ionian Subduction Zone, long-term uplift may be explained as the result of the interseismic continuous elastic deformation due to locking on the plate interface, as proposed for other regions in the Mediterranean realm [128]. However, more investigations are needed to (i) refine the convergence rate, (ii) quantify the "aseismic" amount of subduction, and (iii) take into account deep-seated processes such as mantle upwelling beneath the Calabrian crust [39] if the interseismic continuous elastic deformation is aimed to be studied for this region. Indeed, it is worth noting that the Ionian subduction associated with the process of slab rollback [130] and, in places, the possibility of mantle upwelling beneath the Calabrian continental crust [127,131] identifies its onset synchronously to the uplift of the investigated forearc since the Middle Pleistocene. However, even though geodetic measurements show a spatially homogeneous uplift for the Calabrian-Peloritani forearc [132] with local geodetic-derived subsidence evidence in places likely controlled by interseismic slip on faults on their shear zone [133], as shown in this study (Figure 7), geological and long-term values of "regional" uplift rates still remain unclear and likely spatially changing throughout the forearc. Indeed, different values of geological uplift rates have been mapped beyond fault tips of upperplate active faults where the "local" signal from faults should be zero, such as in NE Sicily (0.89 mm/yr) [63], close to Vibo Valentia in the Capo Vaticano Peninsula (1.75 mm/yr) [9] and this study (0.4-0.52 mm/yr). The cause of this variation should be tested, with more detailed mapping of marine terraces beyond and over fault tips of upper-plate active faults, and future investigations should address whether this variation could be either associated with the locations throughout the forearc and their different distances from the subduction zone or locally differential mantle upwelling process.
Moreover, we stress that the forearc shows prominent intracrustal vertical deformation as also shown in this study from long-term uplift measurements. Quaternary deformation has mostly acted in the form of active normal faulting controlling the topography of this region (Figure 1), with topographic structural highs identified by footwalls of active normal faults [9,63,82,134,135]. Indeed, it has been demonstrated that upper-plate seismogenic sources have been affecting the "regional" uplift signal, implying that these normal faults should be better identified and mapped. For this reason, if the "local" signal from active faults is not considered, this may lead to mistaken subduction interface slip distributions and consequently to a misleading seismic hazard assessment for subduction-related megathrust earthquakes. This would improve the seismic hazard assessment associated with both (i) upper-plate seismic deformation and (ii) subduction-related seismic events.

Co-Seismic Deformation Scenario Associated with the West Crati Fault
We also simulate co-seismic vertical displacement assuming the simplest scenario that the entire~40 km West Crati Fault would rupture during a seismic event ( Figure 10). We performed this by inputting the curvilinear trace of the mapped West Crati Fault into the Coulomb 3.4 software [105], using a relatively new MATLAB code [106]. This code permits us to model vertical displacements arising from earthquakes on faults with variable-strike geometry. Indeed, it is well established that fault geometry influences Coulomb stress transfer [105,106,109] and therefore will influence strain and displacements surrounding the causative fault. We assume an averaged dip angle of 45 • , e.g., [35,36], and an assumed maximum slip at a depth of 3 m at around 9-10 km ( Figure 10, Inset a), for a total coseismic rupture of~40 km long. These parameters are typical for normal faulting-related earthquakes in the Apennines and the Calabrian Forearc [107,[136][137][138] to simulate reliable co-seismic vertical displacements.
We show that the co-seismic footwall uplift region (light-yellow-coloured region defined by isolines of co-seismic uplift) coincides with the long-term footwall uplift deformation mapped by looking at raised marine terraces ( Figure 10 and Figure S11 for the C-D profile trace and signal variations). This is another piece of evidence suggesting that marine terraces have been likely deformed by successive seismic events over the Late Quaternary. Moreover, the maximum co-seismic subsidence coincides with the estimated location of the 1184 seismic event, as well as its modelled magnitude of 6.75. Note that we estimate a maximum co-seismic uplift and subsidence, in the centre of the West Crati Fault, of 0.064 m and −0.810 m (Figure 10, Inset B), respectively, which are not unusual values if compared with other normal faulting-related earthquakes such as the 2009 L'Aquila earthquake (M 6.3) [139] and the 1983 Borah Peak earthquake (M 7.3) [140]. It is important to stress that the maximum expected displacement of 1500 mm (1.5 m) used to estimate T mean in Section 4.1, if a~40 km long fault seismically ruptures, is not the vertical component (throw); it is yet estimated along the displacement D (inset B from Figure 9 for clarification) as suggested by others [68].
We interestingly report that the maximum co-seismic subsidence (dark blue region, Figure 10) coincides with the area where we show the maximum rate of GNSS-derived subsidence within the hanging wall of the fault (Figure 7b).
For this reason, we finally speculate that future works can be made based on an old but simple concept that seismicity rates on brittle faults are controlled by strain rates of viscous deformation at depth on faults' shear zones [141,142]. Indeed, more investigations are needed to test if GNSS measurements on the hanging wall of the West Crati Fault are showing an active shear zone, slipping at its long-term interseismic rate, causing subsidence at the surface, after having ruled out local non-tectonic signals from GNSS measurements; indeed, slip in the upper crust is driven by slip on viscous shear zones in the lower crust, with the viscous flow process on shear zones, as shown by previous investigations [133,141,143].
the 1184 seismic event, as well as its modelled magnitude of 6.75. Note that we estimate a maximum co-seismic uplift and subsidence, in the centre of the West Crati Fault, of 0.064 m and −0.810 m (Figure 10, Inset B), respectively, which are not unusual values if compared with other normal faulting-related earthquakes such as the 2009 L'Aquila earthquake (M 6.3) [139] and the 1983 Borah Peak earthquake (M 7.3) [140]. It is important to stress that the maximum expected displacement of 1500 mm (1.5 m) used to estimate Tmean in Section 4.1, if a ~40 km long fault seismically ruptures, is not the vertical component (throw); it is yet estimated along the displacement D (inset B from Figure 9 for clarification) as suggested by others [68].
We interestingly report that the maximum co-seismic subsidence (dark blue region, Figure 10) coincides with the area where we show the maximum rate of GNSS-derived subsidence within the hanging wall of the fault (Figure 7b). with the associated slip distribution in depth, with a maximum value of 3 m. Co-seismic uplift/subsidence contours produced by our preferred modelled fault in the half-elastic space are shown along the profile A-B, with a maximum co-seismic uplift and subsidence of 0.064 m and 0.810 m, respectively, Inset B. Furthermore, a crustal deformation profile (C-D) is produced comparing the co-seismic and long-term footwall uplift along the strike of the fault ( Figure S11). Note that the estimated magnitude for the 1184 "Valle del Crati" earthquake is 6.75 and matches our modelled magnitude (M 6.8); moreover, the location of this historical seismic event lies where we estimate the maximum co-seismic subsidence.
For this reason, we finally speculate that future works can be made based on an old but simple concept that seismicity rates on brittle faults are controlled by strain rates of viscous deformation at depth on faults' shear zones [141,142]. Indeed, more investigations are needed to test if GNSS measurements on the hanging wall of the West Crati Fault are showing an active shear zone, slipping at its long-term interseismic rate, causing subsid- Figure 10. Modelled co-seismic displacement showing co-seismic uplift and subsidence if the entire fault length of the West Crati Fault (~40 km) would seismically rupture, producing a modelled earthquake of M 6.8. Inset A shows a 3D view of the modelled seismogenic source (the West Crati Fault) with the associated slip distribution in depth, with a maximum value of 3 m. Co-seismic uplift/subsidence contours produced by our preferred modelled fault in the half-elastic space are shown along the profile A-B, with a maximum co-seismic uplift and subsidence of 0.064 m and 0.810 m, respectively, Inset B. Furthermore, a crustal deformation profile (C-D) is produced comparing the co-seismic and long-term footwall uplift along the strike of the fault ( Figure S11). Note that the estimated magnitude for the 1184 "Valle del Crati" earthquake is 6.75 and matches our modelled magnitude (M 6.8); moreover, the location of this historical seismic event lies where we estimate the maximum co-seismic subsidence [29].

Conclusions
In this paper, we use (i) refined GNSS measurements and (ii) the synchronous correlation approach for raised marine terraces modelling to refine rates of crustal deformation at different timescales. This has permitted us calculations of improved rates of uplift through time mapped on the footwall of the West Crati Fault, spanning the Late Quaternary. It has been shown that (i) elevations of marine terraces vary along the strike of the investigated fault, (ii) higher uplift rates are mapped in the centre of the fault, suggesting that the West Crati Fault needs to be considered active, and (iii) older and higher marine terraces are more tilted because of experiencing a longer history of faulting activity.
We show a constant fault-controlled uplift rate over the Late Quaternary, implying that the throw rate on this fault is constant over time, a key observation for the longterm seismic hazard assessment. We also present new refined strain and extension rates throughout northern Calabria of 21 nstrain/yr and~2.3 mm/yr, respectively, across the strike, likely accommodated by the West and East Crati normal faults and the Lakes Fault, suggesting that geodetically derived throw rates may be more the reasonable values to consider for longer timescales. In particular, we qualitatively estimate the throw rate, obtained with a simplified trigonometric approach and assuming an averaged 45 • fault dip angle. We stress that a future development will concern the reanalysis of the data (with further lengthened time series) in a 3D approach and using numerical modelling methods proposed by some in the past [24,144] in order to estimate more quantitively geodetic slip rates across the strike of the Crati valley and northern Calabria. Indeed, we qualitatively estimate a fault throw rate and T mean of~0.8 mm/yr and~1875 years for the West Crati Fault, if an averaged fault dip angle of 45 • is assumed given the computed extension rate from GNSS measurements. Constant rates over the Quaternary of crustal deformation shown in this study are not surprising if we consider that the quasi-E-W oriented ongoing crustal extension is accommodated by a few active faults across the strike. In other words, if more active faults across the strike are accommodating "distributed" strain, this will gradually change their slip rates over time to maintain the regional deformation rate as described by some in the past [32]. Indeed, the Calabrian-Peloritani Forearc is affected by crustal extension accommodated by a more "localized" strain on a few active faults across the strike.
Finally, we stress that it is crucial to consider upper-plate crustal deformation when slip distributions in the subduction interface are estimated using data of uplift to avoid misleading and erroneous subduction-related earthquakes.  Data Availability Statement: Publicly available datasets were analysed in this study. These data can be found here: http://geodesy.unr.edu/ (Nevada Geodetic Laboratory). The results of data that support the findings of this study are available from the corresponding author, M.M., upon reasonable request.