Geomorphology of a recently deglaciated high mountain area in Eastern Anatolia (Turkey)

ABSTRACT We present a geomorphological map of the Cilo Mountain range located in southeast Turkey that illustrates the recent evolution of the second-most glacierized area of Turkey. The map was produced by the manual delineating of landforms from highresolution satellite imagery (Pleiades and Google EarthTM images). Cilo Mountain glaciation played a crucial role in the evolution of the local physical landscape and even today glaciers are significant features. After the termination of Little Ice Age, deglaciation promoted the onset of multiple and intense surface processes; widespread slope instabilities, due to permafrost melting and frost weathering processes, are still the most active processes. Such processes are further promoted by the scarcity of vegetation. The occurrence of several proglacial lakes significantly contributes to the reduction of sediment output from glaciated watersheds, influencing the evolution of proglacial plains and resulting in a variety of landforms that we identified and described in this work.


Introduction
Since the end of the Little Ice Age (LIA), glaciers have retreated worldwide and their terminus positions may have moved upward several kilometres in a relatively short time (Allard et al., 2021;Matthews & Briffa, 2005). The LIA had a global extent: major glacial advances have been registered between the 14th and 17th centuries and a further cold period marked by glacial advance at ca. 1820 and 1855-1860 CE, followed by minor readvances in the 1890s (Holzhauser et al., 2005;Hughes, 2014Hughes, , 2018Le Roy et al., 2017;Matthews & Briffa, 2005;Schimmelpfennig et al., 2012Schimmelpfennig et al., , 2014Styllas et al., 2018). The retreat of glaciers is associated with major (natural and human-induced) climatic changes impacting on the global cryosphere (e.g. Haeberli et al., 2016). Even if the post-LIA retreating phase shows a global trend, a synchronicity in short-timed readvancing phases of alpine glaciers is recognized in the twentieth century, when glaciers regionally experienced minor tongue advances or stillstands in the 1920s, 1970s, 1980s and 1990s (Solomina et al., 2016. In the field, the LIA maximum is commonly recognized by the occurrence of well-preserved moraines, and moraine systems also document post-LIA minor glacier fluctuations (Solomina et al., 2016). The retreat of mountain glaciers also exposes new land surfaces, where new surface processes exist (Carrivick et al., 2018). Most of the latter are still related to the dynamics of glaciers, but they are rapidly associated and/or substituted by water-related, gravity-related, and periglacial processes, shifting into the compound of surface processes defined as paraglacial (Ballantyne, 2002;Knight & Harrison, 2009;Palacios et al., 2021). Such processes affect glacier forelands and surface areas are progressively exposed (e.g. D' Agata et al., 2020). The landforms of glacier forelands are difficult to be represented in a geomorphological map, because of their dynamic nature; they are often characterized by ephemeral duration, and spatial variability (Leigh et al., 2021).
Moreover, deglaciation allows the colonization of organisms and the onset of soil-forming processes, thus promoting the evolution of progressively more complex ecological networks, including feedbacks (Eichel, 2019;Ficetola et al., 2021;Khedim et al., 2021). Consequently, the establishment and development of bacteria, fungi, plant and animal communities and their interactions with geomorphological processes and geological settings of the landscape become crucial in tuning the evolution of recently deglaciated areas (Eichel et al., 2016Franzetti et al., 2020Garavaglia et al., 2010;Moreau et al., 2008;. Several studies have investigated evidence of late-Holocene glaciation in the Mediterranean basin and Eastern Africa (e.g. Hughes, 2014Hughes, , p. 2018Altınay et al., 2020;Gachev, 2020;Allard et al., 2021;Groos et al., 2021), but detailed studies reporting the post-LIA evolution of proglacial areas in the Mediterranean mountain ranges are not available. One of the most interesting and poorly investigated recently deglaciated sectors of the Mediterranean domain is the Cilo Mountain range, located in southeast Turkey. Herein, 29 ice bodies are present, covering 3.82 km 2 , making this region the second-most glacierized area of Turkey after the Mount Ağrı Dağı/Ararat district (Azzoni et al., 2017Sarikaya & Tekeli, 2014).
Cilo Mountain is characterized by several small cirque glaciers and is representative of Mediterranean glaciations, with only small, high elevation and north facing ice bodies. Moreover, the glacierized area around Uludoruk (4136 m asl), the highest peak of the Cilo mountain range, is one of the most interesting high-elevation circum-Mediterranean glaciated areas due to the density and size of the ice bodies present. Geomorphological investigation in the area offers the opportunity of contributing to the reconstruction of the dynamics of Mediterranean glaciers post-LIA deglaciation. For these reasons, we (i) produced the first geomorphological map (Main Map) of the Cilo Mountain, using a remote sensing approach, with particular emphasis on glacial and periglacial landforms related to geomorphic processes typical of the paraglacial stage. Moreover, our map serves as a reference to (ii) reconstruct the characteristics and the evolution of proglacial areas in this specific Mediterranean climatic context.

General setting of the study area
The study area is in the southeast of Turkey (37.3°N, 44°E), in the Hakkari Province at about 15 km from the border with Iraq ( Figure 1). It covers more than 50 km 2 and includes the recently deglaciated areas of the Cilo mountain range that contains three glaciers: the Erinç Glacier, the Avaspi Glacier, and the Uludoruk Glacier. The highest peak of the region, the Uludoruk Dağı (4168 m), is the second highest place in Turkey, after Ağrı Dağı/Ararat (5137 m asl).
The climate of the area is continental. Data from the nearest meteorological station located at Hakkari (20 km from the study area and at 1000 m lower elevation) reports annual total mean precipitation of ca. 744 mm and mean temperature of ca. 9. 5°C (1975-2014) (Satir, 2016).
Considering the geological setting of the area, the Cilo mountain range is located at the southeastern limit of the Taurus Mountain. The area forms an anticlinal system consisting of series of folded arches frequently broken and much deformed, and composed of sedimentary rocks, deposited in the former Tethys Sea (Güldali, 1979). In particular, the geology of the Cilo mountain range consists of Paleozoic and Mesozoic limestones and clastic strata that were deposited along a passive margin on the southern edge of the Arabian microplate (Hanilçi et al., 2020). During the Eocene and the Miocene, basement units were folded and thrusted producing the present day steep mountainous terrain. The study area is characterized by a carbonate bedrock belonging to the Triassic, Jurassic and Cretaceous Cudi Group, and Cretaceous-Eocene limestones belonging to the Şırnak Group.

Previous glaciological and geomorphological studies
The Cilo Mountain range has been investigated by numerous glaciological studies aimed at describing the size and characteristics of ice bodies. The first data were reported by Bobek (1940) who recorded the minimum elevation of local glacier snouts in 1937 at ca. 2600 m asl. Erinç (1952) identified more than 20 glaciers and observed glacier snout altitudes at 2900 m in 1948, suggesting a major frontal retreat of ca. 500 m from the position occupied by the ice according to LIA terminal moraines. Kurter and Sungur (1980) estimated that glaciers of the region covered about 8 km 2 . Kurter (1991) positioned the glacier snout altitude at ca. 3000 m asl. using remote-sensing detection techniques. The increased availability of remote-sensing data led to the compilation of the first atlas of Turkish glaciers (Williams & Ferrigno, 1991), mostly based on a dataset of Landsat imageries between 1975 and 1980. Sarikaya and Tekeli (2014) developed an updated inventory of Turkish glacier identifying 10 glaciers on Cilo Mountain range. Hughes (2014) reported more than 20 glaciers in this area, accounting for two thirds of all the glaciers in Turkey. Yavaşli et al. (2015) investigated the evolution of Turkish glaciation from the 1970s to 2011-2012 and established that in the 1970s glaciers of the Uludoruk area covered 8.5, 5.5 km 2 in the 1980s, 4.1 km 2 in the 1990s, 3.90 km 2 in 2011 and 2013, and 3.51 km 2 in 2011-2012, thus losing almost the 60% of its area in only 40 years and increasing the area of deglaciated ground. Azzoni et al. (2020) in their high-resolution Turkish glacier inventory identified 29 ice bodies in the study area, covering 3.82 km 2 .
Despite the abundance of glaciological data available for the Cilo mountain, very few geomorphological investigations have been performed in the area and most of the available data correspond to the positions and elevations of glacier snouts only.

Methods
Several examples of geomorphological mapping in deglaciated areas have been produced, especially in Alpine (e.g. Bollati et al., 2018;Buter et al., 2020) or high latitude environments (e.g. Everest et al., 2017). Such maps have been realized following different approaches based on field surveys and/or remote sensing analyses (Chandler et al., 2018).
However, the dynamic nature of proglacial areas requires a periodic updating of elaborated maps (Ballantyne, 2002). We used high-resolution and multitemporal remote-sensing data to manually delineate the main geomorphological features of the study area. The Cilo Mountain range is in a semi-arid region, characterized by the complete absence of arboreal vegetation and by a scarce herbaceous cover. This guarantees very good landscape visibility and easily permits the investigation of landforms with satellite imagery in the visible spectrum (Azzoni et al., 2019;Carrivick et al., 2018;Forti et al., 2021). Landforms are classified based on their genetic processes. The legend of the Main Map has been developed on the basis of the guidelines proposed by the Italian Institute for Environmental Protection and Research (ISPRA) (Campobasso et al., 2018) and some modifications have been made to accommodate the methodology in relation of the spatial scale of the map. In this way erosional and depositional landforms are grouped in relation with their genetic process (e.g. gravity-driven landforms are characterized by a red color, glacial ones by a violet color, etc.). Moreover, color intensity refers to the degree of activity of each process (i.e. bright colors refer to active landforms whereas the same color with a darker tonality refer to inactive landforms).
We chose high-resolution Pleiades (spatial resolution: 0.5 m) and Google Earth™ (spatial resolution: ∼2 m) imagery for recognizing and mapping geomorphological units. We used the Pleiades image set TPP1600878788 obtained on 10th of September 2017: this tile was completely cloud-free and with a scarce snow cover, favoring the identification of landforms. Google Earth™ images were collected for September 2006, July 2013, June 2015 and July 2018, with the same sky and snow conditions of the Pleiades frame and were used for checking some areas.
Our geomorphological analyses also took advantage of Digital Elevation Models (DEMs), which allowed better evaluation of both complexity and geometry of the landscape and the computation of 100 m spaced contour lines (vertical interval). In detail, we used ALOS Global Digital Surface Model with a spatial resolution of 30 m and a vertical accuracy of 15 m (2017 data). Contour lines at 100 m derived from this model were used for landform interpretation as well as the related hillshade model.
The bedrock geology and the main structural lineaments have been obtained from the web-GIS of the General Directorate of Mineral Research and Exploration of Turkey (data from Duman et al., 2011) and are simplified using the categories of limestone and radiolarite rocky outcrops Local toponyms were obtained from maps of hiking trails of the Cilo mountain range. Finally, the representation of the glacier surface relies on data elaborated by Azzoni et al. (2020). Although the main geomorphological features are largely controlled by structural and lithological characteristics of the bedrock, structural landforms are difficult to recognize and are mapped from remote-sensed data, and landforms associated with physical-chemical weathering. Consequently, these features have not been detailed in this map, where the major focus is on deglaciation processes. Moreover, the lack of a detailed field survey has prevented an unambiguous description of some deposits formed by the combined action of different processes (glacial, fluvioglacial and gravity-related) that are therefore mapped and classified as stabilized polygenic deposits.

Glaciers and glacial landforms
On Cilo Mountain range we mainly investigated glacial landforms to deepen our knowledge about the characteristics and timing of deglaciation and glacier foreland dynamics.
Considering 2017 data , we found 26 ice bodies, covering 3.36 km 2 , only seven of which have a surface area greater than 0.1 km 2 , whereas two ice bodies, smaller than 0.01 km 2 , are below the general threshold considered to define them as glaciers (Paul et al., 2020). The slight difference between fresh data and those reported in the 2020 inventory of Turkish glaciers  is due to the occurrence of three small ice bodies that fall outside the map. Among the proper ice bodies, seven can be classified as mountain glaciers whereas the other ones are glacierets. The latter are characterized by their small dimensions, wide supraglacial debris coverage and the lack of a clear distinction between accumulation and ablation zones. The largest ice bodies have a wide debris cover and a scarcity of snow cover during the summer season, suggesting a negative mass balance (Zemp et al., 2009). Moreover, most glacier features (e.g. a main northward aspect, lack of crevassed areas, wide debris cover, small accumulation areas) suggest that glaciers are undergoing retreat, as testified by their continuous negative mass balance.
Numerous snow fields are observed in the study area, especially facing northward at the base of steep slopes or in incised valleys. The analysis of multitemporal satellite images from Google Earth TM permits classification as semi-permanent snowfields, due to their persistence during the warmest months of the year. Their occurrence and persistence mostly depend on the high accumulation of winter snow and low summer air temperatures.
Moraine ridges and glacial deposits were identified and mapped. Moraine ridges are clearly visible and well preserved despite active erosional processes on their flank, inducing pseudo-badland landforms typical of such environments (e.g. Bollati et al., 2017;. Based on available data (Hughes, 2014;Yavaşli et al., 2015) and observations on elevation, aspect and size, we can attribute moraine ridges to the last LIA glacier advance.
Moraine ridges were identified in the proximity of present glacierized areas with a maximum distance of ca. 1 km from the present-day glacier terminus (Erinç Glacier) (Figure 2). We have identified 12 moraine systems in correspondence with glaciers no longer present. Four of these are south facing: the positions, elevations and characteristics of moraine ridges are comparable with the systems facing northward. This suggests the existence of former glaciers sustained by substantially different climatic conditions during the LIA. To the best of our knowledge, these are the only LIA moraines of the Mediterranean area that are facing southward (Hughes, 2014). Considering the position of moraine ridges we can roughly estimate that glaciers covered about 6 km 2 in the midnineteenth century that corresponds to the apogee of the LIA.

Periglacial and cryonival landforms
According to Hughes and Woodward (2009), periglacial activity and permafrost occurrence in the area likely played a major role in the evolution of landforms but their effects have not been studied in detail. We noticed a high availability of debris due to intense frost weathering and the lack of plant cover.
Along the north western sector of the study area, three north-facing rock glaciers can be observed; they cover 0.15, 0.13, and 0.03 km 2 (Figure 3). Remotely-sensed images suggest that rock glaciers are lobate with a hummocky surface. At the time of image acquisition, meltwater drainage from the frontal part is observable in correspondence to the snout of each rock glacier. For such reasons, we can propose that identified rock glaciers are still active with widespread presence of buried ice, sufficient to induce internal deformation, creep and down valley movement (Colucci et al., 2016;Cremonese et al., 2011).
Avalanches probably also play an important role in shaping the morphology of this area: massive snow movements can affect debris flow channels during the winter season, however confirmation of these processes can only be done through a field survey.

Gravity-driven landforms
The structural and lithological characteristics of the bedrock, intense frost weathering and the lack of vegetation promote wide slope instability: consequently, few parts of the study area are identified as stabilized surfaces.
Debris flow channels are deeply cut and incise deposits and their occurrence is consistent with the rainfall regime of the region, that is characterized by episodic, high concentrated precipitation events (Satir, 2016).
In the study area, debris flows fans are also common. Debris flow channels and fans, deeply investigated with multi-temporal images for assessing changes in their shape, can be classified as active.
Surface instability is clearly deducible along steep slopes at the head of glaciers. The debuttressing  process affecting rock walls (Ballantyne, 2002), after reduction of ice volume and lowering of the glacier surface, largely contributes to the production of debris that deeply modified the glacier surface characteristics.
In the south eastern sector of the study area, we mapped a major landslide (Figure 4), which is heavily affected by stream water reworking and is partly buried by debris flow fans. The identification of the landslide niche and its track is not possible based on the  satellite images. Such elements can have a role in terms of sediment connectivity, interrupting water flow towards the lower valley (Korup, 2005).

Hydrography and water-related landforms
The hydrography of the Cilo Mountain range is strongly influenced by the topography and structural setting of the area. The hydrographic network of the area is mainly fed by glaciers, semi-permanent snowfields and in part by rock glaciers. All these contributions provide important water flow during the melt season. Short and ephemeral water courses, classified thanks to the multi-temporal images, characterize the steep flanks of the valleys, whereas more organized and braided streams are present at the bottom of valleys. Consequently, fluvial deposits are uncommon and limited to some parts of valley floors. Moreover, topographic characteristics of the area (i.e. steep rock outcrops) prevent the deposition of fluvioglacial sediment that is found only in the proximity of ice bodies.
Local hydrography is also characterized by lakes. We mapped 22 ponds of different sizes. Among these, six are ice-contact lakes (see Figure 2), whose dimensions vary largely according to the season of observation, whereas three are moraine-dammed lakes. The others are bedrock-damned lakes occupying glacial overdeepenings (Otto, 2019). It interesting to note that the effects of linear erosion in proglacial areas seems to be more intense where lakes are absent. The occurrence of (semi)permanent water ponds significantly reduces sediment output from glaciated watersheds and the energy and velocity of meltwater, as reported by Tweed and Carrivick (2015), with a possible impact on the evolution of proglacial plains (Otto, 2019) ( Figure 5).
The dynamic nature of the study area is also highlighted by diffuse and intense rill erosion processes active along proglacial plains and moraines where the scarcity of continuous plant cover prevents stabilization and sediment mobilization becomes significant (Cossart & Fort, 2008). The presence of prominent stream terrace escarpments on the sides of the water course fed by melt waters of the İzbırak Glacier indicates that erosional processes are evident in these areas.

Conclusions
Our first overview on the main geomorphological features of the Cilo Mountain region highlights that a variety of active geomorphological processes shape the local landscape with particular landforms that are presented in the map. Glacial processes prevail at the highest elevations, where recently deglaciated areas are present. At lower altitudes and along the steepest sector of the study area, gravity-related processes, including debris flows, modulate landscape evolution. Along the valley floors, river processes are dominant. In the study area, glacier shrinkage is intense: large terminus retreat since the end of the LIA exposed new land surfaces, where erosional processes are diffuse. Moraine flank instability, rill erosion and fluvioglacial processes shape the proglacial area inside the LIA moraine ridges where indicators of paraglacial readjustment are present.
The case study of the Cilo mountain offers a fresh look on geomorphological processes occurring after the LIA glacial advance in the Mediterranean domain; such processes are widespread in the area but in many cases are overlooked. Understanding recent deglaciation dynamics in such environmental contexts is crucial to forecast the effects of ongoing climate change.
Software Esri ArcGIS 10.2 was used to produce the geomorphological map and specifically the remote-sensing analysis, development, and production of the DEM, drafting of the geomorphological units, and final map layout. Google Earth TM was used for further ground controls of the limits between geomorphological units.

Acknowledgments
The Pleiades data were provided by the European Space Agency (ESA, project: Geomorphological mapping and recent glacier evolution of the Mount Ararat volcanic complex through SPOT and PLEIADES images; ID-32011, project leader RSA). Part of this research was supported by the Italian Ministry of Education, University, and Research (MIUR) through the project "Dipartimenti di Eccellenza 2018-2022" (WP4 -Risorse del Patrimonio Culturale) awarded to the Dipartimento di Scienze della Terra "A. Desio" of the Università degli Studi di Milano. We acknowledge the support for publication fees from the University of Milan through the APC initiative and from grant CTE_-NAZPR21AZERB_01 (to AZ).

Disclosure statement
No potential conflict of interest was reported by the author (s).

Data availability statement
Data available upon request.