Fault Throw and Regional Uplift Histories From Drainage Analysis: Evolution of Southern Italy

Landscapes can record elevation changes caused by multiple tectonic processes. Here, we show how coeval histories of spatially coincident normal faulting and regional uplift can be deconvolved from river networks. We focus on Calabria, a tectonically active region incised by rivers containing knickpoints and knickzones. Marine fauna indicate that Calabria has been uplifted by >1 km since ∼0.8–1.2 Ma, which we used to calibrate parameters in a stream power erosional model. To deconvolve the local and regional uplift contributions to topography, we performed a spatiotemporal inversion of 994 fluvial longitudinal profiles. Uplift rates from fluvial inversion replicate the spatial trend of rates derived from dated Mid‐Late Pleistocene marine terraces, and the magnitude of predicted uplift rates matches the majority of marine terrace uplift rates. We used the predicted uplift history to analyze long‐term fault throw, and combined throw estimates with ratios of footwall uplift to hanging wall subsidence to isolate the nonfault related contribution to uplift. Increases in fault throw rate—which may suggest fault linkage and growth—have been identified on two major faults from fluvial inverse modeling, and total fault throw is consistent with independent estimates. The temporal evolution of nonfault related regional uplift is similar at three locations. Our results may be consistent with toroidal mantle flow generating uplift, perhaps if faulting reduces the strength of the overriding plate. In conclusion, fluvial inverse modeling can be an effective technique to quantify fault array evolution and can deconvolve different sources of uplift that are superimposed in space and time.


Spatial Scales of Uplift and Geomorphic Response
Observational and theoretical studies have demonstrated the influence of tectonic perturbations on the morphology of fluvial networks (e.g., Howard, 1994;Stock & Montgomery, 1999). Rivers are particularly useful for tectonic analysis because, for a particular upstream area, higher uplift rates produce steeper channel slopes (assuming constant sediment cover, precipitation etc.), therefore spatial differences in uplift magnitude may be observed directly from the landscape (e.g., Kirby & Whipple, 2012;Whittaker, 2012). In particular, longitudinal profiles (i.e., plots of channel elevation as a function of downstream distance) usually exhibit a transient response to changes in uplift rate in the form of breaks in slope, known as knickpoints (e.g., Kirby & Whipple, 2012;Whipple & Tucker, 1999). Additionally, river erosion in detachment-limited settings is dominatly an advective process. As the wave of erosion travels upstream through time (assuming erosion rate is linearly proportional to channel slope) the river contains a record of past uplift events (e.g., Loget & Van Den Driessche, 2009;Pritchard et al., 2009;Roberts & White, 2010).
Changes in uplift rate estimated from river profiles have been used to examine causative tectonic processes such as active faulting, fold growth or dynamic topography (e.g., Boulton et al., 2014;Kirby & Whipple, 2001;Roberts & White, 2010;Whittaker & Walker, 2015). Some work has focused on long-wavelength processes by inverting large numbers of river profiles to find continent or island-wide uplift histories (e.g., Czarnota et al., 2014;Fox et al., 2014;Paul et al., 2014;Roberts et al., 2012;Rodríguez Tribaldos et al., 2017), while other studies have investigated smaller scale phenomena (e.g., Goren et al., 2014;Quye-Sawyer et al., 2020). This analysis is the first to quantitatively deconvolve long-wavelength "regional" uplift and short wavelength faulting using river profile inversion.
Geophysical and geomorphological studies suggest that Italy's topography has been generated by active faulting and longer wavelength processes, probably associated with sublithospheric support (e.g., d' Agostino et al., 2001;Faccenna et al., 2014;Faure Walker et al., 2012). However, how different processes have contributed to the formation of topography remains largely unknown, and the rates and magnitudes of each process are poorly quantified throughout the region. The aim of this paper is to investigate faulting and longer wavelength regional uplift in Calabria where geomorphological and archeological observations, and geochronological data, help to constrain landscape evolution over a range of length and time scales Pirrotta et al., 2016;Stanley & Bernasconi, 2012;Westaway, 1993). We use these data alongside 994 river profiles, which cross all major faults in Calabria, and employ a simple stream power relationship to invert their longitudinal profiles for a spatiotemporal uplift history. We show that Calabria's rivers record both regional uplift and changes in fault throw rate.

Geology and Geomorphology of Calabria
The Cretaceous to Eocene collision of the Eurasian and African plates, which resulted in the Alpine and Pyreneean orogenies in Western Europe, caused profound changes to the landscape of the Mediterranean region. The subsequent segmentation of the Alps, accompanied by significant block rotations and magmatism (e.g., Rosenbaum et al., 2002;Savelli, 2002), created positive and negative changes in landscape elevation on geologic and historical time scales (e.g., Antonioli et al., 2009;Braga et al., 2003;Fellin et al., 2005;Ferranti et al., 2008;Scicchitano et al., 2008). However, the extent to which the present-day topography of Southern Italy records crustal stresses, plate flexure, mantle processes, or climate change is poorly understood.
The geology of Calabria reveals the dramatic paleogeographic change of southwest Europe since Late Eocene-Oligocene cessation of Alpine compression. Its basement of granites, gneisses, and schists (Figure 1), which were deformed during the Variscan orogeny, indicate that Calabria was positioned on the Eurasian margin prior to Alpine collision (Rosenbaum et al., 2002;Rossetti et al., 2001Rossetti et al., , 2004. Metamorphosed ophiolites in the Alpine Nappes ( Figure 1) and high pressure-low temperature metamorphism imply that the region was proximal to the subduction front during the closure of Tethys (e.g., Liberi et al., 2006;Pezzino et al., 2008), with localized compression until the Pliocene (Capozzi et al., 2012).
The southern Tyrrhenian Sea has rapidly stretched since the late Miocene separation of Sardinia and Calabria, and ages of dredged oceanic crust reveal episodic oceanic spreading (Rosenbaum & Lister, 2004). Seismic tomography and deep seismicity (35-500 km) indicate that the Tethyan oceanic plate still subducts beneath Calabria (e.g., Chiarabba et al., 2005;Piromallo & Morelli, 2003). An offshore accretionary prism is observed in seismic data from the Ionian Sea (Minelli & Faccenna, 2010). In contrast, active extension is present both onshore Calabria and along its Tyrrhenian coastline, dominatly expressed as a series of NNE-SSW striking normal faults (Figure 2; e.g., Catalano et al., 2008). Numerous historical earthquakes (Figure 3b), many with devastating tsunamis, attest to the recent activity of the majority of these faults (e.g., Catalano et al., 2008;Meschis et al., 2019). This close spatial coupling of compression and extension is also observed further north in the Italian Apennines and is attributed to the roll-back of the cold subducting slab of the Tethyan oceanic plate (Malinverno & Ryan, 1986).
However, despite numerous observations of recent crustal extension, marine terraces, and exposed tidal notches show that much of Calabria has experienced rapid Quaternary uplift (Antonioli et al., 2009). Shear wave anisotropy measurements are consistent with mantle convection around the subducting plate (e.g., Baccheschi et al., 2008;Civello & Margheriti, 2004), which has been recently suggested as the cause of Calabria's long-wavelength uplift Magni et al., 2014), though little work to date has focused on isolating rates of regional uplift from dynamic mantle processes.

Geomorphic Observations of Quaternary Uplift
Early Pleistocene marine terraces reach heights of 1.3 km above sea-level in the Aspromonte region of southern Calabria (Figure 2). These marine terraces, the oldest in the region and the only terraces found in the footwalls of the major faults, are poorly dated with reported ages between 0.58 and 1.8 Ma (e.g., Catalano et al., 2008;Roda-Boluda & Whittaker, 2017;Tortorici et al., 1995). However, a probable age is the Sicilian stage (0.8-1.2 Ma), based primarily on the first appearance of "boreal guests" including Artica islandica and Hyalinea balthica (e.g., Miyauchi et al., 1994). The oldest terraces are easily identified by their well-preserved wave cut platforms flanking the higher relief massifs (Miyauchi et al., 1994;Roda-Boluda & Whittaker, 2017). The widespread nature of these marine terraces demonstrate that the majority of Calabria's topography has probably developed since the Sicilian stage and indicate that much of the region was below sea-level prior to this time. The massifs of weathered Paleozoic crystalline basement (with peaks ≈ 1.8 km above sea-level) are interpreted to have been an archipelago of small islands, surrounded by deltas, and separated by tidal straits, that were exposed above sea-level before the initiation of Pleistocene uplift (e.g., Longhitano, 2011;Longhitano et al., 2012;Rossi et al., 2017;Westaway, 1993).
Last interglacial (MIS 5e) tidal notches and marine terraces are identified along Calabria's coastline due to the presence of Strombus bubonius and other warm-water "Senegalese" fauna, U/Th ages and ammino acid racimization correlation ( Figure 3a and Table 1). Tectonic uplift rates can be calculated from marine terraces using where U is uplift rate, H t is the observed elevation of the marine terrace, S H is the relative sea-level at the time of terrace formation (where positive values of S H denote sea-levels higher than present) and Δt is the time since terrace formation (e.g., Ferranti et al., 2006). The heights of last interglacial terraces are highly variable across Calabria, with the highest uplift rates (>2.5 mm yr −1 for the last 124 ka) observed in footwalls of faults on the Capo Vaticano peninsula (Bianca et al., 2011). Lower uplift rates (0.47-0.89 mm yr −1 for the last 124 ka) exist in the adjacent, and relatively subsiding, hanging wall of the Cittanova fault (Table 1 QUYE-SAWYER ET AL.

10.1029/2020TC006076
3 of 26  Monaco and Tortorici (2000), Catalano et al. (2008), Minelli and Faccenna (2010), and Fiannacca et al. (2015). Topographic contours at 250 m intervals, with spot elevations (gray triangles) of peaks in the Sila and Aspromonte massifs in meters. Active fault traces shown as black lines with ticks on hanging wall. and Figure 3). Terrace heights in the Crotone Basin, >50 km from major active faults, are indicative of consistently low uplift rates (<1 mm yr −1 ;  constraints (see Table 1 and Figure 3: ID 1-5, 38-40), and, as such, different ages have been attributed to the same terrace (e.g., Carobene & Dai Pra, 1990;Westaway, 1993). In general, the superposition of normal faulting and regional uplift complicates terrace correlation across Calabria, and the uneven distribution of uplift constraints can make it difficult to fully quantify fault growth or regional uplift.
Several local studies show that Calabria's rivers have transient longitudinal profiles containing at least one knickpoint (e.g., Pirrotta et al., 2016;Robustelli, 2019;Roda-Boluda & Whittaker, 2017), which also suggests that uplift has varied both spatially and temporally across the region. Catchment averaged erosion rates derived from cosmogenic nuclide concentrations are similarly variable. Erosion rates are generally low at high elevations within the massifs or above fluvial knickpoints (∼0.1 mm yr −1 ), and are higher (up to 1.6 mm yr −1 ) upstream of active faults or in small catchments close to the coast below a major knickpoint (Cyr et al., 2010;Olivetti et al., 2012;Roda-Boluda et al., 2019).

Active Faults in Calabria
Over 130 moderate to large magnitude earthquakes (3.1 ≤ M ≤ 7.1) with well-constrained epicenters have been documented throughout Calabria in the last ca. 1,000 years (Rovida et al., 2016). The wide spatial distribution of their epicenters indicates the presence of many active faults ( Figure 3b). Radiocarbon dating of trenched normal faults and damage to archeological sites provides evidence for Holocene activity of some structures (Cinti et al., 2015;Galli & Bosi, 2002;Galli et al., 2007). Many of Calabria's faults have a clear geomorphic expression that can be mapped from digital elevation models (Figure 2), and the large fault scarps imply that seismicity originated during the Pleistocene (e.g., Catalano et al., 2008;Monaco & Tortorici, 2000). Major faults pertinent to this study will be discussed in detail below.
10.1029/2020TC006076 5 of 26  The NW dipping NE-SW striking Cittanova fault lies entirely onshore, and has the longest fault trace (∼42 km) in Calabria . With fault segments of ∼10 km, the Cittanova fault probably reached its current size though the interaction of a series of en echelon normal faults, whose connecting relay ramps have since been breached (e.g., Fossen & Rotevatn, 2016). This model of fault growth is supported by the presence of knickpoints along tributaries of the Petrace river, which have been interpreted as the geomorphic expression of increases in throw rate (Pirrotta et al., 2016;Roda-Boluda & Whittaker, 2017). Further north, the Serre fault has a similar en echelon morphology and a length of 35 km (e.g., Galli et al., 2008). Along with the Armo fault in the south, they form a linked fault array (Roda-Boluda & Whittaker, 2017), which was probably responsible for the 6.74 ≤ M ≤ 7.1 earthquake sequence in 1783 (Galli & Bosi, 2002).
Published estimates of average throw rate since the onset of faulting for the Cittanova fault lie in the range 0.4 mm yr −1 to 1.4 Whittaker, 2017;Westaway, 1993). Throw rate estimates are similar for the Serre fault, ranging from 0.6-0.7 mm yr −1   Whittaker, 2017). These calculations are based upon an assumed age of the oldest offset marine terrace (Section 1.2.1).
The smaller Scilla, Santa Eufemia, and Reggio Calabria faults lie close to the Messina Strait in the south west of the region, creating a half-graben that is clearly expressed in the topography of the Aspromonte area ( Figure 2b). Synchronous terrace correlation shows that the Vibo fault, on the Tyrrhenian coast of central Calabria, has experienced a throw rate of ∼1 mm yr −1 since 340 ka (Roberts et al., 2013). In the north of Calabria lies the Crati basin, a graben bounded by the West and East Crati faults. Both faults strike approximately N-S and their traces can be mapped at the surface for ∼50 km (Figures 1 and 2). Offset horizons in reflection seismic data indicate an average throw rate for the East Crati fault of ≥0.9 mm yr −1 since 0.7 Ma (Spina et al., 2011). This estimate agrees with an average throw rate of 1.3 0.7 0.5 -mm yr −1 calculated using geomorphic measurements (Roda-Boluda & Whittaker, 2017). Cosmogenic nuclide catchment averaged erosion rates from the footwalls of the Serre-Cittanova-Armo fault array vary along strike, and some erosion rates equal-within error-the throw rates estimated by geomorphic and geologic analyses (Roda-Boluda et al., 2019). On average, however, catchment averaged erosion rates are a factor of two smaller than uplift rates; this discrepancy probably arises because catchments are only partially incised by rivers and may have experienced different amounts of landsliding (Roda-Boluda et al., 2019). Nonetheless, these spatial correlations between erosion rates and throw rates suggest that proxies for surface processes can be used to investigate the timing of active faulting (i.e., we can solve the "inverse problem" of quantifying tectonics from changes in topography).
While the geologic throw and time-averaged displacement rates for the largest faults have been constrained since fault initiation, changes in throw rate have proved more difficult to identify because paleoseismicity can only analyze relatively short time scales compared to geological or geomorphological data (e.g., Galli et al., 2007;Roda-Boluda & Whittaker, 2017). In this paper, we investigate whether fluvial inversion can help to further constrain the temporal history of active faulting in Calabria. In particular, we will focus on the East Crati, Serre, and Cittanova faults. Elevation errors from Ferranti et al. (2006). Uplift rates calculated using Equation 1 assuming MIS 5e occurred during 120-130 ka, with sea-level 3-9 m above that of the present-day . Paleosea-levels and durations of other highstands from Siddall et al. (2007Siddall et al. ( , 2010. Dating method: TL, Thermoluminescence; OSL, Optically stimulated luminescence; SB, Strombus bubonius; SF, Senegalese fauna (not S. bubonius); RC, radiocarbon (calibrated age used); U/Th, U-series or U/Th-series; AM, ammino acid racimization. Assume terrace correlation if no dating method is given. References: (1)

Longitudinal Profile Generation
To extract a fluvial drainage network across Calabria, Esri's steepest descent flow routing algorithms (Flow Direction and Flow Length), were applied to the SRTM 1 arc second (≈30 m spatial resolution) digital elevation model (Stucky de Quay et al., 2017;Tarboton, 1997). An upstream drainage area of 0.32 km 2 is assumed to approximate the threshold for fluvial incision, and cells with this upstream area were systematically sampled to provide the heads of rivers for this study. This technique results in good spatial coverage of the fluvial network and does not bias against rivers of a particular length or stream order assuming that more rivers are extracted from larger catchments. The cumulative number of cells that flow into each catchment (Flow Accumulation) was multiplied by cell resolution (30 × 30 m) to calculate upstream area, A. The morphology of the extracted fluvial drainage network was verified using a combination of aerial photography, published maps (e.g., Pirrotta et al., 2016), and field surveying. The result of longitudinal profile extraction is shown in Figure 4.
Two versions of this river inventory were used for fluvial inverse modeling: The first contains a drainage network across the whole of Calabria, as presented in Figure 4. For the second inventory, we removed all rivers draining the large Crati Basin, where observations of alluviated channels close to the river mouth suggest that a stream power erosion model may be less appropriate. The results of the inverse model from the second river inventory are presented in the supporting information; the differences between the two models are quantified and discussed therein, and in Section 3.1.

Stream Power Erosion Models
Field observations show that many of Calabria's large rivers flow over bedrock with sparse alluvial cover, particularly in the vicinity of the normal faults in the west of the region (e.g., Roda-Boluda & Whittaker, 2017), which suggests that fluvial erosion can be approximated using a detachment-limited model (e.g., Howard, 1994). Erosion rate in a stream power model is parametrized as a function of channel slope, width, and discharge (Howard, 1994). Upstream area, A-measured from digital elevation models-is a useful surrogate for discharge and channel width, which are difficult to quantify over geological time scales. Assuming the rate of elevation change, ∂z/∂t, is the sum of uplift rate, U, and erosion rate, E, a simple version of the stream power model can be expressed as where K is a constant of proportionality often linked to erodibility of the bedrock (e.g., Lague, 2014;Whipple, 2004) and ∂z/∂x is the longitudinal channel slope. Exponents m and n are positive and are usually empirically evaluated.
The exponent, n, determines the dependency of erosion rate on channel gradient and in theory controls the rate of landscape response to perturbation. If n is not equal to 1, the record of tectonic signals can be lost through the formation of shocks and discontinuities (e.g., Harel et al., 2016;Lague, 2014;Pritchard et al., 2009;Royden & Perron, 2013). While theoretical considerations may predict that n > 1-if erosion is controlled by thresholds associated with stochastic weather events for instance (e.g.,  Apennines and Southern Italy have suggested that n ∼ 1 is reasonable in this setting. For instance, the magnitudes and distributions of unit stream power scale predictably with structural and geomorphic measures of footwall uplift and fault throw rate for rivers crossing faults in the Central Apennines of Italy, consistent with n = 1 (Whittaker et al., 2007), while a compilation of catchments crossing active faults in Calabria show that normalized channel steepness index scales linearly with catchment throw rates with no distinct lithological control (cf. Roda-Boluda & Whittaker, 2017; Roda-Boluda et al., 2019) An analysis of knickpoint retreat rates for Italian channels, when hydraulic width scaling is included, also indicates that n = 1 is plausible (Whittaker & Boulton, 2012). Similarly, joint-inversion of drainage networks for uplift rate produced low misfits when n ≈ 1 (e.g., McNab et al., 2018;Paul et al., 2014;Rudge et al., 2015), and a unit stream power model was derived from longitudinal profile morphology of rivers in central Sardinia (Quye-Sawyer et al., 2020). If n ≈ 1, there is a simple, physical relationship between erosion process and channel slope (e.g., Whipple & Tucker, 1999), and the stream power model can therefore be solved using a computationally efficient linearized inversion approach (Glotzbach, 2015;Goren et al., 2014;Rudge et al., 2015). Consequently, we proceed with n = 1, though we acknowledge that the value of this exponent remains contentious, and we therefore return to this assumption in the discussion.
An increase in uplift rate can produce changes in the slope, ∂z/∂x, of longitudinal river profiles known as knickpoints and knickzones. However, an important consideration when interpreting the shape of longitudinal river profiles is the contribution from changes in bedrock competence and discharge. Tensile or compressive rock strength are often used as proxies for bedrock erodibility as a function of lithology (e.g., Sklar & Dietrich, 1998;Roberts & White, 2010;Zondervan et al., 2020). In Calabria, the compressive strength of bedrock along river channels has been recently measured using a Schmidt hammer by Roda-Boluda et al. (2018). These authors found that median Schmidt hammer rebound values were generally low, <35, suggesting that bedrock is weak across a range of lithologies. These observations indicate that lithology probably does not determine the position of Calabria's knickpoints, therefore we may make the simplifying assumption that K is a constant. In addition, if knickpoints are generated by differences in rock strength, we may expect knickpoints to systematically correlate with the position of lithologic transitions (e.g., Wobus et al., 2006). Therefore, we will compare the location of channel slope discontinuities with mapped bedrock geology to evaluate the assumption that changes in lithology do not control the shape of longitudinal profiles.
It is possible for fluvial drainage networks to be modified during glacial periods. A few glacial deposits have been mapped on the highest peaks in the Pollino range, on Mt Sila and in northeastern Calabria (Palmentola et al., 1990). However, since terminal moraines are found >1,400 m above sea-level, and are mainly distributed in an area that lies upstream of the threshold for fluvial incision, we conclude that Pleistocene glaciation had a negligible effect on Calabria's fluvial drainage network (Palmentola et al., 1990).
Mean annual precipitation measured across Calabria indicates that present-day coupling between elevation and precipitation is very weak (D'Arcy & Whittaker, 2014). Moreover, as Calabria has been rapidly uplifted from sea-level during the last ∼1 Myr, it is unlikely that Pleistocene orographic precipitation was more significant than at present (Section 1.2.1). Paleoclimate reconstructions suggest rainfall in Southern Europe did not greatly differ between glacial and interglacial periods (Braconnot et al., 2007). Therefore, climatic changes are unlikely to drive long time period differences in fluvial erosion rate across Calabria, and we will assume that discharge, which controls erosion rate in the stream power model through upstream area, A, does not vary through time to avoid unconstrained model inputs.
The major drainage divide passes through the high relief massifs in central Calabria (Figures 2 and 4), implying that large scale drainage reorganization has not occurred since uplift initiated at ∼1 Ma.
Consequently, we suggest that the majority of observed knickpoints are unlikely to have been generated by drainage divide migration (cf. Willett et al., 2014). Instead, the high number of knickpoints and knickzones across the region, many of which are far from the major drainage divide or upstream of active faults, suggest that fluvial channels are responding to rock uplift at a variety of spatial and temporal scales.

Fluvial Inverse Modeling
We used the joint spatial and temporal fluvial inversion model of Rudge et al. (2015) to predict the cumulative uplift of Calabria since the exposure of the oldest marine terrace at 0.8-1.2 Ma (Section 1.2.1). The advantages of using this type of inverse model include the ability to simultaneously analyze large numbers of river profiles and to calculate uplift rates without the need to pick or classify knickpoints. Moreover, the details of fault location, activity, and linkage history do not need to be established in advance.
The inverse model solves for the spatial distribution of uplift rate on a regular triangular grid that was generated from evenly spaced vertices 10 km apart. A 10 km vertex spacing ensures that at least part of a river exists within the vast majority of grid cells, so the inverse model should be able to resolve recent uplift rates for most of Calabria. This vertex spacing is generally less than the fault separation (Figure 3), and is much smaller than the area believed to be influenced by regional uplift (Section 1.2), therefore our inverse model may capture uplift caused by both the faulting and regional uplift processes that are known to modify Calabria's landscape. By modeling uplift on an arbitrary grid (i.e., without specifying fault position a priori) we can investigate if the inverse model replicates expected geologic behavior such as divergence in uplift rate between the footwall and hanging wall of a mapped fault. Spatial variation in uplift rate was linearly interpolated between vertices using barycentric coordinates. Predicted uplift rate was permitted to vary at 30 evenly distributed time steps.
From Equation 2, the time, τ, for a knickpoint to travel between longitudinal distances x 0 and x 1 can be written as where K is a proxy for bedrock erodibility and A(x) m defines how erosion depends on upstream catchment area, A. Therefore, the predicted elevation, z t , of a river channel as a function of distance, x, can be calculated using where U(x(t), t) is uplift rate as a function of space and time, which is integrated along the time-longitudinal distance path of Equation 3 to derive elevation. For the methods employed in this analysis, Equations 3 and 4 were evaluated using the trapezium rule in order to find the uplift history that produced the observed longitudinal profiles (Rudge et al., 2015).
Equations 3 and 4 show that the stream power incision model (Equation 2) can be linearized such that, z = MU, where elevation and uplift values are denoted by the vectors z and U, respectively. This problem tends to be underetermined (i.e., there are more possible uplift models than can be constrained by fluvial profile observations alone), so the inversion model minimizes where the value of λ s determines damping in space, s. Time at the present-day is denoted by t = 0, and t max is the maximum possible τ for all rivers assuming that a knickpoint can travel from the river mouth to the river head (Rudge et al., 2015). Note that knickpoints can be generated at any position along the river profile using this inverse scheme. Equation 5 was minimized using the nonnegative least-squares Broyden-Fletcher-Goldfarb-Shanno algorithm of Zhu et al. (1997). The initial uplift rate guess for least-squares minimization is U = 0 at all nodes in space and time. A positive uplift rate as a function of space and time was incorporated at subsequent iterations if required to produce a better fit between observed and predicted longitudinal profiles. We assume that Equation 5 is minimized when its value differs by <10 −6 between consecutive iterations. The uplift rate as a function of space and time that minimizes Equation 5 is henceforth known as the best-fitting uplift model.
We followed Parker (1977)'s protocol to seek the smoothest model with the lowest root-mean-squared (rms) misfit, which we will evaluate using independent geologic constraints. In general, inverse models that are highly damped (e.g., λ s ≫ 1) produce smooth uplift with large rms misfit. A very smooth model (large "model norm") might not incorporate short wavelength changes in uplift related to normal faulting and as such would be unsuitable for Calabria. However, models with little damping (e.g., λ s ≪ 1) can over-fit the data and may be fitting noise (e.g., Parker, 1977). We performed a systematic test of model damping, in which λ s was varied between 10 −3 and 10 3 , to find an appropriate value of λ s for this model, and we subsequently evaluate the influence of spatial damping on apparent fault timing (see Section 3.2 and supporting information).
We calculated the root-mean-squared (rms) misfit to evaluate the extent to which river profiles predicted by the best-fitting uplift model correspond to the observed longitudinal profiles. The rms misfit, H, was calculated using where N is the total number of elevation measurements, σ is the error in the observed data, and z o and z t are the elevations of observed and predicted longitudinal river profiles, respectively. The absolute vertical error of SRTM 1 arc second data in high relief regions is ≈16 m (e.g., Mukul et al., 2017), so we set σ to 16 m.
Inverse approaches can systematically test how the exponent of upstream area, m, affects rms misfit and calculated uplift. Most published values of m lie between 0.2 and 1.0, and m = 0.5 is commonly reported for fluvial settings (e.g., Bishop et al., 2005;Howard & Kerby, 1983;Loget & Van Den Driessche, 2009;Quye-Sawyer et al., 2020). Therefore, we repeated the inversion procedure with values of m between 0.1 and 1.0 and assumed that suitable average m values will produce low rms misfits. We also used the inversion modeling to evaluate the average value of bedrock erodibility, K, for Calabria given the time constraints on the age of the upper terrace (∼0.8-1.2 Ma, see Section 1.2.1).
To test the accuracy of our uplift model, we compared uplift rates relative to present-day sea-level calculated from marine terrace heights to those predicted by inverse modeling since interglacials MIS 5 and MIS 7. Although, this comparison was restricted to marine terraces with absolute dating constraints, the analysis still incorporates localities across the region, including the minimum and maximum uplift rates since the last interglacial highstand. For terraces >2 km away from a model vertex, the cumulative uplift map was linearly interpolated so terrace uplift rate and inverse model uplift rate were compared at the same spatial location. As geologic and geomorphic evidence suggests most of Calabria was a submarine environment prior to early Pleistocene time (Section 1.2.1), and to facilitate comparison between the uplift rates predicted by fluvial inversion and marine terrace elevations, we have opted to use sea-level as the most appropriate river base level in this study.

Deconvolution of Normal Faulting and Regional Uplift
Calabria is experiencing simultaneous regional uplift and extensional faulting, which has resulted in some fault hanging walls being uplifted relative to sea-level ( Figure 2). To deconvolve regional uplift and normal faulting, we first extracted cumulative uplift from the best-fitting inverse model at locations in the footwalls and hanging walls of mapped faults to estimate long-term throw rates. We subsequently used ratios of footwall uplift to hanging wall subsidence to estimate regional uplift through time at the same location.
If the oldest terrace (Sicilian stage, 0.8-1.2 Ma) can be correlated across the tip of a fault, being observed in both the footwall and proximal hanging wall, we can calculate regional uplift using where ΔZ h and ΔZ f are changes in the elevation of hanging wall and footwall, respectively. U r , S h , and U f are the rates of regional uplift, hanging wall subsidence and footwall uplift between times t and T (Figure 5a). α is the ratio of hanging wall subsidence to footwall uplift. Substituting Equations 8 and 9 into Equation 7, and rearranging, yields cumulative regional uplift, such that Equation 10 can be applied to the inverse model output at every time step to estimate regional uplift through time (Figures 5b and 5c).

Results and Discussion
The majority of Calabria's rivers contain at least one knickpoint or knickzone between ≈100 and 1,200 m above sea-level (Figures 6 and 7). Although some breaks in channel slope occur on the boundary between different lithologies, such as the knickpoint at 42 m upstream on the Tacina river, knickpoints are not always present on the boundary between these rock types-no knickpoint is present on the boundary between the Oligo-Miocene basins and the igneous basement on the Neto river, for example ( Figure 6). In contrast, most knickpoints reside within a single mapped unit and many are observed upstream of normal faults (Figures 6a, 6b, 6d, and 6f). These observations suggest that Calabria's rivers record uplift that varies in both space and time, and changes in channel slope are not primarily driven by variable bedrock erodibility, in agreement with existing studies at a smaller scale (Glotzbach, 2015;Roda-Boluda & Whittaker, 2017).
For the inverse model, a value of λ s ≈ 1 produces a combination of suitable model roughness (a small "model norm") and low rms misfit for Calabria (Figure 8a), which is similar to the optimal value used in many previous studies (e.g., Conway-Jones et al., 2019; Roberts et al., 2018;Rudge et al., 2015). Therefore, our uplift analysis will initially consider inverse models with λ s = 1. We will subsequently test the influence of the λ s value on the apparent fault timing and throw rates inferred from the fluvial inverse model.
For λ s = 1, the inverse model fits the data poorly if m ≲ 0.3 or m ≳ 0.75, which is consistent with previous inversion studies (Figure 8b; e.g., Roberts et al., 2012). To further constrain the value of m, we compared the elevation of Capo Vaticano's highest terrace, with a mean elevation of 550 m above sea-level, to the cumulative uplift calculated from vertices that intersect the terrace (Figure 8c). The highest terrace of Capo Vaticano is one of the most geographically extensive in the region, extending over several vertices in the model mesh, so represents a good location to test the suitability of the inverse model parameters. We aimed to produce models with similar mean elevation to this terrace that also generated theoretical river profiles with a low rms misfit. These results suggest that m = 0.65 is appropriate for fluvial erosion in Calabria (Figures 8b and 8c).
10.1029/2020TC006076 12 of 26 Figure 5. Graphical representation of procedure to deconvolve regional uplift from results of fluvial inverse modeling.
(a) Fault cross section showing the relative observed uplift, ΔZ f , in fault footwall (circles) and ΔZ h in fault hanging wall (diamonds) for constant regional uplift (blue arrows). U f and S h indicate magnitudes of footwall uplift rate and hanging wall subsidence rate, respectively. Dashed dark blue line represents the magnitude of regional uplift, U r between time t and the present-day t = 0. (b) Uplift as a function time in the footwall and hanging wall (circles and diamonds, respectively) with calculated regional uplift denoted by blue dashes. (c) As (b) but with regional uplift linearly interpolated between time t and the present-day.  Figure 7 shows the best-fitting longitudinal profiles of the eight catchments labeled in Figure 4. The rms misfit, H, is 1.63 for the best-fitting uplift model when m = 0.65 and λ s = 1.0. Although the H value is close to unity, implying that-on average-the inverse model almost replicates the observed longitudinal profiles within error, some rivers have better fits than others. Therefore, we calculated the difference between the observed channel elevation and the channel elevation predicted by the inverse model (the "elevation residual,"  o t i i z z ) as a function of downstream distance for every river. The elevation residuals are normally distributed with a mean of −0.04 m (Figure 9), which suggests that the majority of channel elevations are replicated accurately by the inverse model and elevation is not systematically underpredicted or overpredicted. The standard deviation of the elevation residuals is 26 m, which is the same order of magnitude as the absolute vertical error of the SRTM data set. The largest elevation residuals (∼10 2 m) occur in steep headwaters and across lakes or dams, such as at ≈50 m upstream on the Neto river (Figures 9 and 7).   Figure 4 (plotted as gray lines) and theoretical river profiles (red lines) calculated using uplift history shown in Figure 10.
In general, high residuals are principally a function of model damping, though the accuracy of the SRTM data is also likely to decrease in the steep and narrow topography of Calabria's headwaters (Miliaresis & Paraschou, 2005;Mukul et al., 2017).
Cumulative uplift for the last 20 model time steps is shown in Figure 10a, and uplift rates at each of these time steps are illustrated in Figure 11. We intend to use our fluvial inversion model to analyze the uplift that produced the Pleistocene-Recent marine terraces, therefore the first time step at which the inversion produces uplift is designated an age of 0.8-1.2 Ma (based upon the age of the oldest marine terrace, see Section 1.2). The age range on the maps in Figures 10a and 11 encompasses the uncertainty in the oldest terrace age at all subsequent time steps. An initial uplift time of 0.8-1.2 Ma corresponds to an average bedrock erodibility K = 0.82-1.22 m (1−2m) Myr −1 , noting that an older landscape age would linearly decrease K, and a younger landscape age would linearly increase K because bedrock erodibility is directly proportional to erosion rate according to Equation 2.
Model coverage, whose value depends upon the number of channel measurements between mesh vertices as well as stream power parameters K and m, decreases at earlier model time steps (Figure 10b). This decrease in model coverage occurs because the wave of fluvial erosion continually migrates upstream through time according to the stream power equation, though an uplift event occurring at a place and time when model coverage is >0 should still be recorded on river profiles today (Figure 10b). Conversely, some knickpoints may have reached the river head between the start of uplift and the present-day, so uplift events that produced those knickpoints would not be resolved by the inverse model. Nonetheless, model coverage is >0 over most of Calabria during the last ∼700 ka, which implies that an uplift history can be produced for most of the region at the majority of time steps.
Predicted cumulative uplift from inverse modeling first exceeds 1 km magnitude in the north of Calabria (at ∼300 ka) then in the Aspromonte region ( Figure 10). Uplift of the Serre and Sila Massifs is calculated to occur from 550 ka in the model with initial uplift at 1 Ma, with ≈1 km of uplift prior to 100 ka in the Serre area. A similar pattern of uplift is predicted from the model for the Sila Massif. More than 500 m of cumulative uplift is observed on the Capo Vaticano peninsula by 72-108 ka, and uplift on the east coast of Calabria is typically <500 m throughout the model run. In the hanging wall of the Cittanova fault, and the Crati and Crotone basins, cumulative uplift does not exceed 300 m. Calculated uplift on the footwalls of the Serre and Cittanova faults is initially localized close to the center of modern day fault traces (see 217/325 ka map in Figure 10). Significant cumulative footwall uplift is then observed along a greater extent of the fault array in subsequent time steps.

Evaluation of Fluvial Inversion Results
Given that we have used a simple stream power-based erosion equation to model landscape evolution, are the uplift rates calculated from fluvial inverse modeling comparable to existing uplift rate estimates? The majority of uplift rates calculated from the model replicate, within error, uplift rates derived from Mid-Late Pleistocene terrace heights ( Figure 12 and Table 1). In Figure 12, a range of modeled uplift rates are presented (e.g., 1.6-2.2 mm yr −1 for terrace ID = 10) because these ranges take into account the uncertainty in age of the oldest marine terrace QUYE-SAWYER ET AL.  (i.e., 0.8-1.2 Ma), which was used to calibrate erodibility, K, for the inverse model. The highest inverse model uplift rate of 2.5-3.3 mm yr −1 since MIS 5e coincides with highest observed uplift rate from a terrace on the Capo Vaticano peninsula (Table 1: ID = 14). The smallest uplift rate from the inversion model occurs at the location of one of the lowest last interglacial terraces, near the town of Crotone (Table 1: ID = 34).
The maximum cumulative uplift from the inverse model is 2,077 m (Figure 10: 0 ka panel), situated on a vertex close to the northern drainage divide of the Crati catchment near Monte Pollino (2,248 m). Large magnitudes of uplift (∼1 km) are also predicted at the Sila, Serre, and Aspomonte massifs during the youngest time steps (Figure 10). However, the fluvial inverse model assumes that all topography must be generated between 0.8-1.2 Ma and the present-day, while the massifs probably had preexisting relief of ∼10 2 m in the Sicilian stage in contrast with the majority of Calabria (Section 1.2). This may explain the high modeled uplift rates at the massifs since 100 ka. Uplift at the massifs is unlikely to be added at the start of the model because model coverage is very poor in these locations and at these time steps (Figure 10b).
In addition, we stress that the results presented here are based on the assumption that river erosion in Calabria can be approximated by a detachment-limited stream power model over the last ≈1 Myr. This assumption is probably valid for the majority of Calabria's rivers, especially those in the south of the region that are actively incising across several faults with negligible sedimentation in the uplifting hanging walls (Figure 2). However, some low-lying rivers, such as the those in the large Crati basin, presently contain alluvial channels close to the catchment mouth.
Although we have few constraints on the long-term erosional behavior of these channels, the assumption of detachment-limited erosion may not be appropriate in these areas. Consequently, we removed all rivers within the Crati catchment from the inverse model data as a test to investigate the potential effect of excluding these rivers on our results (supporting information and Figures S1a and S1b). The difference in predicted uplift between the model containing all rivers and the model without the Crati catchment generally does not exceed ±30 m over a 200-300 kyr time interval ( Figure S1c). The uplift difference is usually <±10 m in the areas containing the Cittanova and Serre faults and the dated marine terraces used to compare model and marine terrace uplift rates (Figures S1 and S2). Consequently, we conclude that our results are not materially influenced by the inclusion of the Crati basin in our inverse model (further details provided in the supporting information).
Finally, our analysis also assumes that slope exponent n = 1 in the stream power model. While there is ongoing discussion about the value of this exponent in a number of settings (e.g., Lague, 2014), we are encouraged that we obtain both a low residual misfit between the majority of longitudinal profiles and good spatial replication of uplift rate patterns denoted by Late Pleistocene marine terraces. We therefore suggest that a detachment-limited stream power model with n = 1 and m = 0.65 is appropriate to derive a plausible uplift history for Calabria over the last 1 Myr. We therefore proceed to analyze what the inverse model implies about the magnitude of regional uplift and the evolution of throw rates for Calabria's faults.

Fault Throw and Regional Uplift
The results from the inverse model provide an opportunity to analyze the temporal evolution of throw rate for the Serre and East Crati faults. The throw of these faults can be analyzed using fluvial profiles because the thickness of hanging wall sediment is small (Roda-Boluda & Whittaker, 2017), as expected where hanging walls have experienced significant uplift. For instance, in the Crati basin, reflection seismic and well data indicate that Middle Pleistocene to Recent sediment thickness does not exceed 200 m (Spina et al., 2011). For hanging wall sediment of negligible thickness, the difference in cumulative uplift between footwall and QUYE-SAWYER ET AL.  hanging wall approximates fault throw. Cumulative uplift from the inversion model was extracted from loci 5 km from the Serre and East Crati faults, in directions perpendicular to the fault traces, at locations where the oldest marine terrace is present in both footwall and hanging wall (Figures 13b and 13d). For the Serre fault, the most extensive footwall terraced area occurs on the southern part of the fault, while for the East Crati fault we could extract uplift from the fault center (Figures 13b and 13d). We interpret the divergence of cumulative uplift at these loci as the onset of faulting. We initially discuss the results for λ s = 1 (Figure 13), and results with damping parameter λ s in the range 0.5-5 are included in supporting information.
Divergence in cumulative uplift indicates that movement on the Serre fault began at ∼650 ka if regional uplift initiates at 1 Ma (770 ka if regional uplift begins at 1.2 Ma; 510 ka if regional uplift begins at 0.8 Ma), which is 300 ka before movement on the East Crati fault (Figure 13). This result agrees with QUYE-SAWYER ET AL. asynchronous fault initiation estimated from marine terraces offset by faults elsewhere in Calabria (e.g., Zecchin et al., 2004). The total amount of throw on the Serre and East Crati faults predicted from fluvial inverse modeling (650 and 800 m, respectively) agrees well with stratigraphic observations and measurements of relief (Roda-Boluda & Whittaker, 2017), which also gives confidence to our results.
We acknowledge that spatial damping of uplift rates in the model, determined by the value of λ s , may affect the estimates of both fault initiation time and total throw magnitude. Results where the damping parameter λ s is varied between 0.5 and 5.0 are presented in Figure S3. These results show that reducing λ s to 0.5 would imply fault initiation at 400 Ma for the East Crati Fault and 600 ka for the Serre fault (assuming initial uplift at 1 Ma). Apparent throw estimates for the present-day are ∼100 m larger than the equivalent interpretation QUYE-SAWYER ET AL.

10.1029/2020TC006076
18 of 26 Figure 11. Uplift rates producing the best-fitting fluvial inverse model. Gray circles = inversion model vertices at 10 km vertical and horizontal separation. Age ranges show propagated uncertainty from age of oldest marine terrace (minimum age assumes uplift began at 0.8 Ma; maximum age assumes uplift began at 1.2 Ma). Marine terraces map produced using median uplift rates from independent observations in Table 1 if λ s = 1, but still lie within the range predicted by independent data. Conversely, an increase in λ s decreases the inferred age of fault initiation, and λ s = 5 produces an unrealistically small throw magnitude for the East Crati fault (500 m).
Assuming uplift initiates at 1 Ma, and λ s = 1, average throw rates since the onset of faulting are 1.1 mm yr −1 for the Serre fault and 2.3 mm yr −1 for the East Crati fault (Figure 13), which are broadly consistent with previous estimates. The modeled throw rate on the Serre fault increases markedly at 100 ka (≈120 ka if regional uplift begins at 1.2 Ma; ≈80 ka if regional uplift begins at 0.8 Ma), which probably records the linkage of fault segments as inferred for many fault arrays in the Apennines and elsewhere (e.g., Faure Walker et al., 2009;Hopkins & Dawers, 2015). An increase in throw rate is also apparent in the λ s = 5 and λ s = 0.5 models ( Figure S3).
The outcome of the inverse model (with uplift initiating at 1 Ma) predicts a gradual increase in throw rate since ∼0.3 Ma for the East Crati fault, similar to the calculated throw rate for the Serre fault (≈4 mm yr −1 ) between 120 and 0 ka. The high throw rates predicted by the fluvial inverse model imply that there is a large seismic hazard in the region, and the rates are faster than those predicted from fault scarp trenching (≥0.44 mm yr −1 ) for one strand of the Cittanova fault (Galli & Bosi, 2002). While the throw rates predicted by these methods are significantly different, they are not necessarily incompatible with each other. First, the apparent discrepancy between the inverse model throw rates and the fault trenching throw rates may arise because of temporal earthquake clustering (fault trenching throw rates are averaged over only 25 ka), spatial variation in slip along the fault array-fault trenching rates were obtained near the north tip of the Cittanova fault (Galli & Bosi, 2002)-or the assumptions used to estimate the initial uplift time in the inverse model. Second, Figure 12 shows that uplift rates predicted by the inverse model only replicate uplift rates calculated from marine terrace elevations within a factor of two. When this uncertainty is taken into account, the fault throw rates predicted by inverse modeling are consistent with those in the central Apennines (e.g., Morewood & Roberts, 2000;Roberts & Michetti, 2004 Table 1. Dashed lines denote theoretical 1:1, 2:1, and 1:2 ratios of uplift rates calculated from inverse modeling and marine terraces. Solid line represents linear regression, through the graph origin, between median uplift rates. (b) Map of Mid-Late Pleistocene terraces. Blue circles refer to the locations of MIS 5 and MIS 7 terraces used in this analysis. Locations enclosed in a black circle are all terraces with independent dating constraints (e.g., OSL, biostratigraphic correlation), which includes MIS 3 and Holocene terraces (Table 1). Catchment averaged erosion rates (0.35 mm yr −1 for the southern tip of the Serre fault and 0.32 mm yr −1 for the East Crati fault) are approximately an order of magnitude smaller than our predicted fault throw rates (Roda-Boluda et al., 2019). The large difference between the modeled uplift rates and erosion rates partially arises because the upstream reaches of many rivers have not reached equilibrium with recent uplift rates, so catchment averaged erosion rates may not balance uplift rates across the entire catchment. The difference between uplift rates and measured erosion rates may also reflect the different time scales of investigation. The mean integration time scales of the cosmogenic nuclide erosion rates are 1.7 and 1.9 kyr, respectively (Roda-Boluda et al., 2019), while the fluvial inverse model only solves for uplift at 36/54 kyr time steps (Figures 10  and 11) so cannot capture rapid fluctuations in erosion rate.
For the Serre fault, similar uplift rates are observed in both the modern footwall and hanging wall prior to ∼0.6 Ma, and the inverse model estimates ≈50 m of uplift between 0.6 and 1.0 Ma (Figure 13a). At the location of the East Crati fault, the inverse model predicts ≈150 m of uplift before faulting begins at ∼0.3 Ma (Figure 13b). These results agree with suggestions of regional uplift preceding the onset of normal faulting in Calabria.
We will now use measured ratios of footwall uplift to hanging wall subsidence to calculate regional uplift from ∼1 Ma to the present-day using the methods in Section 2.4. Terraces are present in both the footwall and proximal hanging wall of the Serre and East Crati faults, and the oldest terrace (Sicilian stage, 0.8-1.2 Ma) can be correlated across the tip of the Serre fault, therefore regional uplift can be isolated using Equation 10. For Calabria, published estimates of the ratio of hanging wall subsidence to footwall uplift, α, lie in the range 1-2, with ≈1.6 calculated from observations on the Armo-Cittanova-Serre fault array (Roda-Boluda & Whittaker, 2017).
For values of α within the published range, the total amount of regional uplift calculated from the inverse model lies between 750 and 900 m, and modeled regional uplift rates increase through time for both the Serre and East Crati faults (Figure 13e). The best-fitting uplift model estimates the same magnitude of regional uplift for both faults from 240 ka to the present. Uncertainties in α produce small <10 2 m error bars on estimated regional uplift (Figure 13e, gray shading) so do not greatly affect our interpretations of regional uplift history.
The oldest terrace offset by the Cittanova fault is well-preserved in both the footwall and in the hanging wall of the Piani d'Aspromonte (Figures 2a and 2b). Therefore, it is possible to estimate the total amount of cumulative uplift solely from terrace observations. In this location, footwall elevation is the sum of regional uplift, Cittanova footwall uplift and a small amount of vertical motion from the nearby Santa Eufemia and Scilla faults. To estimate the magnitude of uplift generated by other faults, we subtracted the height of the footwall in the Petrace drainage basin, beyond the northern tip of the Santa Eufemia and Scilla faults, from the height of the footwall at Aspromonte (Figure 2b). For simplicity, we assume that footwall uplift is similar along strike between these two locations (though uplift may have been greatest toward the center of the footwall). Therefore, the magnitude of the regional uplift at the Cittanova fault, for α between 1 and 2, is 620-920 m ( Figure 13). Footwall uplift since the last interglacial is extracted from the inversion model; hanging wall uplift over the same time period is the height of the marine terrace in the Petrace drainage basin ( (e) Calculated regional uplift assuming hanging wall subsidence to footwall uplift ratio α = 1.6 (dashed lines) and α in the range 1-2 (gray shading). Light blue bars: Estimated regional uplift for the Cittanova fault at the present-day and for the last interglacial. (f) Modified version of Figure 5a. Dashed line represents the magnitude of regional uplift without fault movements, assuming α = 1.6.
The similar magnitude and rate of modeled regional uplift indicate that residual (i.e., nonfault related) uplift is broadly uniform across central Calabria between the Cittanova and East Crati faults ( Figure 13). Such similarities are unlikely if apparent observations of residual regional uplift result from the superposition of footwall uplift from multiple large normal faults, and agrees with suggestions of a long-wavelength uplift process operating in the region.
A striking feature of our modeled regional uplift is the increase in regional uplift rate toward the present-day, an increase which has also been suggested from a comparison of Holocene and MIS 5e marine terrace data . Although uplift models derived from fluvial profile inversion cannot definitively identify the cause of landscape change, we can compare the spatial and temporal uplift calculated by the inverse model to uplift patterns predicted by specific geological processes.
Long-wavelength regional uplift of Calabria has been attributed to processes either operating in the sublithospheric mantle, lower crustal flow, or decoupling of the overriding and subducting plates (e.g., Faccenna et al., 2011;Gvirtzman & Nur, 1999;Westaway & Bridgland, 2007;Wortel & Spakman, 2000). For example, Wortel and Spakman (2000) propose that a tear in the subducting slab, which has been imaged using P wave tomography beneath Southern Italy, could generate long-lived regional uplift due to rebound of the overriding plate. The timing of this slab tear probably coincides with the formation of oceanic crust in the Marsili basin between 1.6 and 2.1 Ma (Guillaume et al., 2010;Nicolosi et al., 2006), where oceanic spreading is indicative of an increase in stretching rate after narrowing of the subducting plate. The results from the inverse model suggest that regional uplift rates have increased toward the present-day which, assuming slab tear is complete, would be inconsistent with decreasing uplift rates predicted during rebound of the lithosphere to reach a new equilibrium elevation (e.g., Buiter et al., 2002). However, we cannot rule out a time delay between detachment of the subducting slab and uplift of the overriding plate, which appears to depend on the depth of subduction (Duretz et al., 2011), or an additional, incipient slab tear of smaller magnitude that may be inferred from mantle seismicity (Maesano et al., 2017). Therefore, only multiple episodes of slab tear, or a time delay between slab tear and rebound, would appear to account for the modeled increase in regional uplift rate.
However, toroidal mantle flow around the subducting slab beneath Calabria has also been inferred from shear wave splitting measurements (Civello & Margheriti, 2004) and predicted from seismic tomography, where it correlates well with high topography (Faccenna & Becker, 2010). Toroidal flow may generate continued uplift as long as roll-back operates, though its rate probably changes through time depending on the trench retreat velocity and plate width (Piromallo et al., 2006;Schellart, 2004). Moreover, toroidal flow may degrade the lithospheric thermal boundary layer (Zandt & Humphreys, 2008), which could produce uplift if the mantle lithosphere is thinned more than the crust (e.g., Esedo et al., 2012). While toroidal flow may be responsible for some uplift of the Calabrian Arc, could toroidal mantle flow account for the temporally increasing uplift rate predicted by Figure 13e? Extension of the lithosphere reduces its elastic thickness, which may make the overriding plate more susceptible to deformation caused by asthenospheric flow (e.g., d 'Agostino et al., 2001). Therefore Calabria may become more easily deformed by toroidal mantle flow as faults grow and interact over time, which could result in a temporally increasing regional uplift rate. We hypothesize that if stretching and thinning of the overriding plate has always occurred alongside regional uplift from asthenospheric flow, then the increase in uplift rate predicted by the inverse model could be consistent with ongoing toroidal mantle flow. Results from the inverse model may therefore emphasize the importance of considering geodynamic processes in both lithosphere and asthenosphere, which is often neglected-or difficult to replicate-in numerical or physical models.

Conclusions
We have utilized a spatial and temporal inversion of river longitudinal profiles to calculate uplift of Calabria, Southern Italy. Erosion rates in a stream power model were calibrated using the age of the oldest marine terrace exposed throughout Calabria. Uplift calculated by fluvial inverse modeling is consistent with uplift rates derived from dated last interglacial marine terraces, which indicates that a simple stream power equation can effectively model uplift and erosion in Calabria. Our results are consistent with variable uplift of Calabria since the Early Pleistocene from normal faults and regional processes, predicting 650 and 800 m of total apparent throw on the Serre and East Crati faults, respectively. Fault throw calculated from fluvial inversion is consistent with independent measurements of structural relief, and increases in throw rate are suggestive of fault interaction and linkage. Fluvial inversion therefore is shown to be a useful technique to analyze fault array evolution. Nonfault related (i.e., regional) cumulative uplift superimposed on three of Calabria's major faults is responsible for ≈850 m of uplift, and regional uplift rates appear to have increased toward the present-day. An increase in regional uplift rate may indicate the combined effect of lithospheric weakness and ongoing mantle flow processes.