Abstract

Unsteady base-level fall at river mouths generates knickpoints that migrate as a transient upstream through the drainage network, climbing at the same rate as long as the fluvial erosion process follows a detachment-limited stream power law. Here we demonstrate unsteady and nonuniform rock uplift using knickpoints as geomorphic markers in streams draining the eastern flank of the Peloritani Mountains (northeast Sicily), the footwall of an ∼40-km-long offshore northeast-southwest–oriented normal fault where the uplift is documented by a flight of mapped and dated Pleistocene marine terraces. Using slope-area analysis on the major streams, we project the tops of prominent knickpoints down to the coast, intersecting the marine terraces, thus providing an age for that specific knickpoint and the paleo–longitudinal profile. We model the migration rate of those dated knickpoints to locally solve for the parameters in the detachment-limited stream power law, and apply the results to model the age of other knickpoints with no clear connection to marine terraces. In summary, we demonstrate that the eastern Peloritani Mountains have been nonuniformly uplifted in an along-strike elliptical pattern, consistent with the general model for the footwall of an active normal fault. A calculation of the long-term erosion rate by the volume beneath the dated paleo–longitudinal profiles reveals a tight positive nonlinear relationship with the modeled normalized channel steepness (ksn). Our analysis provides a method for using knickpoints as geomorphic markers in steep, rapidly eroding landscapes that commonly lack datable river terraces.

INTRODUCTION

Active tectonic and geodynamic processes can be quantified by measuring rock and surface deformation using geomorphic and stratigraphic markers (e.g., Willet et al., 2006). River long profiles (Whipple, 2004), and drainage networks (Perron et al., 2012; Willett et al., 2014), for example, are particularly sensitive to base-level fall, creating geomorphic markers such as terraces (reviewed in Pazzaglia, 2013) and knickpoints (Berlin and Anderson, 2007) by unsteady and transient river incision.

Unsteady base-level fall at the mouths of rivers, such as those that flow directly into the sea, generates knickpoints that migrate upstream through the drainage network. These knickpoints will vertically climb all stream channels at the same rate if their migration follows a detachment-limited stream power law, where the power dependency of the slope term is ∼1 and the drainage area is steady (Howard and Kerby, 1983; Seidl and Dietrich, 1992; Seidl, 1993; Howard, 1994; Seidl et al., 1994; Sklar and Dietrich, 1998; Stock and Montgomery, 1999; Whipple and Tucker, 1999; Kirby and Whipple, 2001; Snyder et al., 2000; Dietrich et al., 2003; Zaprowski et al., 2001). If the age of knickpoints can be known or modeled (Berlin and Anderson, 2007; Crosby and Whipple, 2006), they can be used as geomorphic markers, much in the same way that river terraces are used, to document fluvial incision in response to active tectonics.

In this study we investigate the relationship among prominent knickpoints at various elevations and distances from the stream mouth to dated marine terraces in northeastern Sicily (Catalano and De Guidi, 2003). Our goal is to extend the marine terrace ages inland, using the knickpoints as geomorphic markers, to measure the steadiness and uniformity of rock-uplift and erosion rates. We apply this approach for an actively uplifting region in the footwall of a major offshore normal fault. The advantage of doing the study here is that we can control for rock type, rate of uplift at the coast, and mapped structures. We show convincingly that some knickpoints are genetically related to dated marine terraces, and with this information we can model the upstream migration of knickpoints through the landscape. Our overall goal is to use the knickpoints as geomorphic markers that will allow us to investigate nonuniform and unsteady uplift as well as unsteadiness in watershed divides (drainage capture) in this actively deforming setting. The approach and results are generally portable to any geomorphic setting where an age model for traditional geomorphic markers, like terraces, can be extended more broadly using the drainage network and its knickpoints. Furthermore, the geomorphic analysis allows us to add to the emerging nonlinear relationship between long-term erosion rate and normalized channel steepness index (ksn) (Safran et al., 2005; Harkins et al., 2007; Ouimet et al., 2009; DiBiase et al., 2010; DiBiase and Whipple, 2011; Kirby and Whipple, 2012) and it further informs the tectonic setting for this part of the Apennine chain.

GEOLOGIC AND TECTONIC SETTING

The study area takes advantage of correlative river knickpoints and coastal marine terraces that are well preserved along the eastern side of the Peloritani Mountains (Fig. 1), a rugged, deeply incised terrain that is the southern termination of Calabrian arc (Fig. 1). Northeastern Sicily, located along the Africa-Europe convergent plate boundary (Serpelloni et al., 2007) (Fig. 1), represents a segment of the Apenninic-Maghrebian orogen, where the east-southeast–migrating Calabrian forearc (Calabrian arc; Amodio-Morelli et al., 1976; Scandone, 1979) meets the east-west–oriented, southeast-verging Sicilian orogenic belt, currently colliding against the Nubian plate margin (Dewey et al., 1989; Boccaletti et al., 1990; Lentini et al., 1995, 1996) (Fig. 1). The bedrock of northeastern Sicily consists of several stacked tectonic-metamorphic units and their limestone-dolomite, marly limestone, and lesser siliciclastic, highly stretched Mesozoic-Cenozoic sedimentary cover (Lentini and Vezzani, 1975).

The imbricated basement units of the Calabrian arc, including northeastern Sicily, underwent shortening, followed by crustal extension as they rotated into their current positions (Ghisetti and Vezzani, 1982; Lentini et al., 1995; Monaco et al., 1996; Faccenna et al., 2001; Catalano et al., 2008; Minelli and Faccenna, 2010; Olivetti et al., 2010). Low-grade metamorphic rocks consist of Hercynian metapelites, metarenites and phyllites, metabasites, crystalline limestone, and calc-schists (Capo Sant’Andrea unit, Taormina unit, San Marco d’Alunzio unit, Mandanici unit; Lentini et al., 2000) that are structurally overlain by Aspromonte unit high-grade gneiss locally intruded by granite (Lentini and Vezzani, 1975; Lentini et al., 2000) (Fig. 2).

These tectonometamorphic units overlie the Paleogene Tethyan accretionary wedge terranes (Sicilide units in Fig. 2B; Lentini et al., 1994, 2000), draping a wide ramp anticline of the deformed African continental margin at depth. Locally, the metamorphic basement nappes are unconformably overlain by flysch of the syncollisional late Oligocene–early Miocene Capo d’Orlando Formation (Lentini and Vezzani, 1975, 1978; Lentini et al., 1995, 2000; Catalano and Di Stefano, 1996), consisting of arenaceous-conglomeratic facies passing upward and laterally to postorogenic arenaceous-pelite terrigenous units.

Pliocene–Pleistocene east-southeastward migration of the Calabrian arc, retreat of the European upper plate, and roll-back of the subducted Ionian slab, have been accommodated by northwest-southeast–oriented right-lateral strike-slip faults (Finetti et al., 1996; Lentini et al., 1995, 2000), and the extensional opening of the backarc Tyrrhenian Basin (Ghisetti and Vezzani, 1982; Malinverno and Ryan 1986; Faccenna et al., 2001). Beginning in the middle Pleistocene, the southern edge of the Calabrian arc fragmented as a consequence of the crustal stretching, forming a new linked zone of extension across southern Calabria and eastern Sicily (Monaco and Tortorici, 2000; Catalano et al., 2008).

The geodetic (Goes et al., 2004; Hollenstein et al., 2003; D’Agostino and Selvaggi, 2004) (Fig. 1) and seismologic data, such as the seismic swarm of June–September 2011 (see focal mechanism of mainshock, Fig. 1), are consistent with the kinematic semi-independence of northeastern Sicily from the Calabrian sector of the arc as a separate crustal block (e.g., TGRC versus FOSS; Fig. 1). The kinematic independence of northeast Sicily is further supported by morphostructural and mesoscale structural data consistent with reactivation of Pliocene–Pleistocene northwest-southeast–oriented dextral strike-slip faults (Pavano et al., 2012, 2015) (NOTO versus FOSS; Fig. 1). Mantle upwelling processes at the edge of the Ionian slab have been invoked in order to explain geological, structural, and rapid rock uplift characteristics of the Calabrian arc region (Gvirtzman and Nur, 1999; Faccenna et al., 2001, 2014; Catalano et al., 2011).

Northeastern Sicily, along the Ionian side of the Peloritani Mountains, is bound by the north-northeast–south-southwest–oriented, east-southeast–dipping Taormina fault, an active normal fault in the north-northeast–south-southwest– to north-northwest–south-southeast–oriented extensional belt (Monaco and Tortorici, 2000; Catalano et al., 2008) (Figs. 1 and 2A). Well-exposed, deformed, and dated marine terraces are preserved at the coast in the footwall of the Taormina fault (Catalano and De Guidi, 2003; De Guidi et al., 2003). The marine terraces, ranging in age from MIS 7.1 (Marine Oxygen Isotope Stage) to MIS 3.3, indicate that the coast between Serra San Biagio and Messina (see Fig. 2A for locations) is being uplifted nonuniformly with the greatest uplift occurring near Nizza di Sicilia (∼1.7 m/k.y.) (see Fig. 2A for location), which is approximately midway between the more slowly uplifting fault termini at Taormina and Messina (∼1 m/k.y.; Table 1; Catalano and De Guidi, 2003; De Guidi et al., 2002). Accordingly, the local half-ellipse–shaped uplift of the footwall of the Taormina fault (Figs. 3A, 3B) is superimposed on the 1.1 m/k.y. regional uplift (Catalano and Di Stefano, 1997; Catalano and De Guidi, 2003).

GEOMORPHOLOGICAL SETTING

The topography of northeastern Sicily is dominated by the northeast-southwest–trending Peloritani ridge (Fig. 1), which separates a moderately asymmetric, more gentle Tyrrhenian (western) flank from the steeper Ionian (eastern) side (Fig. 2B). There is little orographic precipitation effect in the Peloritani Mountains; both flanks receive an annual rainfall of ∼900 mm (Fig. 2C; see www.osservatorioacque.it). Along the eastern flank of Peloritani Ridge, non-uniform topographic slope and local relief shows that the highest values of these two metrics mainly fall in the central portion of the studied area (Figs. 3C, 3D).

The Ionian-flank rivers are arranged in subparallel elongate ∼120°-oriented watersheds that share the Peloritani ridge divide with the opposing Tyrrhenian flank drainages. There is little field evidence for watershed divide migration in the Peloritani Mountains. There is a general absence of obvious windgaps and barbed tributaries. Some channels locally resembling a barbed pattern are instead probably influenced by northeast-southwest–oriented bedrock structures.

The range is asymmetric in cross section (Fig. 2B) with a widely incised maximum elevation envelope, which is mirrored in the relief, and a mean elevation accordant with the subenvelope. The eastern side shows a more uniform and steeper maximum elevation envelope not as deeply incised as its western counterpart. Between the eastern and western flanks, the central part of the mountains is characterized by a relatively low relief landscape that includes the drainage divide (Fig. 2B). Although Pleistocene paleoclimatic nonuniformity cannot be ruled out and is not known for this part of Sicily, given the relatively uniform rock type and uniform Holocene climate and rainfall, the topographies of the Ionian flank watersheds are consistent with a recent tectonically driven base-level fall forcing.

Marine terraces and associated fluvial fan deltas have been mapped at the mouths of the rivers that drain to the Ionian coast. These terraces are correlated to eustatic highstands of the middle and late Pleistocene, ranging from MIS 7.1 (200 ka) to MIS 3.3 (60 ka), providing a local chronostratigraphic model (De Guidi et al., 2002; Catalano and De Guidi, 2003; Catalano et al., 2003) (Table 1). In summary, these terraces create a coast-parallel antiform trend (Fig. 3A), increasing in elevation from the region near Messina (Briga section in Table 1), where they infer an uplift rate of ∼1.1 m/k.y. toward the south, reaching maximum elevations near Nizza di Sicilia with an inferred uplift rate of ∼1.74 m/k.y. (Table 1). South of Nizza di Sicilia, the entire flight of terraces starts to plunge and converge, reaching the area of the Alcantara River mouth, with an inferred rock uplift rate of ∼0.9 m/k.y. We note that the correlations and rates from the cited studies do not take into account the effects of the glacial isostatic adjustment for local eustatic base level for the Mediterranean region (Lambeck and Purcell, 2005; Peltier, 2004; Stocchi and Spada, 2009), but the correction in our case would be small and is not thought to contribute much to the subsequent uncertainties in our analysis.

LONG PROFILE ANALYSIS

Steady-state rivers show concave-upward power-function long profiles (equilibrium profile; Mackin, 1948; Hack, 1957), traditionally attributed to the well-known inverse discharge-slope relationship necessary to generate the shear stresses that entrain and move sediment (Gilbert, 1877; Baker and Ritter, 1975; Leopold and Bull, 1979; Bull, 1991). Deviation from this condition is commonly indicative of a transient phase of landscape change, represented by knickpoints in the river profile.

Impulsive base-level fall at the mouth of a stream drives fluvial incision into bedrock via a knickpoint that migrates upstream. For settings where rock uplift and rock type throughout the watershed are uniform, the channel reach upstream of the knickpoint remains graded to the higher, paleo–base-level condition. In this case, every knickpoint can be considered as a time-dependent geomorphic marker that migrates through the landscape, lowering the river profile to the new base level. The specific erosion processes, commonly modeled as a detachment-limited stream power dependent relationship (Crosby and Whipple, 2006; Berlin and Anderson, 2007), dictate the rate at which the knickpoint migrates upstream through the bedrock channel.

The processes of fluvial incision into bedrock are rooted in the stream power incision model (Howard and Kerby, 1983), widely applied and accepted in quantitative geomorphology as a basis for long profile and knickpoint migration modeling. For bedrock channels, considering steady, uniform flow and that discharge scales linearly with drainage basin area (A, m2) (Seidl and Dietrich, 1992), the general form of the stream power erosion (E, m/yr) model is 
graphic
where S (dimensionless) is the stream slope and the m and n exponents are positive constants depending on the erosion processes, hydrological features of the stream, and channel geometry (Howard and Kerby, 1983; Howard, 1994; Whipple and Tucker, 1999; Crosby and Whipple, 2006), and K [m(1–2m)/yr] is the detachment-limited erosion coefficient linked to both climatic and lithological factors (Snyder et al., 2000). The time rate of change in the elevation of the river channel (dz/dt) is simply the difference between the rate of rock uplift (U) and the rate of incision, 
graphic
By substitution into Equation 1, with slope (S) expressed as dz/dx, 
graphic
we arrive at a general expression that can be rearranged for the equilibrium case where dz/dt = 0, and, solved for slope, reveals 
graphic
Equation 4 has the same form as Hack’s law (Hack, 1957; Flint 1974; Snyder et al., 2000; Kirby and Whipple, 2001), 
graphic
and 
graphic
where ks (m) is the steepness index that describes the stream gradient and θ is the dimensionless concavity index (m/n), a morphometric parameter that defines the rate of change of slope in the channel gradient as a function of the drainage area growing downstream. Comparing Equations 4 and 5 implies that ks scales with rock uplift and/or rock type (K) (Whipple and Tucker, 1999; Snyder et al., 2000). Given a constant K, several studies now have suggested a linear, or near-linear, relation between ks and U, whereas they have confirmed that θ remains relatively constant for a wide range of tectonic and climatic settings, commonly assuming values ranging between ∼0.3 and 0.6 with a mean value of 0.45 (Snyder et al., 2000; Kirby and Whipple, 2001; Kirby et al., 2003).
For the detachment-limited plucking erosion process case (Hancock et al., 1998), many studies (Whipple et al., 2000) have shown that the n exponent is ∼1, which allows for a simplification of Equation 1, assuming that E = dz/dt (Berlin and Anderson, 2007), to 
graphic
where dx/dt represents the celerity of a headward-migrating knickpoint. This celerity model (Berlin and Anderson, 2007) assumes that the basin area is steady, and does not include variations in channel width, sediment flux, erosion threshold, or the influences of unsteady climate. Nevertheless, based on a dynamic step approach (Crosby and Whipple, 2006) the detachment-limited erosion coefficient (K) and the exponent of drainage area (m) can be solved for using a forced two-parameter search. The best combination of these parameters is set by the lowest misfit between the modeled and the real position of knickpoints celerity using the contributing drainage area. In fact, if dxn represents discrete, constant distance interval covered by a retreating knickpoint along the river profile and An is the corresponding drainage area, changing with distance, the time interval dtn can be obtained by 
graphic
and the total time t required for a given knickpoint to migrate along the river profile from the mouth up to its current position is represented by 
graphic

METHODS

Quantitative topographic metrics are extracted from a digital elevation model (DEM) with a resolution of 10 m, assembled from aerial photogrammetric data (flight ATA 97). We explore the hypsometric integral (HI) of Peloritani Mountains basins that expresses the distribution of relief as a percentage of normalized drainage area. We use the Adb Toolbox open-source software, freely accessible on the Italian Ministry of Environment official Web site (http://www.pcn.minambiente.it/GN/) for this analysis. The corresponding hypsometric integral is computed as elevation-relief ratio (Pike and Wilson, 1971; Mayer, 1990): 
graphic
where Zmean, Zmaximum, and Zminimum are elevation in meters.

Next, we extracted the fluvial network, deriving stream long profile data and drainage basin areas, that we modeled in Excel, SigmaPlot, or using in-house developed FORTRAN scripts. Channel long profiles were smoothed using a loess filter utilizing a local moving linear regression window one-tenth of the entire data set.

We calculate the stream length index (SL; Hack, 1957, 1973), 
graphic
where dz/dx is the channel reach slope and L is the channel distance (m) from the center of that reach to the drainage divide. The SL index, a rapidly calculated complement but not a replacement for slope-area analysis, aids in identifying anomalously steep channel reaches that may be either upstream-migrating transients, or lithologically fixed knickpoints (Miller, 1991; Gardner, 1983; Weissel and Seidl, 1997; Duvall et al., 2004; Frankel et al., 2007). Our study area features two groups of trunk channels, one that mostly or totally flows on low-grade metamorphic rocks and one that mostly or totally traverses high-grade metamorphic rocks. We can distinguish knickpoints that are clearly lithologically controlled and fixed (KPl) from the migrating ones that are not lithologically controlled (KPm) by plotting the high SL indices on a 1:100,000 scale geologic map (Lentini et al., 1998).

In slope-area space we determine ks and θ for the trunk channel reach located directly upstream from each migrating knickpoint. Using a rearranged Equation 5, we are able to reconstruct the paleo–longitudinal profiles of the river (Clark et al., 2006), projecting them downstream from the top of a knickpoint to the coastal terraces (Olivetti et al., 2012; Legrain et al., 2014). Similarly, as a constraint in the knickpoint celerity model, we use three knickpoints, which we recognize along the D’Agrò, Alì, and Briga torrents (Fig. 2C), clearly correlative to the marine terraces assigned in Catalano and De Guidi (2003) to MIS 5.1 (80 ka). Using Equations 7–9, the forced two-parameters search, and dynamic step approach (Crosby and Whipple, 2006; Berlin and Anderson, 2007), we achieve the best fit of the detachment-limited erosion coefficient (K) and the exponent of drainage area (m), that solves for the knickpoint location following 80 k.y. of upstream migration. Thus, we use these K and m values in Equation 9 to forward model the age (t) of any other migrating knickpoint, enabling a comparison to the age of the downstream intersecting marine terrace to our modeled age, as both an independent test of the model and in assembling a population of knickpoints that we can use as dated geomorphic markers.

Using a reference area and concavity (A–θref), where θref is 0.45, which is similar to our measured mean concavity (0.44), and assuming a constant K permits us to explore relationships between a normalized channel steepness index ksn and rock uplift (Snyder et al., 2000; Kirby and Whipple, 2012) using Equations 4 and 5. Also following Snyder et al. (2000), we can and do make the alternative assumption that U is steady, as defined by the terrace uplift rate (Catalano and De Guidi, 2003; Table 1), and explore the variation in K with respect to spatial variations in rock type.

The reconstructed river paleoprofiles are next used to measure rates of river incision (I) and estimate rates of drainage basin–wide long-term erosion rate (Elt). Elt is defined as the ratio of rock volume (m3) below a paleoprofile to the eroded area (m2) defined by the intersection of a horizontally projected plane, with elevation fixed by a paleoprofile, and the DEM topography (Frankel and Pazzaglia, 2005) divided by the age of the paleoprofile. For basins that host more than one migrating knickpoint and paleoprofile (D’Agrò, Allume, Briga, and Galati torrents; Fig. 2C) we choose, as representative rate of erosion for each basin, those calculated only from the youngest, lowest paleo–longitudinal profile (typically the 80 ka profile; Table 2). Similarly, incision rate (I) is calculated as the mean vertical distance of the modeled paleoprofiles to their present river profiles divided by the age of the paleoprofile. In several basins we can take advantage of multiple paleoprofiles to calculate incision incrementally (Ii) below subsequent paleoprofiles, avoiding Sadler effects that are introduced by comparing to the modern channel (Gallen et al., 2015). In synthesizing all data, we plot the long-term erosion rate (Elt) against the mean vertical incision rate (It).

These measures of incision and erosion allow for a direct comparison to the basin’s normalized steepness index (ksn) and the local relief (LR, m). The latter is obtained for each basin by sampling the 10 m DEM by a 100-m-wide search radius, then computing their standard deviation. Covariance of erosion rate and mean relief should be linear, following the well-known Ahnert (1970) relationship, where slopes are below the landsliding threshold (Gallen et al., 2013), a condition that we note occurs, but is not generally true for our landscape.

RESULTS

We present first the basin-scale observations of hypsometry (Fig. 4) and relief before discussing more focused consideration of the long profiles, knickpoint retreat, and rates of incision and erosion (Table 2). All of these metrics are placed in the context of the rock uplift rate, measured from the marine terraces, and rock type, subset as a percentage of low-grade or high-grade metamorphic rocks (Fig. 5A). The hypsometric curves of most basins (Figs. 4A, 4B) have asymmetric sigmoidal shapes that are above the hillslope fluvial equilibrium curve. The hypsometric integral (HI) ranges from 0.35 to 0.5 with a mean of 0.42, except for the segmented curve of the Galati basin that is below the equilibrium curve and has a correspondingly low HI of ∼0.28 (Table 2; Fig. 4A). As a whole, the trend of the hypsometric curves is consistent with catchment topography. As suggested by the clustering of the elevation distribution modes (Fig. 4C), most of basin elevations are higher than the mean of ∼500 m, with ∼13% of the drainage area at ∼600 m. When plotted as a function of drainage basin position along the coast, from south to north, the HI values do not covary with marine terrace uplift rate (Figs. 5A, 5B). In contrast, local relief (LR), with an average value of ∼75 ± 25 m, follows a general elliptical pattern rising from low values of ∼60 m in the south to high values of ∼100 m in the center, before falling back to low values of ∼50 m in the north, a pattern that mimics the marine terrace uplift rate (Figs. 5A, 5B).

From the drainage basin-scale observations we now characterize the river channels, their long profiles, concavity, steepness, and knickpoints. Long profiles have concavity values (θ) ranging from 0.32 to 0.66 (Table 2), in good agreement with theoretical values (0.35–0.6; Whipple and Tucker, 1999) and well within the range published in previous studies (0.3–1.2; Snyder et al., 2000; Kirby and Whipple, 2001; Kirby et al., 2003; Whipple, 2004). The highest θ values (0.57–0.66) tend to cluster in the northern part of the study area (Fig. 5C) coincident with the outcrop of the highest grade metamorphic rocks, whereas to the south the rivers flow on lower grade metamorphic terrains (Fig. 5A), and have correspondingly lower concavities, with average values of ∼0.36 (Table 2; Fig. 5C).

Similarly, the normalized channel steepness (ksn) index using a reference concavity of 0.45 has relatively low values (111–165 m), and mean of 142 m, for the southernmost basins that are underlain by low-grade metamorphic rocks (Table 2; Fig. 5A). Coincident with the outcrop of high-grade metamorphic rocks in the northern basins (Fig. 5A), ksn increases, reaching a maximum value of ∼197 m (Fig. 5C). In general, moving from south to the north ksn mimics the elliptical distribution of rock uplift (Fig. 5A), with maximum values in the central portion of the studied area and minimum values at the two termini (Fig. 5C). When subset by rock type, ksn shows a positive linear relationship with rock uplift (Fig. 5D).

In contrast, the erosion coefficient (K) mirrors the trend for channel steepness (Table 2). In the southern half of the study area K is relatively high with a mean of ∼1.30 ± 0.18 × 10−5 (m0.14/k.y.), that drops abruptly for the northern half of the study area to a mean of ∼0.79 ± 0.27 × 10−5 (m0.14/k.y.) at about the point where the percentage of high-grade metamorphic rocks increases to >20%. There is no clear relation of K to uplift rate (Figs. 5A, 5C).

In detail, the long profiles contain knickpoints that we intend to use to extend the basin-wide and long profile data summarized in Figure 5 upstream into the footwall block of the Taormina fault. We identify 19 knickpoints (Figs. 6 and 7; GSA Data Repository Figs. DR1–DR41), 16 of which we determine to be migrating (16 KPm; Table 3) and only 3 of which are clearly lithologically controlled (3 KPl) by carbonate rock types in the Letojanni torrent (Figs. 6 and 7). These lithological knickpoints also are easily recognizable in the field as waterfalls (inset, Fig. 7). In contrast, in the northern part of the study area there are two basins, the Itala and Giampilieri torrents, that do not show any obvious knickpoints along their profiles.

Modeling of the MIS 5.1 knickpoint (80 ka) location in the D’Agro, Alì, and Briga drainages (Table 3) using the two-parameter forced search for K and m returned values of 1.69 × 10−4 [m(1–2m)/yr] and 0.34, respectively. Forward modeling using these values makes a prediction about the age of other knickpoints in the long profiles and also allows us to compare the model prediction of any knickpoint to its real location (Fig. 8; Table 3). The best fit between modeled and actual knickpoint positions ranges between 1% and the 9% of the stream’s total length (Fig. 8; Fig. DR5).

Similarly, a plot of these modeled knickpoint ages at the downstream projection of their reconstructed long profiles shows remarkable correspondence to the marine terraces (Fig. 9A). With the exception of the oldest knickpoint of Briga torrent, showing a paleomouth projection 40 m higher than its presumed correlative 125 ka paleoshoreline, the distribution of the projected paleothalwegs (triangles in Fig. 9A) very closely mirrors the coastline-parallel convex-upward trend of the marine terraces (lines with error bars in Fig. 9A), within uncertainties of ∼±15 m, comparable to the constant ±10 m (Catalano and De Guidi, 2003) associated with each marine terrace inner edge position. The correspondence is particularly strong for the MIS 5.5 (125 ka) and MIS 5.1 (80 ka) knickpoints, whereas it is more difficult to see a correspondence for the few and sparse data associated with MIS 7.1 (200 ka) and MIS 3.3 (60 ka). There is no marine terrace correlated to the oldest migrating knickpoint recognized along the Savoca torrent, predicted to be to the MIS 9.3 (330 ka) in our model and projected at an elevation of ∼460 m. In summary, we find that the ages of the modeled paleo–longitudinal profile fit very well those deriving from dated marine terraces, plotting essentially on a 1:1 line (Fig. 9B).

Average knickpoint modeled upstream retreat rates scale with basin size and vary from a high of ∼87 m/k.y. for the youngest knickpoint of the D’Agrò torrent to ∼22 m/k.y. for the oldest knickpoint, of the Galati torrent (Fig. 9C). Overall, the mean rate of knickpoint retreat (being fully cognizant that this rate is not constant but rather slowing as the knickpoints migrate upstream) is ∼50 m/k.y., comparable to global published compilations (e.g., Loget and Van Den Driessche, 2009).

With the modeled long profiles, knickpoints, and their correlation downstream to the marine terraces established, paleo–longitudinal profiles can be used to measure long-term rates of incision (I) and erosion (Elt) for 10 of the 13 basins (Fig. 10; Table 2). The Elt values range between 0.22 m/k.y. (Galati torrent) and 0.92 m/k.y. (Briga torrent), and average ∼0.45 m/k.y. for all watersheds. The highest values occur using the 80 ka paleo–river profile for the Briga torrent (0.92 m/k.y.) and Alì torrent (0.66 m/k.y.) (Table 2). In comparison, the rate of vertical incision (I) ranges between ∼0.39 m/k.y. and ∼1.18 m/k.y., with average value of ∼0.69 m/k.y. (Table 2). These incision rates are ∼150% greater than the average Elt (Fig. 10; Table 2) and are not steady with all rivers, generally decelerating in rate toward the present (Fig. 11). We apply uncertainties to these incision and erosion rate values assuming a reasonable uncertainty of ∼2 m in vertical elevation on DEMs, and an uncertainty of the age of the geomorphic markers (marine terraces) of ∼10 k.y. following Equation 1 in Lambeck et al. (2004).

DISCUSSION

We use our results to forward a conceptual model that describes knickpoints as geomorphic markers through an exploration of the relationships among rock uplift, incision, knickpoint retreat, erosion, and rock type. Specifically, we interpret our data to indicate a primary rock uplift (base-level fall) control on these key geomorphic metrics with secondary contributions from rock-type variability. The correspondence between knickpoint elevation and the pattern of uplift inferred from coastal marine terraces (Figs. 8 and 9) indicates that the knickpoints are migrating through a single, largely coherent tectonic block in the footwall of the Taormina normal fault, and that drainage area for the stream network is largely unchanged over late Pleistocene time scales. This interpretation dispenses quickly with one of the key questions we posed at the beginning of our study regarding the uniformity of uplift of the Taormina fault footwall, that at least appears to be a coherent block with little obvious internal deformation.

In the context of this uplifting footwall block, we burrow into the details of interpreting the geomorphic metrics and long profile forms. We deal first with the influence of rock type in trying to qualitatively assess its role in shaping the upstream migration of the knickpoints (Fig. 5). The coast-parallel distribution of the basin hypsometric integral is slightly sensitive to uplift, as evidenced by the curves near the reference curve of the equilibrium between hillslope and fluvial erosion processes. This pseudoequilibrium condition reflects watershed morphology characterized by narrow divides with a general absence of relic low-relief surfaces and basins of the fixed area (Strahler, 1952). In contrast, the influence of rock type on the hypsometry of small catchments (Lifton and Chase, 1992; Hurtrez and Lucazeau, 1999) leads to the largest variation in HI occurring in the Letojanni watershed, where different level of resistant limestone cross the watershed at different elevations. Conversely, the hypsometric integral dependence on basin area and shape (Willgoose and Hancock, 1998; Hurtrez et al., 1999) seems to influence the anomalously low HI of the very small and narrow Galati drainage.

Similarly, and as expected, θ appears to be insensitive to rock uplift. The increases in profile concavity observed in the three northernmost watersheds, draining high-grade metamorphic rocks, are suggestive of a change in bedrock erodibility. In this regard, following Snyder et al. (2000), we note that ksn, calculated using Equation 5, and K, calculated using Equation 6, are generally anticorrelated; ksn increases dramatically for the northern watersheds that drain high-grade metamorphic rocks and K correspondingly falls. The erosion constant (K) is affected by several factors, some of which are independent and others mutually dependent (Snyder et al., 2000; Goswami et al., 2012). For example, nonuniform precipitation can affect K, but we suspect that this is minimal for our study under at least modern climatic conditions (Fig. 2C), a conclusion also reached by Goswami et al. (2012). Channel width, alluvial cover, sediment flux, and debris-flow processes also affect K, but these fluvial system characteristics likely also vary with hillslope processes and soils, which are rooted in the bedrock underlying the watershed. For example, Goswami et al. (2012) did not find a systematic relationship between channel widths and drainage area and argued for local dense fracturing of the bedrock as the underlying cause. We similarly conclude that the sharp change in K at the transition from low-grade to high-grade metamorphic rocks is primarily driven by rock type. South to north, local relief (LR) shows an increasing trend, reaching maximum values at the Fiumedinisi and Alì basins. Toward the north, local relief starts to decline, down to values comparable to those of the southern basins. Essentially, it generally reproduces the rock uplift trend along the coast, rather than showing any influence from the change in rock type.

On one hand, the change in rock type parallel to the coast limits our ability to directly correlate rock uplift with ksn and our estimates for erosion and river incision. On the other hand, it also presents the opportunity to explore the effect of rock type, assuming that the precipitation and fluvial processes discussed here are subordinate factors. Starting in the southern terminus of our study watersheds and continuing north for 15 km we note that there is a corresponding rise in ksn values with rock uplift rates (Figs. 5A, 5C). The normalized steepness index then rises quickly at ∼23 km (Fig. 5C), coincident with the transition to the high-grade metamorphic rocks. Although noisy, the trend in ksn decreases north of 23 km, mimicking the fall in the rate of rock uplift (Fig. 5A). The plot of ksn dependence on the rate of rock uplift (Fig. 5C) collectively has a positive linear trend (Goswami et al., 2012), but the correlation is much stronger when the data are divided into subsets by rock type (Fig. 5D). We conclude that channel steepness in our study area is being driven primarily by rock uplift and is secondarily enhanced by harder rocks. This conclusion indicates that detachment-limited fluvial processes described by Equation 1 and its derivatives is a reasonable model for the Peloritani Mountains landscape and justifies our subsequent knickpoint modeling approach and results.

Using the profile concavity and steepness upstream of a knickpoint, we are able to show a good correspondence between the mapped elevation of marine terraces and projected paleo–longitudinal profiles (Figs. 7 and 9). These results are consistent with unsteady, episodic base-level fall at the coastline, probably linked to the interplay of rock uplift and glacioeustasy. A notable exception to the predicted projected elevation of a paleoprofile and the actual elevation of its corresponding marine terrace occurs for the Briga torrent drainage. Here, the paleo–river mouth associated with the MIS 5.5 marine terraces is predicted to be ∼40 m higher than the corresponding marine terraces inner edge, which would be consistent with a more rapidly migrating knickpoint driven presumably by greater than expected discharge. However, there is no obvious evidence of drainage capture in this watershed, which leaves local tectonic deformation associated with an unmapped fault or a local change in rock type that significantly alters erodibility (K). A decrease in K due to intense local rock fracturing, such as the mechanism described in Goswami et al. (2012), would predict that the Briga knickpoints can migrate headward faster than adjacent watersheds, traveling 250 m upstream and correspondingly rising 20 m higher. However, softer rocks also imply lower channel reach steepness, an observation documented by Goswami et al. (2012); this leaves open the possibility that the rise in knickpoint elevation in the Briga watershed could be caused by local fault-related offset that is gently deforming the Taormina footwall in this particular location.

Other than the Briga torrent drainage, the confidence that we have in successfully predicting knickpoint elevations across the east flank of the Peloritani Mountains opens the opportunity to explore the relationship between the long-term erosion rate (Elt), ksn (Cyr et al., 2010; DiBiase et al., 2010; Fig. 10A), local relief (LR) (Fig. 10B), and incision (I) (Fig. 10C). Our ksn and Elt data directly overlap published data from similar studies that have proposed a general power law relationship between ksn and erosion rate when erosion rate is <∼0.40 m/k.y. (Safran et al., 2005; Harkins et al., 2007; Ouimet et al., 2009; Miller et al., 2011; DiBiase et al., 2010; DiBiase and Whipple, 2011; Kirby and Whipple 2012), following 
graphic
where ksn is normalized channel steepness index, Elt is long-term erosion rate, and Φ is the power law exponent (Fig. 10A). According to this relation, our ksn becomes insensitive to Elt for values >∼0.40 m/k.y. (Fig. 10A). This finding is significant given the diverse ways in which the erosion rates reported in studies cited here are calculated and how they contrast with our approach. The Santo Stefano watershed plots above the rest of our data in Figure 10A. Evidently there are reasons specific to this basin, such as being fully contained in the high-grade metamorphic rocks, that cause the channel to be particularly steep (high ksn).

The long-term erosion rates (Elt) show a complicated dependence on topographic relief in the basins that does not improve when plotted in terms of rock type (Fig. 10B). The well-known Ahnert (1970) linear relationship between erosion rate and local relief holds for creep-dominated soil-mantled landscapes where the maximum erosion rates do not exceed 0.40 m/k.y. (Montgomery and Brandon, 2002). Basin-wide erosion rates in excess of that value seem to be landslide dominated (Roering et al., 1999; Heimsath et al., 2012), and we note that most of our average reconstructed erosion rates are ∼0.40 m/k.y. or greater, consistent with a landslide-dominated erosion process for our basins, an interpretation that is qualitatively supported by historic accounts of landslides as well as visual identification from publicly available satellite imagery.

In a steady or nearly steady-state landscape, our measure for long-term erosion Elt should not only co-vary with the rate of rock uplift and ksn, but also be similar to other measures of landscape change, namely thermochronologic and cosmogenic estimates of erosion and river incision. The two published cosmogenic erosion rates of 0.97 ± 0.11 and 1.44 ± 0.40 m/k.y. (Cyr et al., 2010) and our measured river incision rates that range from ∼0.5 to 1.2 m/k.y. are similar to rock uplift rates measured from the marine terraces. The long-term erosion rates are slower, but still comparable. Apatite fission track and apatite (uranium-thorium)/He results from the Peloritani Mountains suggest post-Pliocene rates of erosion as high as 1.3 m/k.y. for the western part of the range where analyzed samples are buried by Pleistocene deposits (Olivetti et al., 2010). However, for the eastern side of the range for samples collected just south of our study area near Taormina (Fig. 1), post-Pliocene rates of erosion are modeled to be between ∼0.3 and 0.4 m/k.y. (Olivetti et al., 2010), very similar to our mean Elt rates of 0.45 m/k.y. and similar to the long-term erosion rates determined cosmogenically for nearby Calabria (Cyr et al., 2010; Olivetti et al., 2012). With respect to incision rate (I), we find that it mostly scales with, but is consistently more rapid than, Elt (Fig. 10C), indicating an imbalance between the lowering of the trunk streams and the basin-wide erosive response and the growing of topographic relief. This imbalance could mean that our long-term erosion rate estimations are a minimum, given that they are calculated for the lower part of the watershed where the valley bottom has filled with Holocene alluvium. Alternatively, as suggested by the relationship between HI and rock uplift (U), the imbalance could also result from the hillslope lag time response to channel lowering that, generally, also depends on rock type (K), as suggested by the higher discrepancies between Elt and I for basins developed on high-grade metamorphic rocks (Table 2).

We find that the total incision rate (It; m/k.y.), computed for any paleo–longitudinal profile considering as datum the modern channel bed, records rate unsteadiness over time, showing a general decrease for the past 200 k.y (Fig. 11A). This trend is reflected in the distribution of the incremental incision rates (Iri; m/k.y.) and the incremental incision (Ii; m) between subsequent vertically stacked paleo–longitudinal profiles (Figs. 11B, 11C) where the rate calculation does not include the modern channel. Except for the D’Agrò torrent that initially increases in erosion rate between 200 and 125 ka (Fig. 11C), the majority of the results point to slowing incision. Decelerating incision is particularly evident when the rate calculation includes the modern channel (Fig. 11A), but we note the complications that arise in interpreting incision rate data when using a mobile datum like the modern channel (Gallen et al., 2015). In considering only the rates derived from stacked paleoprofiles (Figs. 11B, 11C) and no lag time between rock uplift and fluvial incision, it is possible that the incision data indicate a slowing in the rate of footwall uplift. Alternatively, considering a significant lag time, they indicate a deceleration in incision rate for a footwall that was uplifted impulsively and has now slowed or stopped, with the rivers carving into and relaxing their concavities. In either case, the incision established significant relief in all of the watersheds that are capable of generating and moving large volumes of sediment from valley slopes by landslide and fluvial processes, sediments that, once deposited, transiently and locally have aggraded the valley bottoms.

CONCLUSIONS

Most large knickpoints in the channels draining the east flank of the Peloritani Mountains of northeast Sicily have been created by unsteady glacioeustatic base-level fall at the coast superimposed on the rapidly uplifting footwall block of the Taormina normal fault. Knickpoint celerity and longitudinal profile modeling coupled with a stratigraphic age model derived from uplifted marine terraces define a suite of geomorphic markers that constrain the long-term rates of erosion and river incision as well as the distribution of rock uplift of the footwall block.

A detachment-limited knickpoint celerity model successfully predicts the headward migration of knickpoints from the coast to their current positions in their respective watersheds, documenting fluvial incision in response to active tectonics. In this respect, the knickpoints have ages, and can be correlated from basin to basin like any other geomorphic marker. This correlation reveals that deformation of the Peloritani Mountains follows our expectations for uplift of the footwall of the Taormina fault (Catalano and De Guidi, 2003; Catalano et al., 2008), confirming that it is a coherent, rigid tectonic block, excluding any significant internal deformation. Watersheds for which the celerity model fails to predict knickpoint locations could be explained in terms of channel incision through locally variable rock type and the concomitant variations in K, which are not captured by the constant K input to the model, excluding the involvement of any processes linked to the watershed divides unsteadiness (e.g., fluvial capture).

Moreover, the modeled long profile analysis points out that the trunk streams are almost uniform in their concavity, well within the range of the published values (θ = 0.45), but have normalized steepness indices (ksn) and erodibility (K) that vary in terms of rock uplift and rock type, respectively. Changes in rock type from low-grade metamorphic rocks in the southern two-thirds of the study region juxtaposed with high-grade metamorphic rocks in the northern third of the study area serve to uniformly increase ksn at the expense of K for the harder high-grade metamorphic rocks, but the general impact of tectonic forcing of the uplifting footwall block is still evident in the distribution of ksn. Similar to the conclusion of Goswami et al. (2012), we do not observe variations in K that obviously correspond to the distribution of precipitation. All of these results are consistent with the findings of a similar tectonic and climatic setting in northern California (Snyder et al., 2000).

Base-level falls, as documented by upstream-retreating knickpoints, trigger river entrenchment, steeper slopes, and erosion. Basin-wide erosion rates calculated by the volume of rock removed seaward of a given knickpoint and beneath its paleoprofile projection do not show dependence on mean local relief, regardless of the dominant rock type, arguing for steep slopes that are landslide governed rather than soil mantled and creep dominated. Growing relief and landslide processes are supported by the imbalance between basin-wide long-term erosion rates (∼0.45 m/k.y.) and mean incision rates approaching 1 m/k.y. However, the incision rates appear to be decelerating toward the present, due to either landslides transiently feeding larger amounts of material to the fluvial beds shielding the incision process, or simply reflecting the relaxing stage of the fluvial incision processes into an impulsively uplifted footwall block.

In this regard, our findings indicate that the Ionian side of the northeastern Sicily responded to the tectonic control of a broader extensional belt (Monaco and Tortorici, 2000; Catalano et al., 2008), over the past ∼200 k.y., as a coherently uplifted block. These inferences are consistent with the broader tectonic picture of the southern termination of the Calabrian arc, where northeastern Sicily behaves as an independent crustal domain that, since the middle Pleistocene (600 ka; Catalano and Di Stefano, 1997), uplifts and moves at different rates with respect to the adjacent areas, reactivating previous faults (Pavano et al., 2015).

Overall, the results of our study strongly emphasize the opportunity to employ knickpoints as useful geomorphic markers, documenting fluvial incision in response to rock uplift, with particular applications to those regions lacking in fluvial and/or marine terraces. As an additional applicative outcome, these conclusions could have important impacts on the geologic risk assessment of the region, with particular reference to the region of watersheds crossed by geologically recent retreating knickpoints, where the instability of valley walls could be more impressive, as documented by several catastrophic debris flows in the past decade (e.g., 1 October 2009). In addition, urbanization along the catchment banks has increased almost as quickly.

We thank Lithosphere associate editor Kurt Stuewe, Valerio Olivetti, and an anonymous reviewer for constructive reviews that significantly improved the manuscript. This work represents part of the doctoral dissertation of Pavano, who was supported by the University of Catania to study abroad at Lehigh University for one semester. This paper was funded by the Ministry of Education, Universities and Research of Italy (Progetti di Ricerca di Interesse Nazionale project, Birth and death of oceanic basin: Geodynamic processes from the rifting to the continental collision in the Mediterranean and circum-Mediterranean orogens, grants to G. Capponi and Catalano).

1GSA Data Repository Item 2016317, Figures DR1–DR4: Results of the longitudinal profile analysis of Pagliara, Allume, Fiumedinisi, Alì, Itala, Giampilieri, Briga, Santo Stefano, Galati, and Mili torrents. The main graph displays the change in both the drainage area and the SL index referenced to underlying rock type (G—gneiss; F—phyllite; C—carbonates; COF—Capo d’Orlando flysch). The colors of the modeled paleothalwegs are coordinated to both the slope-area data of the above A-S plots used for the reconstruction and the regressed reaches of long profile upstream from any knickpoint. Table DR1: Modeled versus actual knickpoint location, is available at www.geosociety.org/pubs/ft2016.htm, or on request from editing@geosociety.org.