A glimpse into the northernmost thermo-erosion gullies in Svalbard archipelago and their implications for Arctic cultural heritage

Gully erosion is one of the most destructive geomorphological processes on relatively flat surfaces. This is exacerbated in the Arctic regions, where gullies are referred to as thermo-erosion gullies because of their unique connection to permafrost. As the surface of the permafrost freezes and thaws

Permafrost is ground that has remained at or below 0 degreesC for two or more consecutive years (Schuur et al. 2015). The seasonal temperature fluctuations only induce thaw of the uppermost soil-layer, called the active layer, while the underlying soils and bedrock remain frozen. The current rapidly warming climate induces increasingly thicker active layers, as well as an increase in the temperature of the underlying permafrost, which, if passing 0 degreesC, is no longer permanently frozen (Obu et al., 2019). Within the Arctic and sub-Arctic areas, gullies are known as thermo-erosion gullies. A thermo-erosion gully may initiate when heat transfer from small water tracks penetrating the surrounding soils and increasing the depth of the active layer locally. Below the active layer the ground remains permanently frozen with normally a very high ice content (up to over 90% in some circumstances but often around 75%) in the upper part of the permanently frozen part, called the transition zone (Iwahana et al., 2014). This ice-rich transition zone will, if thawed, release such excess of water that it may initiate both small scale fluvial processes and small slumps or grain collapses, if occurring in any kind of slope. Warmer or longer melting seasons may cause increasingly deeper active layers, thawing the icerich transition layer sediments and causing oversaturated soils. If the surface soils are removed by small scale gullying and/or slumping, new ice-rich soils will be exposed to thawing, initiating a feed-back loop of seasonal gully expansion. The geomorphological process may thus start humbly, for example during an unusual warm spell or prolonged thawing season one year, but as soon as a depression is formed, more gravitational slumping may accelerate the process rapidly. This, together with mechanical erosion and gravitational processes causes general soil subsidence and channel incision through thermo-erosion processes and results in a mixing of soil horizons and continued erosion of adjacent soils (Kokelj and Jorgenson, 2013). Thermo-erosion gullies enlarge both retrogressively upslope and through deepening and widening of the initial incision.
The rapidly increasing trend in temperature and precipitations since the 1970 s has initiated general permafrost degradation in many Arctic regions (Biskaborn et al., 2019), and has strong influences on landforms, ecosystems, geomorphological processes (Berthling et al., 2020), infrastructure (Hjort et al., 2018;Karjalainen et al., 2019), and society (Ford et al., 2021). The warming trend in the Arctic has been twice as rapid as the global average (Bintanja, 2018) and leads to increased and varied thermokarst activity: thawing ice-wedges, degrading pingos, thaw slumps, thermo-erosion gullies, coastal erosion, etc. These processes are also known as cryospheric hazards, which may pose a danger for both society (Ding et al., 2021) and cultural heritage (Ljungqvist et al., 2020). Svalbard is experiencing amplified climate change when compared to the global average (van Pelt et al., 2019), and represents a hot-spot for scientists to focus on: land cover and ice-wedge mapping (Bartsch et al., 2016;Lousada et al., 2018), the surface morphology of fans (De Haas et al., 2015), modelling the velocity of glaciers and ice caps (Strozzi et al., 2017), plastic pollution (Collard et al., 2021), glacier retreat (Kociuba et al., 2021) and coastal erosion (Jaskólski et al., 2018;Zagórski et al., 2020;Nicu et al., 2020;Nicu et al., 2021a).
With a scenario of continued increase in greenhouse gas emissions, climate change impacts for Arctic Norway and Svalbard towards the year 2100 is predicted to result in significant increase of mean annual air temperature (MAAT) and precipitation, and general decrease of snow, glaciers, sea ice, and permafrost across the whole latitude spectrum (Hanssen-Bauer et al., 2017). A few studies have approached thermoerosion gullying in cold regions: Nunavut -Canada (Fortier et al., 2007;Godin and Fortier, 2012;Godin et al., 2014;Perreault et al., 2016Perreault et al., , 2017, Alaska -USA (Tape et al., 2011;Lara et al, 2019), Yamal Peninsula -Russia (Sidorchuk and Matveeva, 2020;Sidorchuk, 2020), Greenland (Pastor et al., 2020;Christensen et al., 2021). Thermo-erosion gullies are of additional interest because through their evolution, they contribute to the release of vast amounts of methane (CH4, a greenhouse gas), stored in the active layer. Therefore, studying thermo-erosion gullying also represents a base to identify source areas where a higher level of methane is being released into the atmosphere (Oberle et al., 2019). With all the efforts being made by scientists, the Arctic region of Svalbard has been ignored in regard to the thermo-erosion gullying process.
The general effects of climate change on global cultural heritages have been acknowledged and many studies tackled this issue (Sesana et al., 2021). While there are previous studies of gully erosion effects on cultural heritage at a global level (Nicu, 2018;Ciampalini et al., 2019;Nicu, 2019;Lombardo et al., 2020), there are no studies that approach the effects of thermo-erosion gullying on Arctic cultural heritage; even though about 180.000 archaeological sites are registered in the Arctic (Hollesen et al., 2018). Climate change-induced processes affecting cultural heritage (Ford et al., 2006) is even less studied in Svalbard; with the only studies being very recent on thaw slumping (climate-induced landslides; see Nicu et al., 2021b), coastal erosion 2021a) and the decay of wooden features (Mattson et al., 2010. In contrast, anthropogenic influence on cultural heritage in Svalbard, such as the Arctic tourism, has been studied in more detail (Hagen et al., 2012;Thuestad et al., 2015;Jaskólski, 2021).
The present study aims to i) create an inventory and have a general overview of thermo-erosion gullying in Nordenskiöld Land (Svalbard); ii) identify the main controlling factors of thermo-erosion gullying; iii) present the monitoring results of small-scale fluvial gullying; iv) analyse the implications of thermo-erosion gullying processes in regard to Arctic cultural heritage sites and a warming climate. All the above points will tackle the mentioned issues for the first time. The results of this study will shed light on the environmental factors controlling thermo-erosion gullying and may be used to evaluate the present state of cultural heritage sites as well as cultural heritage assessment and management. Local authorities and stakeholders in Svalbard are the main end-users of the results and will benefit in prioritising any future mitigation measures in the fast-changing Arctic landscape. Also, this contribution provides the first inventory of thermo-erosion gullying, creating the potential to monitor future development and speed of change, which is of interest not only for Arctic cultural heritage preservation, but for permafrost geomorphology and Arctic climate change in general. Moreover, knowledge of the factors controlling thermo-erosion gullying and their interactions is highly significant for developing effective gully management strategies.

Study area
The study area is located in central Spitsbergen (Fig. 1), which is the largest island of the Svalbard archipelago (governed by Norway and established by the Spitsbergen Treaty from 9 February 1920); Svalbard archipelago, with an area of approximately 61020 km 2 is located ~ 1100 km south of the North Pole and ~ 800 km north of the coast of Norway (Zwoliński et al., 2013).
As a consequence of its location, Svalbard forms one of the most important and strategic terrestrial nodes on Earth, separating the Greenland Sea, the Barents Sea, and the Arctic Ocean (Jaskólski et al., 2018). Svalbard also represents one of the most diverse geological features in the world, where sections representing most of the Earth's history are accessible. Another significant characteristic of Svalbard's geology is the presence of terrestrially exposed sedimentary successions that are rare or do not exist in other places in northern Europe (Elvevold et al., 2007).
Nordenskiöld Land, named after the Finnish-Swedish geologist Nils Adolf Erik Nordenskiöld (Norwegian Polar Institute, 2021a,b), forms the largely ice-free peninsula between Van Mijenfjorden, Isfjorden and Bellsund. Geologically, the peninsula is part of the contact zone of two large structures of the first order: the horst-anticlinorium of the western coast of Spitsbergen and the West Spitsbergen graben like trough. Quaternary deposits are represented by raised marine sediments in the lowlands, glacial and glacio-fluvial deposits, extensive and complex slope deposits, areas of aeolian sediment cover and extensive in situ weathering of bedrock. The landscapes are diverse: from watershed peaks to the landscapes of U-shaped valleys, small valley glaciers and moraines, and coastal plains; peaks, mountains slopes, and moraines are covered by primary and desert-arctic soils with thin herbaceous-mosslichen groups (Demidov et al., 2020). Both sediments and bedrock is influenced by the perennial frost in the ground (permafrost).
The climate of the area is influenced by the warm West Spitsbergen Sea Current and warm and humid air masses from the Atlantic Ocean; as a consequence, the air temperature is higher than would be expected for the Arctic latitude. There has been a considerable increase in the amount of precipitation over the last century; annual precipitation in 1940 was 482 mm and in 2018 was 704 mm. The highest amount of precipitation occurs from October to March (which overlaps with the period of high cyclonic activity), whilst the lowest amount occurs from April to July. Precipitation during the winter is 1.5-2 times higher compared to summer (Demidov et al., 2020). Also, the mean air temperature has increased over Svalbard by 3-5 degreesC between 1971 and 2017, with the greatest change in the winter months (Hanssen-Bauer et al., 2019).
From a hydrological point of view, the largest rivers in Nordenskiöld Land (greater than20 km in length) are Coles, Gren, and Hollendar with an annual discharge of about 40 mil. m 3 . The continuous discharge of rivers lasts about 5 months, and it usually ends by the beginning of October (Romashova et al., 2019).

Arctic cultural heritage
Arctic cultural heritage represents a special type of heritage that needs extra attention and constant protection (Bogolyubova et al., 2019) due to the harsh conditions that are changing every year in a continuously warming climate. The Svalbard Environmental Protection Act states that all remains from before 1946 are automatically protected and considered cultural heritage sites, along with a 100 m buffer area around them. Cultural heritage in Svalbard is dominated by fragile wooden structures (Fig. 2a), and thermokarst (ground surface displacement and erosion caused by permanent degradation of permafrost) and solifluction (slow movement of surface soil on slopes caused by freeze and thaw processes) processes are therefore a real threat.
The threatened CH includes all types of constructions, like buildings ( Out of these, graves are the most common cultural heritage sites in Svalbard (Prestvold, 2008). There are about 8300 cultural heritage items officially registered in Svalbard and Nordenskiöld Land area holds about 10% of them (872 items, Fig. 1).
One of the first initiatives for the protection of Arctic cultural heritage was taken through the initiation of the Nordic Action Plan to Protect the Natural and Cultural Heritage of Arctic -Greenland, Iceland and Svalbard; the plan was aimed at contributing to the realisation of the goals in the Nordic Environmental Strategy and the Arctic Programme for Co-operation. It was approved on August 23, 1999, by the Nordic Environmental Ministers in Iceland (Nielsen, 2006).
Cultural heritage sites in the Arctic express the capability of the people to adapt to the cold climate and survive in tough conditions. Svalbard's cultural heritage has an even greater value when it is regarded as part of unique Svalbard's "haunted landscape" (Kinossian, 2020) and is considered to be internationally valuable (Hacquebord, 2001). Thus, monitoring and maintenance of cultural heritage sites in the Arctic is considered highly relevant (Dahle et al., 2000).

Methodology
In order to build a comprehensive thermo-erosion gully inventory for the study area, the most recent orthophotos (5 × 5 m pixel size) acquired in 2009-2011 from the Web Map Services (WMS) of the NPI (Norwegian Polar Institute/USGS Landsat, 2021) were interpreted. The Gullies were identified morphologically and digitised on-screen as polygons and integrated into a GIS database. The data is stored in WGS_1984_UTM_-Zone_33N coordinate reference system. Different parameters and statistical analyses were applied to understand the most important factors controlling thermo-erosion gullying.
To quantify the spatial distribution of gullies, we generated the point and area densities of gullies (i.e., the magnitude of number and size of gullies per unit area) using a Kernel Density function. To further analyse the size distribution of gullies, we also examined their frequency-area distribution (FAD). This procedure is particularly common in case of I.C. Nicu et al. landslides, being used to identify the probabilities of landslides belonging to specific size bins. Specifically, the FAD of medium and large landslides has been long recognized as distribution exhibiting power-law scaling and the variation between the probabilities of different size bins characterized by the slope of the distribution, called power-law exponent (β; e.g., Malamud et al., 2014). Topography is most likely the main factor controlling the β and thus, specific β ranges are often used in the literature for different environmental settings (e.g., ten Brink et al., 2009, Liucci et al., 2017. In this regard, FADs of landslides are exploited for various purposes including quantitative landslides hazard assessments (Guzzetti et al., 2005) or denudation caused by landslides (Hovius et al., 1997). However, size statistics of gullies has still not become commonplace in the context of gullies, although it could be potentially used to leverage our understanding regarding hazard caused by gully erosion. In this study, to identify β and the power-law fit for the FAD of our gully inventory, we used the code provided by Tanyas et al. (2018) and Clauset et al. (2009), respectively.
In addition to this, we also explored whether locations where gullies are present behave differently from locations where these processes are absent, in terms of terrain properties. Similarly, we assessed whether gullies have diagnostic morphological characteristics that can be retrieved from the polygonal inventory. We also analyse whether gullies originate in specific lithotypes. Lithotypes were extracted from the Norwegian Polar Institute WebGIS service Geological Map, scale 1:250000 (Norwegian Polar Institute, 2021a,b). These analyses completed the digital cartographic assessment and were further complemented by field-based surveys.
For small-scale gully erosion monitoring in Hiorthhamn, topographic surveys were made during field trips in late August/early September in 2019 and 2020. The timing was chosen to capture the end of summer situation, with the deeper active-layer depths and after the initial geomorphological activity during the development of an active layer. The surveys were made with a Trimble S5Series Motorized total station and a Trimble TSC3 controller; detailed point coordinates along each gully limit and thalweg were measured. The data acquired was used to calculate some basic scale (Ltlength, Dhhorizontal distance, Ldlength between gully head and gully mouth, Hheight,) and shape parameters (Cvvertical curvature, Gaaverage gradient, Aarea) ( Table 1) (Ding et al., 2017). The three gullies were specifically chosen for monitoring due to their location in the proximity of fragile cultural heritage.
Cultural heritage data was retrieved from the online Norwegian national heritage database -Askeladden, Riksantikvaren (2021). The study area has 872 cultural heritage sites, which are registered as points; and an additional 100-m buffer area around each one of them, according to Svalbard Environmental Protection Act. This dataset was used to check how many CH sites (along with buffer areas) intersect with a gully, and it was made with the help of the Intersection tool from ArcGIS.

Thermo-erosion gullies data analysis
A total of 810 TEG were digitised. They are represented as points in Fig. 3a. Different types of thermo-erosion gullies were mapped; a few examples are shown in Fig. 3b, c, d, e, f, g. Following the density analysis performed, the highest density (Fig. 4a) of thermo-erosion gullies is found north-east from Barentsburg, along the shore between Heerodden and Kapp Laila and Hollendardalen valley. This area is followed by the three second-densest areas located in Sassendalen (north-eastern part of Nordenskiöld land), on the shores of Van Muydenbukta bay and on the land slip between Kapp Morton and Svartodden (in the south-east). The large majority of the thermo-erosion gullies are located along the coastline. Judging by their area density (Fig. 4b), the largest thermoerosion gullies occur at the highest density point, located at the mouth of Hollendardalen valley.
As for the FAD of thermo-erosion gullies (Fig. 5), a size distribution that is quite similar to the FAD of a landslide inventory were obtained. The power-law exponent (β) of the inventory was calculated as ~ 2.5, which is also similar to the values suggested for subaerial landslide inventories (i.e., β mean = 2.4-2.5; Malamud et al., 2014;Tanyas et al., 2018). Fig. 5 also shows that the FAD diverges from the power-law towards the small gullies, which bears similarities with landslide inventories. We should stress that these similarities should not be attributed to the same physical explanations. For instance, variation from high to low frequencies observed from small to large landslides, on one hand, is associated with the transition from shear resistance controlled by cohesion to friction, respectively (e.g., Bennett et al., 2012). On the other hand, gullies are enlarged as a result of gradual surface processes and therefore, large gullies are, in fact, expanded versions of small ones. However, taking aside the similarities between FADs of landslides and gullies, examining the FADs of more gully inventories would help us to understand if there is a specific ratio (i.e., a specific range for β) between the frequencies of small and large gullies. If this is the case, understanding the factors governing β could lead to developing more robust hazard assessment studies in terms of gully erosion.
The FAD shown above essentially provides an overview on the number of small to large gullies within a study area. But, as the traditional visualization is based on logarithmic scale, the actual distribution of the gully dimensions in the metric system is lost. For this reason, we also provide an overview of the gully area distribution in their original scale and complement this information with two additional morphometric characteristics, these being the maximum distance and the elongation index. These three parameters are shown in Fig. 6 where, as in many other natural processes, the respective distributions appear to be heavy-tailed. This is particularly true for the gully area, with a distribution that features large numbers of small gullies with a surface of around 500 m 2 , and a progressive number of fewer gullies of medium to large sizes, up to a maximum close to 6000 m 2 . Interestingly, the distributions of the two additional parameters computed are still positively skewed but more and more light-tailed. The bulk of the maximum gully distance (Fig. 6, central panel) distribution is approximately 30 m in length and only a few have matured into landforms with lengths greater than 300 m within this timeframe. To make matters more complex, the timeframe has varied since the gullies are often situated on isostatically uplifted previous marine deposits with varying age of uplift over the sea level.  As for the shape of the gullies, most of them appear to have an elongation index of 2.2 (Fig. 6, right panel). This value implies that their shape tends to be rounded and only few cases reach elongation indices above 5, a threshold marking very elongated landforms. The information these three plots convey supports the statement that the genesis of gullies in Nordenskiöld Land is predominantly recent. Before describing the results of the field surveys though, we opted to summarize the distribution of terrain characteristics at locations with or without gully occurrences.

Terrain factors characterising thermo-erosion gullies in Nordenskiöld Land
Results from the terrain factor analysis is shown in Fig. 7 where   Fig. 3. a. Mapped thermo-erosion gullies in the study area represented as point features. A few examples of different types of thermo-erosion gullies on Svalbard; b. Gully cut into old raised beach deposits (Isfjord Radio in the background); c. Gully cut into uplifted beach and marine deposits at the mouth of Linnéelva river close to the CH at Russekeila; d. Gully cut into uplifted beach deposit at Finneset (next to Barentsburg); e, f. Gullies cut into the fluvial fan surface at Hiorthhamn (Adventfjorden) affecting cultural heritage remains from mining industry; g. Gully located in the proximity of Svea cut into raised marine deposits with soliflucted surface. marked differences in the distribution of terrain characteristics between locations labelled as NG (No-Gully) and G (Gully), are indicators of environmental attributes that favour gully occurrences. We stress here that to extract the G information, we have created a buffer of 30 m around each TEG polygon, subtracting the TEG surface itself. This implies that the terrain characteristics we report in Fig. 7 for gully presences are diagnostic of the relatively stable neighbourhood where the gully initiated, which we investigated thanks to a 5 m resolution DEM. For each buffer, we then computed the mean of each terrain parameter, whereas for the NG cases, we used the whole distribution across the study site.
For the Elevation, the bulk of the gullies are in the range between the sea level and 50 m a.s.l., with a few exceptions up to around 200 m a.s.l. Conversely, Eastness and Northness show a similar exposition distribution in both cases, although gullies are slightly more predominant on Eastward and Northward facing slopes. An interesting situation is shown by looking at the Elevation, Slope and Topographic Roughness because these three parameters convey part of the very same spatial signal but also bring something different at the same time. For instance, no gully occurs at elevations greater than 200 m a.s.l., but gullies occur up to a maximum steepness of 37 degrees and generally in a mildly rough terrain.
As for the influence of the terrain curvatures, these appear to be also of less relevance. This is particularly true for the profile curvature (PRC) whereas in the case of planar curvature (PLC), the distribution of PLC values extracted at locations where thermo-erosion gullies occur predominantly shows negative (laterally concave) and zero (linear surface) values. This implies that there are no gullies where PLC is positive (laterally convex).
Morphologically, this is coherent with what we expect from a geomorphological standpoint for landforms that tend to exhibit negative planar curvatures correspond to landscape incisions where gullies  I.C. Nicu et al. naturally initiate and develop. Conversely, positive PLC values are diagnostic of terrains where water diverges rather than converge, limiting the surface incision, a condition for the genesis of thermoerosion gullies. Ultimately, we also tried to investigate whether underlying lithologies predispose the occurrence of thermo-erosion gullies. This attempt is challenging since there is no coherent dataset over unconsolidated sediments or bedrock lithologies on Svalbard, but only a composite, where areas with thick unconsolidated sediments (class MOR, MAD and GFD) are only coarsely defined, but the same area are lacking lithological information. In the same way areas with shale dominated lithology are so physically weathered by the harsh Arctic climate that the residual material have mechanical properties similar to unconsolidated sediments. Fig. 8 highlights that thermo-erosion gullies are predominant in unconsolidated substrata (MAD and GFD). This is logical since unconsolidated sediments have open pore-spaces, which is be the base for high accumulation of internal ice in permafrost regions, and hence vulnerable to thermo-karst processes during thawing seasons. This may also be due to the natural spatial distribution of these parent materials in the lower elevations of the landscape. They mainly occur along valley bottoms and the coastline where the combined action of fluvial and coastal erosion and thermal variations may destabilize the soil leading to the genesis of thermo-erosion gullies.
Aside from unconsolidated sediments, the underlying lithotypes with the highest frequency of gullies are sandstones (SND) and shale-rich parent materials (SMS and CSS). The former lithotype will, as climate induces rapid surface weathering, alter the grain-to-grain bond, tend to produce loose soil columns, which can be mobilized with relative ease as permafrost thaws. Shale-rich bedrock have the added instability of giving rise to clayey soils which are naturally weak to variations in water content. As permafrost thaws is releasing water into the soil; this is combined with the swelling effect of clay minerals to rapidly reduce stability in the sediment.

Monitoring of small-scale gully erosion
The results of the monitoring of three small-scale fluvial gullies are shown in Fig. 9. The gullies are located on the north-eastern shores of Adventfjorden, approximately 3 km north of Longyearbyen (Fig. 9a), at the mouth of Telegrafdalen river. Limits (Fig. 9b) and main characteristics of the three gullies are shown in Table 3.
During the two years monitoring, there were very few changes in the gullies morphology. The length (Lt) varies between 10.7 and 16.65 m for gully number 01 and gully number 02, respectively. It can be observed that the length of the gullies has diminished from 2019 to 2020 monitoring. This is due to the fact that there is a lot of sorted and unsorted material coming down from the fluvial system upstream, filling up the gully depression. This happens during short term episodes of high energy surface overland flow. Also, being close to the coastline, the gullies are under the effect of the tides, which is also transporting sand and driftwood into the gullies; hence, influencing the thalweg length of the gully. Lt has a similar effect on the Ld, which has also decreased ( Table 3). There was no significant gully head advancement for the three gullies, except for gully number 01, which advanced with 0.5 m in one year. (Fig. 9b).
The same tendency can be observed for the H, which has also decreased for all gullies. Cv values are all close to 1, which shows the fact that the gullies longitudinal profile is close to a straight line shape, which is also visible in Fig. 8b. Ga values are low, indicating the fact that the gullies are not very active, and the erosion is not severe. Regarding the A of the gullies, there has been an increase; this is due to the lateral collapses, especially for the gully 01 and 02 (Fig. 9b). These lateral collapses seems to be due to undercutting through high tide/wave induced action in the lower parts of the gully.

Cultural heritage analysis
The majority of Svalbard's CH are located in the coastal area (Nicu et al., 2021b). Following the Intersection analysis, out of 872 CH sites, 44 sites (Supplementary material 1) were found to be vulnerable to future evolution of TEG; out of 44 CH sites, 29 are also exposed to coastal erosion and three to river erosion. They are mainly located near the two main settlements, Longyearbyen and Barents burg, and Colesbukta.
The grave visible in Fig. 10a is located in Finneset, at approximately 1 km SSE from Barentsburg. The description of the grave is as follows "well-defined grave, located between two brooks in a flat area about 20 m from the shoreline; the grave is oriented E-W, with well-preserved cross. The grave is covered with a thick rope, about 40 cm high with sharp edges" (Askeladden, Riksantikvaren, 2021). The description could be more precise (with some additional information observed in the field during the summer of 2021), also referring to the landscape around the cross; the two "brooks" are actually two well-separated active gullies of considerable size. The one in the northern part has a length of about 165 m and the one in the southern part a length of about 155 m. The cross itself is also vulnerable to coastal erosion.
The house remains from Fig. 10b are described as "rectangular house I.C. Nicu et al. remnant, oriented E-W. The house is marked as a burial in the ground with remnants of peat walls, walls of 30 cm wide, 20 cm high and 10 cm in depth. Remains of dry vertical logs of 50-100 cm in height" (Askeladden, Riksantikvaren, 2021). This is located approximately 125 m SSW from the grave, on the south part of the 165 m gully. Also vulnerable to coastal erosion.
The Security Zone (Fig. 10c) has no description; however, this refers to the 100 m buffer zone around a protected cultural heritage site, in this case the remains of the cable car that used to transport coal in the past (Fig. 10f). The security zone is also vulnerable to coastal erosion (Jaskólski et al., 2018).
The remains of the road (Fig. 10d) are described as "the road from Skjaeringa to the preparation plant; it was built in 1967 and was in use until 1985, when a new road was built closer to the sea. The road is built up on the terrain and gravelled. A possible theory of the name may be that it originates from a steep road" (Askeladden, Riksantikvaren, 2021). The new road refers to the present day road that is connecting Longyearbyen to the Longyearbyen airport.
The plane wreck in Adventdalen (Fig. 10e) is described as "remains of wings and cockpit, scattered remains in the area around" (Askeladden, Riksantikvaren, 2021). No further details were provided. However, more details were found about the context of the plane wreck. In the past, planes were landing in Adventdalen (the old airport), about 4 km east of Longyearbyen, until the new airport at Hotellneset was built in 1973 and opened in 1975. The wreck is represented by a German Ju 88, which on June 14, 1942, flying from Banak in north Norway landed on the airstrip in Adventdalen; it was damaged because the ground was soft. Almost two weeks later, it was seriously damaged during an attack by an allied Catalina; following this, it was left behind. Today, the wreck of the German Ju 88 is very close to Adventelva River (Stange, 2019) (as visible in the background of Fig. 10e). Besides the gully that is visible in the foreground of Fig. 10e, the wreck is vulnerable to Adventelva river erosion.

Discussion
There is an immediate need to study the implications of climateinduced geomorphological processes and hazards in the Arctic with the potential to affect cultural heritage sites. One such relevant process is thermo-erosion gullying, which can be postulated to be increasingly active under a warming climate (Godin et al., 2019). Literature has shown that the peak of erosion rate and enlargement of a gully occurs during the first part of its entire lifetime (Sidorchuk, 1999(Sidorchuk, , 2006. Godin et al., 2014 study from Canada showed that active gullying triggered not only a change in drainage pattern but also more efficient drainage of the neighbouring wetlands; moreover, a stabilised state was usually reached between 5 and 10 years after the initial incidence of erosion. It was shown that gullies were often co-existing with thaw slumps, sinkholes and tunnels and that their erosion rates were highly variable over the last four decades ; with mean values of 4.39 m y -1 in length and 61.4 m 2 y -1 in area. In this study, we were not able to provide such long-term analysis, taking into consideration that our main data set comes from only one set of spatial observations. As soon as new high resolution spatial data becomes available (the present ones are from 2009 to 2011), a spatio-temporal analysis will be made. Also, the Godin et al., 2014 erosion rates are not comparable with the ones from the oneyear small-scale gullying monitoring that was made.
The Devon Island study found that: gullies are found an all slopes, regardless of orientation; geology has a significant role in gully initiation, development and morphology, especially geological formations dolostone and limestone. In this study, we find that sandstones are dominant of the solid bedrock lithologies, while in Nordenskiöld Land, slopes oriented towards north, and east are more favourable to gully initiation and development. The distribution of gullies only under 200 m a.s.l. in this study points to the importance of the substrate (sediment) for gully formation. Central Svalbard has experienced an iso-static rebound of up to 65-70 m after the last ice-age (Lønne, 2005) -meaning areas under this elevation often have a marine-originating sediment cover. The lower valleys in this paraglacial landscape are also the sites for both glacial-and glacio-fluvial sediment deposition, while higher areas are exhibiting more bare rocks, slope process deposits and in situ weathering material. Sidorchuk, 2020 study from Russia (Yamal Peninsula) revealed highly uneven distribution of erosion potential level with the maximums on steep banks of the river valleys and in gully heads; moreover, the erosion potential was linked to landforms with steep slopes and unstable vegetation cover. Anthropogenic interventions, such as railways, roads and pipelines were also found to influence the initiation and development of gullies. In our case, gullies can be found on up to 37 degrees slopes and numerous gullies are developed around the two main  Table 2 (here we show only the underlying bedrock for the gullies, leaving other lithotypes unmarked). I.C. Nicu et al. settlements -Longyearbyen and Barentsburg. At this point, however, it cannot be concluded weather direct physical anthropogenic interventions has had a significant contribution to the initiation and development of gullies.
Gullies are landforms that widen and lengthen in time as land degradation takes place (in this case laterally thawing of permafrost). And the very fact that most of the mapped gullies are small, short and rounded in overall shape supports the hypothesis that they still have a large potential to mature. These considerations are valid looking at the study area as a whole and a detailed assessment is required to confirm or reject this hypothesis. Other considerations would be to look and try to detect the time the gullies have been active and try to detect any on-off mode in activity, for example triggered by interchanging warmer and colder periods, especially concerning the lengths of thawing seasons and maybe maximum summer temperatures. This is however difficult since spatial data, aerial and satellite photographs are sparse in time and space over the arctic back in time.
Further research is needed on Svalbard's thermo-erosion gullies, especially regarding the spatio-temporal development and connection to sediment types, to better understand local and regional factors controlling distribution and future development. Understanding the largescale pattern represents a future research direction, along with continuing the small-scale gully erosion monitoring. Another future research direction is towards (multi)hazard predictive modelling of Arctic thermo-erosion gullies and retrogressive thaw slumps, given the logical link between thermal-gully initiation and activity, and arctic climate warming caused by the geomorphological mechanism linked to the thickness of active layers in permafrost regions.
The here presented thermo-erosion gully inventory is valuable when regarded from a climate change perspective, as it can be utilised by scientists from various backgrounds, such as cryosphere research, geomorphology, statistical modelling of climate-induced geomorphological processes, hydrologists, cultural heritage planners, etc. The inventory and maps also represent a new standalone dataset of high importance for  both local, regional and international context and can be used in future spatio-temporal analyses and/or multi-hazard assessment of climateinduced processes (Lombardo et al., 2020;Saha et al., 2021;Javidan et al., 2021). The accumulation and combination of these types of datasets could have a high impact on mitigating climate change effects on pan-Arctic infrastructure (Hjort et al., 2018) and Arctic cultural heritage landscapes (Nicu et al., 2021b;Thuestad et al., 2015Thuestad et al., , 2021. As Arctic cultural heritage landscapes represent a projection of the past, it is in our moral responsibility to try to protect and preserve them.

Conclusions
Nordenskiöld Land, Svalbard represents a typical Arctic permafrost landscape with a multitude of unique cultural heritage sites and remnants, which has not previously been studied for thermo-erosion gullying hazards. We present for the first time a comprehensive dataset of 810 gullies in this area, their general topographical and lithological setting, together with a spatio-statistical analyses comparing them with protected and vulnerable cultural heritages in the same area. We also present some observations and hypothesis into geographically and sedimentological controlling factors of thermo-erosion gully activity and spatial distribution, together with an initial short-term small study monitoring of small-scale gullying. The results so far indicate that thermo-erosion gullying and cultural heritage in Nordenskiöld land coexists on a regional level, although direct erosion and hazard for CH is not yet at any critical scale. The spatial coastal and low-land focus of both thermo-erosion gullies and CH however gives reason for future concern. A rapidly warming Arctic climate would favour increased permafrost active layer thickness and hence potential accelerated thermo-erosion gullying in susceptible areas, which we show overlap with the highest CH densities. Although not conclusive in understanding present and future development of thermo-erosion gullying, this study underpins the relevance and importance of geohazard mapping and research in Arctic settings where CH and present human activity and infrastructure is the most vulnerable to increasingly rapid changes in natural physical and geomorphological conditions and processes.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.