Out of phase Quaternary uplift-rate changes reveal normal fault interaction, implied by deformed marine palaeoshorelines

We have mapped and constrained the timing of tectonically deformed uplifted Late Quaternary palaeoshorelines in the Messina Strait, southern Italy, an area above a subduction zone containing active normal faults. The palaeoshorelines are preserved from up to thirteen Late Quaternary sea-level highstands, providing a record of the deformation over this timescale (~500 ka) for the Messina-Taormina Fault, the Reggio Calabria Fault and the Armo Fault. The palaeoshorelines reveal spatial patterns of uplift through time along the strike of these normal faults, and, given the across strike arrangement of the faults, also reveal how the contribution of each fault to the regional strain-rate progressed through time. The results reveal that the uplift rates mapped within the fault hangingwalls and footwalls were not constant through time, with a marked change in the location of strain accumulation at ~50 ka. The uplift rates, once converted into throw-rates, imply that the three faults comprised similar throw-rates prior to ~50 ka (in the range 0.77 – 0.96 mm/yr), with the Armo and Reggio Calabria faults then switching to lower rates (0.32 mm/yr and 0.33 mm/yr respectively), whilst the Messina-Taormina Fault accelerated to 2.34 mm/yr. The regional extension rate, gained by summing the implied heave rates across the three faults, was maintained through time despite this re-organisation of local strain accumulation at ~50 ka. We explain these out-of-phase fault throw-rate changes during the constant-rate regional extension conditions as due to interactions between these upper plate normal faults. We finally discuss how fault throw-rates changing through time may affect a long-term seismic hazard assessment within active normal fault systems.


Introduction
Knowledge of tectonic uplift rates that extend over multiple seismic cycles has been derived by investigating preserved palaeoshorelines from seismically deforming upper plates above subduction zones such as Western Burma, South America, Japan, southern Italy and Greece (e.g. González-Alfaro et al., 2018;Ott et al., 2019;Roberts et al., 2009;Roberts et al., 2013;Robertson et al., 2019;Saillard et al., 2009;Shikakura, 2014;Shyu et al., 2018). Such studies can provide continuous along-strike profiles of deformation rate of faults responsible for the uplift over different timescales, and across multiple faults. However, such information has rarely been used to constrain how neighbouring faults share the task of accommodating the regional strain-rate over tens to hundreds of thousands of years, providing information on the dynamics of crustal deformation. For northern and central Calabria, Italy, an area located in the upper plate of the Ionian Subduction Zone, observations of knickpoints in rivers crossing faults have been used to suggest that some faults change their slip-rates over the last 300 ka (Quye-sawyer et al., 2021;Roda-Boluda and Whittaker, 2017). This information provides important new insights that hint at the larger scale dynamic processes that occur in such extensional settings above subduction zones. However, these studies did not concentrate on whether fault slip-rate changes were associated with changes in regional extensional rates, and the timing of changes was relatively imprecise as results came from an inversion of river profile data and not absolute dating of deformed horizons. Hence, it is unclear how the slip-rate changes over the Middle to Late Pleistocene can be interpreted in terms of regional dynamics associated with the subduction processes. Such information would be useful, because elsewhere interaction between multiple acrossstrike active normal faults has been shown to be an important process during intraplate crustal extension located within zones of convergence between colliding tectonic plates, and other tectonic settings; the details of the interaction, revealed by data on how the faults experience changes in slip-rate through time, are essential to understand the dynamics governing the overall deformation (Anders and Schlische, 1989;Contreras et al., 2000;Cowie and Roberts, 2001;Cowie and Shipton, 1998;Gupta et al., 1998;Mcleod et al., 2000;Nicol et al., 2006;Roberts et al., 2002). The Messina Strait is located on the extending upper plate of the Ionian Subduction Zone (ISZ; southern Italy) (Malinverno and Ryan, 1986), within an area of active convergence between the African and Eurasian plates, and is considered one of the most tectonically-active regions of the entire Mediterranean (Monaco and Tortorici, 2000;Serpelloni et al., 2010) (Fig. 1). The geodetically constrained crustal extension rate of ~3 mm/yr is accommodated within the Messina Strait by synthetic and antithetic normal faults (Doglioni et al., 2012;Monaco and Tortorici, 2000;Serpelloni et al., 2010) (Fig. 1), which hosted the most powerful recorded earthquakes during the 20th and the 21st centuries in Europe to date (1908 Messina Earthquake -M 7.1; e.g. Aloisi et al., 2013;Meschis et al., 2019;Valensise and Pantosti, 1992). However, Late Quaternary crustal deformation rates associated with these seismically active normal faults are poorly constrained and so the processes that control the dynamics of the deformation are poorly understood. Uplifted Late Quaternary palaeoshorelines suggest Middle to Late Pleistocene deformation Stewart et al., 1997;Valensise and Pantosti, 1992), but rates of activity through time, and the likely effects of interaction on the closely-spaced normal faults are still unclear. Indeed, some have even questioned if these faults are still active (Argnani et al., 2009). Previous studies have shown deformed uplifted Late Quaternary palaeoshorelines mapped along the strike of the offshore Messina-Taormina Fault on its onshore footwall, between Messina town and Taormina town Pavano et al., 2016;Stewart et al., 1997). Likewise, major capable onshore normal faults, such as the Reggio Calabria Fault and the Armo Fault, and their influence on the "regional" uplift throughout the Quaternary on the Calabrian side, have been previously investigated by studying sequences of uplifted Late Quaternary palaeoshorelines both on their footwalls and hangingwalls Miyauchi et al., 1994;Monaco and Tortorici, 2000;Valensise and Pantosti, 1992). Yet, for one of the most seismically active regions in the Mediterranean realm like the fault-bounded Messina Strait, uncertainties remain on (i) how slip-rates on these crustal-scale faults compare between neighbouring structures, (ii) whether the faults interact and hence (iii) what controls the changes in the deformation through time in response to the tectonic forcing.
In this paper we use flights of marine terraces to reconstruct the recent tectonic and geomorphological evolution of southernmost Calabria and northern Sicily. Uplift rates of the footwall and hanging wall blocks of normal faults, several tens of kilometres in length, in the two sides of the Messina Strait (Fig. 1), and the fault throw rates are refined using (i) a synchronous correlation approach driven by age controls from the literature and (ii) new absolute dated age controls. From the modelling, an uplift rate change at 50 ka is identified and we demonstrate synchronous, out-of-phase changes in uplift rate, showing how these relate to slip-rate changes and regional extension rate and hence the regional processes of the subduction system. We finally discuss implications for seismic hazard.

The geological background of the Messina Strait region
The investigated area lies within the uplifting and seismically extending Calabrian-Peloritani forearc, where the origin of the "regional" signal of uplift, associated with either with the Ionian subduction process (Malinverno and Ryan, 1986) or with mantle upwelling Nur, 1999a, 1999b), is still debated, In particular, the Messina Strait, separating Calabria from Sicily, southern Italy, is a downfaulted graben between inward-dipping Quaternary normal faults, deforming a pre-existing fold and thrust belt made of Palaeozoic, Mesozoic and Cenozoic rock formations interposed during Alpine thrusting (Ghisetti, 1984;Malinverno and Ryan, 1986;Fig. 1). This strait is bounded by active normal faults onshore Monaco and Tortorici, 2000;Stewart et al., 1997) and offshore Doglioni et al., 2012;Meschis et al., 2019). In particular, the offshore and E-dipping Messina-Taormina Fault offsets pre-Quaternary basement and has propagated upwards to produce a fault-related anticline on its footwall, from Messina town to Taormina town Meschis et al., 2019;Pavano et al., 2016;Stewart et al., 1997;Fig. 1). This fault has deformed the Late Quaternary-Holocene palaeoshorelines outcropping onshore on the Sicilian side of the Messina Strait De Guidi et al., 2003;Pavano et al., 2016;Stewart et al., 1997). Likewise, W-dipping antithetic normal faults, such as the Armo Fault and the Reggio Calabria Fault, control the coastline topography on the Calabrian side of the Messina Strait, deforming the Palaeozoic crystalline thrust sheets and its clastic Cenozoic cover, with the Middle-Late Pleistocene shallow marine deposits, progressively onlapping on fault planes and suggesting persistent faulting (Ghisetti, 1984;Monaco and Tortorici, 2000).
Existing ages of marine terraces mapped in the Messina Strait are briefly summarized here and in Table 1. In particular, along the mapped palaeoshorelines on the footwall of the Messina-Taormina Fault, shallow marine shells have been dated close to Taormina town, identifying the 125 ka (MIS 5e) palaeoshoreline (Antonioli et al., 2006b). By using the Electron Spin Resonance technique on marine mollusc shells such as Patella and Venerupid, shallow marine deposits have been dated (Table 1). Moreover, eolian non-fossiliferous sediments have been dated north of Reggio Calabria (Fig. 1) by applying a thermoluminescence dating technique, identifying the 50 ka palaeoshoreline (Balescu et al., 1997). Thermoluminescence has also been used to date shallow marine deposits from marine terraces, south of Reggio Calabria, identifying the 125 ka (MIS 5e) marine terrace (Balescu et al., 1997). South of Reggio Calabria and within the hangingwall of the Armo Fault, Persististrombus latus a fossil that is commonly used to date marine oxygen isotope stage 5 deposits in the coasts of Italy, have been collected in situ from within shallow sea deposits and used to identify the 125 ka palaeoshoreline (Aloisi et al., 2013;Antonioli et al., 2006a;Antonioli et al., 2006b;Dumas et al., 2005;Ferranti et al., 2006;Miyauchi et al., 1994).  (Tables 1 and 2). Inner edges are shown with different coloured dots indicating refined ages (Table 3). Inset A, modified from Ghisetti (1984), shows a sketched structural scheme of the Messina Strait, accommodating ongoing crustal extension of ~3 mm/yr (Serpelloni et al., 2010). Inset B shows Quaternary active normal faults within the Calabrian Arc. Zoom in on these locations (dotted black squares) are shown in supplementary materials (Fig. SM 1,SM 2 and SM 3).

Methods
In this section, we present approaches and methodologies used in this work to refine fault-related deformation rates within the crustal extending Messina Strait.

Topographic analysis: DEMs and fieldwork to map palaeoshorelines
Detailed topographic and geomorphological analysis has been carried out using 10 m high resolution DEMs (Tarquini et al., 2012) to create serial topographic profiles along the strikes of the Messina-Taormina Fault, the Armo Fault and the Reggio Calabria Fault ( Fig. 1 and Supplmementary Materials (SM) Figs. SM 1-3). These topographic profiles were chosen where the geomorphology of palaeoshorelines is appropriate due to their preservation from erosion (e.g. Fig. 2) as well as being close to previously mapped locations, that have been used as a guide (Aloisi et al., 2013;Antonioli et al., 2006a;Antonioli et al., 2006b;Balescu et al., 1997;De Guidi et al., 2003). It is important to note that in this study palaeoshorelines are identified as the inner edge of a marine terrace where palaeo mean sea level would stand during a sea level highstand, ruling out distal deposits which may be formed at tens meters of water depth from the inner edge. The topographic profiles cover the lengths of the faults, from their centres towards their lateral tips, but complete coverage along strike has not been achieved as the Reggio-Calabria Fault goes offshore in the SW and palaeoshorelines are unclear and poorly preserved due to erosion in the high topography at the NE-end of the Armo Fault. For the Messina-Taormina Fault along strike coverage has been achieved, and palaeoshoreline elevations have been measured along the entire fault length (Fig. 1). In particular, we have been identifying both erosional and depositional marine terraces and erosional coastal landforms such as truncation surfaces, beach pebbles/cobbles, marine conglomerates, flat-surfaces cut into bedrock by erosion of wave action with evidence of lithophagid borings and millholes (or marine erosion pan) that suggest deposition and re-working in a shoreface environment (Armijo et al., 1996;Ferranti et al., 2006;Miller and Mason, 1994;Robertson et al., 2019;Robertson et al., 2020). They are unconformably overlying the Palaeozoic crystalline thrust sheets and/or Mesozoic/Neogene limestones (Figs. 2 and 3). It is important to highlight that all these features form within the subtidal/intertidal zone a few decimetres to some metres down-dip of the inner edge of the marine terrace. Indeed, gently-sloping surfaces have been interpreted as palaeoshoreface surfaces cut by wave-action and bounded up-dip by fault scarp-like resembling palaeo-sea-cliffs, following previous marine terrace investigations (Armijo et al., 1996;Roberts et al., 2009;Roberts et al., 2013;Robertson et al., 2019Robertson et al., , 2020. Once palaeoshorelines and their associated palaeoshoreface deposits and wave-cut platforms were identified in the field (Fig. 2), their UTM coordinates and elevation were recorded using a hand-held GPS with barometric altimeters with a precision of ±3 m), with the barometer calibrated to sea-level every few hours. Field-based palaeoshoreline elevations were georeferenced onto DEMs to confirm their reliability, following approaches from recent studies Meschis et al., 2020;Roberts et al., 2013;Robertson et al., 2019).

234 U/ 230 Th coral dating
New age controls have been obtained through dating of corals with the 234 U/ 230 Th dating technique. In particular, dated corals allow us to place some constraints on the age for a prominent marine terrace on the footwall of the Reggio Calabria Fault ( Table 2). The corals we collected were from a calcarenite boulder within shallow marine deposits, from a footwall terrace yet mapped in this paper along our Profile 16 and by other authors along the hangingwall cut-off of the Reggio Calabria Fault   (Fig. 1). We then carried out 234 U/ 230 Th dating at the Geochronology and Tracers Facility of the British Geological Survey, Keyworth, UK. Here, we use total dissolution methods as outlined by Crémière et al. (2016), with ratios of the isotopes measured on a Neptune Plus MC-ICP-MS (Multi-Collector Inductively Coupled Plasma Mass Spectrometer). We undertook coral preparation and cleaning, a critical phase before 234 U/ 230 Th analysis, at Birkbeck College. In particular, we isolated millimetre-scale fragments of pure aragonitic corallite, carefully selected and systematically cleaned. We removed as much as possible of any material showing evidence of alteration and/or detrital matrix on the outside of the corallite wall. We did this, following approaches from previous studies (e.g. Meschis et al., 2020; Roberts Table 1 Age controls from literature used to drive our synchronous correlation approach for each investigated fault. Note that elevations palaeoshorelines from this study and previous ones used as age controls may differ because samples may have been collected in a different part of the same terrace (Profiles 12 and 14).  , and refined ages are assigned using a synchronous correlation approach. Field evidence of palaeoshorelines such as lithophagid borings (b) and a Mesozoic limestone-made wave-cut platform overlies by shallow marine conglomerate identifying an unconformity (c) are shown. Location of this profile is also shown in Fig. 1. Robertson et al., 2020), mechanically using a scalpel under a microscope and/or chemically by washing the wall with HCl (10 %) for 5-10 s, followed by thorough rinsing with ultrapure water. Moreover, we separated and removed coral septa from the corallite wall because they are thinner and more prone to diagenetic alteration (Roberts et al., 2009;Robertson et al., 2020). Following other studies (e.g. Meschis et al., 2020), we also examined both bulk corallite wall fragments (2-10 mg) and powdered subsamples obtained using a computer-controlled drill equipped with a 200 μm drill bit in order to avoid any portion of the coral that showed discolouration or other evidence of alteration. Note that samples 1-14 have shown evidence of open system and then discarded from this study yet samples 15-17 allow us to obtain a relative . Inner edge elevations with refined ages are also shown in Table 3.
age control used to refine rates of crustal deformation. Finally, we highlight that our terrace modelling by applying a synchronous correlation approach (described in the next section) is driven by age controls from literature as shown in Table 1.

Synchronous correlation approach for multiple palaeoshorelines and multiple sea-level highstands modelling
This technique has been used within the Mediterranean realm and elsewhere Meschis et al., 2020;Pedoja et al., 2018;Roberts et al., 2009;Roberts et al., 2013;Robertson et al., 2019;De Santis et al., 2021). This technique is built on the concept that sea-level highstands, which are thought to produce palaeoshorelines (Lajoie, 1986), are unequally-spaced in time, implying that for a constant uplift rate through time, palaeoshorelines will be unequally-spaced in elevation (Houghton et al., 2003;Roberts et al., 2009). To use this approach, it is important to obtain at least one age for a palaeoshoreline across the investigated area. The age control drives the iteration of uplift rate to replicate the elevation of the dated palaeoshoreline, but also predicts the elevations of palaeoshorelines from other Quaternary sea-level highstands, for comparison with the elevations of mapped yet un-dated palaeoshorelines. The modelling can involve either a single constant uplift rate or can include multiple changes in uplift rate at times that can be specified by the user. For example, Roberts et al. (2009), studying central Greece, found that a single change in uplift rate through time was necessary to replicate measured palaeoshoreline elevations, whereas Roberts et al. (2013), studying central Calabria (southern Italy), found that no change in uplift rate was needed and a constant uplift rate throughout the Late Quaternary could explain the measured palaeoshorelines. In summary, this approach tests if un-dated palaeoshorelines can be modelled by iterated uplift rates implied by the dated palaeoshorelines. This iteration is enabled by a "Terrace Calculator" built in Excel, available in the literature where data from sea-level highstands are input . The "Terrace Calculator" uses an input uplift rate (UR), which is iterated to calculate 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 this formula: PE = (UR x T) + SLE (e.g. Meschis et al., 2020). More particularly, the uplift for the whole profile is iterated based on the palaeoshoreline with the age control such that the measured elevation matches that predicted by the "Terrace Calculator". Then the measured elevations are matched to the predicted elevations (within a range +/− 10 m) which allows undated palaeoshorelines to be allocated to sea-level highstands.
Note that iterations of uplift rates are driven by reliable and robust age controls from literature shown in Table 1. This approach forces the user to maximise the coefficient of determination (R 2 value) for linear regression analysis through data including the "measured" (or mapped) elevations of all palaeoshoreline elevations on a topographic profile (GIS analysis and field) and the "predicted" (or expected) elevations which identify sea-level highstand elevations implied by iterating uplift rate values driven by age constrains and by means of available sea-level curves for the last ~0.9-1 My (Table SM 1) (Rohling et al., 2014;Siddall et al., 2003).
We discuss how uplift-rates can be used to constrain slip-rates on the faults in the discussion section. However, note that the topographic profiles we study exist in the footwall of the Messina-Taormina Fault and the hangingwalls of the Reggio Calabria and Armo faults; this must be factored into interpretations of the throw-rates and (and slip-rates if the fault dip is known) because an increase in footwall uplift rate implies an increase in slip-rate whereas an increase in hangingwall uplift rate may imply a decrease in slip-rate given a constant background regional uplift rate. Table 3 shows all mapped palaeoshoreline elevations defined in this study with ages constrained by the synchronous correlation carried out in this paper; data cover the regions associated with all 3 faults we have studied across. These elevation data support the plots in Fig. 3, which are shown as examples (and Figs. SM 4-10 for our synchronous correlations applied to all topographic profiles across the 3 investigated faults), with detailed analysis of topographic profiles across the footwall of the Messina-Taormina Fault and the hangingwalls of the Armo and Reggio Calabria faults. The following sections describe the procedures we performed to gather the data in Table 3 and justify their robustness.

Check of the reproducibility of topographic measurements of palaeoshoreline elevations
We conducted a check of the consistency between the two datasets from this work of elevations of palaeoshorelines using (i) a handheld GPS with barometric altimeter in the field (Fig. 2) and (ii) digital mapping from the DEM (Fig. 3). In particular, elevation data from our DEM-based analysis were compared with elevations from our fieldbased palaeoshoreline mapping, also revisiting palaeoshoreline locations from the literature using them as a guide for our fieldwork.

Table 2
Measurements of 234 U/ 230 Th Isotope ratios of coral samples collected on the footwall of the Reggio Calabria Fault (see Fig. 1 Table 3 Mapped inner edge elevations of marine terraces from fieldwork and DEM analysis with refined ages assigned using a synchronous correlation approach are shown. Note that all UTM coordinate are projected in WG84 33 • S grid zone. Note that age controls from literature used for our synchronous correlation modelling are stated by an asterisk (*) and bold-highlighted.
Profile number (sea level highstand referred to Fig. 1  Comparison of the elevations shows a robust correlation with R 2 values >0.99 (Fig. 4). This suggests that: (i) elevations measured elsewhere in DEMs are likely to be reliable and (ii) that our mapping is consistent with mapping data published by other authors in terms of locations and elevations of palaeoshorelines (Antonioli et al., 2006b;De Guidi et al., 2003;Miyauchi et al., 1994), but ages of un-dated marine terraces may need to be reviewed.

Comparison of the elevations of measured palaeoshorelines and the elevations of palaeoshorelines predicted by iterated uplift rate scenarios
We use age controls (Tables 1 and 2) combined with the Terrace Calculator to iterate uplift rates to obtain a correlation between mapped palaeoshoreline elevations and predicted sea level highstand elevations for undated yet mapped palaeoshorelines. This allows us to investigate palaeoshoreline elevations and uplift rates along the strike of the 3 investigated faults. In particular, uplift rates were iterated (see Figs. SM 4-10) and linear regression was used to assess the robustness of correlations between the predicted elevations iteratively calculated given (i) an uplift rate value driven by either new age controls from this study and from literature (Aloisi et al., 2013;Antonioli et al., 2006b;Balescu et al., 1997; Tables 1 and 2) and (ii) fixed values for sea level relative to present-day sea-level for several well-known highstands (Table SM 1; Rohling et al., 2014;Siddall et al., 2003) and those mapped in DEMs . Note that uncertainty associated with sea level curves through the Late Quaternary differ depending on the sea level curve used. For those used in this work the stated uncertainties are 12 m and 6 m (Rohling et al., 2014;Siddall et al., 2003). These values are nearly in the same order of the 10-m DEM resolution and does not affect significantly our refined chronology of marine terraces. Moreover, in a previous study  where similarly a synchronous correlation is applied to study a sequence of uplifted marine terraces deformed by normal faulting activity (with very similar uplift rate values), different sea level curves were tested for sensitivity test to show if differences on final results on derived uplift rates were found. In particular, it is shown that regardless of the sea level curves used (e.g. Siddall et al., 2003;Waelbroeck et al., 2002), very similar results for uplift rates and then tectonic implications were obtained. It is also important to highlight that Siddall et al. (2003) better describes sealevel changes through the Middle Pleistocene to the present for the Mediterranean Sea. Fig. 5 shows correlations between measured and predicted palaeoshorelines elevations using both constant (Fig. 5a, b,) and fluctuating uplift rates through time (Fig. 5d, e, f and Figs. SM 11-18). Initially, we assumed constant uplift rates through time. Correlations between measured and predicted palaeoshoreline elevations produced reasonably good correlations (R 2 values in the range of 0.9677-0.9954; y = 0.9665× to y = 0.903×; Fig. 5a, b and c). However, after a number of iterations, we found that correlations improved if we included changing uplift-rates through time, specifically a change in uplift-rate on all 3 faults at 50 ka (R 2 values >0.9985; y = 0.9697× to y = 1.0167× - Fig. 5d, e and f), with a very robust correlation for all topographic profiles (R 2 > 0.99 - Fig. 5g, h and i). Coefficient of determination, R 2 , for each topographic profile (1− 22) for the fluctuating uplift rate scenario is higher compared to the constant uplift rate scenario (Figs. SM 11-18). We did not enforce this value of 50 ka into the calculations. Instead, this value arose as it produced the best fit for each of the 3 faults, even though they were analysed separately. The fact that the 50 ka value produces the best fits on 3 separate faults in 3 separate calculations suggests that it may well signify a major reorganisation of strain distribution in the rift. The result of our analyses tested this hypothesis by exploring the implications of this possibility. The results of assuming changing rates at 50 ka are detailed in Figs. 5 and 6. Fig. 6a, b and c show spatial variations of uplift along the strikes of the investigated normal faults recorded by the geometry of the palaeoshorelines. The fact that the palaeoshorelines are tilted and/or folded along strike suggests faulting activity spanning the Middle-Late Pleistocene. Furthermore, Fig. 6d, e and f imply that uplift rates vary along the strike of the faults and through time. Note that our interpretation of the uplift rates for the footwall of the Messina-Taormina Fault differs from that of Pavano et al. (2016) who utilise elevation data and uplift rates from ,  and De Guidi et al. (2002). Similar to our interpretation, Pavano et al. (2016) show uplift changing along the coast with minima in the NE and SW near Capo Peloro and Taormina respectively, and a maximum near to Nizza di Sicilia (located in Fig. 1), however they suggest the MIS 5e terrace (~125 ka) is at ~220 m near Nizza di Sicilia whereas we suggest this palaeoshoreline (125 ka -MIS 5e) is at ~170 m elevation (topographic profiles 8 and 9 show the palaeoshoreline elevation between 160 m and 173 m), defining an amplitude of uplift-variation along strike that is bigger than in our interpretation. Moreover, we show that in most cases prominent palaeoshorelines on the sea level curve, mapped elsewhere in the Mediterranean realm such as the 125 ka (MIS 5e), 240 (MIS 7e) and the 340 ka (MIS 9e) (e.g. Cerrone et al., 2021b;Meschis et al., 2020;Roberts et al., 2013;Robertson et al., 2019), identify clear geomorphic inner edges of marine terraces on topographic profiles. More subtle sea level highstands such as the 76.5 ka (MIS 5a) as recognized by others in northern Calabria (Cerrone et al., 2021a), the second peak of the MIS 5e (119 ka) already recognized in SE Sicily  and Apulia region (De Santis et al., 2021) within the Italian territory, and the 410 ka (MIS 11c) are also identified in places across all the 3 faults (in total, we identify palaeoshorelines from 50, 76, 100, 125, 200, 240, 310, 340, 410, 478 ka). In our interpretation, values of uplift rates in the footwall of the Messina-Taormina Fault vary from ~0.5 mm/yr close to the NNE tip zone to ~0.8 mm/yr in the centre of the fault pre-50 ka; post-50 ka uplift rates accelerated from ~1.2 mm/yr close to the NNE tip zone to ~1.9 mm/yr in the centre of the fault (Fig. 6d); this implies an uplift rate increase through time by a factor in the range of 2.375-2.400 (Table 4a). For the hangingwall of the Reggio-Calabria Fault, values for uplift rates for pre-50 ka are 0.9 mm/yr beyond the NNE tip zone and 0.55 mm/yr close to the centre of the fault, with accelerated uplift-rate values after 50 ka of 2 mm/yr beyond the tip zone and 1.33 mm/yr close to the centre of the fault (Fig. 6e); this indicates an uplift-rate increase by a factor in the range of 2.222-2.418 (Table 4a). It is important to note that close to the Scilla region, our uplift rates from Profile 14 for the post-50 ka (2 mm/yr) are in agreement with Late Holocene uplift rates (~2 mm/ yr) proposed by others (Antonioli et al., 2006a;Antonioli et al., 2021;Ferranti et al., 2007). For the hangingwall of the Armo Fault, uplift-rates vary between 0.85 mm/yr close to the southern tip and 0.56 mm/yr close to the centre of the fault prior to ~50 ka, with accelerated uplift rate values after 50 ka of 2.2 mm/yr close to the southern tip and 1.45 mm/yr close to the centre of the fault (Fig. 6f); this indicates an upliftrate increase by a factor in the range of 2.588-2.589 (Table 4a). In summary, the range of uplift-rate increase factors from these 6 different locations (3 fault centres and close to 3 fault tips) is small (2.222-2.589), with an average value of 2.43. This small range, and the fact that all 6 values were derived from independent calculations, suggests this is unlikely to be due to chance, but rather reflects a real change in the distribution of strain that affects all 3 faults in a similar manner. We discuss this further below, but here note that change in fault-controlled uplift rates over the Late Quaternary are not unprecedented; for instance, within the Gulf of Corinth in Greece, tectonically deformed palaeoshorelines mapped within the footwall of the South Alkyonides Fault indicate that uplift rates have varied through time, with an upliftrate change factor of ~3.20 at ~175 ka, suggesting an acceleration of slip-rate for this normal fault (Roberts et al., 2009).
We calculate how much the mapped palaeoshorelines have been tilted along-strike due to displacement gradients along the fault. This tests the hypothesis that the faults have been active throughout the Late Quaternary, because if this is the case older palaeoshorelines will have experienced a longer history of deformation and hence should exhibit higher tilt angles (e.g. Armijo et al., 1996;Meschis et al., 2018;Roberts et al., 2013;Robertson et al., 2019). Results confirm that the topographically higher and older palaeoshorelines present higher tilt angle values, demonstrating that they have experienced a longer faulting history, with progressive active faulting throughout the Middle-Late Pleistocene (Fig. 7). It is important to note that similar along-strike variations of Late Holocene uplift, associated rates and tilting over 10 km along the strike of the Messina-Taormina Fault, close to its southern tip, had been reported by previous studies, suggesting faulting activity (Antonioli et  To gain information on fault slip-rates from uplift-rates it is important to note that the acceleration in uplift at ~50 ka, measured on (i) hangingwalls of Reggio Calabria Fault and Armo Fault and (ii) footwall of the Messina-Taormina Fault, must reflect the interplay between faultrelated vertical motions and vertical motions that are more "regional" in nature either related to the underlying subduction system or mantle upwelling processes through slab window (Gvirtzman and Nur, 1999a;Malinverno and Ryan, 1986;Westaway, 1993). It is widely accepted that there is an important regional uplift in this region, but debate surrounds the precise value, and how it can be separated from more subduction/ mantle upwelling-related regional uplift variations produced by normal faulting (e.g. Meschis et al., 2018;Roberts et al., 2013). However, if the regional uplift stayed constant through time, and affects both the hangingwall and footwalls of the active normal faults, then the following is implied: (i) the footwall uplift on the Messina-Taormina Fault has accelerated at ~50 ka, suggesting that the throw-rate on the fault increased at ~50 ka; (ii) the component of hangingwall subsidence produced by slip on the Reggio Calabria Fault and Armo Fault has decelerated at ~50 ka implying that the throw-rates on these faults decreased at ~50 ka. If our assumption is correct, and regional upliftrates have remained constant through time, this out-of-phase behaviour, where the slip-rates on some faults accelerate accompanied by a synchronous deceleration of slip on others, implies that some sort of fault interaction may be operating.

Field observations that help to convert uplift-rates into throw-rates on the faults
We have collected a number of field observations that allow us to convert uplift-rate histories into throw-rate histories on the faults.
The long-term throw-rate on the Reggio Calabria Fault can be constrained because we have observed in the field a mapped faulted-offset of the 125 ka (MIS 5e) marine terrace along the fault-crossing Profile 16. For the hangingwall of the fault, age constraints in the literature (Balescu et al., 1997 - Table 1) driving our synchronous correlation technique suggest that a prominent terrace composed of marine conglomerates dates to the 5e highstand at ~125 ka. These marine  Fig. SM 11-18. The predicted elevations, representing the synchronously calculated sea-level highstand elevations, indicate a fluctuating uplift rate through time, that has been derived by iterating uplift rates to find the best match to the measured and mapped palaeoshorelines. Note that "measured" elevations represent palaeoshoreline elevations mapped in the 10 m high resolution DEMs. Coefficient of determination, R 2 value, has been used between these two datasets to quantify the best fit for all topographic profiles presented in this paper, with a value >0.99.
conglomeratic deposits are associated with the highest terrace surface along the Profile 16 at ~123 m elevation (Fig. 8). On the footwall cut-off a marine conglomerate deposit exists at ~190 m (Fig. 8), which can be mapped up-dip to a prominent terrace surface. Although we have not achieved a direct age determination of the footwall deposits, we have derived ages for material that is contained within the same terrace, and hence pre-dates this conglomerate, sampled at an elevation of 230 m. We achieved 234 U/ 230 Th age determinations for a death assemblage of detrital corals contained within a boulder of cemented shallow marine sands in this conglomerate of 449 ± 17 ka, 384 ± 11 ka and 480 ± 25 ka ( Table 2). These ages suggest a minimum age of formation of this boulder of 384 ± 11.5 ka (implying the 448 ± 17 ka and 480 ± 25 ka corals are detrital ages for corals inherited from a previous highstand at 384 ± 11 ka). However, the boulder itself is detrital, and its age should be older than the terrace deposits within which it is found. Thus, 384 ± 11 ka is the maximum age of the marine terrace deposit containing the boulder, and the actual age must be younger. It is suggested its age is 125 ka because a very lithologically similar layer of marine conglomerate that caps the marine terrace related to the 125 ka highstand (Balescu et al., 1997;, has been mapped on the hangingwall of the same fault with a base at ~123 m. If the terraces are one and the same, this implies a vertical offset across the fault of ~77 m over 125 ka (Fig. 8). Note that this throw is derived by looking at the same faulted marine conglomeratic deposits in the hangingwall and footwall cut-offs. This would suggest a long-term fault throw-rate averaged over 125 kyrs of 0.62 ± 0.04 mm/yr, in agreement with previous studies which have used offsets of Tyrrhenian deposits and the height of triangular facets (0.6 mm/yr - Monaco and Tortorici, 2000). However, the interpreted increase in uplift-rate in the hangingwall at ~50 ka means that the 125 kyrs average throw-rate should not be used over shorter time intervals. If the uplift-rate change factor is applied (2.43) for the ~77 m of mapped vertical offset (Fig. 8), the throw-rate was 0.80 mm/yr from 125 ka to 50 ka and 0.33 mm/yr from 50 ka to the present day, implying that the throw-rate decreased by a factor 0.41 at ~50 ka ( Fig. 8 and Table SM 2

for throw calculations).
For the Armo Fault, its long-term throw rate can be constrained due to the mapped offset of a marine terrace associated to the 478 ka sea level highstand, a refined age derived by applying a synchronous correlation approach along the fault-crossing Profile 20 (Figs. 9 and 1 for profile location). In particular, we show that the topographically highest and oldest terrace mapped on the hanging wall of the Armo Fault at 355 m dates from 478 ka, according to our synchronous correlations. This terrace, as well as all the lower terraces mapped along the Profile 20 ( Fig. 9), is cut into Plio-Pleistocene bioclastic and siliciclastic marine sandy deposits (Chiarella et al., 2021). A flat terraced surface at the footwall cut-off of the Armo Fault is mapped at ~700 m cut into the Palaeozoic rocks. If this footwall terrace is assumed to be coeval and with the same age of the hangingwall cut-off terrace mapped at 355 m with a synchronously-derived age of 478 ka, this implies a vertical displacement of 345 m, suggesting a throw-rate of 0.72 ± 0.012 mm/yr averaged over 478 kyrs in agreement with other geoscientists that have independently studied this fault through geomorphological investigations (e.g. 0.7 mm/yr -Roda-Boluda and Whittaker, 2017). However, the interpreted decrease in throw-rate at ~50 ka from uplift rates suggests that the 478 kyrs average throw-rate should not be used over shorter temporal windows. If the uplift-rate change factor is applied (2.43) then from 478 ka to 50 ka the throw-rate was 0.77 mm/yr and from 50 ka to present day 0.32 mm/yr for 345 m of offset, implying that the throw-rate decreased by a factor 0.41 at ~50 ka in the hangingwall ( Fig. 9 and Table SM 2 for throw calculations).
For the Messina-Taormina Fault, we do not have direct measurements of faulted offsets of Late Quaternary deposits. However, we derive a throw-rate by subtracting the rate of regional uplift and utilising a value for the ratio of footwall uplift and hangingwall subsidence. In particular, we map the highest uplift rates of 1.9 mm/yr (50 ka to present) in centre of the fault along Profiles 8 and 9 (Fig. 6d). Previous studies have suggested that a "regional" component of uplift, with a rate of ~1 mm/yr for this region (Ferranti et al., 2007;Spampinato et al., 2014;Westaway, 1993), and our observations close to the fault tip, where the effect of faulting would be minimised, approach this value, supporting the ~1 mm/yr "regional" uplift rate. For these reasons, we subtract 1 mm/yr from the "total" mapped post-50 ka uplift rate in the centre of the fault of 1.9 mm/yr (1.9-1 mm/yr = 0.9 mm/yr of footwall uplift) and we then apply a uplift/subsidence ratio of 1:1.6, a value proposed for active normal faults within the Calabrian Arc by previous studies (Quye-sawyer et al., 2021). This implies a hangingwall subsidence of 1.44 mm/yr that, if summed with 0.9 mm/yr of footwall uplift, implies a throw-rate of 2.34 mm/yr for the last 50 kyrs (0.9 mm/yr footwall uplift plus 1.44 mm/yr hangingwall subsidence). If the uplift rate change factor is applied (2.43), then a throw-rate of 0.96 mm/yr before 50 ka is suggested.
For Calabria, Quye-Sawyer et al. (2021) and Roda-Boluda and Whittaker (2017) suggest an increase in throw-rate on the Serre, East Crati and Cittanova faults (Fig. 1) between ~100-300 ka, using observations of knickpoints along rivers crossing these faults and incising into marine terraces, attributing these changes to the interaction of these faults during fault growth and linkage. The questions that arise are (a) whether these slip-rate changes are related to a change in rate across the whole subduction system or due to a change in the internal organisation of strain of the upper plate of the subduction system, and (b) how the changes in rate influence seismic hazard. We discuss a comparison between the rates we have measured, and regional extension rates measured with GPS and then discuss seismic hazard (Fig. 10).

The relationship between fluctuating slip-rates on individual faults and the regional extension rate
Constraints on the regional extension rate are available from GPS observations and modelling based on those observations. Serpelloni et al. (2010) observed a maximum extensional strain-rate of ~65 nanostrains/yr and an extension rate of 3 mm/yr across the Messina Strait (Fig. 10). Inversion of these observations is attempted to resolve slip-rate on faults and the bootstrap analysis of model uncertainties. It shows that for a modelled fault that "dips gently SE-ward" located in the Messina Strait with dimensions and location similar to the Messina-Taormina fault considered herein, optimal values of 3.5 + 2.0/− 1.3 and 1.6 + 0.3/− 0.2 mm/yr for the dip-slip and strike-slip components are found, respectively, locked above 7.6 + 4.6/− 2.9 km depth. We show the implied heave-rate through time for these inversion results as an inset in Fig. 10d.
In order to compare our results with those from Serpelloni et al. (2010), we have summed the implied throw-rates on the faults through time and converted these into heave-rates (extension rates) by assuming fault dips ranging from 40 to 70 o ( Fig. 10d and Table SM 3 for calculations). Although the throw-rates on the faults decrease at 50 ka by a factor of 0.41 on the Reggio Calabria Fault and on the Armo Fault, alongside an increase by a factor of 2.43 at 50 ka on Messina-Taormina Fault, the combined effect of these out-of-phase throw-rate histories is to maintain a relatively constant heave-rate (see before and after the red dashed line in Fig. 10d). These results are not inconsistent with those from Serpelloni et al. (2010), although we prefer higher values for fault dip constrained with field measurement from outcropping fault planes for the Armo and Reggio-Calbaria faults, and modelling of coseismic effects from the 1908 earthquake . This suggests that these slip-rate changes are not related to a change in rate across the whole subduction system, but rather due to a change in the internal organisation of strain of the upper plate of the subduction system. It is worth noting that these faults are closely spaced across-strike, with a distance <10-15 km at the base of the seismogenic layer, where large earthquakes likely nucleate. This distance is much less than half a fault length (~29 km for the Messina-Taormina Fault which is 58 km in length), a distance within fault interaction is expected. This is consistent with the conclusions from Quye-Sawyer et al. (2021) and Roda-Boluda and Whittaker (2017) who suggest that slip-rate changes in Calabria in the Late Quaternary are due to interactions between faults. Where one fault increases its slip-rate other faults across strike decrease their sliprates to maintain constant regional strain-rates, and this has been reported in a number of other locations (the North Sea, central Italy, central Greece, Baja California, New Zealand; (Cowie et al., 2005(Cowie et al., , 2017Cowie and Roberts, 2001;Gupta et al., 1998;Nicol et al., 2006;Roberts et al., 2002;Roberts and Michetti, 2004), and in models of activity rates on faults (Cowie et al., 2012;Cowie, 1998;Cowie and Roberts, 2001;Mildon et al., 2019;Sgambato et al., 2020).
Our observations suggest summed heave-rates in the range of 0.92-3.56 mm/yr across faults that are assumed to dip in the range of 40-70 o , but note that in detail it is challenging to reconcile this with the heave rates of 0.75-4.76 mm/yr derived through inversion of GPS data by Serpelloni et al. (2010) because they suggest a 30 + 1.1/− 0.7 • SE-  Serpelloni et al. (2010). Possible explanations of this disparity are unclear, but may include factors such as (a) that Serpelloni et al. (2010) only included a single fault in their inversion, whilst the palaeoshoreline data we present show that at least 3 faults have been active during the Late Quaternary, (b) that rates measured over interseismic timescales where elastic processes dominate (e.g. decadal timescales from GPS), may differ from rates measured over multiple seismic cycles (10 4 -10 5 year timescales from deformed Quaternary palaeoshorelines) and (c) that the nearby Mt. Etna volcano has undergone changes in its chemistry due to re-organisation of the volcanic plumbing system in the Late Quaternary (plume-arc transition recorded by Etnean primitive melt inclusions from 125 kato present; Schiano et al., 2001), or even at ~60 ka (Barreca et al., 2018), and these changes may be associated with changes in the spatial distribution of strain within the overall subduction system. More work is needed to reconcile rates measured over different timescales. Serpelloni et al. (2010) do show that changes in slip-rate of several millimetres per year can be produced by changes in the locking depths of faults in the overall subduction system, and hence how faults interact, and this may be a way forward towards reconciling slip-rates measured over different timescales. However, for now we conclude that the changes in slip-rates we have constrained are broadly consistent with those in Serpelloni et al. (2010), but our results suggest that rates of late Quaternary slip appear to have been out-of-phase on faults located across strike from each other, Fig. 7. Tilt angle values calculated for each mapped palaeoshoreline shown in Fig. 6a-c, showing that older palaeoshorelines have higher tilt angles, suggesting that they have experienced a longer history of differential uplift, and that differential uplift has been ongoing progressively during the Late Quaternary. Note that for the Messina-Taormina Fault, a flex point is not recorded because we were not able to map palaeoshorelines belonging to the 50 ka sea level highstand. Inset A shows a no-scaled cartoon sketch of the tilt angle of a palaeoshoreline along the strike of a fault.
Note that values of tilt angle for each investigated palaeoshoreline have been calculated as a tan − 1 of a gradient "m" of a straight-line equation (y = mx), as proposed by previous studies Meschis et al., 2020;Robertson et al., 2019). and the summed heave rate is relatively constant through time, suggesting that the slip-rate changes are not related to a change in rate across the whole subduction system, but instead due to a change in the internal organisation of strain of the upper plate of the subduction system. Finally, we speculate, that this apparent decoupling between subduction rates and slip rates on faults in the upper plate may have implications in terms of maximum magnitude for shallow crustal compressional earthquakes along the Ionian subduction zone. Based on the Italian seismic catalogue, likely the most complete catalogue worldwide, we have no record of strong subduction earthquakes in the Calabrian arc (e.g. Maesano et al., 2017;Sgroi et al., 2021). Ionian subduction is relatively aseismic at least in the first 30-60 km. This fits very well with the seismotectonic interpretations proposed in the literature (e.g. Nur, 1999a, 1999b), showing that the asthenosphere is shallow and just beneath the Calabrian arc; and describing thermal conditions in the crust so that the coupling between the subducting plate and the Calabrian arc is very poor. We speculate that this puts a limit (Mw ca. 6.0) to the maximum compressive earthquake that can occur in the shallow crust (first 30-60 km) of the Calabrian arc. In this scenario, strong shallow crustal earthquakes are limited to the upper plate, while aseismic creep dominates along the subduction interface. We finally stress that our study suggests that more  Table 2 (Samples 15, 16 and 17) and the significance of the ages described in the text. Note that elevations of the same marine conglomeratic deposits are mapped in the hangingwall cut-off and footwall cut-off, implying a throw of 77 m.
investigations are needed, combining geodetic and geologic measurements, to better refine rates of crustal deformation within upper plates above subduction zones at different timescales.

Fault interaction and its impact on seismic hazard assessment within the Messina Strait
In relation to the seismic hazard affecting the Messina Strait, the newly-refined long-term fault throw-rates in this paper, and the associated T mean values (or also known as Earthquake Recurrence Interval), play a crucial role for seismic hazard calculations. For example, Earthquake Recurrence Intervals (or T mean ) for a given earthquake magnitude are shorter for higher fault slip-rates (e.g. Cowie and Roberts, 2001). We are aware that in some cases rates of deformation are averaged over longer and different timescales (e.g. Martín-Banda et al., 2021;Roberts et al., 2013;Robertson et al., 2019) compared to this study because it is simply difficult to obtain more detailed faulting activity through time. However, herein, using a throw-rate on the Reggio Calabria Fault for instance over 125 kyr would provide a misleadingly short value for T mean . This is because the throw-rate decreased by a factor of 0.41 at 50 ka (throw-rate over 125 ka = 0.62 mm/yr; assuming 1 m slip per event implies a T mean of 1612 years; throw-rate over 50 ka = 0.33 mm/yr; assuming 1 m slip per event implies a T mean of 3030 years). Similarly, for the Armo Fault, whose throw-rate has decreased by a factor of 0.41 at 50 ka (throw-rate of 0.32 mm/yr), this would imply a T mean of 3125 years assuming a 1 m slip per event, rather than a T mean of 1388 years for throw-rate of 0.72 mm/yr over the last 478 kyrs. For the Messina-Taormina Fault, we estimate a T mean of 427 years assuming 1 m slip per event or 855 years assuming 2 mm slip per event for a throw-rate of 2.34 mm/yr over the last 50 ka. These values are summarized in Table 5. Thus, in this paper, we emphasise that it is vital to assess not just the throw-rate, but also the variability in throw-rates through time on a fault so that an appropriate throw-rate can be used in seismic hazard calculations. However, to do this, one must identify the correct ages of geological markers used in defining slip-rates.
Overall, our results produce a pattern of uplift through time that differs from previous authors (e.g. De Guidi et al., 2002;Pavano et al., 2016). In particular, our interpretations suggest a major re-organisation of the activity rates on the three main faults in the region at ~50 ka that has not been noted by previous authors. Changes in activity rate have been reported from further north on the Cittanova, Serre and East Crati faults (Quye-Sawyer et al., 2021;Roda-Boluda and Whittaker, 2017). Finally, it is important to note that recently authors have claimed a possible different seismogenic source from the Messina-Taormina Fault for the 1908 Messina Earthquake (Barreca et al., 2021). In particular, they propose an offshore NNE-oriented fault, named "W-Fault", with its onshore prolongation bending towards southern Calabria, for a total length of ~34 km (Fig. 15 from Barreca et al., 2021), close to our Profiles 14 and 15 on the hangingwall of the Reggio Calabria Fault. However, note that this geometry and orientation would not explain the geomorphology of wellmapped and tectonically deformed marine terraces in this paper and by others on the footwall of the Messina-Taormina Fault. Indeed, where lower uplift rates are mapped close to Messina town in this paper and already proposed in other studies (e.g. Pavano et al., 2016) coincides with the centre of the "W-Fault" proposed by Barreca et al., 2021, and thus where a maximum long-term footwallrelated uplift rate should be expected. Furthermore, on its onshore Fig. 9. Topographic profile and field sketch showing mapped uplifted palaeoshorelines on the hangingwall of the Armo Fault cut into Plio-Pleistocene marine deposits. It is also shown the Armo Fault offsetting the synchronously-derived 478-aged marine terrace, suggesting faulting activity over the Late Quaternary. Bottom photo shows some of the palaeoshorelines mapped on the hangingwall of the Armo Fault along the Profile 20. Profile location is shown in Fig. 1.  Fig. 10. (a) Throw-rates through time are shown for the Messina-Taormina Fault, Reggio Calabria Fault and Armo Fault. Changing throw rates through time have been calculated from faulted coeval marine terraces for the Reggio Calabria Fault and Armo Fault. In (b) and (c) a tectonic sketch with refined throw-rate through time, pre and post 50 ka, with horizonal extension rates from Serpelloni et al. (2010). Note that photos are from Google EarthTM (Maxar Technologies, Landsat). In (d), values of throw-rates-converted heave rates are shown, iterating different fault dip-angles for all 3 investigated faults, compared to those proposed by Serpelloni et al. (2010) using geodetic measurements. In (d), we assume that fault dip angles for all 3 investigated faults are the same because we have no information for those values. The contrast of changing throw-rates of individual faults with constant heave-rates when summed all the faults suggests that faults are interacting to maintain the regional extension rate.
prolongation there is no either co-seismic evidence or mapped fault scarp likely produced by successive earthquakes; instead, the Reggio Calabria Fault is well-known and previously mapped by several studies (Aloisi et al., 2013;Basili et al., 2008;INGV -DISS Working Group, 2018;Monaco and Tortorici, 2000), deforming the marine terraces outcropping along the Calabrian coastline of the Strait.

Conclusion
A new uplift rate scenario has been presented within the tectonically extending Messina Strait, southern Italy, on the upper plate of the Ionian Subduction zone by refining ages of tectonically deformed Late Quaternary marine terraces. Palaeoshorelines investigated on the Sicilian coast, lying on the footwall of the offshore Messina-Taormina Fault, show along-strike deformed geometry, implying faulting activity over the Late Quaternary. Similarly, long-term faulting activity is implied for the Reggio Calabria Fault and the Armo Fault by showing along-strike deformed geometries of palaeoshorelines. Changes of uplift rates along-strike of the investigated faults confirm: (i) tectonic subsidence of the hangingwalls, which is presumably counteracting the "regional" uplift signal possibly associated either with the Ionian subduction process (Malinverno and Ryan, 1986;Roberts et al., 2013) or with mantle upwelling (Gvirtzman and Nur, 1999b), and (ii) the footwall tectonic uplift. The changes in uplift rates are determined to be correlated with some faults speeding up and others slowing down at around 50 ka. We stress that this new scenario of temporally changing fault throw-rates has critical tectonic implications and effect the long-term seismic hazard approach either for the Messina Strait region and worldwide.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.

Table 5
Values of Earthquake Recurrence Interval (T mean ) are summarized showing the variability if pre-50 ka rates and/or long-term rates are used instead of post-50 ka rates. Note that throw-rates for the Messina-Taormina Fault are for the first time estimated in this paper, implying that no long-term rates are available from literature.
Fault name (slip per event in m) T mean using 125 ka long-term throwrate (years) T mean using pre-50 ka throw-rate (years) T mean using post-50 ka throw-rate (years)