The Uplift of an Early Stage Collisional Plateau Unraveled by Fluvial Network Analysis and River Longitudinal Profile Inversion: The Case of the Eastern Anatolian Plateau

Orogenic plateaus of continental collision zones exhibit landforms and fluvial networks that retain first‐order information on their topographic evolution and vertical growth. The inversion of river longitudinal profiles allows to reconstruct the base level fall history of plateaus, supporting the study of landscape evolution in the frame of geodynamic models. The Eastern Anatolian Plateau (EAP) of the Arabia‐Eurasia collision zone is a plateau at an early stage of development. It stands at ∼2,000 m, presents endorheic basins, and is drained by three main river networks. Seismic data indicate a thinned lithospheric mantle that explains the late Cenozoic volcanic activity. Despite the number of studies on the EAP uplift, its history is still debated. In this study we investigated the EAP hydrography and topography, and we inverted the longitudinal profile of one of the main rivers: the Arax River. The results describe a high‐standing, low‐relief plateau drained by a hydrography, controlled by active tectonics. Longitudinal profiles and χ‐plots illustrate rivers in disequilibrium with channel steepness increasing downstream. The elevation of marine deposits indicates a surface uplift of ∼2,000 m, ∼500 m of which are of residual topography. This upheaval occurred by two increases: the first one at 10–11 Ma with the opening of a slab window and the arrival of a mantle flow from Arabia and the second one at ∼5 Ma with the continued inflow coupled with isostasy. Our results describe the early stage of collisional plateau development, distinguishing the contribution of deep processes and isostasy to the topographic growth.

of its recent surface uplift history.In this respect, the EAP is exceptionally favorable because, although it is still partially internally drained, portions of it have been integrated into exorheic drainage systems.Fluvial incision, together with active tectonics and volcanism, is the reason why the EAP has a more rugged topography (38% of hillslope with a steepness <5°) than an internally drained plateau like the Tibetan one, where almost 60% of the slopes are <5° and the remaining steeper slopes are associated with recent and active tectonics (Babault &Van der Driessen, 2013).So, the fluvial network of EAP, integrating into part of it, retains first order information on the modes of its fluvial capture as well as on the interplay between climate and tectonic processes, which are crucial to reconstruct its topographic evolution.Moreover, recent methodological advances (Clementucci et al., 2023;Fox et al., 2014Fox et al., , 2020;;Gallen, 2018;Goren et al., 2014;Pazzaglia & Fisher, 2022;Roberts et al., 2012;Smith et al., 2022) that allow to invert river longitudinal profiles to reconstruct the base level fall history can provide information on the rate of uplift that are difficult to obtain by different approaches.This information can be used to integrate the evidence of surface uplift recorded by geologic and geomorphic data with the geodynamic models that explain the impact of deep processes on the topographic growth of the plateau.
We studied the EAP by analyzing the general shape of topography (filtering it at different wavelengths) and the geometry of the hydrography (river longitudinal profiles, χ-plots, drainage divide mobility), and coupling them with geologic data.In particular, using geological maps and DEMs, we derived the elevation of the surface that marks the change from a marine to non-marine depositional environment to get minimum estimates of the mean surface uplift.Finally, we inverted the longitudinal profile of the Arax River (Arax R.), one of the main streams draining a large portion of the plateau to reconstruct its base level fall history.This approach provided, for the first time, information concerning timing and uplift rates of the EAP in the context of a progressive rise of the topography and allowed linking surface processes with geodynamic models.Thus, we propose an example of the evolution of a continent-continent collisional plateau at an early stage of development, evidencing how deep processes together with isostacy contribute to its formation and uplift.

Geological Overview
The Eastern Anatolia Plateau (EAP) is a dome-shaped geomorphic feature of the Arabia-Eurasia collision zone that represents the lateral continuation of the central Anatolian Plateau (Schildgen et al., 2014) and the Iranian Plateau to the west (Ballato et al., 2017).The EAP is characterized by a rolling low-relief topography standing at a mean elevation of ∼2,000 m (Figure 2a).Many volcanoes rise from this surface, some of which reach almost 5,000 m.To the north and to the south, the topography decreases down to the lowlands of the Kura basin, the Euphrates and Tigris Rivers valley and to the Black Sea quite abruptly, as evidenced by the high values of slope (Figure 2b).The plateau holds depressions of different size, the larger of which are the Ararat basin, that is drained by the Arax R., and the internally drained Lake Urmia.Indeed, several of these depressions are endorheic, including the lakes Van and Sevan, both located in areas where the topography reaches the highest elevation (Figure 2).
The high-standing plateau is drained by three main fluvial networks (Figure 2b): the Ҫoruh River flowing to the Black Sea, the Kura-Arax drainage system to the Caspian Sea, and the Euphrates-Tigris drainage system to the Persian Gulf.The hydrography does not show a typical radial pattern, but the rivers integrate into the plateau in a disorganized way, locally following tectonic structures (see Figure S1 in Supporting Information S1).The presence of internally drained basins contributes to consolidate the appearance of a drainage system in a transient state of organization.The fluvial network scarcely incised the plateau except for its external slopes.
The EAP consists of a mosaic of plates that were accreted to the southern Eurasian margin during the progressive closure of different branches of the Neotethys ocean (e.g., Barrier et al., 2018;Okay & Tüysüz, 1999).The closure of the northernmost branch of the Neotethys occurred from the late Cretaceous to the Paleocene-early Eocene and led to the development of the Izmir-Ankara-Erzinçan-Sevan-Akera suture zone (e.g., Barrier et al., 2018;Okay & Tüysüz, 1999).The continental collision between the southern Eurasian Margin (Pontides) and the Tauride-Anatolide-Armenia block was diachronous and oceanic subduction in eastern Anatolia may have continued until the Miocene (Gürer and van Hinsbergen, 2019).This former plate boundary is marked by several ophiolite belts (Rolland et al., 2020 and references therein) extending from western Anatolia to eastern Iran and separates the southern Eurasian margin from the Tauride-Anatolide-Armenia block (Figure 1).Oceanic subduction and continental collision were associated with upper plate deformation and erosional exhumation as documented by structural, stratigraphic and low-temperature thermochronological data from the Eastern Pontides 10.1029/2022TC007737 3 of 31 (Robertson et al., 2013), the Adjara Trialeti fold and thrust belt (Gusmeo et al., 2021) and the lesser Caucasus (Cavazza et al., 2019).Currently, these mountain ranges constitute the northern margin of the Eastern Anatolia Plateau (Figure 2a).
After plate suturing, approximately from the early-middle Eocene, a regional phase of subsidence associated with extensional tectonics and volcanism occurred from central-western Anatolia to southeast Iran (e.g., Okay et al., 2020;Verdel et al., 2011;Vincent et al., 2005).This extensional phase has been interpreted to reflect rollback processes along a north-dipping subduction zone located to the south of the Tauride-Anatolide-Armenia block (i.e., along the southern branch of the Neotethys ocean; e.g., Vincent et al., 2005;Verdel et al., 2011), possi bly in association with slab folding at depth (e.g., Boutoux et al., 2021).Alternatively, this phase of magmatism, at least for the Eastern and Central Pontides, may represent a post-collisional event associated with lithospheric thinning during lithospheric mantle delamination (i.e., the convective removal of thick lithospheric root following the collision of Eurasia with the Tauride-Anatolide-Armenia block (e.g., Göçmengil et al., 2019).Finally, early to late Eocene volcaniclastic and turbidite sediments were also deposited in the Sivas Basin (westernmost sector of Eastern Anatolia), although they appear to reflect deposition in a foredeep of a growing orogen localized along the northern margin of the Anatolide-Tauride block (Legeay et al., 2019 and references therein).
The southern margin of the EAP includes the Bitlis-Pütürge massifs, a tectonic unit that experienced high-pressure-low-temperature metamorphism in the latest Cretaceous.This massif has been interpreted to represent either the vestige of a microplate located to the north of the Arabian plate (Oberhänsli et al., 2013) or the northern margin of the Arabia indenter that experienced metamorphism during the regional phase of late Cretaceous obduction (Şengör & Yılmaz, 1981).
The following closure of the southern branch of the Neotethys ocean along the Bitlis-Zagros suture zone led to the amalgamation of Arabia to the composite southern Eurasian margin (Figure 1).The timing of this recent collision is highly debated with estimates ranging from the latest Eocene-earliest Oligocene (possibly within a southeastward younging trend, see Darin & Umhoefer, 2022 for a recent review) to the early-middle Miocene (e.g., Cavazza et al., 2018;Hüsing et al., 2009;Okay et al., 2010).Late Eocene-Oligocene, deformation processes were spatially uniform across most of the upper plate from the northern margin of the collision zone (Lesser Caucasus and western Greater Caucasus) to the plateau interior, and did not produce significant exhumation (Albino et al., 2014;Forte et al., 2014;Vincent et al., 2016).Moreover, from the latest Eocene, plate convergence was accommodated by oroclinal bending processes that led to the progressive acquisition of the current arcuate geometry of the southern European margin from the Talesh Mountains of western Iran to the Eastern Pontides (Van der Boon et al., 2018).
Interestingly, an early-middle Miocene regional phase of subsidence occurred as documented by shallow-water marine limestones that can be found from SE Iran to NE Turkey (e.g., Okay et al., 2020).The mechanisms responsible for such a widespread subsidence are poorly understood and may involve slab steepening (Bottrill et al., 2012) or deep subduction dynamics (Boutoux et al., 2021).Available ages indicate that the plateau started to be emerged after ∼15 Ma (e.g., Okay et al., 2020 and references therein).From the middle-late Miocene, the margins of the plateau experienced accelerated exhumation, possibly in response to a more advanced stage of the Arabia-Eurasia collision (Ballato et al., 2011;Cavazza et al., 2018Cavazza et al., , 2019;;Darin & Umhoefer, 2022;Gusmeo et al., 2021;Okay et al., 2010).Thermal histories from low-temperature thermochronological data (mostly apatite fission track and apatite U-Th/He) indicate that accelerated exhumation along the southern plateau margin started sometime between 18 and 14 Ma in the Bitlis mountains (Cavazza et al., 2018) and from 14 Ma in the Zagros (Kurdistan region of Iraq; Koshnaw et al., 2020).Similarly, the northern plateau margin experienced accelerated exhumation from 18 to 14 Ma and possibly later in the Eastern Pontides (Albino et al., 2014), the Adjara Trialeti fold and thrust belt (Gusmeo et al., 2021), the Lesser Caucasus (Cavazza et al., 2019) and the Talesh Mountains (Chu et al., 2021;Madanipour et al., 2017).The magnitude of erosional exhumation along the southern plateau margin is greater than recorded along the northern margin (cooling from 200 to 150°C vs. 150-100°C, respectively) suggesting up to 4 km of difference in the thickness of eroded material (5-7 km vs. 3-5 km, respectively).Subsequently, the orogenic wedge expanded outward, incorporating new sectors of the Arabian (Cavazza et al., 2018) and the southern Eurasian foreland (Kura Basin; Gusmeo et al., 2021).Similarly, deformation propagated northward in Greater Caucasus, which experienced enhanced exhumation and topographic growth from ∼5 Ma, possibly after the closure of a relict Mesozoic back-arc basin located between the Greater and the Lesser Caucasus (Cowgill et al., 2016 and references therein).Overall, this pattern of exhumation indicates that most of the plate convergence that occurred during the last 18-15 Ma was accommodated along the margins of the plateau.
Collisional deformation appears to be almost coeval with a regional phase of voluminous volcanism that led to the widespread deposition of volcanics starting from 11 Ma (e.g., Keskin, 2003).These volcanics cover a large portion of the plateau (Figure 3) and mask older stratigraphic and tectonic contacts.Geochemical studies indicate that such a magma was originated by melting of a metasomatized and thickened Anatolian subcontinental lithospheric mantle that was destabilized by delamination processes after slab breakoff (e.g., Rabayrol et al., 2019 and references therein).The occurrence of volcanism is considered to be coeval with the development a system of lithospheric scale, dextral (North Anatolian Fault) and sinistral (East Anatolian Fault) strike slip faults that haves accommodated the westward escape of the Anatolian microplate from 13 to 11 Ma (e.g., Ballato et al., 2018;Racano et al., 2023;Şengör et al., 2005).The two faults intersect at the Karliova triple junction (KTJ in Figure 2a; see also Figure S1 in Supporting Information S1) along the western side of the plateau.The coincidence of all these events led some authors to suggest that there should be a genetic link between deep seated-processes (break-off, asthenospheric upwelling, large-scale mantle convection, mantle delamination, Anatolia extrusion), and large-scale surface uplift (e.g., Faccenna et al., 2006Faccenna et al., , 2013;;Keskin, 2003;Memis et al., 2020;Schildgen et al., 2014;Şengör et al., 2003).
Finally, GPS data indicate that at the longitude of Eastern Anatolia, the Arabian plate is moving northward to north-northwestward at ca. 20 mm/yr (Ahadov & Jin, 2017;Karakhanyan et al., 2013;Reilinger et al., 2006).To a first approximation, most of the plate convergence is accommodated through shortening along the southern (Bitlis and western continuation of the Zagros fold and thrust belt) and northern (Kura Basin and Greater Caucasus) margins of the collision zone (Ahadov & Jin, 2017;Reilinger et al., 2006).The plateau interior, instead, presents a complex pattern of deformation with a system of NW-SE to NNW-SSE right-lateral strike-slip faults accommodating a few mm/yr of slip, which can be locally associated with an extensional or a compressional component of less than 1 mm/yr (Karakhanyan et al., 2013).

Methods
To reconstruct the landscape evolution of the study area and in particular its uplift history, we analyzed the pattern of both topography and hydrography at the regional and local scale.We also generated paleolongitudinal profiles of the Kura and Arax Rivers and we inverted the longitudinal profile of the Arax R. to get more information on the modes and the rates of the surface uplift.

Map of the Elevation of the Change From Marine to Continental Deposition (Surface Uplift Map)
Using published geologic maps (Geological Map of Turkey scale 1:250.000,Geological Map of Georgia scale 1:500.000,Geological Map of Iran scale 1:1.000.000,Geological Map of Iraq scale 1.000.000,Geological Map of Arzebaijan scale 1:500.000,Geological Map of Armenia scale 1:500.000)and available scientific articles (Cavazza et al., 2018(Cavazza et al., , 2019;;Gelisli & Maden, 2006;Hüsing et al., 2009;Kergaravat et al., 2016;Okay et al., 2020;Popov et al., 2004Popov et al., , 2006;;Reichenbacher et al., 2011;Rolland et al., 2020;Varol et al., 2016;Şen et al., 2011), the stratigraphic boundary between the youngest marine deposits and the superseding continental sediments was 10.1029/2022TC007737 6 of 31 mapped along the entire study area (Figure S2 in Supporting Information S1).The present-day elevations of the obtained points were subsequently interpolated using the natural neighbor algorithm (Sibson, 1981) to obtain a surface uplift map with respect to present-day sea level.It is important to underline that this interpolation is local, considering the closest subset of data points, and does not generate unrealistic peaks or pits, ridges or valleys.The resulting surface passes through all the points (Figure S2 in Supporting Information S1), but it is smooth where they are not present.The available data are not uniformly spread throughout the entire study area affecting locally the geometry of the resulting surface (see Figure S2 in Supporting Information S1).

Filtered Topography at Different Wavelengths
Regional-scale morphologies of the EAP can be related to regional tectonics as well as to mantle dynamics.To distinguish and investigate them, we produced three maps by smoothing the raw elevation data (ETOPO1 DEM) in the frequency domain by a circular low-pass filter.The same method has been previously used in areas such as Yellowstone and the Colorado Plateau (Roy et al., 2009;Wegmann et al., 2007), Eastern Africa (Sembroni et al., 2016(Sembroni et al., , 2021)), the Apennines and Carpathians (D'Agostino & McKenzie, 1999;Faccenna et al., 2011;Molin et al., 2004Molin et al., , 2012) ) to isolate mantle-related topographic features.The size of these topographic signals is depend ent on the intensity, scale and depth of the mantle processes that generate them (Braun, 2010).For this reason, it is important to choose the appropriate filter wavelengths (Molin et al., 2012).Here, we chose three cut-off wavelengths: 50, 150, and 300 km.The first one represents a sort of threshold between the influence of crustal and mantle processes (Molin et al., 2004(Molin et al., , 2012;;Wegmann et al., 2007).The wavelengths of 150 and 300 km allowed to filter out the topographic signals of tectonic depressions and fluvial wide valleys and to catch the possible topographic component related to deeper mantle processes.

Topographic Swath Profile
A swath profile is a topographic profile displaying the pattern of the maximum, minimum and mean elevation within a rectangular area (observation window).It is generated by sampling the elevation data along strips perpendicular to the main direction of the swath (Isacks, 1992;Telbisz et al., 2013).The algebraic difference between the maximum and the minimum elevations at each point of the profile quantifies the local relief that, in a non-glaciated area, mimics the fluvial incision: in active tectonics landscape, regions of high local relief often correspond with areas of incision in response to uplift (Molin et al., 2004).The results are then condensed into a single plot allowing the description and quantification of general topographic features in an observation window (Isacks, 1992;Molin et al., 2004;Ponza et al., 2010).In this work, a SW-NE trending swath profile has been extracted from the ETOPO1 DEM by tracing 25 equally spaced topographic profiles into a 100 km wide window.The elevation data along each profile have been sampled every 2 km.
To better describe the topographic pattern, the plot of the swath profile was combined with three curves extracted from the topography filtered at 50, 150, 300 km wavelengths along the middle line of the observation window.
In the EAP, to identify the influence of tectonics on rivers at a regional and local scale, we extracted from the 90 m resolution SRTM DEM the longitudinal profiles of 33 rivers from the major drainage basins of the area, identified the major knickpoints, and calculated θ and k s from the log-log S versus A plot by the Stream Profiler tool (http://www.geomorphtools.org;Wobus et al., 2006).This tool removes spikes from the raw profiles with a moving smoothing window 250 m-in-length and it normalizes the steepness index (k sn ) by a reference concavity θ ref = 0.45.This normalization is needed because the two indices co-vary and to allow the comparison of longitudinal profiles with different drainage areas (Wobus et al., 2006).Here we chose θ ref = 0.45 since it appears to be most selected value in tectonically active areas (Kirby & Whipple, 2012;Leonard et al., 2023;Wobus et al., 2006) and falls within our θ estimates.
In the stream power model, river incision (E) into bedrock can be expressed as a power low relation between the upstream drainage area, A, and the local gradient, S (Howard, 1994;Whipple & Tucker, 2002): where K is rock erodibility and the exponents m and n are positive constants that depend on basin hydrology, channel geometry, lithology, and erosion processes (Howard, 1994;Whipple & Tucker, 1999).
Equation 2 can be solved for S as Considering steady-state conditions (dz/dt = E = U = 0), the integration of this equation from the base level (x b ) to an observation point (x) predicts the elevation of a channel profile (Perron & Royden, 2013): To use Equation 4to represent the river profile with units of length on both axes, a reference drainage area A 0 is introduced: In the χ plot, the transformed longitudinal profiles appear relatively linearized with the slope describing the steepness index.This representation allows a more precise knickpoints and channel slope analysis and a more direct comparison among different profiles.

Divide Stability
In order to quantify how changes in river basin drainage area have affected the EAP landscape evolution, we performed a drainage divide stability analysis among the major watersheds of the area.
The migration of drainage divide was first studied by Gilbert (1877) who proposed that the topographic asymmetry across a divide relates to different erosion rates on its sides.Such a difference forces the divide to move toward the side with a lower hillslope gradient and erosion rate.Lately, Willett et al. (2014) introduced a new method, which allows to evaluate the divide mobility through the parameter χ.Differences in χ values at the channel heads along adjacent catchments are indicative of disequilibrium conditions, which should be associated with the motion of the drainage divide toward the channels with higher χ.Instead, channels from adjacent catchments with similar χ are considered to be stable (Willett et al., 2014).Despite the quick and easy applicability of the χ method, several concerns have been raised (Whipple et al., 2017).To solve these problems, Whipple et al. (2017) proposed a set of three metrics (Gilbert metrics) to be compared with the χ analysis: the across divide difference in channel elevation at a reference drainage area, the mean headwater hillslope gradient, and the mean headwater local relief.
The analysis of the drainage divides stability of the EAP was carried out following Whipple et al. (2017) using the MatLab tool (http://github.com/amforte/DivideTools)elaborated by Forte and Whipple (2018).Considering the dependence of χ on outlet elevation, we cropped the 90 m resolution SRTM DEM to an elevation of 400 m, which approximately coincides with the bedrock-alluvial transition all over the study area.Finally, to visualize the result of the divide mobility analysis, we extracted the map of mean along channel local relief which represents the most efficient parameter to quantify divide stability (Forte & Whipple, 2018;Whipple et al., 2017).

Paleolongitudinal Profile
Assuming that the concave reach of a river profile upstream of a knickpoint was in equilibrium with previous base level conditions (Gallen et al., 2013), we reconstructed the shape of the ancient profile of the Kura and Arax Rivers.In particular, we performed the linear regressions on log-transformed slope versus area data above major knickpoints and, using a modified version of Equation 1, we reconstructed the paleo-longitudinal profile of the rivers: the channel upstream of a knickpoint is projected downstream to the previous outlet (Clark et al., 2006;Legrain et al., 2014;Olivetti et al., 2012;Pavano et al., 2016;Sembroni et al., 2020;Siravo et al., 2021).Excluding substantial Quaternary tilting of the upstream portion, the difference in elevation between the present and the reconstructed outlets can be considered as a minimum value of the amount of surface uplift since the onset of the erosional wave due to a fall in base level and presently represented by the migrating knickpoint (Kirby & Whipple, 2012 and references therein;Gallen et al., 2013).A similar result could be also obtained from a χ plot projecting the extrapolation of the steepness of a river segment to the elevation axis and estimating the amount of surface uplift from the difference between projected and present outlet elevation (see Clementucci et al., 2023;Fox et al., 2020;Smith et al., 2022).

River Longitudinal Profile Inversion
To reconstruct the base level history of the study area we inverted the longitudinal profile of the Arax R. starting from the elevation data of SRTM DEM (90 m of resolution) and following Pazzaglia and Fisher (2022).
As written above (Section 3.4), the value of k s is normalized by a reference concavity (θ ref ) to obtain k sn .Assuming detachment-limited condition, n = ∼1 (Tucker & Whipple, 2002;Whipple, 2004) and θ ref = 0.45, k sn can be expressed as This assumption is questionable since many rivers are not purely detachment-limited channel but are mixed (transition in time and/or in space from detachment-limited to transport-limited conditions).Moreover, channel morphology as well as longitudinal profile parameters are not diagnostic of one of the two conditions (Whipple & Tucker, 2002).In this respect, knickpoints are the landscape features that are typical of detachment-limited systems, since transport-limited ones are characterized by gradual changes in channel gradient (Whipple, 2004;Whipple & Tucker, 2002).In general, the simplification n = 1 does not strongly impact on the general results of the river longitudinal profile inversion (Gallen, 2018;Goren et al., 2014;Pazzaglia & Fisher, 2022;Pritchard et al., 2009;Racano et al., 2023;Roberts et al., 2012).
In case of a knickpoint migrating from the river outlet (x b ) to a given point (x) along a river channel (erosional wave), the response time of the system (τ) is expressed by combining Equations 1 and 6 and substituting dz/dx for S (Gallen, 2018;Pazzaglia & Fisher, 2022;Whipple & Tucker, 1999): Moreover, the elevation of the channel can be predicted by Equation 4 (Perron & Royden, 2013).In this way τ, χ, and z(x) provide the basis for a linear inversion of channels in order to unravel the base level fall history of a given area (U) (Goren et al., 2014;Pazzaglia & Fisher, 2022;Pritchard et al., 2009;Roberts & White, 2010).Some studies have reconstructed it with the assumption of uniform or nearly uniform uplift across the study area, a base level fall occurring at the mouth of the stream, and an erodibility coefficient spatially variable and dependent on bedrock lithology (e.g.Gallen, 2018;Pazzaglia & Fisher, 2022).Other similar inverse modelings consider K spatially uniform and solve for a non-uniform uplift (Goren et al., 2014;Pritchard et al., 2009;Roberts & White, 2010).
In the model used in this study (Gallen, 2018;Pazzaglia & Fisher, 2022), the spatial variability of the erodibility coefficient K is accomplished using the vectorized Lithologic Map of the World (Hartmann & Moosdorf, 2012a, 2012b; https://www.dropbox.com/s/9vuowtebp9f1iud/LiMW_GIS%202015.gdb.zip?dl=0; in Figure 3 is displayed the part of this map regarding the EAP) in order to define 13 domains of common rock erodibility (Table 1) where the mean value of k sn is calculated.
By Equation 7, the K values are calculated for each erodibility domain starting from the best fit k sn and the average erosion rates (in our model we used the values obtained from 10 Be concentrations in alluvium sampled along the Arax R. main channel; see location in Figure 2b; Kaveh Firouz, 2018).So, the K values are not predetermined, but reflect primary and secondary (like weathering, joints) properties of the rocks that define the susceptibility to erosion.The uncertainties associated with k sn and averaged erosion rate are calculated by a Monte Carlo routine with 1,000 independent simulations (Gallen, 2018;Pazzaglia & Fisher, 2022).The K values, associated to specific rock type and mapped onto the fluvial network, result in a τ-z plot.This rock-types specific τ values and the channel elevations are used in a linear inversion model (Goren et al., 2014;Pazzaglia & Fisher, 2022) to determine the base level fall history.For more details see Gallen (2018) and Pazzaglia and Fisher (2022).

Map of the Elevation of the Change From Marine to Continental Deposition (Surface Uplift Map)
To evaluate the surface uplift of the EAP, we generated a map that interpolates the current position of the stratigraphic boundary between the youngest marine sediments and subaerial deposits cropping out on the plateau (Figure 4; see also Figure S2 in Supporting Information S1).The age of these deposits ranges from ∼40 Ma to the north to 13-14 Ma in the rest of the plateau (Popov et al., 2004 and reference therein;Okay et al., 2020 and references therein), thus the inferred surface uplift pattern does not a have a temporal meaning but just a geometric one.Although most of these sediments were deposited in shallow-water marine systems and the stratigraphic boundary marks the change to a subaerial sedimentation, the surface uplift estimates contain uncertainties related to the sea level variations.In the last 35 Ma, the Parathetys fluctuated between ±200 m respect to the present sea level (Popov et al., 2010).So, our estimates of surface uplift are affected by this amount of uncertainty.The resulting map describes a general surface uplift of more than 2,000 m with respect to the present sea level.The highest values (∼3,000 m) are centered in the L. Van region, but they decrease rapidly along the steep Bitlis flank down to the Euphrates-Tigris lowland.The topographic depression enclosing L. Urmia disappears in the surface uplift map while lower values of surface uplift are observed in the Ararat basin (Figure 4).The local highs (e.g.L. Van area) and lows (e.g.Ararat basin) are surrounded by tectonic structures.

Filtered Topography at Different Wavelengths
To discern the contribution of crustal and subcrustal processes to the construction of topography at regional scale, we filtered the ETOPO1 elevation data at the wavelength of 50, 150, and 300 km (Figure 5).The filtered topography at 50 km describes the plateau as fragmented into local highs and lows bounded by tectonic structures (Figure 5a).Among the topographic highs (L.Sevan area, Lesser Caucasus, Eastern Pontides, L. Van area, Iranian Plateau), the highest values are just east of L. Van, in the very same area where a low velocity anomaly at a depth of 50 km is located (Medved et al., 2021).The filtered topography at 50 km presents also two main depressions in correspondence with the Ararat basin and L. Urmia.At a wavelength of 150 km, the filtered topography shows two highs, one centered in the EAP and one in the Iranian Plateau, whereas at a wavelength of 300 km, this latter disappears (Figures 5b and 5c).
In Figures 5d and 5e, we report the residual (i.e.isostatically anomalous) topography modified from Faccenna and Becker (2020).It has been estimated by subtracting the Airy compensated models of crustal structure (in other words, the isostatic compensated component of topography) from the actual topography.The crustal density in Figure 5d is constant, whereas it is variable in Figure 5e (see Faccenna & Becker, 2020 for a detailed explanation).Both maps show a positive signal over the EAP with the highest values (500-750 m) centered in the area of L. Van.Higher values are also present in correspondence with Eastern Taurus, Eastern Pontides, and the north-western portion of the Lesser Caucasus.The main negative signal is located in the Kura basin (Figures 5d  and 5e).These results are in agreement with the values calculated by Şengül Uluocak et al. (2021).

Topographic Swath Profile
We extracted a NNE-SSW oriented swath profile perpendicular to most of the compressive and strike slip structures (Figure 2c).The patterns of maximum, mean, and minimum elevation from the Bitlis to the Lesser Caucasus   1.000.000,Geological Map of Iraq scale 1.000.000,Geological Map of Arzebaijan scale 1:500.000,Geological Map of Armenia scale 1:500.000)and from scientific articles (Cavazza et al., 2018(Cavazza et al., , 2019;;Gelisli & Maden, 2006;Hüsing et al., 2009;Kergaravat et al., 2016;Okay et al., 2020;Popov et al., 2004Popov et al., , 2006;;Reichenbacher et al., 2011;Rolland et al., 2020;Varol et al., 2016;Şen et al., 2011).Inside the plateau, the location of the highest values of uplift are just north of L. Van while local minimum values are in correspondence with the Ararat pull-apart basin.again 3,000 m, but drops down to 800 m in the Ararat basin, increases again to almost 3,000 m in the Lesser Caucasus and finally decreases to 300-400 m in the Kura basin, which is delimited by the Greater Caucasus to the NNE.
In the swath profile plot (Figure 2c), we included three curves extracted from the surfaces of Figures 5a-5c

River Longitudinal Profiles and χ-Plots
To study the drainage evolution of EAP we extracted the longitudinal profiles (see Figures S3-S19 in Supporting Information S1) and the χ-plots of 33 rivers (see location in Figure 2b).To quantify the general shape of the longitudinal profiles, we also computed the concavity (θ) and normalized steepness (k sn ) indices (Table 2).Overall, the analyzed profiles are characterized by low/very low values of concavity (<0.40) with few exceptions (channels 7, 16, 22, 23, 26, 27, 33) that vary between 0.41 and 0.49 (Table 2).Most of the steepness indices are below the value of 100 and are mostly located in the plateau interior, including the one relative to the L. Van tributaries.Values above 100, instead, are observed in rivers clustering along the plateau margin such as the Lesser Caucasus (channels 7, 8, 9), the Eastern Pontides (channels 14, 15, 16), and the Bitlis Mts/L.Van area (channels 4, 28) (Table 2, Figure 2b).
The rivers draining the northern sector of the plateau (Kura, Arax, Ҫoruh drainage basins, channels from 1 to 16; Figure 2b) have concave up longitudinal profiles with knickpoints and knickzones at ∼1,500 m, upstream of which the rivers flow across the high-standing plateau (Figures 6a-6c).Here minor knickpoints are grouped around ∼2,500 (channels 1, 3, 4, 5, 13, 14, 15) and ∼3,000 m (channel 7 and 8) of elevation.Channel 10 has a knickpoint at ∼1,300 m where the stream passes across the Pambak-Sevan-Sunik Fault changing abruptly direction from NE-SW to NW-SE (Figure S1 in Supporting Information S1; Karakhanian et al., 2004).Channel 4 has a knickpoint at ∼900 that corresponds with the location of the intersection between the valley and the strike-slip Maku Fault (Figure S1 in Supporting Information S1; Karakhanian et al., 2004).Just below ∼700 m, the Arax R. (channel 3) longitudinal profile has a wide knickzone that corresponds upstream with the southern termination of the active strike-slip Nakhichevan Fault, one of the structures that originated the pull-apart Ararat basin, and downstream with an incised narrow valley (Figure S1 in Supporting Information S1; Karakhanian et al., 2004) (Figures 2 and 6b).In this area the active transpressive Pambak-Sevan-Sunik fault system passes across the Arax R. valley (Figure S1 in Supporting Information S1; Philip et al., 2001;Karakhanian et al., 2004).The mentioned convexities are evident in the χ-plots (Figures 7a-7c) where they mark a change in Residual topography maps (see text for explanation) modified from Faccenna and Becker (2020).Crustal density is constant in (d) and variable from CRUST1 model in (e).The values of density used for crust, lithosphere, and asthenosphere are 2,781 kg/m 3 , 3,300 kg/m 3 , and 3,226 kg/m 3 respectively.The EAP is characterized by a mean positive residual of ∼500 m with the highest (500-750 m) values located between the Van and Urmia lakes.the steepness of the channels.Actually, the river segments can be gathered into three groups: the ones above 2,500-3,000 m that show a low steepness (trend 1), the ones between ∼1,500 and 2,500-3,000 m, that are steeper (trend 2), and the ones below ∼1,500 m, that are the steepest ones (trend 3) (Figures 7a-7c).Note the knickzone at ∼700 m of channel 3 (Arax R.) that disrupts trend 3 into two segments by a decrease in χ (Figure 7b).
Among the rivers that drain the southern sector of the plateau (Figure 2b), the profiles of the Euphrates R. and its tributaries are concave-up but characterized by several knickpoints and knickzones (Figure 6d).Some of them correspond with small tectonics basins, the intersection with faults including the Northern Anatolian Fault (Figure S1 in Supporting Information S1; Gürbüz & Şaroğlu, 2019) or mark the entrance into the high-standing plateau at ∼2,500 m.In channels 17 (western Euphrates R.) and 18 (eastern Euphrates R. or Murat R.) several knickpoints are generated by dams and the relative reservoirs, but the one just below 1,500 m marks the Northern Anatolian Fault.The χ-plot (Figure 7d) shows that the entire channel 17 and the downstream portion of its tributaries have the same steepness (trend 4).The upstream segments of the tributaries have lower values of χ and often have the same steepness of trend 1 above ∼2,500 m of elevation.
The profiles of the Tigris R. (channel 22) and its tributaries show different shapes (Figure 6e).The Tigris R. profile is concave-up but characterized by a main knickpoint at ∼1,300 m of elevation located where the river passes across the Eastern Anatolian Fault.During the 1874 and 1875 earthquakes the 4 m dip slip component of the Palu-Lake Hazar segment of the Eastern Anatolian Fault beheaded the Tigris R. stopping the flow between the lake and the downstream fluvial network (Figure S1 in Supporting Information S1; Cetin et al., 2003).The other three profiles are steeper (especially channels 27 and 28; Table 2) with wide knickzones that correspond with several parallel thrusts and dextral strike slip faults of the Bitlis-Zagros Suture Zone (Figure S1 in Supporting Information S1; Copley & Jackson, 2006;Şengör et al., 2008;Sağlam Selҫ;uk, 2016;Sanҫar, 2021).Channel 27 profile includes just two rectilinear segments, it is shorter than the other tributaries, and is sourced from lower elevation (Figures 2  and 6e).This peculiarity could be related to the history of the L. Van that experienced variation in water volume and lake level causing the opening of the system to the southwest (Kuzucouglu et al., 2010;Stockhecke et al., 2014;Tomonaga et al., 2017).Here a windgap is located along the divide between the L. Van drainage and channel 27 suggesting that this river has been a lake outlet (Figure 2b).The χ-plot (Figure 7e) shows that the upper segment of channel 28 has a similar steepness of trend 4 above ∼1,700 m.Channel 28 has a steeper segment downstream with lower values in χ, showing a pattern that is typical of drainage area gaining (Willet et al., 2014).Conversely, the channel 22 (Tigris R.) is less steep than trend 4 (Figure 7e) suggesting a possible loss of drainage area.
The longitudinal profiles of the tributaries of L. Van (Figure 6f) are concave up and just the one of channel 25 shows a main knickpoint at ∼2,000 m.This is located where the river, after meandering across a marshy alluvial plain, flows into a narrow valley: here there is a step between the Erciş and the Dorutay right strike faults (Figure S1 in Supporting Information S1; McEnzie et al., 2016 and reference therein;Sağlam Selҫ;uk, 2016;Cukur et al., 2017).Channel 26 has the same steepness of trend 4 above ∼1,700 m, similarly to channel 28 that sources in the same area (Figures 2 and 7g).
The tributaries of L. Urmia have similar longitudinal profiles (Figure 6g): just few of them show a main knickpoint.In channels 29 and 31, it corresponds to a dam and the relative reservoir.In channel 33 it stands at ∼1,500 m, where the change from a wide alluvial plain upstream to a narrow and steep valley downstream occurs.This river steep segment ends where it intercepts the Northern Tabriz Fault (Figure S1 in Supporting Information S1; Karakhanian et al., 2004;Solaymani Azad et al., 2015).
In the χ-plot, L. Urmia and L. Van tributaries shows similar steepness that apparently mimic Trend 4 (Figures 7f  and 7g).

Divide Stability
We created a map displaying the mean local relief along the channels draining the study area (Figure 8) since this parameter is considered particularly efficient in quantifying divide stability (Forte & Whipple, 2018;Whipple et al., 2017)   the Pontides, the Lesser Caucasus) or where they are still trying to integrate into it (e.g., the Tigris R. tributaries draining the Bitlis).In the case of Arax R., high values are located downstream and upstream of the Ararat basin.

Paleolongitudinal Profiles
We reconstructed the paleolongitudinal profiles of the Kura R. (channel 1) and of its tributary Arax R. (channel 3) taking into consideration the segments that characterized their present longitudinal profiles (Figure 9).We have chosen the Arax R. because its long profile has been inverted (see the following Section 4.7) and the Kura R. as a term of comparison.
The Kura R. has a profile composed of two concave-up segments separated by a knickpoint standing at ∼1,700 m (Figure 9a).Reconstructing the upstream reach, we defined an outlet i.e., 1,000 m above the present one.So, the low relief upland standing above ∼1,500 m was affected by a base level lowering of ∼1,000 m.
The profile of the Arax R. includes three concave-up segments separated by two knickpoints at ∼700 and ∼1,500 m of elevation.The results of the modeling of the two upstream ones indicate that the Arax R. was apparently affected by two base level lowering of 500 m each.It is interesting to remember that (a) the pull-apart Ararat basin, that is coincident with the middle profile segment, was fully subaerial at least since 6.1 Ma (Popov et al., 2004); (b) the knickzone at ∼700 m corresponds with the Pambak-Sevan-Sunik fault and the Sardarapat-Nakhichevan fault.

River Longitudinal Profile Inversion
We inverted the longitudinal profile of the Arax R (channel 3) to reconstruct the base level fall history of the study area.As outlet we chose the Lesser Caucasus mountain front for several reasons: (a) although at present the Arax is a tributary of the Kura R., in the past they were separated (Hoogendoorn et al., 2005;Jorissen et al., 2020;Mamedov, 1997) Moreover, the location of the chosen outlet at the Lesser Caucasus mountain front could be considered constant in the last 12 Ma, whereas the Kura Basin has been alternatively submerged and emerged several times (Popov et al., 2004;Forte & Cowgill, 2013;Jorissen et al., 2020 and reference therein).
The result of the inverse modeling of the Arax R. basin (Figure 10 and Tables S1-S6 in Supporting Information S1) shows that the uplift rates inferred from Monte Carlo simulations (gray lines) have a variability ≥±0.1 mm/yr for each time step for ages older than 12.25 and 10.75 Ma in case of non-uniform and uniform rock-type, respectively.This wide range in rates is possibly related to the large variation in steepness index of the low order streams (i.e., of stream segments close to the divide).For this reason, only the younger portions of both reconstructions can be considered reliable.Moreover, according to Popov et al. (2004), between 14 and 13 Ma the Arax R. basin area was almost completely submerged and represented a shallow shelf marine environment apart from the most upstream portion.
The model result for non-uniform erodibility (see Table 1 for K values) indicates that between 12.25 and 10 Ma the uplift rate is ∼0.1 mm/yr, but from 10 to 5.25 Ma it increases at an almost constant rate up to more than 0.2 mm.In the elevation versus τ diagram, the knickpoint standing at 1,500 m corresponds to a response time of ∼5.5 Ma (Figure 10b).Then the uplift rate slightly decreases until 3 Ma, when it abruptly reaches the maximum value of 0.3 mm/yr at 2.75 Ma.It is noteworthy that between 5.5 and 2.75 Ma the Monte Carlo simulations show a wider variation in uplift rate.After a new decrease in uplift rate down to a minimum of 0.17 mm/yr at around 1.5 Ma, a new increase characterizes the last 0.25 Ma.
The values of k sn on which the modeling is based are affected by variations in rock-type erodibility and uplift rates.So, the non-uniform erodibility of rock types affects the base level history, increasing the variability in uplift rate and decelerating its decrease toward the present because of the high erodible rocks outcropping close to the river outlet (Pazzaglia & Fisher, 2022).So, we also modeled the Arax R. profile accounting for a uniform erodibility (mean values of K = 2.17 • 10 −6 ± 2.3 • 10 −7 m 0.1 yr −1 ), just to reduce the variability in uplift rate and to get a general picture (see Tables S7-S12 in Supporting Information S1).The result is characterized by fewer oscillations and by a general increase in uplift rate from ∼9 to ∼1.5 Ma.In detail, it is possible to discern three main time spans: (a) between 12 and 9 Ma the uplift rate is constantly around 0.05 mm/yr; (b) from 9 to 5 Ma it grows up to 0.2 m/ yr; (c) after a pause at around 4-5 Ma, it increases rapidly up to 0.35 mm/yr at ∼1 Ma and then decays (Figure 10).This pattern characterized by initial constant uplift rate, a first increase, a pause, a second more marked increase, and a final decrease is recognizable also in the result of model for non-uniform erodibility.The time intervals and the uplift rates are slightly different, but the patterns of the recent (last 12 Ma) uplift history are very similar.The results of both models are consistent with the ones from a similar approach by McNab et al. (2018).Although this modeling is more general and applied to the entire Anatolia, the uplift history for the EAP is very similar except for the last 5 Ma because it also includes the Kura and the Euphrates rivers basins.These two rivers, as evidenced by the χ-plots, have a different shape with respect to the Arax R., and hence are recording a different evolution.

Discussion
In this study the analysis of topography and hydrography coupled with geologic data and a river inversion model allowed us to reconstruct the base level fall history of EAP.The plateau has a symmetrical topography standing at a mean elevation of ∼2,000 m, steep and highly dissected margins and a low relief upland with both internally and externally drained basins (Figures 2, 5, and 8).Available low-temperature thermochronological data suggest the occurrence of a widespread phase of moderate late Eocene to early Oligocene exhumation (Albino et al., 2014).
From the middle-late Miocene, however, a renewed phase of exhumation has started along the northern and southern plateau margins where the topographic relief is greater (Cavazza et al., 2018(Cavazza et al., , 2019;;Chu et al., 2021;Gusmeo et al., 2021;Madanipour et al., 2017).Importantly, this recent exhumation event is not recorded in the low relief areas of the plateau interior (Albino et al., 2014).Stratigraphic evidence summarized in Popov et al. (2004) show that the EAP was close to the sea level during most of its history, being affected by several marine transgression-regression cycles until the late middle Miocene (12-11 Ma) when it became completely emerged.The map of the elevation of the change from marine to continental deposition (Figure 4) documents ∼2,000 m of surface uplift throughout the plateau.The age of the marine sediment varies from ∼40 Ma to the north to 13-14 Ma in the rest of the plateau (Popov et al., 2004;Okay et al., 2020 and references therein) indicating that the emergence and possibly the uplift were diachronous.These observations suggest that the onset of late Cenozoic surface uplift is almost coeval with the onset of exhumation along the plateau margins and of widespread magmatism (Rabayrol et al., 2019).Moreover, the occurrence of ∼500 m of positive residual topography (Figure 5) over a mean elevation of ∼2,000 m suggests a dynamic contribution of deep-mantle processes during the buildup of the EAP.
The general shape of the topographic dome centered around L. Van (Figures 2c and 5), however, has been disrupted by tectonic structures.Both Figures 4 and 5a show that local topographic highs and lows are bordered (black line) and its tributaries (gray lines).(b) τ-z plot of the Arax R. and of its tributaries.The color bar is scaled according to the rock domains from Table 1).(c) τ-U plot displaying the modeled rate of base level fall with non-uniform rock types (K values from Table 1).The vertical dashed line indicates the response time (12.25 Myr) when the Monte Carlo simulation have a variation in U ≤ 0.1 mm/yr.(d) τ-U plot displaying the rate of base level fall with a uniform rock type (K = 2.17 • 10 −6 ± 2.3 • 10 −7 m 0.1 yr −1 ).The vertical dashed line indicates the response time (10.75Myr) when the Monte Carlo simulations have a variation in U ≤ 0.1 mm/yr.In both (c and d) two horizontal arrows mark the time interval of the slowing emergence of the plateau and the onset of magmatism.
by faults.The highest values of steepness indices of the river longitudinal profiles are clustered in the uplifting blocks of these faults (i.e., in the Lesser Caucasus, in the Eastern Pontides, and in the Bitlis Mts/L.Van area but excluding the L. Van tributaries; Table 2, Figures 2b, 4, and 5a), confirming the occurrence of localized high uplift rates.Thus, tectonic deformation has also played an important role in shaping the EAP.
The fingerprint of deep-seated and shallower tectonic processes is retained in the landforms and the land surface hydrography.Although the fluvial network is integrating in the upland landscape of the plateau, deep fluvial incision in response to uplift is still confined along the plateau flanks where exhumation is greater (Figures 2 and 8).
The entrance of rivers into the high-standing surface is marked by a knickpoint at ∼1,500 m in the northern sector (Kura, Arax, Ҫoruh basins), but is absent or not evident to the south ( Figures 2, 6, and 7).The presence of large endorheic basins (Van, Sevan, Urmia lakes) as well as the migration of the drainage divide toward the plateau interior (Figure 8) record the transient state of the fluvial network organization, which is partially hampered by the active tectonics that beheaded and deviated rivers, generating steps in their channels.Indeed, many river longitudinal profiles have knickpoints and knickzones in correspondence with tectonic structures (Figure 6).An example is the wide knickzone below ∼700 m in the Arax R. (channel 3) longitudinal profile: the southern termination of the active strike-slip Nakhichevan Fault (Karakhanian et al., 2004), keeps the knickzone fixed in the landscape down-throwing continuously the floor of the pull-apart basin and blocking any possible erosional wave generated by a base level lowering (Figures 2 and 6b).
Interestingly, in the χ-plots, the river profiles (Figure 7) show different trends that allow to group the rivers draining the northern sector (the Kura, Arax, Ҫoruh basins characterized by three common trends) from the ones that drain the southern one (the Euphrates and Tigris basins characterized by a fourth different trend) and from the endorheic basins (Van and Urmia lakes).The four trends suggest a difference in the landscape evolution of the EAP.To investigate this latter, we modeled the uplift history of the Arax R. alone since erosion rate data are available only for the Arax and Kura basins (Kaveh Firouz, 2018).Unfortunately, it is impossible to model the uplift history of the Kura R. because some of its tributaries drain the Greater Caucasus and the result would mix the history of this mountain chain with the EAP one.
To facilitate the reconstruction of the EAP landscape evolution, we discuss separately the northern and southern sectors of the plateau.

The Northern Sector (Kura, Arax, Ҫoruh Drainage Systems)
The integration of the stream network analysis of rivers with the reconstruction of the base level fall history from the inversion of the Arax longitudinal profile (the results obtained accounting non-uniform and uniform rock-types) provides a scheme of the landscape evolution in the northern sector of the EAP.In particular, the χ-plots show three segments with uniform slope like the inverse modeling results of the Arax R., which indicate three main events (Figures 7 and 10).
1. Channels 7 and 8, that drain the Lesser Caucasus, present a river segment at elevation higher than 3,000 m.Similarly, channels 1, 3, 4, 5, 13, 14, 15 have a low steep segment, but at elevation higher than ∼2,500 m.This common high-standing reach (Figures 6 and 7) is interpreted to record the capture of an ancient landscape that developed close to the base level after the first emergence of the northern EAP.The rest of the northern sector of EAP was fully emerged since the early Miocene (20.5-19 Ma;Popov et al., 2004).In the early middle Miocene (16-15 Ma) a transgression flooded the lowlands that remained submerged in the Middle Miocene (14-13 Ma); only in the late Middle Miocene (12-11 Ma) the northern EAP was almost completely emerged (Popov et al., 2004).This configuration fits well with the modeling of the Arax R. profile: considering that the model is not reliable for response time >12 Ma, the inversion results indicate that before ∼10 Ma, the basin was affected by a low, almost constant uplift rate (∼0.1 mm/yr) that probably allowed the development of a low relief landscape, now standing at elevation higher than 2,500 m (Figure 10).It is noteworthy that most of the right tributaries of the Arax R. whose headwaters are along the divide with the Euphrates, Van and Urmia basins as well as the most upstream portion of the Arax R. itself have erosion rates equal or lower than 0.1 mm/yr (Kaveh Firouz, 2018).These low values fit well with a low-relief landscape developed close to base level.2. Most of the rivers of the Kura, Arax, Ҫoruh basins have a segment characterized by a similar steepness, usually delimited downstream by a knickpoint at ∼1,500 m, and steeper than Trend 1, named Trend 2 (Figures 6  and 7).This trend records a second ancient landscape, characterized by low local relief and developed close to base level that was subjected to an uplift of ∼1,000 m, as suggested by the paleo-longitudinal profile of the Kura and Arax rivers (Figure 9).The steeper Trend 2 apparently agrees with the progressive increase in uplift rate (up to more than 0.2 mm/yr) from the modeling of the Arax R. profile between 10 and 5 Ma.It is noteworthy that processes like lithospheric mantle delamination, slab break-off and upwelling of hot mantle materials are dated back to ∼10-15 Ma (Keskin, 2003(Keskin, , 2007;;Şengör et al., 2003;Faccenna et al., 2013 and references therein) as well as the onset of most of the volcanic activity, which started at ∼11 Ma (Keskin, 2003;Rabayrol et al., 2019).Moreover, at the same age an acceleration in exhumation rate has been observed in the Lesser Caucasus, Eastern Pontides and Talesh Mts (Cavazza et al., 2019;Chu et al., 2021;Gusmeo et al., 2021;Madanipour et al., 2017).3. The χ-plots of all the rivers of the Kura, Arax, Ҫoruh basins show a much steeper segment (named Trend 3) downstream of the knickpoint at 1,500 m (Figures 6 and 7).This result indicates that a regional uplift event affected the entire northern portion of the EAP in more recent time.The plot of uplift history of the Arax R (Figure 10) displays that after 5 Ma the uplift rate remained constant for a short time interval and then accelerated again until ∼2 Ma before decaying.Moreover, the knickpoint at ∼1,500 m has a response time of ∼5 Ma (Figure 10b).This event could not be related to a drop of ∼600-1,500 m of the Caspian Sea level from ∼5.5 to 3.2 Ma (Forte & Cowgill, 2013) as suggested for the Alborz mountains of north Iran (Ballato et al., 2015).The reason is that the Black Sea experienced much smaller oscillations between 0 and -100 m (Loget et al., 2006;Popov et al., 2010;Abdullayev et al., 2012;Forte & Cowgill, 2013 and references therein;Tari et al., 2015 and references therein), but the Ҫoruh R. and its tributaries have Trend 3 in their χ-plots like the rivers of the Kura and Arax drainage systems that flow into the Caspian Sea.It is interesting that both paleolongitudinal profiles of the Kura and Arax rivers reconstructed from their segment upstream of the knickpoint at ∼1,500 m indicate a drop in elevation of ∼1,000 m suggesting the occurrence of a similar process (Figure 9).In both basins the mean erosion rate is ∼0.2 mm/yr (Kaveh Firouz, 2018): if we hypothesize a similar uplift rate, the incision of ∼1,000 m will occur in 5 Myr, a result similar to the one from the inversion of the Arax R. profile.
The general pattern described above (i.e., an acceleration followed by a decay in base level fall rate), has been interpreted as a transient pulse in uplift rate driven by dynamic mantle processes (Pazzaglia & Fisher, 2022).In our case, this pattern is disrupted by an abrupt increase in uplift rate starting at ∼5 Ma.

The Southern Sector (Euphrates and Tigris Drainage Systems)
In the southern sector of the EAP, the presence of a large number of dams makes impossible to model the base level lowering history for the Euphrates and Tigris rivers basins.Therefore, the evolution of the southern sector of the EAP could be reconstructed just coupling our data with literature information.
The χ-plot of the Euphrates R. basin (Figure 7d) shows that the whole course of channel 17 (western Euphrates) and the downstream portion of its tributaries have a steepness (Trend 4) that is different from the ones found to the north (i.e., lower than Trend 3, but higher than Trends 1 and 2).Upstream, the higher values of χ, lower steepness, and the presence of Trend 1 indicate that these rivers are losing portion of their drainage area at elevations higher than 2,000-2,500 m by capture events.Indeed, the upstream portion of the drainage basin has been emerged since the middle Miocene (16-15 Ma; Popov et al., 2004).The emerged area progressively enlarged, and the coastline migrated to the south, probably defining Trend 4. It is noteworthy that this trend, although it could be apparently detected also in some channels of the Tigris, L. Van and L. Urmia drainage basins, is probably typical of the Euphrates drainage basin that is located at the periphery of the dome-like topography of EAP.
The χ-plot of the Tigris R. basin (Figure 7e) shows that, apart from channel 22 (Tigris R.), the river profiles have a steeper segment upstream of the confluence with the Tigris R.This feature could be at least partially linked to the exhumation history of their drainage area.The Euphrates and Tigris rivers pass across the westernmost portion of the Bitlis Mts, where the mean topography presents lower elevation and where exhumation is younger than 20 Ma (Cavazza et al., 2018).Channels 23, 27, 28 drain the portion of the Bitlis Mts that were exhumed after 15 Ma (Albino et al., 2014;Cavazza et al., 2018) and are characterized by a much rougher and steeper topography.
Some ages of fluvial terraces of both the Euphrates and Tigris rivers (e.g., Bridgland et al., 2017;Demir et al., 2008Demir et al., , 2009;;Westaway et al., 2009;and references therein;Stow et al., 2020) allow to calculate rough estimates of the incision rates that range from ∼0.03 mm/yr in the Arabian Plate to 0.07 mm/yr in the southern EAP over the last ∼2 Ma.These values appear to be much lower than the 10 Be-derived estimates of ∼0.2 mm/yr (Arax and Kura rivers basins) averaged over shorter time scales for the northern EAP (Kaveh Firouz, 2018).Apart from the ones relative to the Arabian foreland, the low values in the southern EAP could be due again to their peripheral location with respect to the center of the dome (Figure 5).
In summary Trend 4 possibly records the progressive emergence since 16-15 Ma of the southwestern EAP coupled with the influence of the slowly uplifting foreland (Euphrates-Tigris valley) on the downstream portion of fluvial network that induced the migration of the coastline to the south.Trend 4 is not present in the southeastern portion of EAP because the center of the dome (L.Van area) is closer, and the topography is rougher making the rivers steeper.
In conclusion, despite some differences among the southern and northern portion of the EAP, the uplift history of the plateau is mainly characterized by (a) a slow emergence that allowed the development of low relief landscapes that now is at elevation higher than ∼1,500 m; (b) a first increase in uplift rate around 10 Ma; (c) an acceleration in uplift rate (at least since 5 Ma) that probably generated most of the present topography.

The Uplift of the EAP in the Arabia-Eurasia Collision Zone
Our data show that, on ∼2,000 m of elevation of the EAP, ∼1,500 m are isostatically compensated and ∼500 m are positive residual topography.Moreover, the combination of river longitudinal profile analysis and inversion allows to define three main steps in the uplift history of the EAP and hence provide the basis for answering a first-order question: in the framework of the Arabia-Eurasia collision, what processes have driven this peculiar upheaval pattern of the plateau uplift?In general terms, uplift may be attributed to the isostatic response to crustal thickening or to a dynamic mantle origin (i.e., to a transient feature generated by mantle convection).The presence of a low velocity anomaly, positive free air gravity anomaly, large wavelength uplift, deep source volcanism, and positive residual topography are considered markers of a dynamic component of topography.This means that at least a fraction of the elevation of the EAP should be dynamically supported (Faccenna et al., 2013;Goğus & Pysklywec, 2008;Keskin, 2003;Kounoudis et al., 2020;Şengör et al., 2003;Şengül Uluocak et al., 2021).In particular, three classes of models have been proposed so far.
• The first one considers the uplift resulting from the removal of the mantle lithosphere in the form of mantle delamination or in the form of Rayleigh-Taylor instabilities.Being colder and hence denser than the underlying asthenosphere, the removal of the mantle lithosphere results first in a subsidence and then in a net gain in elevation (Goğus & Pysklywec, 2008;Keskin, 2003;Kounoudis et al., 2020).Surface wave inversions (Şengör et al., 2003) and joint inversion of receiver functions and surface waves (Gök et al., 2007) support this class of models suggesting a reduced or absent mantle lithosphere beneath the EAP.• The second class of models considers a break in the slab favored by the protracted continental collision as the main cause of the uplift (Bottrill et al., 2012;Faccenna et al., 2006;Keskin, 2003;Schildgen et al., 2014).This is in good agreement with the lack of high velocity anomaly (Govers & Fichtner, 2016;Kounoudis et al., 2020;Piromallo & Morelli, 2003) and with a change in volcanic signature around 10 Ma (Faccenna et al., 2013 and references therein).• The third class of models considers the plateau uplift resulting from a mantle asthenosphere support, perhaps related to the presence of hotter-than-average mantle, that may derive from the Arabian plate (Faccenna et al., 2013;Keskin, 2007;Şengül Uluocak et al., 2021).
Those three mechanisms can be also considered as part of a single evolution because the inflow of hotter mantle material in the EAP can occur only after the slab break-off (Faccenna et al., 2013).
Our results can be integrated into the geodynamic models proposed by Keskin (2003Keskin ( , 2007) ) and Faccenna et al. (2013).In particular, the combination of geological data with our modeled uplift rates suggests that before ∼10-11 Ma the study area was characterized by low uplift rate that allowed the progressive emergence of the EAP (Popov et al., 2004) and the development of a low relief landscape now at elevations above 2,500-3,000 m, that is drained by low steepness channels (Trend 1; Figures 2, 7, 8 and 10).At this time (around 15-16 Ma) the EAP was shortening and thickening over the subducting slab (Figure 11a).In particular, the onset of fault-related exhumation occurred in the Bitlis, in the Eastern Pontides, and in the Lesser Caucasus (Cavazza et al., 2019 and reference therein).
Our analysis shows a two-step, probably progressive increase in uplift rate at 10 and then at 5 Ma.The first step, around 10-11 Ma, can be related to the episode of slab break-off (Faccenna et al., 2006;Schildgen et al., 2014), as attested by the change in volcanic signature with Arabian mantle imprinting in Karacadağ first and later in the Anatolia volcanism (Figure 11b; Keskin, 2003;Faccenna et al., 2013;Agostini et al., 2021).The surface effect of this process can be twofold.The uplift related to slab break off on one side (Duretz et al., 2011), and the inflow of hotter asthenosphere flowing northward from the Arabian Plate through the slab window.Together those processes may be responsible for the acceleration of surface uplift rate (Figure 10) inducing the headward propagation of an erosion wave and the development of a landscape with channels steeper (trend 2) than the upstream segments (Figures 2, 7, 8, and 10).
Around 4-5 Ma, the slab window continued to widen, and the inflow from the Arabia mantle reached the central and the western Anatolia (Faccenna et al., 2013), providing further supplementary support to the EAP's topography.Moreover, the presence of hot asthenosphere induced the thermal erosion of the lithosphere and the heating of the crust, strongly reducing the lithosphere strength and favoring its deformation (Keskin, 2003(Keskin, , 2007;;Lanari et al., 2023).This may result in a further increase of the surface uplift and rate, with a combination of static and dynamic contribution (Figures 10 and 11c).The ongoing shortening, the erosion of the lithosphere, and the dynamic support of the mantle flow grew the topography, generating the upstream propagation of a third wave of incision and a steepening of river channels (Trend 3; Figures 5a  and 7a-7c).To the west, where the plate convergence is mostly accommodated by strike slip tectonics (north and east Anatolian faults zones), the sole dynamic support of the mantle flow results in a lower topography and less steep channels (Trend 4).
The dynamic support of the mantle progressively contributed to the surface uplift of EAP and is still responsible for about ∼500 m of the present-day positive residual topography.The inflow of the Arabian mantle may have reached the southern margin of the Central Anatolian Plateau during the Pleistocene as suggested by the two-to three-fold increase in surface uplift rates at ∼1.6 Ma (Schildgen et al., 2012) following the onset of uplift sometime after 8 Ma (Cosentino et al., 2012).There, however, an additional peak in uplift rates occurred around the middle Pleistocene most likely during crustal rebound following slab tearing and break-off (Racano et al., 2021).

Conclusion
The Eastern Anatolian Plateau (EAP) is considered a continent-continent collision plateau at the early stage of development (Hatzfeld & Molnar, 2010;Keskin, 2007).With respect to the more evolved Tibetan Plateau, the main differences are related to the lesser penetration of the Arabian Plate into Eurasia than India, but the processes that generated the high-standing topography are similar (Hatzfeld & Molnar, 2010).So, to better understand the uplift history of an early stage continent-continent collision plateau, we quantified the surface uplift of EAP, focusing on the geometry of topography and hydrography and inverting the longitudinal profile of the Arax R. The results of our analysis indicate.
• The plateau has dome shaped topography (Figures 2 and 5) characterized by a high-standing low relief surface that experienced ∼2,000 m of surface uplift and that reaches maximum values of ∼3,000 m in the area of L. Van (Figure 4).On this ∼2,000 m of elevation, ∼1,500 m are isostatically compensated and ∼500 m are supported by dynamic processes.• Tectonics structures, both compressive and strike slip with oblique component, generated local topographic highs and lows that include endorheic basins (Figures 2-5).• The fluvial network pattern indicates a disorganization: rivers are trying to integrate into the high-standing plateau by deeply incising its flanks, capturing progressively the internally drained portion of it, and pushing the drainage divide toward the plateau interior (Figures 2 and 8).In this context, the river pattern is locally controlled by tectonic structures.The irregular shape of the longitudinal profiles (Figure 6) confirms the transient state of disequilibrium of hydrography and the χ-plots (Figure 7) describe a progressive channel steepening apart for those rivers that drain the western peripherical portion of the plateau.• The inversion of the Arax R. longitudinal profile (Figure 10), supported by the reconstruction of the Kura and Arax rivers paleolongitudinal profiles, as well as by the χ-plots (Figures 7 and 9), allowed to describe the base level lowering history of the EAP, evaluating the timing and rates.In detail, the results describe a three-phases pattern characterized by a progressive increase in uplift rate.• By integrating the reconstructed uplift history of EAP with available geodynamic models, we propose that before ∼10 Ma the study area was characterized by a low, almost constant uplift rate of ∼0.1 mm/yr.At 10-11 Ma, the opening of the slab window concomitant with the arrival of the mantle flow from beneath Arabia supported dynamically the topography of the EAP making the uplift rate to increase up to ∼0.2 mm/ yr.At ∼5 Ma, the continued inflow from Arabia together with the isostatic response to the ongoing shortening and crustal thickening started after 15 Ma at both northern and southern margin of EAP induced a further increase in uplift rate that exceeded 0.3 mm/yr.
In conclusion, the EAP represents an example of a collisional plateau at an early stage of development, initially slowly emergent, but successively affected by a first acceleration in uplift rate induced by mantle processes and by a second stronger increase due to the addition of the isostatic contribution of crustal shortening and thickening to the dynamic support.

Figure 1 .
Figure 1.Tectonic sketch of the study area and surroundings (modified from Cavazza et al. (2019)).The rectangle indicates the location of the Eastern Anatolia Plateau (EAP) as part of the Anatolian-Iranian Plateau.

Figure 2 .
Figure 2. (a) Topography (from SRTM DEM at 90 m of resolution) of the study area including the main tectonic structures.The white lines are drainage divides whereas the white rectangle is the observation window of the swath profiles of (c).The label KTJ marks the location of the Karliova Triple Junction.(b) Slope map of the study area showing the high-standing plateau characterized by low values of hillslope and the surrounding flanks by high values.The black lines are the study rivers named by a number, the white dashed line indicates the boundary of the Eastern Anatolia Plateau (EAP) that merges to the Iranian Plateau to the SE, the red star marks the location of a wind gap.The white dots along the Arax R. indicate the sampling locations for the estimation of the erosion rates, whereas the side numbers indicate the results in mm/yr (from Kaveh Firouz (2018)).(c) NNE-SSW oriented swath profile across the EAP.The three colored curves are extracted from the filtered topography of Figures 5a-5c (see the text for explanation) along the middle line of the observation window of the swath profile.

Figure 3 .
Figure 3. Geological map of the study area modified from the Lithologic Map of the World(Hartmann & Moosdorf, 2012a, 2012b).The main tectonic structures are reported in black.
are very similar, showing an almost symmetrical dome-like topography.From the lowlands of the Euphrates and Tigris rivers valley (around 300 m) the topography increases rapidly up to 3,000 m in the Bitlis Mts.This pattern is abruptly interrupted by the internally drained depression (∼1,700 m) of the L. Van.It is interesting to note that (a) the SW sector of the lake is ∼450 m deep(Cukur et al., 2017), indicating a variation in elevation of almost 1,800 m from the top of the Bitlis Mts to the bottom of the lake; (b) the internally drained basin (Figure2) is located in a topographic high.To the NNE, beyond the basin, the maximum topography reaches

Figure 4 .
Figure 4. Surface uplift map elaborated by the interpolation of the elevation of outcrops characterized by the change from marine to continental deposition.The data derive from published geologic maps (Geological Map of Turkey scale 1:250.000,Geological Map of Georgia scale 1:500.000,Geological Map of Iran scale 1:1.000.000,Geological Map of Iraq scale 1.000.000,Geological Map of Arzebaijan scale 1:500.000,Geological Map of Armenia scale 1:500.000)and from scientific articles(Cavazza et al., 2018(Cavazza et al., , 2019;;Gelisli & Maden, 2006;Hüsing et al., 2009;Kergaravat et al., 2016;Okay et al., 2020;Popov et al., 2004Popov et al., , 2006;;Reichenbacher et al., 2011;Rolland et al., 2020;Varol et al., 2016;Şen et al., 2011).Inside the plateau, the location of the highest values of uplift are just north of L. Van while local minimum values are in correspondence with the Ararat pull-apart basin.
. The one relative to the topography filtered at a wavelength of 50 km well describes the crustal component of topography, showing the dome extending from the Euphrates and Tigris lowlands to the south to the Kura basin to the north and evidencing the highs (e.g., the Bitlis Mts or the Lesser Caucasus) and lows (the L. Van depression and the Ararat basin) that characterized the plateau.At longer wavelengths (150 and 300 km), the crustal component disappears and the EAP topography takes the shape of a dome centered in the L. Van area.

Figure 5 .
Figure 5. (a) Filtered topography at a wavelength of 50 km, a threshold between crustal and deep-seated processes.Local highs and lows are bounded by tectonic structures (in black).(b) Filtered topography at a wavelength of 150 km showing the dome morphology of the Eastern Anatolia Plateau (EAP) separated from the Iranian Plateau by a local low around Lake Urmia.(c) Filtered topography at a wavelength of 300 km describing the dome centered in the L. Van area.(d and e)Residual topography maps (see text for explanation) modified fromFaccenna and Becker (2020).Crustal density is constant in (d) and variable from CRUST1 model in (e).The values of density used for crust, lithosphere, and asthenosphere are 2,781 kg/m 3 , 3,300 kg/m 3 , and 3,226 kg/m 3 respectively.The EAP is characterized by a mean positive residual of ∼500 m with the highest (500-750 m) values located between the Van and Urmia lakes.
and locally the fluvial incision.The comparison of the along channel local relief values across the main divides allows to understand the occurrence of ongoing drainage divide migration or stability in the study area.Inside the plateau, the values of the along channel local relief appear to be very similar across the divides indicating stability.Different values are located uphill on the plateau flanks, suggesting a migration of the external divides toward the plateau interior.Indeed, the highest values of local relief, indicating a deeper incision, are where the rivers pass across the plateau flanks (e.g., the Euphrates to the E and the rivers that drain

Figure 7 .
Figure 7. (a) χ-plots of the rivers of the Kura drainage basins.Note that some rivers are characterized by an upper segment above an elevation higher than ∼2,500 m (Trend 1, orange dashed line), a middle segment between ∼2,500 and ∼1,500 m (Trend 2, light blue dashed line) and a lower and steeper segment below ∼1,500 m (Trend 3, red dashed line).(b) χ-plots of the rivers of the Arax drainage basin.Note the presence of Trends 1, 2, 3. (c) χ-plots of the rivers of the Ҫoruh drainage basin showing the Trends 1, 2, 3. (d) χ-plots of the rivers of the Euphrates drainage basin.Note that most of the rivers display a similar channel steepness (Trend 4, light green dashed line).(e) χ-plots of the rivers of the Tigris basin.Apparently, they do not show a common trend, apart the upstream segment of channel 28 that resembles Trend 4. (f) χ-plots of the tributaries of the Lake Urmia drainage basin.(g) χ-plots of the tributaries of the Lake Van drainage.

Figure 8 .
Figure 8. Map of the mean local relief along the channels draining the study area, generated following Forte and Whipple (2018).The relief values across the drainage divides (black lines) indicate the possible direction of migration.In detail, the divide moves toward local lower values (cold colors) as indicated by the white arrows.In the study area, the highest values are located along the Eastern Anatolia Plateau (EAP) flanks where fluvial incision is deeper and the drainage divides migrate toward the plateau interior.There, the internal divides appear to be stable.

Figure 9 .
Figure9.Reconstructed paleolongitudinal profiles of the Kura and Arax Rivers starting from the upstream segments of their present longitudinal profiles (see text for detailed explanation).Due to the irregularity of the segments the uncertainty is quite high, but in both Kura and Arax Rivers, the modeling of the segment above the knickpoint at ∼1,500 m results in a base level lowering at the outlet of ∼1,000 m.
; (b) in the Kura basin, the Arax R.changed its location as indicated by abandoned meanders and meander scars apparently shifting to the north; (c) the Kura basin, a subbasin of the Southern Caspian Sea (Bagirov & Lerche, 1999), has been flooded several times during high-stands in Caspian base level (Forte & Cowgill, 2013; Jorissen et al., 2020 and reference therein); (d) the high erodibility of the soft sediments that fill the Kura basin would influence the model results impacting the values of uplift rate (Pazzaglia & Fisher, 2022).

Figure 10 .
Figure 10.Result of the inversion of the longitudinal profile of the Arax R. reconstructing the base level fall history of the study area.(a) χ-z plots of the Arax R.(black line) and its tributaries (gray lines).(b) τ-z plot of the Arax R. and of its tributaries.The color bar is scaled according to the rock domains from Table1).(c) τ-U plot displaying the modeled rate of base level fall with non-uniform rock types (K values from Table1).The vertical dashed line indicates the response time (12.25 Myr) when the Monte Carlo simulation have a variation in U ≤ 0.1 mm/yr.(d) τ-U plot displaying the rate of base level fall with a uniform rock type (K = 2.17 • 10 −6 ± 2.3 • 10 −7 m 0.1 yr −1 ).The vertical dashed line indicates the response time (10.75Myr) when the Monte Carlo simulations have a variation in U ≤ 0.1 mm/yr.In both (c and d) two horizontal arrows mark the time interval of the slowing emergence of the plateau and the onset of magmatism.

Figure 11 .
Figure 11.Cartoon illustrating the uplift of the Eastern Anatolia Plateau (EAP) in the frame of the Arabia-Eurasia collision zone from ∼16 Ma to Present.Note that the red triangles indicate volcanoes, the red arrows the flow of hot asthenosphere.(a) The EAP was shortening and thickening over the subducting slab and sparse volcanism is fed by magma from a mantle source modified by subduction components.A low uplift rate induced the progressive subaerial exposure of EAP and the development of a landscape (Trend 1, Figures 2, 7, 8 and 10).(b)The slab started to break-off, allowing the penetration of hot asthenosphere from Arabia as attested by the spread magmatism with a change in volcanic signature.Both these processes (slab break-off and inflow of hotter asthenosphere) resulted in an increase in surface uplift rate that induced the headward propagation of an erosion wave (Trend 2; Figures2, 7, 8, and 10).(c) The slab window continued to widen, and the inflow of Arabia mantle continued to support the EAP's topography.The thermal erosion of the lithosphere and the heating of the crust reduced the lithosphere strength and favored its deformation resulting in a further acceleration in surface uplift rate and a consequent generation of a third erosion wave (Trend 3; Figures5a, 7a-7c, and 10).The elaboration of the three cross sections is based onAngus et al. (2006),Vanacore et al. (2013),Ogden and Bastow (2022),Darin and Umhoefer (2022) for modern and past crustal thickness, onAngus et al. (2006),Nikogosian et al. (2018),Agostini et al. (2021),Reid et al. (2019),Rabayrol et al. (2019) for modern and past magmatism and magma source, on Ismail-Zadeh et al. (2020),Biryol et al. (2011) for modern mantle structure.

Table 2
Values of the Concavity, the Steepness and the Normalized Steepness Indices of the Studied Rivers and Relative Errors Figure6.Plots of the longitudinal profiles of the rivers draining the study area.The red lines are the main channels (5: Kura R., 7: Arax R., 17: Ҫoruh R.; 21: western Euphrates R.; 26: Tigris R.) and the black ones are tributaries.The red arrows indicate the location of knickpoints corresponding to the main tectonic structures.Note that many knickpoints and knickzones characterize the river longitudinal profiles, despite the upward concavity of all of them.