Isolating Lithologic Versus Tectonic Signals of River Profiles to Test Orogenic Models for the Eastern and Southeastern Carpathians

Fluvial morphology is affected by a wide range of forcing factors, which can be external, such as faulting and changes in climate, or internal, such as variations in rock hardness or degree of fracturing. It is a challenge to separate internal and external forcing factors when they are co‐located or occur coevally. Failure to account for both factors leads to potential misinterpretations. For example, steepening of channel network due to lithologic contrasts could be misinterpreted to be a function of increased tectonic displacements. These misinterpretations are enhanced over large areas, where landscape properties needed to calculate channel steepness (e.g., channel concavity) can vary significantly in space. In this study, we investigate relative channel steepness over the Eastern Carpathians, where it has been proposed that active rock uplift in the Southeastern Carpathians (SEC) gives way N‐ and NW‐wards to ca. 8 Myrs of post‐orogenic quiescence. We develop a technique to quantify relative channel steepness, the relative steepness index, based on a wide range of concavities, and show that the main signal shows an increase in relative steepness index from east to west across the range. Rock hardness measurements and geological studies suggest this difference is driven by lithology. When we isolate channel steepness by lithology to test for ongoing rock uplift along the range, we find steeper channels in the south of the study area compared to the same units in the North. This supports interpretations from longer timescale geological data that active rock uplift is fastest in the southern SEC.

Studies aiming to link topography with tectonics have focused on the main erosive engine of non-glaciated landscapes: the river system (e.g., Ahnert, 1970;Goren, 2016;Hack, 1960;Kirby & Whipple, 2012;Schoenbohm et al., 2004;Seagren & Schoenbohm, 2019;Willett et al., 2014). Amongst quantitative tools developed to describe fluvial morphology, channel steepness, or its normalized equivalent integrating discharge, has been perhaps most widely used. With the reasonable assumption that surface motions directly alter the gradient of channel networks, contrasts in steepness have been interpreted as direct (steepening at fault contacts) or indirect (transient migration of steepening) signs of tectonic activity (e.g., Kirby & Whipple, 2012). However, a variety of different forcings can affect channel steepness resulting in similar morphological expressions (Whipple et al., 2013, and references therein). As one example among many, climatic forcing can affect channel steepness by controlling river discharge via runoff (e.g., Adams et al., 2020). The coexistence of different forcings affecting channel morphology in a similar way raises potential ambiguity regarding the nature of some signals. In this contribution, we focus on the role of lithology on channel steepness in tectonically active landscapes. Where softer rocks give way downstream to harder rocks, a steadily eroding channel will steepen (e.g., Bernard et al., 2019;Forte et al., 2016;Perne et al., 2017;Yanites et al., 2017). Critically, fault displacements commonly juxtapose different rock types, resulting in uncertainty about whether different channel steepnesses on either side of a fault are a function of different uplift rates, rock strength, or both. This common feature of geologically heterogeneous landscapes generates mixed signals in the river network, resulting in ambiguity in interpreting the main forcing controlling the steepening (e.g., Strong et al., 2019).
Here, we attempt to disentangle the role of tectonics from the role of lithology in a tectonically active and lithologically heterogeneous landscape. We focus on the Eastern and Southeastern Carpathians (SEC), where extracting the spatial distribution of active tectonic motions from river profiles is confounded by lithologic contrasts. We use a combination of (a) topographic analysis to extract channel steepness from digital elevation models (DEMs) and (b) field observations and measurements to constrain rock strength for the main lithologies. We then trace lithological units laterally from regions where active tectonics are thought to play a role, northward to where the range has been inactive for several millions of years. Through this exercise, we isolate the signal of active rock uplift on the river profiles from the role of lithology, and hence test tectonic models for the region.

River Long Profiles
Scaling between channel steepness and discharge, or its proxy drainage area, has been qualitatively suggested and observed for over a century: "In general we may say that, if all else is equal, declivity bears an inverse relation to quantity of water" (p. 114 of Gilbert (1877)). In the mid-1950s, Hack (1957) and Morisawa (1962) quantified this qualitative observation, describing a systematic relationship between drainage area and channel gradient. These studies led to the formulation by Morisawa (1962) and later Flint (1974) of a power law describing the commonly observed decrease of channel gradient with increasing drainage area: where S is the river gradient (  dz S dx where z is the elevation and x the distance along the channel); s k the steepness index representing the overall gradient of a river system, a single river or one of its reaches; A the drainage area; and  the concavity index dictating the rate at which channel gradient declines downstream. In order to compare different rivers over one or several networks,  is commonly fixed to a reference value, frequently denoted  ref , in order to extract comparable steepness index values (i.e., normalized to the same value of the concavity index). s k is then referred to as sn k , the normalized channel steepness.
Calculating s k (or sn k ) and determining  (or  ref ) has been traditionally done by applying linear regressions of  ( ) ( ) log S log A plots, where the gradient is - and the intercept s k (e.g., Flint, 1974;Kirby & Whipple, 2012;. However, slope-area plots suffer from significant limitations, mainly linked to the inherently noisy nature of channel gradient derived from DEMs (e.g., . It requires the use of averaging methods (e.g., binning by drainage area and averaging the slope) to exploit the data, inevitably resulting in data loss. An alternative method has been developed to mitigate the effects of GAILLETON ET AL.  Royden et al., 2000). This consists of integrating Equation 1 over the distance of the channel: where b x is the local base level chosen for the analysis (e.g., a basin outlet or fixed elevation; we refer to Forte and Whipple (2018) for a sensitivity analysis on base level choice) and 0 A , a reference drainage area, which is introduced to ensure the units of  do not depend on  . From this equation, Royden et al. (2000) defined a longitudinal coordinate  as: Any point of the channel can be defined using  such as: The  approach normalizes the river profile to a  ref and provides an alternative method to explore S-A relationships. If 0 A is set to a value of unity in Equation 3, then the gradient of  -elevation is equal to s k (e.g., .  has been widely used in various geomorphological studies linking channel morphology to surface processes, to investigate the evolution of drainage divides (e.g., Fan et al., 2018;Forte & Whipple, 2018;Giachetta & Willett, 2018;Guerit et al., 2018;Seagren & Schoenbohm, 2019;Willett et al., 2014) or to derive topographic metrics to describe river networks (e.g., Gailleton et al., 2019;Hergarten et al., 2016;Wang et al., 2017).

Channel Steepness, Tectonics and Lithology
s k has been widely used as a proxy for geomorphological processes. Compilations of detrital cosmogenic nuclide concentrations (e.g., Be 10 ), used to quantify average erosion rates for a given river catchment area (e.g., Bierman & Steig, 1996;Lal, 1991), have demonstrated a direct positive correlation between erosion rate and s k (e.g., Codilean et al., 2018;DiBiase et al., 2010;Harel et al., 2016;Kirby & Whipple, 2012;Mandal et al., 2015;Scherler et al., 2014). This is a direct quantification of early hypotheses that steeper channels should tend to erode more rapidly (e.g., de Lapparent, 1907;Gilbert, 1877). Changes in erosion rates can result from tectonic or climatic forcings, enabling the use of s k to study tectonic or climatic evolution over large areas.
In tectonically active landscapes, changes in s k have been interpreted as a proxy for differential tectonic activity. Wobus, Whipple, and Hodges (2006) linked a sharp increase in channel steepness of the Marsyandi River as it crossed the region of the Main Central Thrust of the central Himalaya to a rock uplift signal related to the tectonic structure, using other proxies of erosion rates to support this hypothesis. This relationship between rock uplift and s k has been thoroughly explored in a range of settings (e.g., Lavé & Avouac, 2001;Seagren & Schoenbohm, 2019;. Previous studies using both topographic data (e.g., Kirby & Whipple, 2012) and numerical models (e.g., Eizenhöfer et al., 2019) have highlighted potential explanations for large breaks in channel steepness. In both these studies, concentrated relative uplift could be caused by deep structures (e.g., midcrustal ramps) under the mountain belt. s k has also been interpreted as an indirect expression of base-level change resulting from tectonics (e.g., Hurst et al., 2019;Ouimet et al., 2009;Steer et al., 2019;  As tectonics, climate and stream piracy can affect channel steepness by inducing external forcings to the river channels, intrinsic forcings (e.g., fractures, weathering, lithology) will also affect s k . Amongst these intrinsic forcings, the effect of differential lithology on fluvial morphology has been a recent focus of geomorphological studies (e.g., Bernard et al., 2019;Campforts et al., 2019;Forte et al., 2016;Kirby et al., 2003;Peifer Bezerra, 2018;Seagren & Schoenbohm, 2019;Strong et al., 2019;Thaler & Covington, 2016;Yanites et al., 2017). Rivers flowing over harder rocks tend to have steeper channels and affect the overall landscape morphology (e.g., Forte et al., 2016;Tucker & Slingerland, 1996;Yanites et al., 2017). This effect is linked to the fact that harder lithologies are more difficult to erode, forcing the channel to steepen to maintain a constant erosion rate.
Studies of entire mountain ranges (e.g., Bernard et al., 2019;Duvall, 2004;Gabet, 2019) have demonstrated the important effect of lithology on channel steepness in syn-to post-orogenic settings, with a positive correlation between sn k and rock strength appearing to be the controlling forcing on landscape morphology in non-glaciated areas. Careful acknowledgment of lithological heterogeneities still permits the interpretation of climatic and tectonic signals from river morphology (e.g., Campforts et al., 2019;Kirby et al., 2003), but can also confuse the signal (e.g., Strong et al., 2019) and potentially lead to misinterpretation. In this study, we focus on cases where contrasts in the erodibility of rock are co-located with possible contrasts in rock uplift. In that case, the origin of channel steepening remains difficult to interpret.

The Orogenic and Geomorphological Evolution of the Eastern and Southeastern Carpathians
The Carpathians are an arcuate mountain range located in the eastern continuation of the Alpine orogenic belt (Figure 1). Previous studies have shown that the overall Carpathian structure formed in response to the Triassic to Tertiary opening and closure of two oceanic realms by subduction and continental collision (details in Csontos & Vörös, 2004;Mațenco, 2017;Săndulescu, 1988;Schmid et al., 2019). In a plate tectonics GAILLETON ET AL.   (2017)). Note the different references used for Prut/ Dniestr, Siret and Focani base-levels. EC = Eastern Carpathians, SEC = South-Eastern Carpathians, NEC = North-Eastern Carpathians, P-T = Post-Tectonic cover (sensu post Late Miocene Collision), CF = Convolute Flysches and C-S = Ceahlău-Severin. Note the post-tectonic cover is not displayed on this figure for clarity purposes, as it would cover the whole map, but is displayed in the rest of the manuscript's figures with the blue color. The main frontal thrust is displayed in black where reaching the surface and gray where buried below the sediments of the Focani basin. scenario, the studied area of the Eastern and SEC is made up by two basement-bearing continental mega-units in an upper plate position, the European (sensu largo) continental foreland in a lower plate position, and a thin-skinned thrust and fold belt deformed at or near their subduction contact (Figures 1 and 2). The European foreland is furthermore overlain by a foredeep that locally reaches 13 km in depth in the area of the Focani Basin ( Figure 2, Tărăpoancă et al., 2003).

Tectonic Evolution
The Middle Jurassic opening of the Alpine Tethys was followed by the Cretaceous-Miocene closure of its Pienides-Magura and Ceahlău-Severin branches (Figure 1, Plasȋenka, 2018;Săndulescu, 1988;Schmid et al., 2008). The closure scraped off sediments deposited over the subducting ocean and its eastern passive continental margin by forming a thin-skinned system of thrust sheets, grouped in nappes emplaced in a foreland-breaking sequence from the Cretaceous (Ceahlău), late Oligocene to Early Miocene (Convolute Flysch, Audia/Macla), middle Miocene (Tarcau, Marginal Folds), to late middle Miocene to Early late Miocene (Subcarpathian) times (Figures 1 and 2). The thin-skinned deformation took place until around 9-8 Ma when the main crustal subduction zone was locked by the continental collision (Mațenco, 2017;Schmid et al., 2008, and references therein). Low temperature thermochronology studies, primarily apatite fission tracks and apatite U-Th/He, have shown that the thin-skinned accretion was associated with gradual exhumation. Exhumation of up to 6 km took place at average rates of below 1 mm/yr and peaked between 13 and 8 Ma during the Miocene collision (Gröger et al., 2008;Merten et al., 2010;Necea, 2010;Necea et al., 2021;Sanders et al., 1999). The exhumation was spatially distributed throughout the thin-skinned nappes with higher values in their center (around the Tarcau and Marginal Folds nappes in Figure 2). Similar exhumation rates were also interpreted in the northern part of the Eastern Carpathians (EC) during two periods of exhumation, one more rapid between 12 and 5 Ma and another after 5 Ma. In this area, the exhumation history is interpreted to be driven by the erosion of a thickened wedge after the cessation of shortening at 12-11 Ma, associated either with slab break-off or with the end of subduction (Andreucci et al., 2015).
While tectonic activity remained minor elsewhere, a further deformation episode took place after 8 Ma in the area of the SEC. The formation of high-angle thick-skinned reverse faults truncating both the basement and the overlying thin-skinned thrust belt at depth created a crustal root presently located beneath GAILLETON ET AL.   Necea (2010). Note that (a) and (b) share the same x-axis as distance along the cross-section. the external parts of the thrust belt (Figure 2), as proven by seismic, gravity and magnetic studies (e.g., Bocin et al., 2005Bocin et al., , 2009Hauser et al., 2007). This deformation was associated with gradually accelerating exhumation at values between 1.5-5 mm/yr in the external part of the orogenic wedge, located above the thick-skinned reverse faults Necea, 2010). This presently active deformation was also coeval with subsidence in the foreland at values of 1-3 mm/yr, which created the overall synclinal geometry of the Focani Basin ( Figure 2, Leever et al., 2006;Mațenco et al., 2007;Tărăpoancă et al., 2003). It was also coeval with smaller amounts of subsidence in the order of hundreds of meters, creating the shallow Braov and Târgu Secuiesc intramontane basins, which covered most of the internal part of the orogenic wedge and its Dacia basement (Figure 1). These differential vertical motions are thought to be related to an asthenospheric circuit driven by the sinking Vrancea slab, still (barely) attached to the overlying lithosphere in the final stages of slab detachment (Ismail-Zadeh et al., 2012;Martin & Wenzel, 2006;Mațenco et al., 2016). The post-8 Ma tectonic structures of the SEC, deformation along thick skinned reverse faults and the larger underlying mantle circuit, are presently active, as demonstrated by the large intermediate mantle (70-220 km) seismicity of the Vrancea slab, the moderate seismicity of the overlying crust (Bocin et al., 2009;Ismail-Zadeh et al., 2012;Oncescu & Bonjer, 1997;Radulian et al., 2000), and GPS movements reaching up to 7 mm/yr (Schmitt et al., 2007;van der Hoeven et al., 2005), together with interpretations from studies of the mantle structure, anisotropy and attenuation (Bokelmann & Rodler, 2014;Ivan, 2007;Martin & Wenzel, 2006;Popa et al., 2005Popa et al., , 2008Popa et al., , 2008Russo et al., 2005).

Lithology and Geomorphology
The Eastern and SEC show a large diversity of mostly clastic, but also carbonatic lithologies across the orogenic strike, which maintains a remarkable continuity in the same tectonic units over hundreds of kilometers along its strike. The Cretaceous -Paleogene sedimentation is generally dominated by a deep-water mixture between pelagic and dominantly turbiditic ("flysch") sedimentation, with shallower shelf to alluvial coarse sediments deposited in forearc basins over the accretionary wedge during peak tectonic moments (such as the Albian Ceahlău conglomerates), well described in numerous regional or local studies (e.g., Belayouni et al., 2009;Melinte-Dobrinescu et al., 2008;Miclău et al., 2009;Olariu et al., 2014;Săndulescu, Krautner, et al., 1981;Săndulescu, tefănescu, et al., 1981;Roban et al., 2017). A gradual transition toward a regressive basin fill ("molasse") and coarser deposition took place during the Miocene continental collision in the more external Marginal Folds and Subcarpathian nappes, while the foredeep contains a middle Miocene -Pleistocene transition from shallow-water marine and lacustrine sedimentation dominated by an orbitally forced cyclicity to a deltaic and alluvial continental sedimentation (e.g., Jipa & Olariu, 2013;Săndulescu, tefănescu, et al., 1981;Stoica et al., 2013;Vasiliev et al., 2004).
Geomorphological studies available in the Eastern and SEC (Rădoane et al., 2017 and references therein) are in general agreement with the tectonic scenario described above. These studies have inferred that the EC have a general topography that mirrors the decay of an older (Miocene) orogenic buildup, with longitudinal river profiles trending toward an equilibrium, and sediments generated dominantly by river channel erosion. In contrast, the SEC have a young and actively changing topography, shown by a significant disequilibrium in longitudinal river profiles, sediments generated dominantly by recycling landslides, rapid uplift observed in geomorphic markers such as terraces, migration of knickpoints, water divides, and possible piracy events derived from  profiles (see also Bălteanu et al., 2010;Cristea, 2014Cristea, , 2015Leever, 2007;Necea et al., 2005Necea et al., , 2013Rădoane et al., 2003;ter Borgh, 2013). These studies also suggested that recent tectonics may have shifted the presently observed main water divide separating rivers draining to the European foreland from those draining to the Transylvanian hinterland and the middle of the thin-skinned wedge in the central part of the SEC (compare maps in Figure 1). Furthermore, the tectonically induced differential vertical movements may have triggered a general drainage re-organization with rivers being deflected toward the center of the Focani Basin (Fielitz & Seghedi, 2005 and references therein). While all these indications point toward a differentiation in the Eastern and SEC between the erosion of older tectonic relief and a topography controlled by active tectonics, respectively, the mechanisms responsible for the significant variability observed locally are less understood. For instance, structural and geomorphological studies have suggested that the Pleistocene to recent uplift of the SEC has migrated eastwards through time toward the Focani Basin ( Figure 2, Molin et al., 2012;Necea et al., 2005Necea et al., , 2013, qualitatively interpreted as an effect of the Vrancea slab steepening and retreating in the same direction (e.g., Mațenco et al., 2007).
On this first order pattern, the locally observed influence of lithological strength contrasts on the surface morphology and heterogeneities in normalized channel steepness (Cristea, 2015;Rădoane et al., 2017) still has to be quantified.
The potential influence of climate on fluvial morphology in the study area was assessed by Necea et al. (2013), particularly for river incision during the Pleistocene. During the early Pleistocene, the Black Sea experienced sea level drop of up to 150 m (Winguth et al., 2000). However, the river network in the SEC was only indirectly connected to the Black Sea via the Danube river, crossing different faults and basins, and Necea et al. (2013) demonstrated that the wavelength of incision recorded by the terrace system is not correlated to the sea level drop of the Black Sea. Furthermore, during the Middle Pleistocene, local glaciation occurred, especially toward the internal Carpathians where NW to WNW winds induced a wetter climate compared to the SE Carpathians. However, most of the glaciation occurred in the Southern Carpathians (e.g., the Fagaras) and to a lesser extent in the Northern Carpathians, neither of which are in our study areas.
In summary, all previous studies have suggested that the fluvial morphology is controlled by local and regional tectonics modulated by lithological variations. We build on these studies by applying our fluvial geomorphometry and channel steepness analysis at the scale of the entire Eastern and SEC for rivers draining into the European hinterland. Furthermore, we explore the consistency of channel steepness variations across ranges of concavity indices constrained in the field area. We delimit the area into three regions controlled by different base levels ( Figure 1): (a) the Focani Basin area, which aggregates rivers draining into the SEC foreland basin, (b) the Siret base level, aggregating rivers into the foreland basin along the entire chain, and (c) the Prut base level and the associated drainage system, which is used as a reference area located far into the European foreland that is, not directly linked to Carpathians mountain building processes. Our analysis specifically excludes the southern-most termination of the SEC (the Ialomita catchment) with a Danube river base level (Figure 1), as this is affected by significant strike-slip to transpressive deformation and recent salt diapirism (Mațenco & Bertotti, 2000). In the same area, our analysis also excludes the comparatively smaller internal part of the orogenic wedge that drains into the Transylvanian hinterland.

Rock Strength
We apply a semi-qualitative approach to estimate rock strength. First, the extent of the tecto-lithologic units is estimated using the compilation of 1:50,000, 1:200,000 and 1:500,000 geological maps (published by the Geological Institute of Romania) and the studies of Mațenco et al. (2010) and Mațenco (2017). The Ukrainian section of the map has been completed and extrapolated using the extent of tectonic units in Andreucci et al. (2015), with some spatial approximations and unit grouping to match nomenclature in the different datasets. We also acknowledge that local lithostratigraphic variation can occur within each tectonic unit, however we prioritize N-S continuity in our unit definition in order to track changes along the mountain range over adding localized sub-units. We take account of potential internal changes in rock strength within each unit using Mațenco and Bertotti (2000), which compiles detailed local stratigraphic information (e.g., Joja et al., 1968;Săndulescu, 1984) and allows us to determine if potential local anomalies are due to internal lithologic contrasts. The chosen grouping allows us to (a) follow the continuous northward evolution of channel steepness along similar units, and (b) encompass large-scale signals.
We then measure the uniaxial compressive strength of the rock through the study area. Schmidt hammer measurements were carried out in the field on rock outcrops, where we focused on fresh rock surfaces. The Schmidt hammer, type N in this study, records a "rebound value" between 10 and 100 where higher values denote high elastic strength of the rock. We also record the outcrops where the rock was too soft to be tested, that is, where the Schmidt hammer did not encounter enough resistance from the rock to return a measurement. The rebound value can be converted to compressive strength using a chart provided and calibrated with the equipment used in the field.

Digital Elevation Model, Preprocessing, River Network and Climate
We use the publicly available ALOS World 3D 30 (AW3D30) meter resolution topographic data set for the study (Tadono et al., 2016). It has been shown to better capture accurate channel elevations than 30 m SRTM data and, in some cases, 12 m TanDEM-X topographic data (Boulton & Stokes, 2018;Mudd, 2020;Schwanghart & Scherler, 2017).
The raw DEM has some internal depressions, which spuriously stop flow routing in the DEM and therefore break the drainage area accumulation. Different solutions to filling such depressions exist, but we chose to use a carving algorithm (Lindsay, 2016). Filling algorithms tend to affect an area upstream of numerical dams or depressions, and we wish to minimize the number of pixels affected by pre-processing.
However, a preliminary step is required as AW3D30 contains a small number of pit artifacts. These can be tens of meters deep and, based on inspection of satellite imagery, appear to be correlated with reflective surfaces (the AW3D30 data set is generated from multispectral imagery). Although their area is small enough to not significantly affect the river extraction, these artifacts affect the carving algorithm by forcing unrealistic trenches to drain them. We therefore use a localized filling algorithm on these pits prior to the carving to minimize DEM corrections while ensuring realistic flow routing. Details about the process are available in the supporting materials.
Drainage area and flow direction is extracted using a D8 algorithm (O'Callaghan & Mark, 1984), and we extract the channel network using a drainage area threshold of 450,000 2 m for all basins draining to the topographic mountain front in the study area (Romanian South-Eastern and EC). This threshold has been selected to reach a drainage density (of the extracted channel network) sufficient to provide continuous sn k data through the field area while excluding low-drainage area rivers potentially affected by different processes (Whipple et al., 2013, and references therein).
To check the potential influence of climatic effects, we calculate a separate data set with a drainage area weighted by precipitation data. We apply the methodology of Babault et al. (2018) and Harries et al. (2020), using Fick and Hijmans (2017)'s data set for precipitation patterns. Because our study focuses on disentangling the signals of tectonics from lithology, we detail this process in the supporting materials.

k sn Extraction
As shown in Section 2.1, sn k can be represented as the gradient of  -elevation profiles (  dz d ). To calculate these, we must first make some decisions about how to calculate the  coordinate: the choice of b x (i.e., the point at which the numerical integration starts), reference drainage area ( 0 A ) and the reference concavity of the overall river network ( ref ). We set  0 1 A so that the gradient in  -elevation space is equal to sn k . As demonstrated by Forte and Whipple (2018), the choice of b x affects the value of  , and is important to justify when interpreting the absolute value of  . However, we only use the gradient  dz d in this contribution, the value of which is not affected by the choice of b x . In our case, b x only sets the geographical outlet of the drainage network extracted. We therefore arbitrarily fix b x at the topographic mountain front in order to analyze all rivers draining the eastern foreland basins.

Concavity Index
We take particular care when selecting the concavity index, as only sn k values extracted with a same reference concavity ( ref ) can be relevantly compared. Following Mudd et al. (2018); Niemann et al. (2001); Perron and  and Wobus, Crosby, and Whipple (2006), if the correct concavity index is selected, tributaries and the main stem channel should be co-linear, even in transient landscapes. We use a set of algorithms described in Mudd et al. (2018) and Hergarten et al. (2016), aiming to maximize the co-linearity of  -elevation space for each watershed, which is then selected as the most likely value of  ref for that watershed. Uncertainty around that best fit is also calculated by calculating best fit for sub-sets of connected rivers within each watershed

. Segmentation of  -Elevation Profiles
Once  ref has been determined, sn k can be calculated using the gradient of elevation as a function of  . Direct, pixel-by-pixel determination of sn k is sensitive to inherent DEM noise and would require the use of some form of post-processing (e.g., a moving average window) to exploit the results. Such a method would smooth over discontinuities such as knickpoints. Instead, we employ the algorithm described in Mudd et al. (2014), which applies a statistical method to select the most likely combination of linear segments in  -elevation. These linear segments are predicted by the theoretical work of Royden and Perron (2013).
The Mudd et al. (2014) algorithm first selects a user-defined number of adjacent river nodes, referred to as tg n . The algorithm calculates all the combinations of segments composed of a minimum amount of nodes and calculates best-fit metrics for each combination of segments. A good fit to the data is balanced against too high a number of segments (i.e., over fitting) using the Akaike information criterion (Akaike, 1974). Each segment describes a section of river profile as: where   sn M k if  has been calculated with  0 1 A , and  b represents the intercept of each segment. To make sure that small-scale noise does not affect the results, the algorithm repeats this segmentation a user-defined amount of times following a Monte Carlo scheme. At each iteration, one node every sk n is randomly skipped in order to select different subsets of nodes in the Monte Carlo process. The sn k value for each node is the mean value of all the segment slopes involved in the calculation. 1, but others need to be carefully justified as their choice will affect the segmentation process. The size of the segments is a particularly important factor to consider: It will determine the scale represented by sn k variations extracted with the algorithm. The segment size is determined by the number of nodes targeted by each algorithm iteration ( tg n ) and the number of nodes skipped at each Monte Carlo iteration ( sk n ). Higher values for these parameters will tend to generate larger segments, thereby averaging longer river reaches, whereas smaller values will generate smaller segments representing small-scale features. The effects of varying these parameters have been explored in detail by Gailleton et al. (2019).

Relative Steepness Index
As shown in the previous sections, calculating sn k depends on a number of parameters which affect (a) the absolute value of sn k and (b) the scale it represents via the relative size of segments in the profiles. Two populations of sn k , for example, from different watersheds, are directly comparable only if the metric has been calculated with the same parameters (e.g., Hurst et al., 2019;Kirby & Whipple, 2012).
Different values of  , for example, will generate different orders of magnitude of sn k . Large areas, such as entire mountain ranges, will naturally have spatial variations in the concavity index (e.g., Chen et al., 2019;Seagren & Schoenbohm, 2019). In this study, we propose circumventing this limitation by (a) calculating sn k for a wide range of parameters in order to ensure the systematicity of sn k contrasts despite heterogeneity in both  ref and scales and (b) comparing cross-parameter results with a relative steepness index.
The relative steepness index aims to provide a non-dimensional way to compare channel steepness. We chose a statistical approach to allow us to assess which parts of the landscape are showing the steepest sn k independently from its absolute magnitude. To calculate a relative steepness index, we use a statistical metric called the modified z-score (Iglewicz & Hoaglin, 1993), which we denote with i M . i M represents the statistical distribution of a population and allows us to quantify how it varies in space. In statistics, the z-score of a value within a population represents how far the value is from the mean in number of standard deviations. It therefore suits cases where the population can be described as normally distributed around the mean. The modified z-score is a nonparametric version of the z-score, based on the median rather than the mean, and suits our data set better as sn k values are not expected to be normally distributed-particularly in a transient and heterogeneous environment where many locally steepened reaches form groups of outliers.
where , i j M is the modified z-score for pixel i and parameter value combination j. Note that we test a range of  ref for the entire area. The constant 0.6745 is introduced to follow the commonly used standardization of Iglewicz and Hoaglin (1993). Each pixel has a relative steepness index for a given parameter combination is equivalent to the median and higher and lower values denote respectively higher and lower samples compared to the overall population. This method is widely used to detect outliers in large datasets (e.g., Giustacchini et al., 2017). Because all values of , i j M are normalized to the median values and MADs of each parameter value combination, we can use these to compare relative steepness indices amongst sn k data with different parameter values. We therefore refer to the , i j M data as the "relative channel steepness" in all our figures, with values greater than zero representing parts of the channel network that have steepness greater than the median, and values less than zero representing parts of the channel network that are gentler than the median sn k values.
Because of natural spatial heterogeneity in  ref and other parameters, the relative steepness index method can theoretically also be affected by distortions due to not-optimal parameterization (e.g., for sets of parameters j with non optimal  ). However, there is no perfect solutions to bypass this limitation because this approach requires the whole landscape to be considered for each j. An alternative method, potentially safer, would be to only compare sub-basins displaying similar concavities. However, this would only be possible if there are enough rivers, for a single optimal concavity, outcropping everywhere in the field area for each tecto-lithologic unit. This is not the case in our study area. We suggest, however, that the relative steepness index method displays signals that are consistent through the different sets of parameters j, and therefore highlight the signals that can be safely interpreted. We explore and illustrate this aspect in the supporting material.
Another possibility to achieve comparability across ranges of  ref is to adjust the gradient of  -z profiles by changing 0 A in Equation 3 until one reaches similar order of magnitudes in  and therefore sn k . However, the relationship between the gradient of  -z and sn k would be obscured and the choice of 0 A for each  ref more difficult to justify. We therefore opted for a more common and standardized statistical tool to base our method on, as the modified z-score has been commonly utilized in a wide range of applications (e.g., Giustacchini et al., 2017).

Rock Strength
We collected a total of 347 rock strength measurements across the tectonic units in the SEC. The results are quantified using two different metrics: (a) the rebound values (medians and quartiles for each tectonic unit excluding the non-responsive data points) and (b) the proportion of non-responsive measurements for each tectonic unit (Figure 3).

Journal of Geophysical Research: Earth Surface
These results are consistent with qualitative field observations. The first group shows more resistant lithofacies and crops out more frequently in the landscape than the second, which shows fewer, thinner and sparser resistant layers. The results show several trends: Across all studied basins, we find that the concavity indices have a median of  0.35 0.10 (red square in Figure 4) for our study area. In the SEC (basins with northing values ranging from 5,000 to 5,100 km, see Figure 1), the range of values is narrower than in the EC (basins with northing values greater than 5,100 km). The smaller basins within the South-Eastern Carpathians, mainly draining the frontal units (Focani Basin), tend to show higher concavity indices than larger basins. Concavity indices in the EC are more heterogeneous than in other parts of the study area. On the basis of these data, we chose the range 0.2-0.6 for investigating the relative distribution of sn k through our landscape, as it includes all the most likely values in individual basins (excluding two outliers) and most of the interquartile values ( Figure 4).

Relative Channel Steepness
We calculated   Figure 1). The data points represent the median rebounds values, and the error bars the first and third quartiles, respectively. The proportion of non-responsive points is also displayed, as an indirect proxy for the proportion of weak rocks within each unit.
In addition, we tested the potential effect of climate forcing on the relative steepness index in our field area to ensure its effect can be safely ignored. The results for the latter are available in the supporting material. Figure 5 shows the relative steepness index as a function of the northing coordinate. This provides an overview of channel steepness in regards to the different areas of differential tectonics suggested in Section 3. The data is noisy and does not show an obvious N-S trend. There is a sharp increase in the relative steepness index between northing values of 5,000 km and 5,030 km, which may be linked to the bending of the mountain range and incorporation of a few unrepresentative data points in the extreme south of the Buzau watershed (Figure 1). Three regions host steep channels compared to the rest of the landscape: (a) the Focani Basin area (northing 5,000-5,080 kilometers, HS1 in Figures 5 and 6), (b) in the heart of the EC (northing 5,125-5,240 kilometers, HS2 in Figures 5 and 6, note that in Figure 5, it is expressed by a larger spread of relative steepness), and (c) a less prominent steep area in the Northeastern Carpathians from 5,340 km. These three areas are connected by two regions of lower relative steepness index.

Regional Distribution of Relative Steepness Index
The absence of a monotonic N-S trend is also expressed in a map view (Figure 6) where the medians of all the relative steepness indices suggest a compartmentalized data set. A clear N-S mid-range linear feature sharply separates an eastern region of lower steepness and a western region of higher steepness. This main break in relative steepness index is labeled Main Break in Steepness (MBiS) in Figure 6. The sharpness of the contact is less clear south of 5,160 km. Other less clearly expressed trends can be observed with this map view: (a) Within the western region of high steepness, high patches stand out, particularly at kilometers 5,030 (HS1) and 5,130. (b) Within that same region, localized patches of low values express the presence of high-elevation low-gradient (HELG) valleys in the Buzau, Trotus, Bistrita and Prut watersheds (labeled HELG in Figure 6).    (Figure 1) (c) followed by a sharp decrease until kilometer 5,300 (i.e., the northern part of the Siret base level) and (d) high values in the northernmost area, linked to Prut and Dniestr base level, north to km 5,300. Finally, the basement rocks of the Dacia units locally impose patches of high relative steepness index in the EC where these rocks are largely exposed. Figure 7 also highlights multiple notable behaviors differing from a northward monotonic decay as one moves away from the active vertical motions of the SEC. (a) Although the Subcarpathian nappe has its highest values in the Focani area, it also displays a local peak north of km 5,200, denoting a greater proportion of steeper channels within the Subcarpathian nappe in the area. Note that the exposed surface of this unit decreases northward (Figure 5

Spurious Tectonic Signals
A prominent break in relative steepness index can be seen in Figure 6 (labeled MBiS). It extends along the entire N-S axis of the study area, and is located east of the western drainage divide of the watersheds in the study area. Section 2.2 highlighted tectonics as a common forcing generating similar features. In the Carpathians, recent tectonic activity is concentrated in the southeastern bend of the mountain range (see Section 3). The break in relative steepness index observed in Figure 6 extends far beyond the region where deformation is inferred from other independent proxies, and could be used as an argument for extrapolating recent tectonic activity to the North. However, our rock strength data (Figure 3), combined with apparent tectonic inactivity north of the South-Eastern Carpathians (as suggested by other studies, see Section 3), point to lithology as the main driver of the break in relative steepness index. This is concentrated where the evaporite-rich and highly fractured rocks of the Subcarpathians and sandstone-rich Tarcau and Marginal Fold units lie in contact (e.g., Bernard et al., 2019;Yanites et al., 2017). This highlights the danger of extracting channel metrics at large scale without taking local lithological context into account. This contrast is also prominent in more local data set (see Figures 8 and 9).
This line of reasoning also suggests lithology as a control on more local relative steepness index contrasts, for example: (a) The patch of high relative steepness indices at the top of the Bistrita watershed, described in Section 5.3. Its boundaries correspond to the mostly magmatic rocks of the Dacia basement units and the volcanic rocks linked to Neogene volcanism (dark gray vs. gray in Figure 8). (b) Very sharp and significant drop of relative steepness index (Figures 6 and 7) occurs within the Tarcau nappe around northing kilometers 5,200 to 5,250 (see Figures 7  and 8). Local litho-stratigraphic data (Mațenco & Bertotti, 2000) highlights that this also corresponds to a lithological change from coarsegrained resistant sandstones in the Bistrita valley to finer-grained, often shaly turbidites in the Moldova valley (see Figure 6). It additionally overlaps nearly perfectly with the drainage divide between the Bistrita and Moldova watersheds ( Figure 6); this represents another possible expression of lithologic forcing by "pinning" drainage divides on resistant rocks (e.g., Bernard et al., 2019;Seagren & Schoenbohm, 2019). (c) Low steepness values are observed at the highest, westernmost part of the Bistrita watershed (see Figure 8), corresponding to a switch from the resistant metamorphic rocks of the Dacia basement to softer sedimentary rocks belonging to the Transylvanian Basin (Mațenco, 2017). Figures 8 and 9 illustrate how global and local lithologic forcings can generate relative steepness index contrasts which can potentially lead to spurious tectonic interpretations if lithology is ignored. The localized sharp contrast in relative steepness index in Figure 9 is spatially associated with faults. But it is also superimposed on a smoother, larger-scale tectonic gradient toward the SE and it would be easy to confuse the two. The same morphological signature is also displayed in Figure 8, this time normal to the mountain range. However, despite striking similarities in their fluvial expression, this signal is solely lithologic. Quantifying the relative importance of tectonics versus lithology is not straightforward and is discussed in Section 6.4.

Integration of Relative Steepness Index in the Tectonic Model
Knowing that lithology can influence the patterns of relative steepness index, we must then consider strategies for extracting tectonic signals from lithologically complex terrain (see  Within litho-tectonic units at the eastern edge of the range, that is, the whole area eastern to the main break of steepness (MBiS in Figure 6), we find higher values of relative steepness index in the South Eastern Carpathians (HS1 area of Subcarpathians and Post-Tectonic units in Figure 7, also visible in the close-up of the Buzau catchment in Figure 9). This suggests that there is a tectonic signal of increasing rock uplift rates from north to south in the frontal units, consistent with what has been suggested by structural and exhumation studies (Mațenco, 2017, and references therein). This pattern is particularly clear for the Marginal Folds, the Subcarpathians and the Focani basin/Post-Tectonics units, with all showing a monotonic northward decrease in relative steepness index. When looking at steepness patterns over the entire mountain range, the changes in steepness from different lithologies are greater than the N-S trends within the frontal litho-tectonic units, highlighting how tectonic patterns may be masked by lithologic contrasts in rock erodibility.
Previous studies (see Section 3) also suggest an eastward gradient in vertical motions within the SEC by reactivating deep faults that do not reach the surface. Our data is compatible with this interpretation. Given that the Subcarpathians and the post-tectonic cover are both made of soft rocks, a wholly lithologic signal, would place the mountain front at the Marginal Folds/Subcarpathians boundary. A wholly tectonic signal would show steeper channels in the frontal units than the internal units. The sharpness of the break in relative steepness index demonstrates that the most prominent signal is due to lithology, although the GAILLETON ET AL.  somewhat obscured tectonic signal is still expressed in the topography (see Figure 9). An expression of this mixed signal is the 1,000 m high mountain at the front of the Putna valley made of very soft sedimentary rocks, which can only be explained with recent surface uplift in this area.
Patterns of relative steepness index within the Tarcau unit are more ambiguous than the others. Here, the relative steepness index does not show a gradual decrease northward like other units. It sustains higher values northern than other units before a sharp drop. Given the fact that the Tarcau    rocks in the SEC thin-skinned sedimentary succession and also contain significant northward changes of lithology (see Figure 8), we suggest that the lithologic forcing overprints the tectonic one in this unit.

Non Lithologic Low-Gradient Area Within the South-Eastern Carpathians
Although rock hardness measurements in the South-Eastern Carpathians do not suggest significant lithologic contrasts between the Ceahlău-Severin, Audia/Macla and Tarcau units (Figure 3), the upper parts of the Buzau basin show low values of relative steepness index (see Figure 9). We explain this different behavior using local data from these units. (a) Thermochronometers from Merten et al. (2010) have suggested an older and lower magnitude exhumation of these units through time in the South-East Carpathians (in the Buzau watershed), which can be related to long-wavelength of exhumation related to slab-retreat type of processes (e.g., Mațenco, 2017;Picotti & Pazzaglia, 2008). (b) Several authors (Fielitz & Seghedi, 2005;Necea, 2010;ter Borgh, 2013) suggest a drainage reorganization to explain these high elevation low-gradient valleys. Our data set is consistent with these previous observations, showing steep "aggressive" (sensu Willett et al. (2014)) rivers in the Buzau watershed juxtaposed with an upstream low gradient, diffusive landscape. These values are at odds with the regional pattern of tectonic activity, that is, high tectonic activity in the South-Eastern Carpathians and post-collisional decay in the EC, and bias the global distribution of relative steepness index (Figures 5 and 6) by reducing the regional values.
As tectonics is a common driver for drainage divide reorganization (e.g., Giachetta & Willett, 2018;Seagren & Schoenbohm, 2019;Willett et al., 2014), the two can be linked here. Figure 9 summarizes the local signals observed within the Buzau watershed, illustrating the diversity of local expression of relative steepness index.

Disentangling Tectonics From Lithology
In our study area, separating the effect of lithology and tectonics on channel steepness (or its relative equivalent) is complicated by both their heterogeneity and the fact that lithologic contrasts are also co-located with tectonic contrasts (for example, faults often juxtapose rock types). The isolation of each component is further obscured by the variations in concavity index. Such quantitative interpretations require us to make assumptions, and so it is important to emphasize that undertaking this exercise is somewhat exploratory.
First, we select two units of interest: the Marginal Folds and the Subcarpathians (Figure 1). They capture the main features of relative steepness index contrasts: (a) They outcrop in the active area of the SEC with a proposed tectonic gradient toward the foreland and in the inactive North; (b) the Marginal Fold units are composed of notably harder lithologies than the Subcarpathians (Figure 3). They are also exempt from internal lithological changes and stream piracy, unlike the Tarcau unit for example, (Figures 8 and 9). The relative steepness index of both units for each set of parameters j binned by north and south is displayed in Figure 10a.
The northern area has been tectonically inactive for several million years (see Section 3). On this basis we assume that the contrast in channel steepness (and its relative equivalent) between the two units is the expression of the lithological component of the channel steepening. This allows the isolation of three tectonic signals, relative to lithologic effects ( Figure 10). When we turn our attention to the southern area, where there is a greater difference in RSI between the Subcarpathians and the Marginal folds relative to the north, and we attribute this increase to tectonics (Figure 10a, middle column). Both units exhibit higher RSI values in the south. This suggests a regional tectonic signal affecting both lithologic units.
We note that the change in RSI between the north and south is not the same in two units (Figure 10a, right column). This could mean that the two units are uplifting at different rates in the south, or it could mean that the regional tectonic signal results in the channel steepness index that varies by rock unit. Structural evidence suggest the former is more likely, as it points to greater uplift in the Marginal folds (see Section 3 and Figure 2). We can quantify this potential difference in uplift by calculating the ratio between the change in RSI for the two units between the north and south. That is, we calculate RSI between the Marginal Folds and the Subcarpathians units respectively for the North and for the South. We then calculate the ratio between these two RSI values, with the results shown in Figure 10b. This ratio aims to "normalize" the GAILLETON ET AL.

10.1029/2020JF005970
17 of 26 RSI for the lithologic effect, because we assume there is no tectonic contrast in the dormant north. This suggests that the southern relative steepness indices are roughly 1.06 and 1.33 times that of their northern counterparts, confirming steeper channels in the south in response to active tectonics.
The relative steepness index offers a level of abstraction from sn k which allows cross- ref comparison. But this level of abstraction also restricts further interpretations: It is rather difficult to link a relative steepness index to more instinctive inferred landscape properties such as erosion rates. To achieve this link, further assumptions are required. As a general rule, it has been suggested that sn k is linked to erosion with the following formula (e.g., Whipple et al., 2017): where E is the erosion rate,  is a coefficient linked to local factors such as climate, lithology or fractures, and t an exponent linked to thresholds and hydrologic conditions -ultimately controlling the non-linearity of the relation between erosion and channel steepness. Because we lack the data to constrain these, we use Equation 8 within the detachment limited formulation of the stream power law (Howard & Kerby, 1983): where K is the rock erodibility, which we assume to be directly linked to rock strength in our case as climate influence can be neglected (see Supplemental Material), and the n and m exponents represent the erosion process through the sensitivity of erosion to changes in channel slope and the proportion of precipitation going into runoff. In this frame of reference,  ref m n  / , and Equation 9 can be rearranged to be a function of sn k : We therefore apply a methodology similar to the relative steepness index in order to cancel the lithologic effect, with n values ranging from 0.5 to 2 to represent different behaviors (e.g., Tucker & Slingerland, 1996). In order to extract "lithologically agnostic" erosion rates we manipulate ratios of erosion rates inferred using Equation 10, which is similar to an approach used in other studies (Resentini et al., 2017). The main problem we are facing is the unknown value of K, the erodibility factor. We apply two methods to assess erodibility contrasts, which bypasses the need to constrain K values directly. First we calculate ratios of inferred erosion rates within a single lithologic unit, assuming rock strength does not vary within a rock unit (see Equation (11) for example). Second, if for an area we can assume spatially constant erosion rates across multiple tecto-lithologic units, for example, in a tectonically dormant part of the range, we assume a ratio in n sn k represents a ratio of the K values amongst different rock units.
In our study area, we calculate several ratios (a) between the Marginal Folds and Subcarpathians units in the northern area to isolate the lithologic component between the two units and normalize the ratio of erosion rates in the southern area, and (b) for each of these units between the North and the South to directly infer the increase in erosion rates. This approach is similar to our approach using RSI, but in this case our inferred erosion rates can be directly related to tectonics.
First, we can explore how the inferred erosion rates differ in the north and south for a given lithologic unit. In this case we assume the erodibility for a given unit, U K (the subscript U denotes a given tecto-lithologic unit) is the same in the north and south. We can then then calculate the ratio in the inferred erosion rates between the north and south for a given unit ( , where the subscript U is for a unit, the subscript S is for the southern area and the subscript N is for the northern area. From Equation 11 we see that the ratio of erosion rates between north and south for a given unit reduces to the ratio in their channel steepness index ( sn k ) values raised to the power n.  Figure 11. In both these units the southern area is always inferred to be eroding faster, regardless of the choice of n. We can also explore east-west trends in the southern, more tectonically active area. This is more of a challenge because the lithology varies from east to west. We focus again on the Subcarpathians and the Marginal Folded units. We first make an assumption that in the northern, tectonically inactive area, the erosion rates of these two units are the same: ) should be equal to the inverse ratio of the respective channel steepness indices: where the subscript S denotes the southern area. Figure 11 shows these ration, and the violin is generated by taking the median sn k values over a range of parameter values and constructing the ratios from these resulting median sn k values. Combined, these figures show lower erosion rates in the Marginal Folds than in the Subcarpathians in the South (black violins in Figure 11). Although the range of inferred erosion rates is wide, it is consistent with the observation that the tectonic activity is migrating toward the SE in this area (see Section 3). For both the Subcarpathians and the Marginal Folds, the inferred erosion rates are higher in the south. In Figure 11 the ratio of E is the ratio of inferred erosion from North to South, so for n = 2 a ratio of 0.2 reflects and increase in inferred erosion in the south by a factor of 5. The Marginal Folds also have an increase in inferred erosion rates from north to south, but there is less of an increase relative to the Subcarpathians (Figure 11).
We wish to emphasize once more that these are speculative observations relying on many assumptions. The quantification of spatial variations in erosion rates showed in Figure 11 assumes that the detachment-limited model is a good approximation, which is neither systematic nor straightforward to constrain with wide ranges of possible values for m, n and K (e.g., Harel et al., 2016;Lague, 2014). In addition, the detachment-limited model does not explicitly account sediment fluxes while field observations and satellite im-GAILLETON ET AL. for the E-W plot and Equation 11 for the N-S plots (pink are Subcarpathians and red are the Marginal Folded units) for a range of common n values. A ratio of unity in each plot would indicate the quadrants (either N-S or E-W) were inferred to be eroding at the same rate. For E-W plots, ratios lower than unity mean higher inferred erosion rates in the east (that is, within the Subcarpathians) and in the N-S plots ratios lower than unity mean higher erosion rates in the south. agery suggest some sections of the rivers in the Carpathians are potentially transport-limited. Finally, there is no direct measurement of erosion rates to validate or help constrain our results. However, this first-order exploration represents the first encouraging steps toward quantifying the erosion response and paves the way for future studies exploring these topics with landscape evolution models and cosmogenic-derived erosion rates.

Conclusions
Detecting tectonic signals from channel steepness can be challenging due to lithologic heterogeneity; a common feature of mountain ranges that potentially overprints any tectonic signal and potentially alters their morphological expression. Additionally, exploring channel steepness across a wide geographical range will almost inevitably encompass basins with differing concavities, which can cloud interpretation of channel steepness. In this study, we successfully unravel tectonic and lithologic signals from channel steepness in the Eastern and SEC, a range showing different lithologic and tectonic gradients across multiple scales.
We find that the concavity index, which affects normalized steepness values ( sn k ), varies between approximately 0.2 and 0.6 in the Eastern and SEC. Choosing a single reference concavity might result in misleading sn k values. We therefore developed a method for calculating a relative steepness index that can be applied across basins with different concavities using a modified z-score method that takes into account the nonnormal distribution of channel steepness values across all catchments.
The first order values of relative steepness index across the range show a large contrast between the gentle eastern front of the range and steep areas near the drainage divide. The N-S trends in uplift rates proposed by previous studies are not obviously reflected in the relative steepness data at this scale. However, when we group steepness by litho-tectonic units, the tectonic signal is successfully isolated and highlighted (most prominently on the Post-Tectonic cover, Subcarpathians and Marginal Folds units).
We collected rock hardness data across litho-tectonic units and found that the hardness can be broadly grouped into soft (Post-Tectonic Cover and Subcarpathians) and hard units (the others). This grouping is reflected in the relative channel steepness data: For example, the main break in slope in the Carpathians represents the spatial transition between these groups of units.
Separating the relative steepness by litho-tectonic units, a N-S spatial pattern appears. In the units at the mountain front, this pattern is most clear: Relative steepness is highest in the part of the mountain range where thermochronometers have recorded the highest long-term exhumation rates. In addition, steepness data confirms the migration of the surface uplift pattern toward the east, where thermochronometers show unreset ages and cannot be used to estimate exhumation. Without accounting for lithology, this tectonic signal would have been entirely masked by differences in rock hardness. Spatial trends in the harder rocks toward the peaks of the range show more localized patterns: For example, high-elevation low-gradient valleys expressing localized stream piracy and lithologic variations within hard units explaining other less prominent contrasts in relative steepness indices.
Evaluation of variable rock uplift from channel steepness measurements on the scale of an entire mountain range is challenged by the variability in rock strength and the concavity of the channel profile. Through characterization of channel concavities and independent measures of rock strength, it is possible to isolate for the role of tectonics versus lithology.

Data Availability Statement
The Schmidt Hammer Data utilized for this research is available through Gailleton et al. (2020). The DEMs used for this study have been provided by AW3D of the Japan Aerospace Exploration Agency, available through Tadono et al. (2016). The topographic analysis have been processed with lsdtopytools (Gailleton & Mudd, 2021