Variable Quaternary Uplift Along the Southern Margin of the Central Anatolian Plateau Inferred From Modeling Marine Terrace Sequences

The southern margin of the Central Anatolian Plateau (CAP) records a strong uplift phase after the early Middle Pleistocene, which has been related to the slab break‐off of the subducting Arabian plate beneath the Anatolian microplate. During the last 450 kyr the area underwent an uplift phase at a mean rate of ~3.2 m/kyr, as suggested by Middle Pleistocene marine sediments exposed at ~1,500 m above sea level. These values are significantly higher than the 1.0–1.5 m/kyr estimated since the Late Pleistocene, suggesting temporal variations in uplift rate. To estimate changes in uplift rate during the Pleistocene we studied the marine terraces along the CAP southern margin, mapping the remnants of the platforms and their associated deposits in the field, and used the TerraceM software to identify the position and elevation of associated shoreline angles. We used shoreline angles and the timing of Quaternary marine sedimentation as constrains for a Landscape Evolution Model that simulates wave erosion of an uplifting coast. We applied random optimization algorithms and minimization statistics to find the input parameters that better reproduce the morphology of CAP marine terraces. The best‐fitting uplift rate history suggests a significative increase from 1.9 to 3.5 m/kyr between 500 and 200 kyr, followed by an abrupt decrease to 1.4 m/kyr until the present. Our results agree with slab break‐off models, which suggest a strong uplift pulse during slab rupture followed by a smoother decrease.


Introduction
Orogenic plateaus, defined by widespread low-relief landscapes, strongly elevated in comparison with their surrounding areas, represent singular features of the main orogenic belts. These low-relief elevated areas are responsible for both local and global climate changes. The uplift of the Himalayan-Tibetan Plateau, in the Eurasian orogen, and the Puna Plateau, in the Andean Cordillera, has been related to major climatic changes caused by the growth of these topographic barriers (e.g., Ehlers & Poulsen, 2009;Gregory-Wodzicki, 2000;Harris, 2006;Hartley, 2003;Lenters & Cook, 1997;Molnar et al., 1993;Ruddiman & Kutzbach, 1989;Strecker et al., 2007;Zhisheng et al., 2001). Therefore, defining the timing and rates of plateau growth is crucial for understanding the driving mechanisms.
In Asia Minor, the Central Anatolian region is defined by a high topography area known as the Central Anatolian Plateau (CAP). It marks the western portion of the largest collisional belt on Earth, the Turkish-Iranian-Caucasian-Himalayan-Tibetan orogen, resulting from collision between the Arabian and Indian plates with Eurasia, which in eastern Anatolia started at ca. 20 Ma along the Bitlis-Zagros thrust zone (e.g., Ballato et al., 2011;Okay et al., 2010). along the CAP southern margin documented by its deep erosion Faranda et al., 2013;Radeff et al., 2016Radeff et al., , 2017 have revealed two strong uplift pulses in the late Messinian and middle Pleistocene. Uplift rate estimates range between~4.1 m/kyr for the first pulse and 3.2-3.4 m/kyr for the second . Afterwards, in the late Pleistocene-Holocene, the CAP southern margin has been characterized by continuous uplift rates around 1-1.5 m/kyr Serpelloni et al., 2013).
These two uplift phases at the CAP southern margin, separated by a long period of subsidence (ca. 5 Myr) agree with thermo-mechanical models that suggest a rapid decrease of surface uplift after a slab break off as a consequence of continental collision (Duretz & Gerya, 2013). Lithospheric processes, such as crustal delamination and slab break off, have been invoked to explain these uplift rate variations (e.g., Schildgen et al., 2014). However, despite these estimations on the uplift history of the CAP southern margin, the detailed evolution of Quaternary uplift has yet to be disclosed. The main information on the Quaternary uplift is from the Gülnar area, where middle Pleistocene deep marine deposits are exposed at 1,100 m elevation, suggests an uplift of~1,500 m in the last 450 kyr (~3 m/kyr; . In addition, at Silifke 14 C ages of Holocene sea level markers  indicate lower rates (ca. 1-1.5 m/kyr), which are still high when placed in a global context (Pedoja et al., 2014). However, the causes of the decrease in uplift rates and when it happened are even still open questions. Along these lines, the staircase morphology of the coastal area south of Gülnar offers the possibility to obtain detailed uplift rate estimates to bring light on these questions.
In coastal areas, like the Mediterranean coast of the CAP southern margin, geomorphic features such as marine terraces represent the surface expression of combined tectonics and eustatic sea-level variations, providing detailed information about surface deformation during the Quaternary (e.g., Anderson et al., 1999;Bloom & Yonekura, 1985;Cosentino & Gliozzi, 1988;Kelsey & Bockheim, 1994;Lajoie, 1986). Tectonic uplift is the main factor responsible for the preservation of marine terraces, forming staircase sequences in areas of rapid uplift rates (higher than 0.2 m/kyr, Pedoja et al., 2014). During periods of sea-level lowstand, marine platforms are abandoned, and new platforms can develop during the next sea-level rise (Anderson et al., 1999;Bloom et al., 1974). At lower uplift rates, polygenetic marine terraces known as rasa surfaces (Melnick, 2016;Paskoff, 1977;Pedoja et al., 2014;Regard et al., 2010;Westaway, 1993) can develop. These kinds of terraces are characterized by wide coastal plains bounded by steep cliffs, related to reoccupation during the successive sea-level highstands. Low uplift rates can also produce a total overprinting of older terraces, as demonstrated in some areas of the Mediterranean such Crete (Robertson et al., 2019) and southern Italy (Meschis et al., 2018).
The objective of this study is elucidating the morphostructural and morphostratigraphical evolution of the margin and estimating the temporal evolution of surface uplift during the last 0.5 Ma. To accomplish this objective, we carried out detailed mapping marine terrace sequences along the CAP southern margin, estimating uplift rates using novel numerical approaches. Our results suggest a dramatic but transient increase in uplift suggesting that slab break-off processes play a significant control on vertical movements along the CAP.
The modern topography of the CAP southern margin has been associated to crustal-scale deformations and lithospheric processes Schildgen et al., 2012Schildgen et al., , 2014. In the CAP northern margin the uplift has been related to the strain accumulation along the North Anatolian Fault (Yildirim et al., 2011). The CAP southern margin experienced two strong uplift phases occurred in the late Miocene Tectonics Radeff et al., 2017;Schildgen et al., 2012) and in the middle Pleistocene . These authors correlated the geological observations of these two uplift phases to the lithospheric delamination of the old subducting plate and its break-off, as suggested by mantle tomography and seismic studies (Biryol et al., 2011;Gans et al., 2009;Imprescia et al., 2012;Kalyoncuoğlu et al., 2011;Portner et al., 2018), which indicate the absence of subducting lithospheric slab east of Cyprus. Furthermore, the plateau morphology also shows a notorious asymmetry along the north-south direction and northward tilt that is reflected in the predominant northward direction of the plateau drainage system Schildgen et al., 2014).

Geological Setting of the CAP Southern Margin
The Central Taurides belt, which delimitates the CAP at the southern margin, is characterized by oceanic and continental units associated to shortening during the late Mesozoic-Eocene orogeny ( Figure 2). This shortening phase was associated to the closure of the Neotethys Ocean and the collision between the Taurus carbonate platforms and the crystalline complex of the Central Anatolia (Dixon & Robertson, 1894;Sengör & Yilmaz, 1981). This area records a post-orogenic extension phase occurred during the late Eocene-early Miocene (Akay et al., 1985;Aksu et al., 1992;Görür, 1992;Robertson & Grasso, 1995;Williams & Unlugenç, 1992), associated with the development of suprasutural sedimentary basins as the Mut-Ermenek Basin. These basins record a first depositional cycle of continental units, mostly characterized by basal coarse clastic fluvial and lacustrine beds (Lower Oligocene) overlain by fine-grained lacustrine deposits (upper Oligocene-lower Miocene). These continental units unconformably overlie on the deformed Taurides bedrock and are generally tilted and faulted by the continuing extension .
During the lower-to-upper Miocene, a marine sedimentary succession was deposited in the Mut-Ermenek Basin, comprising shallow-to-deep water carbonates (Bassant et al., 2005;Cipollari et al., 2013) unconformably lying onto the upper Oligocene-lower Miocene continental units (Fakirca Fm) or directly resting onto the highly deformed Taurides Units. The uppermost part of this marine succession was dated late Tortonian Cosentino et al., 2012), currently exposed at elevations of 2 km.
Recently, lower-to-middle Pleistocene deep-marine deposits (Sarıkavak Formation) were recognized at the CAP southern margin, around the city of Gülnar, at elevations between 1,000 and 1,500 m Öğretmen, Ferezza, et al., 2018). In the Gülnar area, the Pleistocene deep-marine deposits unconformably cover the Miocene limestones of the Mut Fm.  and Öğretmen, Ferezza, et al. (2018) dated the youngest marine deposits of the Gülnar area to the middle Pleistocene (Ionian, ca. 450 ka). Considering the present-day elevation and the paleodepth of these marine deposits,  estimated uplift rates around 3.2-3.4 m/kyr. This rapid uplift phase, which affected the CAP southern margin in the last 450 ka, was preceded by a short uplift phase (4.1 m/kyr) that affected the margin during the late Messinian (5.45-5.33 Ma, from Radeff et al., 2016Radeff et al., , 2017. Other recent studies based on paleoclimatic analysis and the study of offshore seismic sections fix the main topographic building of the plateau at the late Miocene (Meijers et al., 2018) or prior to the Pliocene (Fernández-Blanco et al., 2019).

Background of Marine Terraces
Marine terraces are essential geomorphic features that may provide crucial information about paleoclimate and surface deformation rates in coastal areas, in particular regarding vertical tectonic movements (Cosentino & Gliozzi, 1988;Dupré, 1984;Limber & Murray, 2011;Marquardt et al., 2004;Saillard et al., 2009;Trenhaile, 2002). Marine terraces are formed mostly during sea level highstand by the interaction between wave erosion and surface uplift. During highstands, wave erosion drives cliff retreat generating a rocky-shore platform etched on the landscape (Anderson et al., 1999;Bradley, 1957;Lajoie, 1986;Muhs et al., 1990). Rocky-shore platforms are characterized by smooth and low sloping surfaces often covered by shallow marine deposits (Anderson et al., 1999;Bowles & Cowgill, 2012;Merritts et al., 1991). Depending on the thickness of the sedimentary successions they can be divided into wave built and rocky-shore marine terraces. The wave-built marine terraces are gently sloping landforms formed by thick accumulation of shallow marine sediments (Jara-Muñoz & , whereas the rocky-shore marine terraces are landforms generated by wave erosion, comprising only a thin sedimentary cover. Marine

10.1029/2019TC005921
Tectonics terraces are usually abandoned during the following sea level lowstand period, if vertical displacements are faster and/or the next sea level highstand is lower than the previous one the platform can be preserved higher and safe from wave erosion.

Geomorphic Mapping and Field Observation
Fieldwork was focused in the southern part of the CAP margin, between Aydincik, Silifke, and Gülnar area ( Figure 3), characterized by well-preserved marine terraces. We documented the geomorphic characteristics  (Senel, 2002;Ulu, 2002).
of marine terraces and their associated deposits; in particular, we focused on discriminating between continental and marine terraced surfaces. We considered as marine terraces those surfaces characterized by a staircase morphology and characterized by the presence of marine sediments, fossils, and the exposure of ancient shorelines markers, such as lithodomes holes in the outcropping bedrock.

Morphometric Analysis of Marine Terraces
New advances of high-resolution topography have allowed observing geomorphic features like marine terraces with outstanding level of detail, allowing obtaining accurate measurements of terrace elevations that can be used to reconstruct past relative sea-level positions and estimating uplift rates. Morphometric analysis of marine terraces was carried out using 10 m resolution TanDEM-X digital elevation models developed by the German Aerospace Center and based on the twin imaging radar satellite Terra SAR-X. These data were validated using the differential GPS measurements before the morphometric analysis of marine terraces.

Surface Classification Method
The surface classification model (SCM) provides spatial estimation on the distribution of terraced surfaces and is based on the slope and roughness of the topography (Bowles & Cowgill, 2012). It is defined by the equation: where SLP is the topographic slope and the roughness, defined as the standard deviation of the slope calculated within a 3 * 3 roving window ( and Tol (Ionian) sections ; focus on the remnants of younger levels (b and c); marine terrace elevation histogram derived from the SCM (d); and shoreline-angle elevation plot (e).
roughness are normalized and linearly combined to generate the SCM isolating flat and smooth platform surfaces from the rough and steep areas that represent erosional features such as valleys and cliffs. Afterwards, false positive classifications such as valley floors, low dipping structural features, and fluvial terraces were manually removed.
We defined a maximum threshold of slope and roughness in the area. These values were used to normalize the distributions of the morphometric parameters as presented in equation 1. After isolating the flat surfaces, the SCM was intersected with the topography to obtain the elevation distributions of marine terraces. The different terrace levels were defined by differentiating elevation bands from the histograms and corroborated with our field observations to generate a map with the distribution of marine terraces in the study area.

Measuring Marine Terrace Elevations
To measure marine terrace elevation we use the shoreline angle, a geomorphic marker considered as an indicator of past sea-level positions that can be directly correlated with sea-level highstands (Ashby et al., 1987;Bradley & Griggs, 1976;Lajoie, 1986;Scott & Pinter, 2003). The shoreline angle is commonly covered by colluvial sediments and may be estimated from the intersection between the terrace paleo-platform and the paleo-cliff. To measure shoreline angles we use the software TerraceM (Jara-Muñoz et al., 2016, 2019) that uses swath profiles and linear regressions to estimate the shoreline angle position. We use the maximum distribution of elevations from swath profiles to reconstruct the original terrace morphology, devoid of the effects of river incision. Then linear regressions are applied upon selected segments of the paleo-platform and paleo-cliff and extrapolated, locating the shoreline angle at their intersection. Vertical errors are based on the extrapolation of the 2σ ranges from the linear regressions. For this analysis we used 170 swath profiles with a maximum width of 500 m considering the mean spacing of channels and gullies crossing through terrace risers and the resolution of the DEM. The profiles were oriented perpendicularly to paleo-cliffs.

Estimating Uplift Rates
The shoreline angle elevations can be directly correlated with sea-level highstands and used to estimate uplift rates using the equation 2 (e.g., Berryman et al., 1989;Weber, 1990): where E is the terrace elevation and e the sea level at the moment of marine terrace formation and T the age of the terrace (Lajoie, 1986). Nevertheless, uplift rates that include a correction for eustatic sea level might be biased by the uncertainties in the chosen sea-level curve (e.g., de Gelder et al., 2020;Pedoja et al., 2014;Yildirim et al., 2013). Beyond this simplistic approach to calculate uplift rates, temporal variations in uplift rates might produce cumulative vertical displacements that may be difficult to estimate using equation 2, as they depend on the vertical position of the markers before the change of uplift take place. In this regard, the application of LEM may help to account for complex uplift scenarios.

Landscape Evolution Model (LEM)
One of the main uncertainties in estimating uplift rates is the lack of adequate chronological constrain of marine terraces, which we address by using an LEM to model marine terrace ages. In the study area, the only known age comes from the higher terrace level that comprises the upper part of the Sarıkavak Formation dated as younger than 467 ka . However, a single age is not enough to determine if the uplift rate has been steady or variable through time. LEMs have been largely used to analyze the development of coastal landscapes under variation of sea level and uplift, simulating various processes for generating and degrading marine terraces including sediment transport and deposition, cliff retreat and diffusion, fluvial incision, and coral growth among others (e.g., Anderson et al., 1999;Hanks et al., 1984;Jara-Muñoz et al., 2017;Melnick, 2016;Nakamura & Nakamori, 2007;Refice et al., 2012;Shikakura, 2014;Storms & Swift, 2003;Thébaudeau et al., 2013). Based on the wave erosion and dissipation model of Anderson et al. (1999), our LEM simulates the formation of bedrock carved marine terraces under an oscillating sea level and constant or time-varying uplift rate. Recently, Melnick (2016) applied this model to estimate uplift rates along the central Andean coast.
The LEM used in this study is based on the equation of wave energy dissipation proposed by Sunumura (1992) and modified by Anderson et al. (1999), which considers the sea-bed erosion rate as a linear function of the rate of energy dissipation against the sea bed, as presented in the equation 3: 10.1029/2019TC005921 where dz/dt is the vertical sea bed erosion rate, dE/dt is the energy dissipation rate due to bathymetric drag, (dE/dt) 0 is the rate of wave energy dissipation in shallow waters, h is water depth of the local sea floor, and h* is the water depth at which the dissipation rate is 1/e in shallow waters.
The wave base, h wb , is the depth at which dissipation rate is essentially zero and is set as h* = h wb /4. Connecting this to climate signal, the dissipation scale h* is of the order of the wavelength. β is a constant that indicates the efficiency factor relating energy dissipation rate to erosion rate, here considered as 1.
For the simplest case of a planar shelf of slope Θ, the energy available in the wavefield to drive the cliff retreat ΔE is estimated as follows: where V is the component of wavespeed normal to coast.
The horizontal length scale is as follows: corresponding to the length of shelf over which most of the wave energy is dissipated. The expression indicates that regions with extensive shallow shelves should allow less energy to the cliff to drive erosion.
The rate of cliff retreat into the landmass, dx/dt, is assumed to be linearly proportional to the wave energy remaining at the instantaneous shoreline: where β cliff is a constant that folds together the fraction of the energy that is available for the erosion and the resistance of the local rock to the erosion. E cliff is defined by the equation: where E 0 is the original wave energy in the farfield wave climate.
The model input parameters are the slope of the shelf (assumed to be sub-planar), the uplift rate, the erosion rate, the wave height, the diffusion coefficient for the exposed cliff, and a sea-level history.

Comparing Terrace Morphology and LEM
To compare LEM results and the topography of marine terrace sequences, we simulate different scenarios of uplift and erosion rates in order to fit the elevation of measured and modeled shoreline angle elevations. We use a Random Optimization (RO) method that allows finding the best-fitting uplift history that describes the distribution of shoreline angle elevations. The RO uses random LEM uplift rates to generate sets of models, which are compared with measured terraces. This process is repeated 3 times, decreasing and restraining the number and range of random inputs until the fit between model and observations is minimized enough to reduce the uncertainty of our estimates.
We first selected three areas (Figure 8) where shorelines angle were well represented and at similar elevations. We test time-variable uplift rates obtained by bibliography and randomly created. The random pattern of variable uplift rates was generated by (1) fixing the time from 500 to 50 ka before present, (2) starting the model before the first highstand after the end of marine deposition (MIS 11c- , and (3) stopping the simulation immediately after the last highstand (MIS 3) before the present sea level rise in order to preserve the younger marine terrace levels, constrained by the Holocene uplift rate of 1-1.5 m/kyr estimated at Aydıncık-Silifke by Cosentino et al. (2016). We use three arbitrary time steps for variable uplift, between 90 and 110 ka (t2), 170 and 250 ka (t3), and 340 and 420 ka (t4). For each time step we define an arbitrary uplift range varying from 1 to 6 m/kyr based on previous mean estimates for the middle and upper Pleistocene . For the variable uplift scenario, we first created a series of combinations 800 of uplift patterns to feed the LEM; these different patterns were compared with the mean elevations of terrace levels in the three selected boxes to find the range

10.1029/2019TC005921
Tectonics of uplift patterns that best fit with the elevations of shoreline angles. We repeat the RO three times, including 600 uplift combinations in the second test and 400 in the third optimization in order to improve the fit between model and observations.

Marine Terrace and Shoreline Angle Morphology
We recognized seven levels of surfaces with a staircase morphology extending from the topmost part of the Gülnar area down to the Mediterranean coastline. We defined a maximum threshold of 17°and 12 for slope and roughness specifically for the study area, respectively, based on the morphometric characteristics of valleys, gullies, and other erosional features, to map the planar and sub-planar surfaces, then and we used a threefold criteria to define these surfaces as marine terraces: (1) geomorphic criteria: characteristic staircase morphology; (2) morphometric criteria: Due to the temporally pulsed nature of marine terrace formation, the platform surfaces should concentrate into discrete elevation bands if they are marine terraces (Bowles & Cowgill, 2012) (Figure 4d); (3) paleontologic criteria: the presence of remnants of marine organisms representing markers of past sea level (e.g., lithodomes). Based on the presence of these three elements we interpreted these surfaces as remnants of marine terraces (Figure 3).
Because of the prominent erosion, in particular for the middle levels, the terraces outcrop is scattered and forming partly dismantled remnants, making difficult direct mapping or survey. The high middle to late Pleistocene uplift rates  of the CAP southern margin (i.e., high erosion rate), together with lithological variability, seem to be the reason for the bad preservation of marine terraces. Therefore, we mapped shoreline angles in selected areas of best terrace preservation ( Figure 3). Our analysis revealed seven levels of marine terraces recognized in the field and in the digital elevation model. The distribution of marine terraces follows as SW-NE orientation, the width of the terrace levels decrease from top to bottom; for instance, the lower T7 level is represented by small remnants of abrasion surfaces set on the Tauride units near the city of Aydıncık, with an elevation of about 50 m a.s.l. (Figure 3b).
The uppermost level (T1) is the most extended marine surface and is placed on the top of the CAP southern margin (Figures 5a and 5b), with a mean elevation of 1,242 m a.s.l. (Figure 4). This level is characterized by a subhorizontal surface gently tilted eastward of about 0.5°, along a distance of 50 km, ranging from 1,400 m a.s.l. in the westernmost mapped area of the plateau margin top to 1,050 m a.s.l. east the city of Gülnar (Figure 3a). This level is etched on the Quaternary marly-calcarenitic units of the Sarıkavak Formation, in particular around the city of Gülnar, and occasionally on the shallow-water upper Miocene limestones of the Mut Formation Öğretmen, Ferezza, et al., 2018) (Figure 5a).
The boundary between the two formations is marked by sea level markers, as lithodome holes (Figure 5c), which indicate the presence of an ancient coastline (Figure 3a), similar as interpreted by . Considering the 1,500 m of elevation of the upper part of the Sarıkavak Formation, and the end of the sedimentation that occurred around 457 ka Öğretmen, Ferezza, et al., 2018), we relate this terrace level to the first sea level high-stand occurred after the end of marine sedimentation (MIS11-ca. 425 ka).
T2 (Figures 3 and 5a-5d) is located at a mean elevation of 722 m a.s.l. (Figure 4) and is etched on the Miocene limestones of the Mut Formation. The boundary between T1 and T2 surfaces is marked by a steep change in slope in the area of Akpinar Yaylasi-Bozągaç-Yenice (Figure 5a), which indicates the presence of an ancient and partly dismantled cliff; here the shoreline angles of T2 (SrA2) are located at 870 m a.s.l (Figure 4). The top of the surface shows no evidence of karstification and is characterized by a thin sediment coverture of pebbles; furthermore, fragments of ostreae, cardithes, echinoids, and pectinidae were found on top of the surface. Nevertheless, it is not yet clear if these fossils have been reworked from the upper Miocene limestones (Mut Formation) or if they belong to the marine incursion responsible for the T2 marine terrace formation.

Tectonics
The shoreline angle mapping for T1 and T2 is affected by strong errors because of the erosion processes that interested the paleo-platforms and the paleocliffs. In particular the mean elevation of the T1 shoreline angles results lower than the mean elevation of the T1 terrace ( Figure 4) principally because the TanDEM-X tile used for the swath profile analysis does not cover the western boundaries of the mapped paleosurface (Figure 3a), where mean elevation is higher. However, the position of the shoreline angles fits well with the inner boundaries of T1 and, in particular, T2, where the paleocliff is well preserved (Figure 5d), and the linear regression places the mean elevation of the shoreline angles lower than the platform elevation ( Figure 4).
From T3 to T6 ( Figure 3) the surfaces become progressively narrow, with an elongated shape and a trend parallel to the present coast (Figure 6a). The terraces are discontinuous, and west of Silifke, for a distance of about 15 km, it is possible to appreciate the staircase sequence of the terraces from the upper to the lower levels (Figures 3, 6a and 6b). The terraces are etched on the Miocene units and the deformed Tauride bedrock, lacking of terrace deposits at the top. The lowest level T7 is represented by very small and narrow

10.1029/2019TC005921
Tectonics platforms between 25 and 40 m a.s.l. These remnants are not included due to the scale of the map, but the position of the shoreline angles is in accordance with the parallel-to-coast direction (Figures 3 and 6).
Westward, the terraces are represented by poorly preserved platform remnants (Figures 7a-7c), carved on the deformed Tauride Units. The paleosurfaces cut the bedding of the bedrock (Figure 7b) and lack of sediments or fossils at the top. In the coastal areas of Aydıncık and Yenikaş (Figures 3c and 7a), T7 remnants are well-preserved, represented by narrow platforms parallel to the coast lacking of deposits atop.

Shoreline Angle Mapping and Marine Terrace Evolution Modeling
To estimate the Quaternary uplift trend we developed an LEM that simulate the evolution of marine terraces in the CAP southern margin. Because of the bad preservation of all terrace platforms, we used the mean elevation of shoreline angles in three sample areas (Boxes 1-3 in Figures 8a and 8b) to compare with the LEM. We first performed a sensitivity experiment to exploring the contribution of initial erosion rate (E 0 ), coastal geometry, the sea-level curve, and uplift rate (Ur) on surface morphology. We compared E 0 values between 0.1 and 0.9 m/yr (see sensitivity test in supporting information). Considering that the variability of synthetic staircase morphologies for different E 0 is minimum, we used an arbitrary E 0 of 0.3 m/yr. We calculated the mean coast slope (10°) using swath profiles. The wave height estimation is based on satellite altimetry timeseries (the altimeter products were produced by Ssalto/Duacs and distributed by Aviso+, with support from Cnes https://www.aviso.altimetry.fr); we used the 95% tail of the distribution of wave heigths (2.7 m) that represent waves capable to produce enhanced erosion and cliff retreat.
We run the LEM using the sea-level curves of Bintanja et al. (2005) and Rohling et al. (2014) because they cover the analyzed time span and are based on northern hemisphere proxies, used the random-generated uplift sets (Figure 8c), and tested different E 0 , to obtain the best combination of parameters that replicates the evolution of marine terraces at the CAP southern margin (see supporting information).
For each model we calculated the Root Mean Square deviation (RMS); we arbitrarily selected the five uplift trends that generate shoreline angles with the lower RMS and the higher number of shoreline angle fitted with the mapped ones (Figure 9).

Tectonics
The results obtained by the sea-level curve of Bintanja et al. (2005) generates the higher number of fitted shoreline angles and the lower values of RMS ( Figure 9); just for the Box 2 the RMS obtained for the curve of Rohling et al. (2014) is lower, but it is lower also for the number of fitted shoreline angles. Despite the difference between the two simulations, the general trend obtained for each analyzed area is very similar and describes a strong increase in uplift occurred from 500 to 200 ka, with uplift rates that accelerate from 1.7-2.0 to 3.4-3.8 m/kyr, followed by an abrupt and fast decrease until a range of 1.6-1.2 m/kyr at present days ( Figure 10).

Platform Development and Morphostructural Setting of the CAP Southern Margin
Despite the scarcity of sedimentary elements that suggest a direct correlation between some of the CAP southern margin paleosurfaces to marine terraces, their morphostratigraphical relationships (i.e., staircase

Uplift Trend and Tectonic Evolution
As evidenced by , the reconstruction at the CAP southern margin (Gülnar area) of a middle Pleistocene marine shoreline at 1,500 m elevation suggests a strong uplift phase in the last 450 kyr, with mean uplift rates of 3.2-3.4 m/kyr. Nevertheless, along the Mediterranean coastal area of southern Anatolia, markers of ancient sea level have been dated through AMS 14 C dating and U-Th chronology suggesting late Holocene uplift rates around 1-1.2 m/kyr , significantly lower than the middle Pleistocene uplift rates estimated at Gülnar.
Along the CAP southern margin, where the timing of marine paleo-platforms is poorly constrained, the use of an LEM allowed to estimate an uplift history during the last 500 kyr. Mapping shoreline angles to be used as elevation benchmarks and starting from the first highstand (MIS 11, ca. 425 ka) after the youngest deposits of the Sarıkavak Formation found in the Gülnar area (<467 ka), which we assume as the highstand responsible for the origin of the T1 paleosurface, we were able to find the best uplift trend that generates models of marine terraces where the shoreline angle elevations fit with the mapped ones (Figures 11a and 11b).
The first important result of our analysis is the compatibility between the obtained uplift trend and the mean uplift rates estimated by previous authors for the study area by biostratigraphic studies on the Quaternary marine deposits Öğretmen, Ferezza, et al., 2018). The results of the LEM describe bell-shaped uplift trends that increase from 1.84-1.95 m/kyr at 500 ka to 3.58-3.74 m/kyr at 200 ka and then decrease until 1.37-1.35 m/kyr for the present days ( Figure 11b). These data provide a high-resolution uplift trend for the CAP southern margin during the last 500 kyr, matching with estimates by  for the middle Pleistocene marine deposits and by Cosentino et al. (2016) for the late Holocene between Aydıncık and Ayaş, and relates the terrace formation to the main sea-level high stands occurred in the last 500 ky.
This high-resolution uplift trend for the Pleistocene phase of uplift at the CAP southern margin is at odds with recent reconstructions that place the uplift phase prior to the Pliocene (Fernández-Blanco et al., 2019) or during the late Miocene (Meijers et al., 2018). Both these reconstructions do not consider the presence of middle Pleistocene marine deposits up to 1,500 m of elevation in the Gülnar area Öğretmen, Ferezza, et al., 2018). The Pleistocene high-resolution uplift trend of the CAP southern margin, as provided in this paper, is consistent with the surface uplift signals of thermo-mechanical models for experiments of subduction-collision systems (Duretz & Gerya, 2013). In those models, after ca. 15 Myr of continental crust subduction, slab detachment takes place and triggers the sharpest topographic and short-lived surface uplift signal of the experiment, with uplift rates from 1 to 5 m/kyr (Duretz & Gerya, 2013). After the slab detachment, the model shows a rapid decrease of surface uplift, with average uplift rates on the order of 0.2-1 m/kyr. Therefore, the Pleistocene high-resolution uplift trend for the CAP southern margin, presented in this work, is in agreement with the more recent slab break-off model proposed for the Quaternary uplift history of the CAP southern margin  and is at odds with model that consider the plateau margin developed prior to the Pliocene by shortening led by Eurasia-Africa compression (Fernández-Blanco et al., 2019). Following the slab break-off model, and according to our results, the slab detachment of the subducting African-Arabian plate occurred in a time window between 400 and 150 ka, with a peak uplift rate of 3.6-3.7 m/kyr at around 200 ka.

Conclusions
We estimated an uplift history for the CAP southern margin during the past~500 kyr using marine terraces and an LEM. By considering the age of the youngest marine deposits at Gülnar , the staircase morphology of the terraced surfaces reaching 1,400 m, and an LEM that reproduced the evolution of marine terraces, we reconstruct a morphostratigraphic model for the marine terraces at the CAP southern margin. We correlate terrace sequences with the sea-level highstands of the past 500 kyr. Our results bridge caps in previous studies at the Holocene and middle Pleistocene timescales allowing to reconstruct a detailed history of uplift rates along the CAP southern margin, finding that, 1. The CAP southern margin is characterized by a staircase morphology defined by the presence of at least seven marine terraces (T1-T7), from 1,400 m down to 30 m of elevation, due to the interplay between 10.1029/2019TC005921 Tectonics sea-level changes and tectonic uplift. All of them developed after the deposition of the youngest marine deposits (late middle Pleistocene) recently recognized in the Gülnar area Öğretmen, Ferezza, et al., 2018), which were subsequently uplifted at 1,400-1,500 m of elevation. 2. The LEM allows us to correlate the T1-T7 marine terraces that were recognized in the field and through the analysis of a 12 m resolution TanDEM-X digital elevation model, with highstands that occurred from MIS 11 (ca. 425 ka) to MIS 3 (ca. 50 ka). 3. Applying these morphostratigraphical results and the mean elevations of the shoreline angle of each of the T1-T7 marine terraces, we were able to develop an LEM that provided a high-resolution Quaternary uplift pattern showing how the uplift rates varied from a maximum of 3.4-3.8 m/kyr, until 1.6-1.1 m/kyr for the present days. 4. The resulting uplift history that affected the CAP southern margin during the late Quaternary describes a bell-shaped temporal evolution, characterized by a strong increase and decrease of uplift rates in a very short time window (<500 kyr). This result is consistent with the peak of uplift rate that thermo-mechanical models of continental collision show in correspondence with the slab break-off of the subducting plate (Duretz & Gerya, 2013). This observation strengthens further the relationship between the late Quaternary uplift of the CAP southern margin and the geodynamics of the plate boundaries at the convergence between Africa, Arabia, and Eurasia, which underwent continental collision and possibly slab break-off. 5. Following this geodynamical model to explain the late Quaternary uplift of the CAP southern margin, the uplift pattern from the LEM analysis suggests as the break-off of the subducting slab occurred very recently, in a time window between 400 and 150 ka. 6. The LEM results indicate as the present days topography at the CAP southern margin has been acquired in the last 500 kyr and not during the late Miocene-Pliocene as suggested by Meijers et al. (2018) and Fernández-Blanco et al. (2019), with the maximum uplift rates between 350 and 150 ka. 7. The pattern of the late Quaternary uplift rates for the CAP southern margin and Cyprus, as evidenced by this work and that from Palamakumbura et al. (2016) for the Kyrenia Range in Cyprus, shows a faster uplift for the CAP southern margin with respect to the Cilicia Basin and Cyprus. This differential uplift has been accommodated by the late Quaternary activity of the major fault zones affecting the CAP southern margin (e.g., Ecemiş Fault Zone, Kozan Fault Zone, and Kirkkavak Fault Zone).

Data Availability Statement
TanDEM-X data were provided by the German Aerospace Center (DLR) under TanDEM-X Science Proposal IDEM_GEOL0119 via the DLR science portal (https://tandemx-science.dlr.de/cgi-bin/wcm.pl?page=TDM-Proposal-Submission-Procedure).