**Cultural Heritage and Natural Disasters**

Editors

**Ionut Cristi Nicu Alin Mihu-Pintilie Erich Nau**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Ionut Cristi Nicu Norwegian Institute for Cultural Heritage Research (NIKU) Norway Alin Mihu-Pintilie Alexandru Ioan Cuza University of Iasi (UAIC) Romania Erich Nau Norwegian Institute for Cultural Heritage Research (NIKU) Norway

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Sustainability* (ISSN 2071-1050) (available at: https://www.mdpi.com/journal/sustainability/ special issues/CHND).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-1078-1 (Hbk) ISBN 978-3-0365-1079-8 (PDF)**

Cover image courtesy of Ionut Cristi Nicu.

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editors**

**Ionut Cristi Nicu** (Dr.) is a geographer with a highly interdisciplinary background, showing special interest in the evaluation and monitoring of natural hazards (landslides, gully erosion, coastal erosion) with direct effects on immovable cultural heritage. A few keywords that can better describe his research background and interests are soil erosion; natural hazards; landscape archaeology; geomorphology; cultural heritage; remote sensing; GIS; protection of cultural heritage; spatial analysis; predictive modelling. Currently, he is working at the High North Department within the Norwegian Institute for Cultural Heritage Research in Tromsø. His current focus is on the effects of climate change on immovable cultural heritage in Arctic areas (Svalbard), participating in several projects focusing on this issue. His future interests in Svalbard are related to understanding the effects of climate change on permafrost landscape dynamics and future risks for Arctic cultural heritage in an integrated approach for multi-hazard evaluation.

**Alin Mihu-Pintilie** (Dr.) is a geographer interested in advanced spatial analysis using GIS, remote sensing and geophysical techniques for bathymetric and terrestrial mapping, geoarchaeology, underwater archaeology and for reconstructing the alluvial paleo-landscapes. He has been a specialist in continental and oceanic sea diving since 2008 (NAUI Worldwide Patent), has been part of numerous geographical and archaeological expeditions in different parts of Eurasia (Danube Delta, Volga Delta, Black Sea, Baltic Sea, Southern Siberia) and the North Atlantic (Iceland, Norway). In the Special Issue "Cultural Heritage and Natural Disasters" for which he is a Guest Editor , he focused on the following topics: - Integration of remote sensing and GIS for cultural heritage assessment; - Finding the new tools and methods to generate new data to continue possible monitoring processes;


**Erich Nau** (Mag. phil.) is an Austrian archaeologist employed by the Norwegian Institute for Cultural Heritage Research (NIKU). His main expertise and interests lie within archaeological–geophysical prospection and the 3D documentation of cultural heritage. In his current position at NIKU, he is devoted to the further development of a motorized large-scale ground penetrating radar and in particular, its consequent application within archaeological registration projects as a standardized tool. He has led several large commercial projects applying GPR for archaeological registration prior to major infrastructure expansions in Norway. His current focus of research is that of studying the effects of seasonal and environmental factors on the visibility of archaeological structures in GPR data within the "Vestfold Monitoring Project". Within his second field of interest, he is working on the application of state-of-the-art 3D documentation technology such as laser scanning and image-based modelling for cultural heritage documentation and research. The focus in his recent work has been on visualizing large amounts of 3D data and the interactive publication to a wider public. Among other things, he was responsible for the documentation project at Sveagruva on the High Arctic Archipelago of Svalbard. He is a regular participant and presenter at national and international conferences and a member of the International Society of Archaeological Prospection.

## **Preface to "Cultural Heritage and Natural Disasters"**

In recent decades, the intensity and increasing numbers of natural disasters have started to affect immovable cultural heritage around the world. This is important due to the high complexity and impossibility of recovering any kind of data from a cultural heritage asset destroyed by a natural disaster. This SI has brought to the attention of the readers a few significant studies that approach the issue of natural disasters and cultural heritage. Notably, the SI has published papers on topics of interest which were proposed at the beginning of the SI launch. In total, six papers were published as a result of applying different and new methods in the field of cultural heritage. They are spread over different parts of the globe, as follows: Romania (two studies), Svalbard and north Norway (two studies), Spain (one study), Canada (one study).

Three of the papers propose new methods to be applied in the field of cultural heritage: a new method for archaeological predictive modelling (APM), a method that employs image-based modelling to measure the erosion at archaeological sites, and a study that compares the performance of airborne laser scanning and remote sensing in Arctic landscapes. The other three studies employ methods and techniques that are well known and used in the field of cultural heritage; the use of old maps and new orthophotos, combined with field surveys were used to assess the erosion rates of coastal cultural heritage located in Svalbard. A 3D laser scanner was used to document and graphically analyze the pavilions muqarnas at the Court of the Lions in the Alhambra in Granada (a UNESCO World Heritage Site). The usefulness and applicability of satellite-based remote sensing was once again successfully applied to document, map and monitor archaeological sites in urban landscapes (Apulum, Alba Iulia, Romania). More efforts should be taken at an international level to develop new methods in the field of cultural heritage and/or to apply methods that are already used in studies at the border between cultural heritage and natural disasters, as cultural heritage represents the legacy of our ancestors.

> **Ionut Cristi Nicu, Alin Mihu-Pintilie, Erich Nau** *Editors*

## *Article* **Archaeological Surveying of Subarctic and Arctic Landscapes: Comparing the Performance of Airborne Laser Scanning and Remote Sensing Image Data**

**Alma Elizabeth Thuestad 1,\*, Ole Risbøl 2, Jan Ingolf Kleppe 3, Stine Barlindhaug <sup>4</sup> and Elin Rose Myrvoll <sup>4</sup>**


**Abstract:** What can remote sensing contribute to archaeological surveying in subarctic and arctic landscapes? The pros and cons of remote sensing data vary as do areas of utilization and methodological approaches. We assessed the applicability of remote sensing for archaeological surveying of northern landscapes using airborne laser scanning (LiDAR) and satellite and aerial images to map archaeological features as a basis for (a) assessing the pros and cons of the different approaches and (b) assessing the potential detection rate of remote sensing. Interpretation of images and a LiDARbased bare-earth digital terrain model (DTM) was based on visual analyses aided by processing and visualizing techniques. 368 features were identified in the aerial images, 437 in the satellite images and 1186 in the DTM. LiDAR yielded the better result, especially for hunting pits. Image data proved suitable for dwellings and settlement sites. Feature characteristics proved a key factor for detectability, both in LiDAR and image data. This study has shown that LiDAR and remote sensing image data are highly applicable for archaeological surveying in northern landscapes. It showed that a multi-sensor approach contributes to high detection rates. Our results have improved the inventory of archaeological sites in a non-destructive and minimally invasive manner.

**Keywords:** cultural heritage; LiDAR; satellite image; aerial image; High North

#### **1. Introduction**

Remote sensing has become a very valuable resource for modern archaeology as a useful tool for finding and documenting cultural heritage sites. Airborne laser scanning (LiDAR hereafter) and aerial and satellite image data are well established and highly applicable approaches to find and study sites in a non-destructive and minimally invasive manner. Archaeological application of remote sensing data is, however, rarely a straightforward process, and considerable efforts have been put into developing and improving methodological approaches. The pros and cons of different remote sensing data vary and, accordingly, so do areas of utilization and methodological approaches applicable for analysing the data. Other complicating factors are the multitude and variety of archaeological sites and the wide range of natural environments in which they are located.

The northernmost part of Norway has been surveyed for archaeological sites by researchers and management authorities since the early 1920s [1–3]. Survey efforts henceforth focused primarily on coastal and riverine areas as well as other relatively easily accessible areas. Extensive areas, especially in the interior, are thus yet to be systematically surveyed. Consequently, there is a significant disparity in the number of known archaeological sites in inland and coastal areas, as can be seen in Figure 1. Survey projects conducted in inland

**Citation:** Thuestad, A.E.; Risbøl, O.; Kleppe, J.I.; Barlindhaug, S.; Myrvoll, E.R. Archaeological Surveying of Subarctic and Arctic Landscapes: Comparing the Performance of Airborne Laser Scanning and Remote Sensing Image Data. *Sustainability* **2021**, *13*, 1917. https://doi.org/ 10.3390/su13041917

Academic Editor: Alin Mihu-Pintilie

Received: 15 December 2020 Accepted: 4 February 2021 Published: 10 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

areas in recent decades have highlighted the fact that this distribution pattern is not representative of prehistoric activity in this region [4]. There is thus a need to utilize data and methodological approaches that allow remote and hard-to-get-to areas to be surveyed in a low-impact and cost-efficient manner.

A characteristic of cultural heritage in the subarctic and arctic landscapes of Northern Norway is that traces of past human activity are often directly visible on the surface. Sparse vegetation cover, slow regrowth, a relative absence of modern infrastructure and low agricultural activity are factors that have allowed for preservation of sites from the pioneer settlement up until modern times. Remote sensing utilizing LiDAR data and satellite and aerial images has a great potential for improving the inventory of known archaeological sites throughout the region.

The aim of this study is to assess and compare the applicability of LiDAR and remote sensing image data for archaeological surveying in the subarctic and arctic landscapes of Northern Norway. High-resolution satellite images, digital aerial photographs and LiDAR data are used to identify and map potential archaeological features as a basis for (a) assessing the pros and cons of the different approaches, and (b) assessing the potential detection rate of remote sensing in northern landscapes.

**Figure 1.** Cultural heritage sites listed in the national cultural heritage database Askeladden in December 2020, here marked in amber, in and around the Varanger Peninsula in northern- and easternmost Norway. The 25 km<sup>2</sup> large study area, here outlined in red, is located on an isthmus between a river and the inner part of the fjord just south of the Varanger Peninsula. Image ©Norwegian Mapping Authority.

Remote sensing is not a new concept in archaeology [5]. In Norway, systematic use of remote sensing for archaeological purposes began in the early 1960s in conjunction with a large national state-funded survey project conducted between 1960 and 1991. Aerial photographs were then used as base maps for delineating cultural heritage and as a basis for producing Norwegian Public Land Use Maps. The actual surveys were ground-based and conducted through traditional means [6].

Remote sensing platforms, having grown increasingly sophisticated, now offer significantly finer spatial and spectral resolutions than previously available, coincidentally providing data that are progressively more useful for archaeological purposes. Initial examples of archaeological applications of satellite images in Norway used Landsat images with 30 m resolution [7]. Image data are highly useful for monitoring vegetation changes and have proven suitable for assessing subsequent impact on cultural heritage. Landsat data (Landsat 5 Thematic Mapper (TM) and 7 Enhanced Thematic Mapper (ETM+)) with a spatial resolution of 28.5 m were used to show impact on archaeological features through regrowth and reforestation brought on by farm abandonment in Northern Norway [7,8]. In Northern Norway, Quickbird-2 image data (available from DigitalGlobe, a subsidiary of Maxar Technologies, Longmont, CO, USA) with 0.6 m (panchromatic image) and 2.4 m (multispectral image) resolution enabled detection of maintained and abandoned reindeer pens, milking pens and dwellings linked to Sámi reindeer husbandry [9].

Archaeological use of remote sensing image data encompasses quantitative and qualitative techniques and often a combination thereof. Visual analyses, alone or in combination with image enhancement techniques or semi-automatic procedures, have proven useful for detecting surface and subsurface archaeological features, especially large-scale features [10]. Image processing techniques are a valuable contribution as they serve to enhance the interpretability of data by enhancing spatial patterns or local anomalies linked to past human activities [11–13]. Effective automatic procedures have proven challenging for archaeological purposes, but semi-automatic approaches and enhancement techniques work quite well, keeping in mind that their performance has proven somewhat "site specific" and "feature specific" [12,14]. In Norway, a semi-automatic detection methodology was developed and implemented through the CultSearcher software developed by the Norwegian Computing Center (Oslo, Norway) in the early 2000s [15,16]. It was initially based on satellite image data, but the focus has since been directed more towards semi-automatic detection based on LiDAR data [17]. The system detects potential archaeological features, and accordingly, verification is dependent on field surveys. The use of remote sensing image data has proven most successful when archaeological sites are located in landscapes with a relatively uniform topography and vegetation cover. For instance, features located in arable land or marl landscapes appear as anomalies in an otherwise uniform image [16,18].

LiDAR has been available to the archaeological community since the beginning of the millennium [19]. It is an active method which collects data from the ground using laser pulses, rendering it possible to collect very detailed information about the earth surface, which can be used to create high resolution 3D representations of the topography. An important aspect is the ability to classify the collected data in ground and off-ground data, enabling the elimination of laser points bouncing off vegetation, the top of buildings, etcetera and leaving the user with a "bare-earth" model. LiDAR has, due to this principal advantage, gained interest among archaeologists involved in mapping cultural heritage [20]. A main objective of many early studies was to test the suitability of LiDAR in different landscapes and with regard to the type of archaeology specific for the region worked in. LiDAR has proven effective for identifying regularly shaped archaeological features such as burial mounds, charcoal kilns, hunting pits, etcetera [21,22]. Further work has concentrated on elucidating the detection success of mapping cultural heritage within a given area [22,23]. Another field of interest has been different visualising techniques aimed at improving the detection success [24–27], followed by studies comparing different visualising techniques. The results generally point towards the conclusion that one will benefit from using various techniques on the same data set [28–30].

#### **2. Study Area and Cultural Historical Background**

#### *2.1. Study Area*

The 25 km2 large study area is located on an isthmus in the easternmost part of Northern Norway. The area encompasses inland mountainous and lowland landscapes as well as coastal landscapes. According to land resource maps [31] the area outlined in Figure 1 encompasses forested areas (47.3%), mountainous and other areas without vegetation (17.7%), marshes and moorland (32.6%), agricultural areas (1.1%), freshwater rivers and lakes (1%) and ocean (0.3%).

The area is located in a region with a rich cultural heritage going back to the pioneer settlement following the retreat of the ice cap at end of the last Ice Age. The area was selected precisely because it encompasses a large number of cultural heritage sites and features typical for the region, and because it is an area with a very high potential for so far unknown or unlisted sites. The latter is especially true for the inland part of the area. A primary focus of this study is therefore on the inland areas, which are also less disturbed by modern activity than the coastal areas.

#### *2.2. Cultural Historical Background*

Human presence in Northern Norway extends back around 12,000 years. As ice melted and withdrew after the last glacial period, archaeological evidence shows that ice-free areas were quite rapidly inhabited [32,33]. The older prehistory of the region is divided into three periods characterized by different technologies and social organization: The Early Stone Age (10,000–4500 BC), the Late Stone Age (4500–1800 BC) and the Early Metal Age (1800 BC–BC). Latter prehistory is designated the Iron Age (0–1100 AD) and the Medieval period (1100–1600 AD). Throughout prehistory, well into historic times, nomadic or semipermanent settlement based on hunter-gathering, fishing and reindeer herding have been predominant [34].

The inlands of the region where the study area is located is especially known for numerous hunting pit systems reminiscent of early reindeer hunting, and the coastal areas for the multitude settlement traces ranging from the Stone Age to historic times.

Numerically speaking, in the northernmost part of Norway, hunting pits are the most common category of archaeological feature listed in Askeladden, the national cultural heritage database [35]. Hunting pits most often appear as up to 1.5 m deep and 2–5 m wide circular or oval depressions, normally surrounded by a low wall of soil (see Figure 2 for examples). Although disputed due to methodological challenges when dating hunting pits, excavations have yielded dates covering a time span of almost 4000 years from the Late Stone Age to AD 1600 [36–39]. Pits are found throughout Northern Norway, but numbers are highest in inner fjord areas and inland areas along major rivers [36,40,41]. Other features commonly linked to wild reindeer hunting are hunting blinds and leading fences. The high numbers and widespread distribution are indicative of the importance of reindeer hunting from the Stone Age to Sámi subsistence and economy [42,43].

Settlement sites dating to the Late Stone Age and Early Metal Age are often identified through the still visible remains of dwelling structures (see Figure 3 for examples). Sites may encompass anywhere from one dwelling to several tens. The dwellings, which may be slightly dug into the ground, are circular or rectangular structures ranging in size from around 10 m<sup>2</sup> to 70 m2. Sámi sites, the oldest going back around 2500 years, commonly encompass camp sites and dwelling remains in the form of circular and rectangular tur houses.

**Figure 2.** Hunting pits for wild reindeer located in the study area during field work in 2008. Photo ©Ole Risbøl, NTNU.

**Figure 3.** Settlement site with dwellings from the Late Stone Age located in the study area during field work. Photo ©Elin Rose Myrvoll, Sámediggi.

#### **3. Materials and Methods**

*3.1. Satellite and Aerial Image Data*

Image data used in this study encompass WorldView 2 satellite data (DigitalGlobe, Longmont, CO, USA) and aerial images made available by the Norwegian Mapping Authority [44]. Interpretation was based on visual analyses aided by image processing techniques. WorldView-2 images (Table 1) and digital aerial photographs (Table 2) were analysed (see below) using ENVI 5.1 software (Exelis Visual Information Solutions, Inc., Broomfield, CO, USA, a subsidiary of Harris Corporation) and ArcGIS 10.4.1 software (ESRI®ArcMap™, Redlands, CA, USA).

**Satellite Image Data Sensor WorldView-2** Product Type: Ortho Ready Standard 2A Product Option: Bundle-multispectral (4 bands) and panchromatic image Acquisition Data (Y-M-D) 2010-07-26 Cloud Cover (%) 0 Wavelength (μm) Panchromatic: 0.450–0.800 Band 1 (Blue): 0.450–0.510 Band 2 (Green): 0.510–0.580 Band 3 (Red): 0.630–0.690 Band 4 (Near-infrared1): 0.770–0.895 Panchromatic Resolution 0.5 m Multispectral Resolution 2 m

**Table 1.** Metadata WorldView-2 image data.

**Table 2.** Metadata digital aerial images.


<sup>1</sup> The study area is covered by part coverages acquired in 2008 and 2010.

The Gram–Schmidt pan-sharpening procedure using nearest neighbour resampling was utilized to spatially enhance the lower resolution satellite multispectral images, ensuring that all analysed images had a spatial resolution of 0.5 m. Interpretation was aided by image enhancement functions available in ENVI, which served to improve interpretability by making subtle features more visible and increasing distinction between various features. Linear contrast enhancement was used, the histogram manipulated to increase contrast between potential features and the surrounds. Visual interpretation of the satellite images was based on both panchromatic and pan-sharpened multispectral images, which were analysed separately to assess whether the results differed significantly. Prior to the final analysis of the multispectral images, a variety of band combinations was tested to ascertain contrast. The final analyses were based on the following sequences: (a) bands 3, 2, 1 (Red-Green-Blue) and (b) bands 4, 3, 2 (Near-infrared1-Red-Green). The first sequence consists of the three primary colour bands red, green and blue (hereafter RGB) producing a so-called "true colour" image that closely resembles what would be observed by human eyes. The second sequence encompassing the near infrared, red and green band allows for vegetation and variations in vegetation linked to type and condition to be readily detected. The aerial images (see Figure 4 for example) were only available as RGB or "true colour".

Detected potential archaeological features were manually delineated at a 1:2000 resolution.

**Figure 4.** Two images of the same section of the study area; on the left, showing hunting pits in the digital terrain model (DTM) and, on the right, hunting pits in the aerial image. Image: ©Norwegian Mapping Authority.

#### *3.2. LiDAR Data*

Airborne laser scanning (see Table 3 for metadata) covering just over 30 km2 resulted in a total of 154 million collected points, of which 86 million were classified as ground points and 68 million as off-ground points. The average point density varied from 5 to 8 points per m2 with approximately 3 ground points per m2. The dataset was divided into 92 tiles covering 700 × 700 m each. A 25 km<sup>2</sup> large section of this dataset was used as a basis for the analyses in this study. The dataset was analysed and interpreted for the first time in 2007/2008 [23,45] and later used for new analyses carried out by Troms and Finnmark County during the period 2010–2020.

**Table 3.** Metadata airborne laser scanning (LiDAR).


In the 2007/2008 study, a bare-earth digital terrain model (DTM hereafter) generated from the ground points constituted the base for analysis and interpretations. Generating the DTM as well as analysis and interpretations of the consecutive hill-shaded model was carried out using a 3D modelling software called Quick Terrain Modeler (hereafter QTM) (Applied Imagery, Chevy Chase, MD, USA) designed to easily handle and enhance LiDAR generated terrain models in a 3D environment [46]. QTM facilitates real-time manipulation of large amounts of 3D data, allowing you to surf the model and manipulate the light angle and direction as well as exaggerate the elevation (the z-value) and make digital crosssections of detected anomalies etcetera. Using this approach, the DTM was put through visual analysis as a basis for detecting and identifying anomalies thought to be potential archaeological features. The results from the 2007/2008 study were the basis for the present study but were supplemented by further studies in 2010–2020 (see Figure 5).

In this period, a new model based on the same dataset with a 25 cm resolution was built. In addition to traditional hillshade analysis, additional methods were used both for identifying previously unidentified anomalies and structures and refining the visualisation and on-screen identification of these. Methods used were multi-hillshade, slope, skyview factor and openness negative/positive as well as blends of these [27,47,48]. The LiDAR data was used to further investigate the detection success by implementing two recently developed visualizing techniques, "Local Relief Model" [26] and "Sky-View-Factor" [24], as additional approaches to analysing and interpreting the data set. The rebuilt DTM was analysed using the visualization techniques "Sky-View-Factor" as well as "Openness positive" and "Openness negative" [25,28] to further explore recent developments in interpretation of DTMs derived from LiDAR. "Sky-view-factor" and "Openness negative/positive" results were both analysed separately as well as overlain with 25–30% opacity. The latter implementation is in line with Kokalj and Hesse [27].

#### *3.3. Ground Surveys*

Ground surveys were conducted in 2008 in order to verify whether the detected features were in fact cultural heritage, and, if so, to verify the interpretation. They were point-to-point surveys encompassing a selection of the features identified in the LiDAR data, they were not general archaeological surveys of the area. In order to avoid biased data, a random approach was used where every tenth anomaly was systematically investigated as part of the field verification. When the desk-based interpretation was carried out, all anomalies were assigned a unique id. A randomised selection was then made, where all anomalies with 9 as its last figure in the id were chosen: 9, 19, 29, 39, 49 and so on. These were then looked up in the field.

Supplementary field work has been carried out in connection with the second study of the LiDAR data set. The ground truthing efforts were mainly focused on hunting pits and the inland parts of the study area. Other feature types and sites in coastal areas were, to a lesser degree, included to provide a broader set of data as a basis for assessing rate of "right interpretation". A less extensive survey focused on verifying potential dwellings near the coast in the eastern part of the scene. Data from the ground-based surveys were used to assess the rate of "right interpretation", not only for the LiDAR data but also the image data. The assessment of the results from analysing the image data was based on overlapping features, that is, point features placed within 10 m of each other in both survey data and analysis results and visually interpreted to be the same feature.

#### **4. Results**

#### *4.1. Potential Archaeological Features Identified in the Data*

The results from analysing the LiDAR data and satellite and aerial images vary from 368 features identified as potential cultural heritage based on the aerial images, to 437 and 1186 in the satellite images and LiDAR data, respectively, as shown in Table 4. Moreover, 347 and 380 features were detected in the panchromatic and multispectral satellite images, respectively. The higher number of features detected in the multispectral image reflect the higher number of dwellings detected and identified in this image set. As some features were only detected in either the panchromatic and multispectral image, after merging results and removing double entries, visual analyses of the satellite image data resulted in a total of 437 potential archaeological features. Of the 327 hunting pits identified in the two image sets, there is a 73% overlap. Although roughly the same number of pits was detected in the two image sets, 89 were only detected in either the panchromatic or the multispectral image. Study 1 of the LiDAR data resulted in the identification of 893 potential archaeological features. The analysis in LiDAR study 2, which focused on the western part of the scene, resulted in 980 features.

**Table 4.** Potential archaeological features detected in the LiDAR and image data.


<sup>1</sup> Split interpretation; one feature was identified as a structure in the panchromatic image and a dwelling in the multispectral image.

The features detected in all data sets are mostly interpreted as pits, specifically hunting pits for wild reindeer. As can be derived from the results presented in Table 4, hunting pits constitute 75% of the features identified in the panchromatic and multispectral satellite images and 94 % of the features in the aerial images. In the initial LiDAR study, 86% of the features were thought to be hunting pits, while in the subsequent study 100% of the identified anomalies were defined as hunting pits. These pits were detected in high numbers in the western and inland areas of the study area, and many were grouped together and organized in rows forming larger systems as seen in Figure 6.

A number of dwellings were detected in low-lying or coastal areas in the eastern part of the study area. This category encompasses features dating to the Stone Age and Early Metal Age as well as remains of Sámi turf houses, most of which probably are from the 1800 and 1900s. Although it is generally possible to differentiate dwelling sites of varying age through placement and characteristics apparent in the image data, the results listed in Table 4 show the total number to ensure comparability with results based on the LiDAR data. However, it is clear from the aerial and, especially, the satellite images (see Figure 9) that the vegetation in the immediate surrounds of several features interpreted as turf houses has been impacted by human activity over time, possibly indicating long-term settlement. In close proximity to and in conjunction with dwellings, a number of features interpreted as pits and structures of uncertain function were detected. The remaining features listed in Table 4 were identified as stone rings and a mound.

**Figure 6.** A section of the study area showing hunting pits detected in the LiDAR data and the aerial and satellite images. Image: ©Norwegian Mapping Authority.

#### *4.2. Verification of the Potential Aarchaeological Features Identified in the Remote Sensing Data Sets*

The results from ground surveys of selected features are presented in Table 5. A total of 139 features were surveyed and identified as either an archaeological feature or a natural occurrence. The surveys were aimed at verifying or "ground-truthing" the results from the analyses but concurrently led to the discovery of archaeological features that the interpretations of the image and LiDAR data failed to identify. Altogether, 40 features, encompassing 20 meat depots, 15 dwellings, three hunting pits and two rectangular structures, were found and mapped during the fieldwork. As the survey was intended to verify the interpretations from the analyses, these new discoveries were "incidental".

Of the surveyed features identified in the initial study of the LiDAR data, 73% were verified as archaeological features while 27% turned out to be natural features. As much as 55% of the features interpreted as dwellings turned out to be natural features, which constitutes 88% of all misinterpreted features. In comparison, 96% of the 52 field verified hunting pits were correctly interpreted during the analyses and only two turned out to be natural pits. In the second study of the LiDAR data, the ground truthing of hunting pits showed a success rate of 100%. It should also be noted that the methods used allowed a very clear picture of the respective features. The blending of methods, as suggested by Kokalj and Somrak [47], was particularly effective, more so than manually layering results on top of each other. This included identification of individual features of the hunting pits, in particular the interior shape of the bottom of the pits.


**Table 5.** Verification of features detected in the LiDAR and image data.

<sup>1</sup> Split interpretation; one feature was identified as a structure in the panchromatic image and a dwelling in the multispectral image; <sup>2</sup> One feature was interpreted as dwelling during the ground surveys.

> Verification of results from the image analyses showed that 94% of the features surveyed in both aerial and satellite images were archaeological features. Of the dwellings, 100% were archaeological features, while the verification rate for hunting pits detected in the satellite and aerial images was 100% and 96% respectively.

> When returning to the desk to study the models based on the LiDAR data again after the field verification, quite a few of the omitted archaeological features were actually identifiable as such when we were aware of their presence. That was the case with two of the three hunting pits and approximately half of the dwellings that were found during the fieldwork. Some of the surveyed meat depots were detectable as recessions in screes, but these are very hard to distinguish from natural recessions in the same screes.

> The detection success was only marginally affected by applying the two additional visualising techniques, "Local Relief Model" and "Sky-View-Factor". The gain was primarily their ability to enhance the visualisation of already identified features. The analysis related to the second LiDAR study where "Openness negative" and "Openness positive" were used as supplemental visualization techniques demonstrated the potential of the method to successfully identify features in areas with steep relief and large differences in slope, as also shown by Doneus [25]. The combination of these methods with "Sky-View-Factor" was very effective for detecting anomalies. However, detection was not dependent on these visualisations as more than half of the detected hunting pits were also visible using hillshading. A consideration is the fact that the reanalysis focused on the western half of the LiDAR dataset. As there is a partially overlapping dataset in this area, data from the second dataset was rebuilt and compared to the results of the dataset under discussion. The results of this analysis did not differ in terms of detection rate. This did, however, demonstrate that the combination of methods used in the reanalysis is robust. The results from the analysis of the overlapping dataset are not included here.

#### *4.3. Listed Cultural Heritage in the Study Area*

Askeladden, the Norwegian national cultural heritage database, currently contains information on 50 cultural heritage sites encompassing over 400 individual features within the study area [35]. The data have been gathered through traditional means, that is, ground based surveys. Many of the features, especially hunting pits that are a part of larger systems and dwellings located in the large coastal settlement sites, were not individually mapped when the sites were originally surveyed. As an example, one system (Askeladden Id 67178) located in the study area encompassing 55 pits is only mapped as a polygon in the heritage database [35]. The 55 pits are not individually mapped. Comparing our results with the archaeological data available from Askeladden, a significant number of previously unlisted or unknown features were detected, as shown in Figure 7. An important distinction, however, is that the sites and features listed in Askeladden are verified cultural heritage. Until verification, the features listed in Table 4 should be considered potential archaeological sites and features.

**Figure 7.** Sites listed in the national cultural heritage database Askeladden prior to this study are significantly outnumbered by the features detected through the LiDAR data and remote sensing images. Image: ©Norwegian Mapping Authority.

Field surveys (Table 5) showed that the conducted analyses are quite reliable for hunting pits for both remote sensing image data and LiDAR data. Based on the rate of right interpretation indicated in Table 5, it is also possible to extrapolate a theoretical figure for the number of features. For LiDAR study 1, a 96% rate of right interpretation indicates that the number of hunting pits within the study area is 740 and the number of dwellings 53, but that there are no mounds. For LiDAR study 2, the rate is 100%. Further validating the detection rate and consequently the applicability of LiDAR in this landscape, the results from analysing the LiDAR data were a starting point for subsequent surveys conducted by cultural heritage management authorities at Troms and Finnmark County. As shown in Figure 8, there are a far higher number of features listed in Askeladden in 2020 than in 2015. The newly listed features are mostly verified hunting pits.

**Figure 8.** Surveys following the LiDAR studies resulted in a significant increase from 2015 to 2020 in the number of sites and features listed in Askeladden. Image: ©Norwegian Mapping Authority.

When it comes to dwellings as well as the more indefinable structure and pit categorizations, the results are less impressive. Less than half of the verified features identified as dwellings in the DTM were correctly identified (Table 5). On the other hand, the results for the image data were 100%.

#### **5. Discussion**

#### *5.1. Assessing the Pros and Cons of LiDAR Data, Aerial and High-Resolution Satellite Images for Cultural Heritage Surveying*

A starting point for this discussion is the combined results from the analyses listed in Tables 4 and 5. According to these, the DTM yielded a significantly higher number of potential archaeological features compared to both satellite and aerial images. This is especially apparent when it comes to hunting pits. Although the picture is not entirely unequivocal, a substantial number of the hunting pits clearly apparent in the DTM are vague or even non-detectable in the image data as shown in Figure 4. As the suitability of LiDAR regarding hunting pits has previously been proven [23,45], this result is not entirely surprising.

Hunting pits generally form large systems, some of which can consist of hundreds of individual pits. Not only were a larger number of individual pits within systems identified, some smaller systems were only detected through the DTM. Subsequent ground surveys showed that the majority of the detected pits were correctly interpreted. As hunting pits generally are quite uniform circular pits of a regular size and depth, their characteristics mean that LiDAR is highly applicable for detecting this type of archaeological feature. As commented above, the methods used in the analysis of the rebuilt DTM went beyond mere detection. Individual features of the hunting pits could be detected, both the interior shape of the bottom of the hunting pits, as well as features relating to the mounds surrounding them.

The same feature characteristics also proved important when analysing the satellite and aerial images covering the same area. Apart from the regularities in feature characteristics, when hunting pits occur several together in a linear system, as they usually do, they can stand out quite clearly in the images (see Figure 4 for example). Single pits or systems where the pits are located at a distance from each other proved more difficult to identify.

The image analyses and results, especially the distribution of identified pits, indicate that detectability varies with topography and vegetation cover, detectability in this regard is site dependent. Forest and slopes reduce detectability, while hunting pits located in more open and flat landscapes are easier to identify. Vegetation is a critical factor that adversely affects the detection rate when analysing remote sensing images. The option of eliminating vegetation from a dataset by a filtering process is to a great extent the most prominent advantage of using LiDAR. In our case, this is very evident when it comes to detecting hunting pits. The inland parts of the study area are partly covered by birch trees and willow thicket, and the vegetation obscured the visual analyses of the satellite and aerial images.

As Table 5 shows, the ground surveys indicate a high rate of correct interpretation for hunting pits identified in either the LiDAR or image data. The results were not as conclusive regarding potential dwellings. A total of 118, 77 and 10 dwellings were detected in the DTM, the satellite and aerial image data, respectively. However, the rate of correct interpretation is seemingly significantly lower for dwellings than for hunting pits as the field surveys showed 45% of the verified dwellings detected in the DTM to be correctly interpreted. The image data provided better results in the sense that a higher percentage of the verified features proved to be correctly identified. However, considering what is known from the descriptions of the Stone Age settlement sites [35] within the study area, a rather low percentage of the remaining dwellings was detected through the analyses. Looking at the results, the characteristics of the detected dwellings are a key factor for detectability, both for LiDAR and image data. What is thought to be remainders of walls are generally apparent in either dataset. In addition, quite a few of the Stone Age houses are partly dug down and now appear as circular or rectangular depressions which may or may not be surrounded by a low wall.

Apart from feature characteristics such as size, shape and depth/height, discernible differences in the vegetation cover aided detectability when analysing the image data. The results from analysing the panchromatic and multispectral satellite images showed some differences in detection rate. Although the multispectral image provided a slightly higher number of potential features than the panchromatic image, the most notable differences are related to type of feature. The multispectral image provided somewhat better results regarding potential dwellings, pits and other structures linked to settlement sites. Dwelling remains constitute 19% of the features in the multispectral image and 13% of the features identified in the panchromatic image (Table 4). The difference is especially notable when it comes to features located in arable land near the coast. Of the 77 identified dwellings, 39 were detected in only one image. A total of 327 hunting pits were identified in the two image sets. Although nearly the same number of pits was detected in the two image sets and a substantial number of these pits overlap, 89 were only detected in either the panchromatic or the multispectral image.

Running two band combinations, one including the NIR1-band, proved advantageous for detecting and identifying dwellings. As previously mentioned, the NIR1-R-G bands allows for vegetation and variations in vegetation linked to type and condition to be detected. The band combination proved useful for detecting areas where the vegetation had been impacted by human activity, where features such as dwelling remains subsequently were detected. Vegetation and variations in the vegetation cover proved useful for detecting potential features and for determining where to look for archaeological features. This was especially true for the more recent traces of Sámi activity and settlement. Sámi settlement areas may show up quite clearly, as shown in Figure 9 where a large rectangular structure interpreted as a potential turf house of about 12.5 × 6 m was detected in a clearing. It is worth noting that, although this structure was very apparent in the multispectral image and quite easily detectable in the panchromatic image, as can be seen in Figure 9, it was not detected when interpreting the DTM. As this site was not among those surveyed during the field work, it is unclear whether feature characteristics play a role in detectability. As for hunting pits, running the RGB band combination alone proved adequate. The past activities affiliated with reindeer hunting have not left an easily detectable imprint on the vegetation in and around the pits.

Of the detected dwellings, 38 are part of extensive Stone Age settlements located near the coast. For these dwellings, vegetation differences generally seem less distinctive than for the more recent Sámi dwelling remains. The exception was the walls of features located in what is now arable grassland, which appeared somewhat highlighted compared to the surrounding area in the images. The vegetation on and around features which are not located in what is or has recently been agricultural areas does not appear to similarly "enhance" heritage features. Regardless, some dwellings were quite apparent in the images, most of which were found in coastal areas or near rivers in more inland low-lying areas. The disadvantages of a dense vegetation cover when analysing image data were less apparent as most dwelling remains are located in coastal areas, in what is now agricultural areas or cleared settlement areas. Along the coast of northernmost Norway, forest and shrubs are generally sparse or absent, leaving the landscape open and covered by grass and heath. This together with the fact that dwellings that are several thousand years old are still visible on the surface argues for the viability of image data in these areas.

**Figure 9.** The remains of a rectangular turf house is located within the area outlined in red. The dwelling is detectable in the panchromatic and multispectral satellite images. The panchromatic image is on the top right, while the multispectral image is shown on the bottom row with the bands NIR1-R-G on the left and R-G-B on the right. In the LiDAR data on the top left, the remains are not identifiable as a turf house. Note, however, the highly visible hunting pits on the left in the figure. Images: ©Digital Globe.

#### *5.2. The Applicability of LiDAR and Remote Sensing Image Data for Archaeological Surveying in the Subarctic and Arctic Landscapes of Northern Norway*

This study has shown that LiDAR and aerial and satellite images are highly applicable and useful for archaeological surveying in northern subarctic and arctic landscapes. The results have contributed to improve the inventory of known archaeological sites and features throughout the region. At the same time, the analyses and results also show that there are some limitations.

The listed cultural heritage sites within the study area represent various aspects of hunting and gathering, settlement and subsistence as well as religious beliefs and practices from the Early Stone Age into modern times. The national database [35] lists several hunting pit systems (some of which consist of up to around 50 pits), meat depots, aggregations of up to 50 house sites dating to the Stone Age, Stone Age artefact scatters (surface finds), Sámi turf houses, hearths (arran) (of which many are related to reindeer husbandry) and Sámi sacrificial sites. In comparison, our study yielded a limited number of cultural heritage types. Our analyses provided data on mainly two types of archaeological features, hunting pits for reindeer and dwellings, the last category encompassing features dating to the Stone Age and Early Metal Age as well as more recent centuries. To some extent this reflects local cultural history and existing cultural heritage data. The high number of detected hunting pits reflects the fact that hunting pit systems are known to be numerous in this area and is indicative of the importance of reindeer hunting for subsistence and economy in prehistoric and historic times. However, our results did not provide a general representation of cultural heritage sites in the area.

Our study showed remote sensing to be ill-suited for many of the feature types known to be found in this region. Hunting pits and dwellings are both archaeological features that are, comparatively speaking, quite sizable and of uniform design. A range of features linked to Sámi resource and landscape use, are in many cases anything but monumental and often blend in well with the landscape. The most visible remains of Sámi reindeer husbandry camp sites are often the stones delineating the hearths. Meat depots, which are often located in screes, appear as one depression among many and were indistinguishable from the surrounding area in both the DTM and image data. Stone Age activity and settlement areas have in many cases been located through surface artefact scatters which, so far, are undetectable by LiDAR and aerial and satellite remote sensing image data. However, as sites in this area often encompass a range of features, the results may contribute to narrow down the search area and provide a useful indicator of where to conduct further studies and ground surveys. A good example are the point-to-point surveys to selected detected features which also resulted in the discovery of another 40 archaeological features that the analyses of the image and LiDAR data failed to identify. Another example is the dwellings detected in coastal areas where extensive Stone Age settlements are known to be located. The fact that utilising LiDAR and remote sensing image data is shown to be applicable for narrowing down and consequently prioritising search areas is an important argument in favour of integrating such data and methodological approaches in archaeological surveying as they may contribute towards remote and hard-to-get-to areas being surveyed in a lowimpact and cost-efficient manner. Interpretation of images and a DTM based on visual analyses aided by processing and visualizing techniques is time-consuming but compared to the often much more time-consuming and costly field work, well worth the effort in these areas.

As mentioned above, the national database lists several hunting pit systems of varying size within the study area. Earlier surveys in the area generally mapped the extent of the systems and provided an approximate number of pits, but individual features were seldom mapped. This is also the case for several of the Stone Age settlement sites in the area. As the remote sensing data have been used to map individual features, they provide a valuable addition to already listed sites, especially for hunting pit systems where the rate of right interpretation is high. Remote sensing methods can thus add reliable and valuable detailed information of individual features within an archaeological site.

Visual interpretation, as used in this study, is a qualitative and seemingly simple process. Its main advantage compared to digital data processing is the fact that it can be carried out also when object features are not easily distinguishable [12]. However, in order to recognize patterns and characteristics, to extract meaningful data through visual analyses, features have to be identified and allocated into known categories. The practical and realistic use-value of LiDAR data and remote sensing image data is dependent on local conditions pertaining to surface conditions such as topography and vegetation cover, the characteristics of archaeological features and site characteristics such as size, depth/height, uniformity, building materials, placement and, not least, what we expect to find within a given area. Familiarity with cultural history, feature type characteristics and localization factors are of importance when conducting the analyses. As an example, hunting pits for reindeer and remains of Sami reindeer husbandry are within expectations for this area, consequently the data is searched for features with characteristics in line with what is known of archaeological features like hunting pits. Hunting pits generally occur in systems so once one or two potential pits are detected, one will generally detect more in the vicinity. If the interpreter is unaware that large-scale hunting of reindeer was based on hunting pit systems, a pit may very well be dismissed as a pit.

Current knowledge and understanding of the area's tangible cultural history is primarily based on ground-based surveys. Our study resulted in a high number of potential archaeological features far exceeding current knowledge and existing listed sites in the national cultural heritage database. Our analyses of aerial and satellite images as well as LiDAR data have shown potential, especially when combining results from the two sets of

remote sensing data. The advantages of a multi-sensor approach can be instrumental in improving current knowledge. In our case, the LiDAR data provided better results for the more vegetated inland areas, while the satellite images proved useful for detecting sites based on human inflected changes in the vegetation.

#### **6. Conclusions**

What can remote sensing contribute to archaeological surveying in subarctic and arctic landscapes? We have proven remote sensing utilizing LiDAR data and image data to be highly applicable for archaeological surveying in these northern landscapes. This study has led to the identification of 368 cultural heritage features based on aerial images, 437 on satellite images and 1186 based on a DTM covering the study area. Furthermore, ground surveys following the studies have resulted in a significant increase in the number of sites and features listed in the national cultural heritage database, Askeladden. We have shown that a multi-sensor approach contributes to high detection rates and thus can be instrumental in improving current knowledge. At the same time, our results showed that there are limitations. Feature characteristics, topography and vegetation cover proved key factors for detectability. LiDAR and remote sensing image data were shown to be ill-suited for some of the feature types common throughout the region. Even so, the results proved a useful indicator of where to focus attention and further survey activities.

We have contributed to improving the inventory of known archaeological sites in a subarctic and arctic landscape in a non-destructive and minimally invasive manner. A better overview of the archaeological remains in a landscape is crucial for a qualified cultural historical understanding of how the landscape was utilized by humans in the past. Furthermore, it is a precondition for optimizing cultural resource management and for sustainable cultural heritage management. Apart from the value of remote sensing for improving knowledge and understanding about the past and for efficient improvement of data quantity and quality, there is the additional value of remote sensing as a non-invasive means of mapping and monitoring cultural heritage. Considering the vast and hard-to-getto areas, especially in the north of Norway and other northern areas, the substantial costs and potential environmental impact associated with traditional archaeological surveying are a significant inducement to the current focus on further exploring the potential of remote sensing technologies. To further improve the applicability of remote sensing, a fruitful avenue for further research would be to focus on methodological development aimed at detecting and identifying a wider range of the cultural heritage sites and features commonly found throughout the region.

**Author Contributions:** Conceptualization, A.E.T.; methodology, A.E.T., O.R. and J.I.K.; validation, A.E.T., O.R. and J.I.K.; E.R.M.; formal analysis, A.E.T., O.R. and J.I.K.; investigation, A.E.T., O.R. and J.I.K.; writing—original draft preparation, A.E.T., O.R. and J.I.K., S.B.; writing—review and editing, A.E.T.; visualization, A.E.T.; supervision, A.E.T.; project administration, O.R., A.E.T.; funding acquisition, O.R. and E.R.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by The Research Council of Norway, Bruk av avansert teknologi for å kartlegge, forstå, konservere og forvalte kulturarven, grant number 208439 of which this study was a part.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author, the second and third co-author.

**Acknowledgments:** The authors would like to thank colleagues at NIKU's High North Department for constructive feedback. Three anonymous reviewers are kindly acknowledged for helpful comments which contributed to improve this paper.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


### *Article* **Remote Sensing for Cultural Heritage Assessment and Monitoring: The Case Study of Alba Iulia**

**Cristian Moise 1,\*, Iulia Dana Negula 2, Cristina Elena Mihalache 1, Andi Mihai Lazar 2, Andreea Luminita Dedulescu 2, Gabriel Tiberiu Rustoiu 3, Ioan Constantin Inel <sup>3</sup> and Alexandru Badea <sup>2</sup>**


**Abstract:** In recent times, satellite-based remote sensing has a growing role in archaeology and inherently in the cultural heritage management process. This paper demonstrates the potential and usefulness of satellite imagery for the documentation, mapping, monitoring, and in-depth analysis of cultural heritage and the archaeological sites located in urban landscapes. The study focuses on the assessment and monitoring of Alba Iulia, which is one of the Romanian cities with the richest historical past. Multitemporal analysis was performed to identify the land use/land cover changes that might contribute to an increased cultural heritage vulnerability to natural disasters. A special emphasis was dedicated to the assessment of the built-up area growth and consequently of the urbanization trend over a large time interval (30 years). Next, the urbanization and urban area expansion impact was further analyzed by concentrating on the urban heat island within Alba Iulia city and Alba Iulia Fortress (located in the center of the city). As temperature change represents a key element of climate change, the temperature trend within the same temporal framework and its impact on cultural heritage were determined. In the end, with regard to the cultural heritage condition assessment, the research was complemented with an assessment of the urban ground and individual building stability, using persistent scatterer interferometry. The results contribute to the detailed depiction of the cultural heritage site in such a manner that the site is monitored over an extensive timeframe, its current state of conservation is accurately determined, and the future trends can be identified. In conclusion, the present study offers reliable results regarding the main factors that might endanger the cultural heritage site as a basis for future preservation measures.

**Keywords:** remote sensing; Earth observation; satellite imagery; multi-temporal analysis; urban heat island; persistent scatterer interferometry; long-term monitoring; cultural heritage assessment; Alba Iulia (Apulum)

#### **1. Introduction**

The increasing contribution of remote sensing to the archaeology and cultural heritage domain is favored by the variety of satellite sensors, data policies that enable free and open access to satellite imagery, and the use of combined techniques in order to maximize the spectrum of complementary information that can be generated through the exploitation of satellite data for a more accurate and extensive analysis of the archaeological sites and their surroundings [1].

Using very high resolution (VHR) multispectral satellite imagery, Megarry et al. [2] promoted the use of automatic methods (i.e., machine-learning and classification techniques) that strengthen the archaeological prospection process. Likewise, using VHR

**Citation:** Moise, C.; Dana Negula, I.; Mihalache, C.E.; Lazar, A.M.; Dedulescu, A.L.; Rustoiu, G.T.; Inel, I.C.; Badea, A. Remote Sensing for Cultural Heritage Assessment and Monitoring: The Case Study of Alba Iulia. *Sustainability* **2021**, *13*, 1406. https://doi.org/10.3390/su13031406

Academic Editor: Ionut Cristi Nicu Received: 15 December 2020 Accepted: 25 January 2021 Published: 29 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Synthetic Aperture Radar (SAR), Balz et al. [3] determined that remote sensing imagery represents a useful tool for the analysis of the archaeological landscape. Steward et al. [4] used both VHR SAR and multispectral satellite data and assessed their suitability for archaeological prospection and monitoring. The capabilities of satellite remote sensing beyond site detection were proven by Borie et al. [5] in a study that focused on the analysis of specific environmental resources. Malinverni et al. [6] demonstrated the advantage of using remote sensing for the documentation of an archaeological complex that further served as support for perseveration and enhancement efforts. Satellite remote sensing plays a crucial role also in the condition assessment [7] and detection of threats to heritage sites [8,9]. An overview of the innovative products that can be generated using remote sensing data and techniques was presented by Bassier et al. [10].

In this context, the present study expands the knowledge related to the use of satellite remote sensing for archaeological research by focusing on a cultural heritage site located in an urban environment. Satellite imagery enables a better understanding of the environmental changes, thus contributing to the development of enhanced urban management methodologies [11]. Related to the very challenging characteristics of the urban environment, Gamba et al. [12] demonstrated the importance of the combination of different types of data with regard to spatial resolution or diverse polarizations, multi-resolution, multi-scale, and multi-temporal techniques, whereas Wellmann et al. [11] identified that the classification of satellite imagery is one of the most used remote sensing techniques for urban area analysis. The study has three main pillars, namely, (1) the analysis of the land use/land cover (LULC) changes within the administrative boundary that contains the cultural heritage site, with a special emphasis on urbanization and urban growth from 1988 to 2018; (2) the evaluation of the temperature trend in correlation with the urbanization process as well as of the urban heat island (UHI) effect between 1988 and 2019; and (3) the assessment of the ground and individual building stability in order to have a more profound insight into the cultural heritage site condition from 2018 to 2020. Hence, the main goal of the research is represented by the assessment and long-term monitoring of cultural heritage as well as the analysis of cultural heritage vulnerability to natural disasters, for sustainable management measures. The research integrates multitemporal analysis and supervised classification using the maximum likelihood algorithm, based on Landsat and Sentinel-2 images; in addition, Copernicus land cover data sets (i.e., CORINE Land Cover, Urban Atlas) complement the study. The land surface temperature (LST), brightness temperature at sensor, and UHI analysis for thermal stress estimation encompass Landsat data and ground-truth temperature measurements for the validation of the results. Finally, the ground and building stability analysis was performed based on a large series of Sentinel-1 images and using the persistent scatterer interferometry technique (PS-InSAR).

Regarding the specific methods and techniques that are exploited in this study, the scientific literature contains plentiful relevant publications focused on cultural heritage sites, for example: the identification and surveying of land-cover changes [13], the analysis of the spatiotemporal land use/land cover evolution based on Landsat imagery [14], the use of supervised classification for the generation of land cover products based on multispectral satellite images with different spatial resolution [15], the impact of urbanization [16], the use of a CORINE Land Cover (CLC) dataset for the analysis of the urban sprawl [17], and the application of PS-InSAR as a standalone technique or in combination with other remote sensing approaches for the investigation of ground and structural stability [18–22]. The "Copernicus services in support to Cultural Heritage" report, published by the European Commission (EC), refers to LST and the Urban Atlas (the core service products provided by the Copernicus Land Monitoring Service—CLMS) as data useful for cultural heritage management [23]. According to the report, LST enables the detection of thermal anomalies, whereas the Urban Atlas contributes to the mapping of surrounding infrastructure in cultural heritage sites. Some of the research results obtained in this study (e.g., the UHI output) might represent examples of additional products suitable for cultural heritage applications.

The study concentrates on the monitoring and condition assessment of the Alba Iulia Fortress (Figure 1), located in the heart of Alba Iulia city in Romania. For a comprehensive analysis, the research includes the monitoring of the Alba Iulia administrative territory, which has an area of about 102 km<sup>2</sup> (Figure 2).

**Figure 1.** Orthophoto image of Alba Iulia City (© image source: Romanian National Agency for Cadaster and Real Estate Advertising (ANCPI) geoportal—https://geoportal.ancpi.ro/). The sevenpoint star-shaped medieval citadel represents a unique Transylvanian landmark.

**Figure 2.** Location of Alba Iulia city in the central-western part of Romania (© image data source: OpenStreetMap—OSM, geo-spatial.org, Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM) courtesy of the United States Geological Survey (USGS)).

Amongst the factors that might affect a heritage site are urban high rise/urban sprawl, changes to skyline, commercial and industrial development, the use of transport infrastructure, air pollution, environmental or biological factors (e.g., wind, humidity, temperature, dust, water, micro-organisms), climate change (e.g. temperature rise), natural disasters (e.g., floods, earthquakes, landslides), and human activities [24].

Some of these factors have an impact on the Alba Iulia cultural heritage site, hence their early identification and assessment strengthen the cultural heritage safeguarding and conservation measures. For example, the assessment of the urban development trend and its consequences (e.g., increased vehicle traffic, elevated air pollution levels) on the cultural heritage is of high importance for the management authorities since both the authorized and the illegal urban expansion might endanger the ancient Roman city. An accurate temperature trend estimation and delineation of the urban heat islands is a critical parameter for the management authorities, considering that for the optimal preservation of the cultural heritage a lower temperature is preferred. In Alba Iulia, a higher temperature correlated with higher humidity was observed in numerous archaeological sites within the fortress, thus affecting the buried and the ground structures. In addition, taking into account that some of the Alba Iulia's cultural heritage buildings have visible displacements, a remote sensing analysis of the buildings' stability provides the local authorities an accurate overview that might represent the basis for further detailed investigations. At present, the natural disaster risk is minor in the Alba Iulia administrative territorial unit. Catastrophic floods occurred in 1851, 1857, 1870, 1887, 1893, 1970, and 1975 until embankment works were performed on the Mures, River. However, a prevailing problem is that the eastern and southern parts of Alba Iulia city are still affected by high groundwater levels [25]. In order to address the aforementioned aspects, the study was built on three main pillars, each of them capitalizing on different remote sensing technologies and methods, as described in Section 2.3. Methodology.

#### **2. Materials and Methods**

#### *2.1. Description of the Test Site: Alba Iulia (Ancient Apulum)*

The earliest traces of human presence on Alba Iulia's territory are dated to the Neolithic [26], and the most representative site of this period also gives its name to the "Lumea Nouă" culture [27]. Alba Iulia had several major administrative roles throughout history: Apulum (the ancient name of Alba Iulia) was the capital of Roman Dacia [28], later the city was for a short period the capital of the Voivodeship of Transylvania [29,30], then the capital of the Principality of Transylvania [31], the capital of all Romanian medieval states in 1600 [32], the bastion of Austrian rule in Transylvania, the capital of the Great Unification from 1 December, 1918 (which became the national day of Romania), and the coronation capital [33]. All of these administrative functions had an impact on the architectural heritage that is largely preserved today.

In Roman times, Alba Iulia was the place of a legionary fort [34], and also of two cities with the same name that were granted colonial status (the highest a city from Roman Dacia could get); this was a unique situation in the province of Dacia, being also rarely encountered in the entire Roman Empire [28]. Several elements of the former fort of the Legio XIII Gemina are still preserved and can be visited: parts of the enclosure walls [34], the gate (*porta*) *principalis dextra* [35], the *via principalis*, and some of the *principia*.

The Roman fort was most probably built sometime between AD 107 and 108. The structure has a typical rectangular shape and the dimensions of 480 × 440 m, based on the measurements made during the investigations from 2011 [36].

Archaeological investigations conducted in 2011–2012 partially unearthed the *via principalis* and *principia* (Figure 3). A sector of *via principalis* with A length of 7.20 m and A width of 7.30 m was identified on that occasion. The road consisted of gravel bedding, which was probably covered with sandstone slabs. A drain made of bricks, with a depth of 0.30 m, was set under the slabs. The central lane of *via principalis* was reserved for carts, being delimited by two parallel rows of limestone slabs, and the side lanes were used by pedestrians.

**Figure 3.** Alba Iulia: (**a**) *via principalis* and (**b**) *principia* (© image source: the personal photo archive of G. Rustoiu).

The fort's *principia* was partially investigated, allowing the identification of its main components, while the unexcavated ones were graphically reconstructed. Archaeological excavations identified parts of the *atrium*, *armamentarium*, *basilica*, and *tabularium.* The entire *principia* most likely measured 65 × 100/110 m, whereas the investigated area of the inner courtyard—the *atrium*—measured 50.85 × 38.34 m, though the reconstructed total dimensions were 50 × 50 m [36]. The *atrium* was initially covered with a layer of sand and was paved later with stone stabs. The *basilica*, measuring 65 × 30 m, was identified on the western side of the principia, alongside the *aedes*, a central room paved with bricks, which measured 11.22 × 9m[36]. The *aedes* and its nearby rooms were delimited by the portico. The *armamentarium* was partially unearthed on the southern side of the *principia*, one of the rooms having a hypocaust. The entrance to the *praetorium* was on the eastern side, where four bases built of stones and bricks were identified [36].

Following the conquest of Transylvania by the Habsburg Empire and the peace of Karlovitz, Alba Iulia lost its political role but benefited from the construction of the bastion fortress of the Vauban type [37], which is still preserved almost completely today. The bastion fortress has the shape of a star with seven corners, defined by the presence of seven bastions. Aside from them, the fortress also includes a central fort, the former Roman fort, and the medieval fortress, ravelins, banquettes, counterguards, tenailles, and two rows of concentric moats with counterscarps [37]. The structure was primarily built between 1715 and 1738 to serve as the main bastion of the Austrian rule in Transylvania. However, the fortress lost its military importance during construction due to the evolution of European politics, so it was invested with a mainly administrative function that allowed the inclusion of many architectural elements of Baroque origin: four of the six gates (Figure 4a) and the bastions [38], which gave it a particular character among the other bastion fortresses of the Vauban type from Europe. The fortress (the part preserved today) was built largely following the plans of the Italian architect Giovanni Morando Visconti. The Austrian Empire left its mark on the fortress area by building a number of civilian and military institutions and public monuments, one of the most representative being the Military Club—now the Union Hall. It was built between 1898 and 1900 for the local officers. Its ballroom was the place where representatives of the Romanians from Transylvania, Banat, and Partium decided on the unification with Romania on 1 December, 1918 [39]. Another representative monument from the fortress at Alba Iulia is the Coronation Cathedral (Figure 4b) which was built for the coronation of King Ferdinand I and Queen Maria as sovereigns of Greater Romania on 15 October, 1922 [40]. The cathedral was built in 1921–1922 in the Neo-Brâncoveanu style [41]. The plan has the shape of an inscribed Greek cross measuring 53 × 18 m, with a height of 45 m [42]. The church is set into a pavilion enclosure whose western entrance is defined by a bell tower with a height of 58 m [43].

**Figure 4.** Alba Iulia: (**a**) the third gate of the fortress and (**b**) Coronation Cathedral (© image source: personal photo archive of G. Rustoiu).

From the urban planning perspective, Alba Iulia had two distinct areas, namely, a fortified area that included the bastion fortress of the Vauban type and an open urban area that was developed in Roman times within two different locations that evolved into two attested settlements, specifically, Colonia Aurelia Apulensis and Colonia Nova Apulensis. The open urban area was developed mainly around the fortress, on the western, southern, and northern sides, almost entirely at approximately 20 m above the floodplain of the Mures, River [44]. Starting in the 18th century, the draining works performed on the eastern side of the fortress enabled the reconfiguration of the open urban area. Consequently, during the 18th and 19th centuries, the open urban area was relocated mainly to the eastern, southeastern, and northern sides of the fortress [45]. Except for the buildings located in the protected area corresponding to the fortress, the rest of the Alba Iulia buildings were built starting in the middle of the 18th century and predominantly between the end of the 19th century and the beginning of the 20th century [45]. Nowadays, the eastern and southern parts of the city are affected by the high groundwater levels that are caused by the initial marshland area that was drained [25].

#### *2.2. Satellite Imagery and Additional Data Used in the Study*

The case studies exploited satellite data acquired by different missions (i.e., Landsat, Sentinel-1, and Sentinel-2), Copernicus land cover datasets (i.e., CLC, Urban Atlas, and Urban Atlas Change), and ground-truth measurements (i.e., temperature at meteorological ground stations). The Landsat and Sentinel-1/-2 images were obtained free of charge from Earth Explorer data portal [46] and the Copernicus Data Access Hub [47], respectively.

In the first case study (i.e., the multitemporal analysis of the land cover changes), 8 Landsat (30 m spatial resolution) and Sentinel-2 (10 m spatial resolution) images acquired between 1988 and 2018 were used. These spatial resolutions were adequate for the identification of the patterns and the evaluation of the land use/land cover at a regional scale [48].

For the second case study (i.e., the thermal stress and UHI analysis), 31 Landsat images collected between 1988 and 2019 were analyzed (one image was selected for each year of the investigated time interval). The climatic conditions of Romania, which is located in the eastern part of Central Europe, between 43◦ and 48◦ northern latitude and 20◦ and 29◦ eastern longitude, cause August to be the hottest month of the year. For the years when cloud cover was present in the images from August, images from July and September were used in order to fill the gaps. The thermal band 6 of Landsat 5/7 and band 10 of Landsat 8 were used within the study. During the Level 1C processing step, all the thermal bands were resampled to 30 m, thus facilitating the processing of long-term time series.

Within the third case study (i.e., the ground and structural stability analysis), 37 Sentinel-1A (S1A) and 37 Sentinel-1B (S1B) C-band SAR images collected from 2018 to 2020 were exploited. The data were acquired using the Interferometric Wide Swath (IW) mode, in Single Look Complex (SLC) format with VV polarization and a repeat cycle of 12 days, from descending (S1A) and ascending (S1B) orbits.

Additional land cover datasets provided by CLMS [49] were correlated with the multitemporal analysis, thus strengthening the results' validation. CLC provides consistent and thematically detailed information on land cover and land cover changes across Europe [50,51]. The Urban Atlas provides reliable, inter-comparable, high-resolution land use and land cover data for European urban areas [52]. In this study, CLC datasets corresponding to the years 2000, 2006, 2012, and 2018 and Urban Atlas Change datasets covering the intervals of 2006–2012 and 2012–2018 were used.

For the validation of the thermal stress study and UHI analysis, temperature data recorded by the meteorological ground station located in the center of Alba Iulia were employed, courtesy of the National Meteorological Administration.

#### *2.3. Methodology*

2.3.1. Multitemporal Analysis of the Land Cover Changes through Supervised Classification

The dynamic modifications of the land cover can be identified and monitored using satellite imagery [14]. This analysis aims to identify the risk factors of land cover changes that can affect the cultural heritage site. Based on multitemporal satellite data, the case study presents the evolution of the land cover within the study area over a time frame of 30 years. For a comprehensive analysis of a cultural heritage site, Bai et al. [14] demonstrated that the monitoring of its surrounding area is essential. Therefore, the study area was represented by the Alba Iulia territorial administrative unit, which includes the following component localities: Alba Iulia, Bărăbant, , Mices,ti, Oarda, and Pâclis, a.

Supervised classification of the multitemporal images was performed for the identification of the land cover changes, as this approach is suitable for the monitoring and management of the built-up areas [53]. The processing of the Landsat and Sentinel-2 images was performed using the ArcGIS Software. Considering the human activities and the development of functional areas, five main classes were defined and analyzed, namely: Arable land, Grassland, Forest, Built-up area, and Water course. The land cover maps were obtained through the supervised classification of the time series composed of the 8 satellite images, using the maximum likelihood method. The main steps of the supervised classification include the subset of the area of interest, the definition of the training samples, and the generation of the classification results.

The maximum likelihood algorithm represents the most common technique used in the supervised classification process [54,55], which involves the spectral information of each pixel to generate the land cover map. Maximum likelihood is a parametric statistical method that uses training zones defined by the analyst for the classification of the entire scene. The distribution of the training data within each class sample in the multidimensional space is assumed to be normal (Gaussian). Using probability density functions, each pixel of the image is allocated to the class to which it has the highest probability of being a member [56].

The accuracy assessment of the supervised classification represents a significant step in the classification process. A confusion matrix to determine the overall accuracy of the classification was employed. This calculates the accuracy of the classification by comparing the pixels of a reference image considered the ground truth with the classified pixels [57]. The pixels corresponding to all possible correlations between ground truth and the classified results are in the cells of the tables. The columns represent the true classes, whereas the rows correspond to the predicted classes. The correctly identified pixels are distributed on

the major diagonal of the matrix and the non-diagonal cells represent classification error values. The overall accuracy was calculated for each image using the following formula:

$$\text{Overall accuracy} = \frac{\text{the sum of the correctly identified pixels}}{\text{the total number of pixels}} \times 100\tag{1}$$

Next, the land use/land cover changes were analyzed (also including the supplementary land cover datasets), with a special emphasis on the evolution of the urban extension, as the process might affect the cultural heritage site and cause significant changes in temperature values as well [14].

#### 2.3.2. Evaluation of the Thermal Stress Using UHI Analysis

Urbanization creates modifications in the local low atmosphere because of the changes in energy exchange rates. UHI is described as the difference in temperature between the impervious areas and the surrounding rural areas. Urban development creates large areas with low albedo values that store thermal radiation. The phenomena can be most clearly noticed at nighttime when the temperatures in the impervious areas remain higher than in the surrounding areas [58]. Increasing temperatures can have a negative impact on the well-being of citizens, especially in areas where heat waves are frequent. UHI can increase the negative impacts of heat waves, threatening citizens with health issues. The goal of this case study was to check whether in the past decades urban growth combined with global climate change contributed to UHI temperature rise and to identify the areas of Alba Iulia that are subject to the most significant changes.

It is predicted that global climate change will increase the frequency and extension of heat waves around the globe [59]. Combined with the rapid expansion of urban areas, the effects will be more intense than ever before. Living in such conditions will require more resources in order to maintain regular day-to-day activities, mainly to keep indoor habitation spaces at normal temperatures.

UHI can be determined using the thermal infrared band of optical remote sensing satellites [60]. Compared to the in-situ measurements, in which case a limited number of areas can be monitored, satellite remote sensing gives the possibility of obtaining precise measurements over large areas with minimum costs. The availability of the Landsat archive enables the long-term identification of the temperature trend.

After calculating the brightness temperature at sensor, the next step in the assessment of UHI is represented by the delineation of the city limits using Landsat Natural Color Composite images for the years that recorded the biggest changes in the urban footprint. After analyzing the satellite images for the entire investigated period, five representative years were selected: 1994, 2000, 2006, 2010, and 2018. The urban footprints were extracted in vector format (polygon type) using ArcGIS.

For the city limit corresponding to the year 2018, a centroid was created and used as the base for generating a buffer of 7.5 km that enclosed the city (Figure 5). The dimension of the buffer was chosen to be 7.5 km in order for the resulting vector to be larger than the extremities of the city. After the buffer creation, a total of 100 random sample points were created inside the clipped buffer, but outside the city limits. The resulting points were used to calculate the mean temperature outside the city limits for each of the 31 Landsat satellite images exploited in this study. In order to calculate the UHI, the values previously calculated were subtracted from the land surface temperature product.

The goal of the approach was to verify whether there was a correspondence between the expansion of the city limits and its UHI, and temperatures inside the historically important Alba Iulia Fortress. In order to check this hypothesis, another mask was created inside the Alba Iulia city limits to delineate Alba Iulia Fortress. In order to extract and compare the mean values of the temperatures corresponding both to Alba Iulia city and Alba Iulia Fortress, the Zonal Statistics function from the ArcGIS 10.4 software was used. The Pearson correlation coefficient was computed for the statistical analysis of the results.

Next, the validation step was performed using in situ temperature data provided by the National Meteorological Administration.

**Figure 5.** Urban Heat Island assessment workflow.

#### 2.3.3. Assessment of the Ground and Structural Stability Based on PS-InSAR

In this case study, PS-InSAR was considered the most suitable technique for the monitoring of Alba Iulia, based on multitemporal SAR data. PS-InSAR enables the estimation of the displacement velocity at the millimeter level and yields the best results when applied to built-up areas.

The interferometric processing of the SAR images was accomplished using the SARscape tool of the ENVI software. The PS-InSAR technique implies the identification of single coherent pixels [61], named Persistent Scatterers (PSs), for the calculation of their displacement velocity within a large series of SAR data. The free access to Sentinel-1 imagery was the key element of this case study. For a higher accuracy of the displacement estimation and lower atmospheric effects, the stack of SAR images was subject to intensive computation steps that are described hereunder.

Firstly, the Sentinel-1 SLC images were imported and the satellite positions were corrected using the Precise Orbit Determination (POD)—Precise Orbit Ephemerides provided by European Space Agency (ESA). Next, the area of interest was delimitated (i.e., the Alba Iulia administrative boundary) and the satellite images were cropped accordingly. Afterward, an analysis of the spatial and temporal baselines was carried out for the identification of the pairs that exceeded the thresholds and of the most suitable pairs that had the lowest baselines [62].

The displacement in the Line of Sight (LOS) of the satellite could be measured using the phase difference of two SAR acquisitions. The phase difference is affected by various factors that contribute to the phase noise, namely, the geometric decorrelation, the temporal decorrelation, and the atmospheric disturbance [63]. To reduce the phase noise, the SAR images were co-registered into interferograms and multilooked, thus obtaining a better radiometric resolution [63].

The software estimates the displacement velocity in two mandatory steps called the first and second inversion. The first step inversion achieves the residual height and displacement velocity [64]. The radar signal characterized by high coherence was identified in order to analyze its phase history. The selected PS candidates had stable temporal behavior [65]. In the second step inversion, the noise due to the atmospheric component is calculated and subtracted from the phase [63].

The last step is represented by the georeferencing process. The output is a file containing the PSs symbolized as points. In this step, a coherence threshold of 0.65 was applied, considering that coherence ranges from 0 to 1, starting from low data quality to high quality [66]. The LOS and the vertical displacements were generated and projected into the WGS84 cartographic system. Due to the looking angle, the satellite observes and measures the vertical displacements more accurately than the horizontal ones [62]. If the PSs are

moving towards the satellite, the displacement values are positive, and if the movement is away from the satellite the values are negative.

The LOS measurements were transformed into the vertical displacements using the following equation [67]:

$$
\Delta s = \frac{\Delta R}{\cos \theta (inc)} \tag{2}
$$

where Δ*s* is the surface displacement in the vertical direction, Δ*R* is the LOS displacement, and *θ*(*inc*) is the incidence angle.

The great advantage of using satellite remote sensing data in general and the PS-InSAR technique in particular is that investigations can be conducted over large areas, the temporal sampling is suitable for many applications, and accurate results can be immediately obtained. In comparison, the ground techniques that can be used to obtain the same outcomes as those resulting from the PS-InSAR require a more considerable amount of time and effort. The PS-InSAR technique allows for the identification of the potential vulnerable areas and proved to be a crucial tool for a more effective cultural heritage management.

#### **3. Results**

#### *3.1. Case Study (1): Multitemporal Analysis of the Land Cover Changes through Supervised Classification*

The overall accuracy results are represented as a percentage in Table 1, where 100% represents an impeccable classification. The accuracy assessment indicates an overall classification accuracy of 94.59%, which represents a high correspondence between the supervised classification results and the reference image.

The time series supervised classification facilitates the identification of significant urbanization expansion (Figures 6–9), a process that was presented as a threat to the cultural heritage by Agapiou et al. [16] and Drahor [68]. Over the years, Alba Iulia had an impressive evolution of both the industrial sector from 1970–1990 and the development of residential areas after 1990 and especially after 2000, in the northwest of the city and in the area of Pâclis, a to the southwest [69].



**Figure 6.** Land cover classification of satellite imagery from 1988 and 1992.

**Figure 7.** Land cover classification of satellite imagery from 1996 and 2000.

**Figure 8.** Land cover classification of satellite imagery from 2004 and 2006.

**Figure 9.** Land cover classification of satellite imagery from 2012 and 2018.

For the investigated period, the main differences were identified in the class of the built-up area that recorded an increase from 9.7 square kilometers in 1988 to 11.7 square kilometers in 2018. When comparing the first four years of analysis (1988, 1992, 1996, and 2000), the changes in built-up areas were not very visible. The high urban expansion started with the year of 2000, especially in the peri-urban area where residential districts were developed. In close connection with the increase of the residential area, service spaces and commercial, production, or storage areas were also developed [69]. The most remarkable increased percentage value from the Built-up area class was recorded between 2000 and 2004, with 4.1% (Figure 10).

**Figure 10.** Percentage distribution of the land cover classes (investigated time frame: 1988–2018).

To quantify the built-up growth, the area covered by each class was calculated after the transformation of the raster resulting from the supervised classification into polygons. Overall, the built-up area had an ascending trend. These trends were confirmed by the values of CLC data, shown in Figure 11. CLC vector data acquired at Level 3 from the Copernicus Land Monitoring Service was reclassified according to the CLC nomenclature to Level 2 in the following classes: Arable land, Urban fabric, Forest, Pastures, and Inland waters. An overall change map of the area was created based on the Urban Atlas Change data (Figure 12).

The environmental changes in the Alba Iulia territorial administrative unit were directly related to human activities. The risk represented by human pressure can be monitored and managed for safeguarding the cultural heritage area. Remote sensing data are indispensable in this dynamic process of urbanization because an early update of information enables the continuous and comprehensive monitoring and control of urban expansion [70].

The information related to land cover changes and their evolution is essential for a sustainable management of the cultural heritage area and its surroundings. In addition, it represents the basis for the anticipation of human actions and consequently the reduction of the risk that can affect the cultural heritage area.

The results of the multitemporal analysis represent relevant information for the Alba Iulia territorial administrative unit over a long period of time. The General Urban Plan of Alba Iulia encompasses not only the fortress or the city, but also the surrounding areas. Consequently, the corresponding archaeological area has a greater extent. As mentioned before, the monitoring of the surrounding areas and the identification of the land cover changes provide information about the factors that may affect the cultural heritage site in

the future. Specifically, urban expansion also causes increased archaeological intervention risks to the terrains that are subject to a change in land use. These outcomes may be used in the forthcoming development of management strategies in order to protect and conserve Alba Iulia Fortress and the valuable archaeological sites.

**Figure 11.** Evolution of the land cover classes according to the supervised classification and CORINE Land Cover (in 2000, 2006, 2012, and 2018).

**Figure 12.** Distribution of urban changes in the Alba Iulia territorial administrative unit between 2006 and 2018.

It could be argued that the results do not have the required level of detail for the scale of the cultural heritage site (i.e., Alba Iulia Fortress). Indeed, the limitation of this use case is represented by the fact that the small-sized archaeological features could not be monitored given the spatial resolution of the satellite images used in this study. In addition, the spatial resolution hindered the identification of the illegal constructions erected within the cultural heritage protected area. This information is of critical importance since these structures destroy the cultural heritage site. However, the study focuses on and proves that

the monitoring of the surroundings is essential, as the expansion of Alba Iulia city has an undisputable impact on the fortress located right in its center. The results are satisfactory and represent pieces of relevant information that reach their full potential in a broader context, together with other elements that are important to the process of cultural heritage management.

#### *3.2. Case Study (2): Evaluation of the Thermal Stress Using UHI Analysis*

The delineation of the city limits corresponding to the years 1994, 2000, 2006, 2010, and 2018 is presented in Figure 13.

**Figure 13.** Evolution of Alba Iulia city limits over the investigated time frame (1994–2018); background satellite image: Sentinel-2 (© European Space Agency).

During the investigated years, the area of Alba Iulia city doubled in size. In consequence, the UHI footprint expanded accordingly.

The spatial distribution of the UHI over Alba Iulia city and Alba Iulia Fortress during the study period is illustrated in Figure 14. The results were obtained by calculating the average value of every UHI product generated for the study period. Variations could be observed within the limits of the city, with peaks in areas corresponding to industrial and commercial districts with values as much as 8 ◦C higher than the surrounding areas. In the case of Alba Iulia Fortress, the fortification walls acted as a barrier, especially due to the covering vegetation layer. The temperatures measured in the area of the fortification walls averaged 1 ◦C lower than the surrounding areas. Two hotspots caused by the presence of large commercial and residential areas with tall apartment blocks were present to the north of the fortress. The data also shows that, in terms of habitation, the areas less affected by thermal discomfort were the outskirts, with the exception of the northeastern area, which is highly industrialized.

**Figure 14.** Average temperatures of the UHI during the study period (1988–2019); background satellite image: Sentinel-2 (© European Space Agency).

The mean temperature values generated from the processing of the Landsat time series are presented in Figure 15 for Alba Iulia city and Figure 16 for Alba Iulia Fortress. By examining the two charts representing the temperatures both for Alba Iulia city and Alba Iulia Fortress, it can be observed that the temperatures had a rising trend within the study period (about +2 ◦C for the city and +1.5 ◦C for the fortress).

**Figure 15.** Temperature trend inside Alba Iulia city.

**Figure 16.** Temperature trend inside Alba Iulia Fortress.

The two resulting data series representing the evolution of temperatures in the study area were statistically analyzed using R programming language with the help of RStudio, an easy-to-use, open-source, integrated development environment (IDE) software. For this statistical validation, the Pearson linear correlation coefficient was calculated in order to investigate whether there was a correlation between the rise of temperatures in the city and in the fortress. By adding the two Excel files with the temperatures for the city and for the fortress, a graph representing the two data series was plotted. Simultaneously, the "r" coefficient was calculated. The results showed a high positive correlation [71] between the two data series, with a value of 0.83 for the Pearson "r" coefficient (Figure 17).

**Figure 17.** Correlation between temperatures within Alba Iulia city and Alba Iulia Fortress.

The validation of the results was performed by correlating the satellite remote sensing outputs with the temperature data measured by the meteorological ground station located in the center of Alba Iulia city, southeast of Alba Iulia Fortress. Practically, the correlation between the temperature data recorded by the National Meteorological Administration's station and the temperature values of each pixel of the satellite image time series that overlapped with the position of the ground station was statistically analyzed (Figure 18). The results were satisfactory, showing a moderate positive correlation [71] between the two data sets with a Pearson "r" coefficient of 0.61.

**Figure 18.** Correlation between the Landsat brightness temperatures and the in situ measurements (1988–2019).

It could be argued that since the results of the study represent the LST brightness temperature at sensor and not the actual temperature of the land surface, validating them based on the in-situ temperature data is unfeasible. However, the study aimed to identify the temperature trend over the years when the city expanded considerably and corroborate it with the trend of the temperatures recorded by the local weather station over the same time interval.

The UHI analysis has major relevance for preservation measures [72]. It is well known that vegetation induces lower temperatures that are beneficial for a better preservation of cultural heritage by "avoiding overheating on surfaces, reducing thermal stress and also damage due to rising damp and saline efflorescence" [73]. Higher temperatures affect underground and aboveground structures. In addition, higher temperatures lead to higher humidity, as was noticed in numerous archaeological sites within Alba Iulia Fortress. For Alba Iulia's cultural heritage, the optimal temperature is a maximum value of 22 ◦C and the optimal humidity is in the range of 40–60%. In the *principia* the humidity could reach 90%, thus causing the decomposition of the structures. The UHI analysis provides a detailed overview of the heat islands and this information can be directly used by the cultural management authorities to apply corrective actions. For example, the vegetation areas that cover the fortification walls can be extended in order to maintain lower temperatures. Likewise, the areas that are affected by higher temperatures require the planting of vegetation or the extension of the protection areas. In correlation with the degree of physical structure degradation, the planted vegetation should be grass or other types of vegetation that would not further deteriorate the cultural heritage (for example, the roots of trees alter the structures, causing disintegration).

#### *3.3. Case Study (3): Assessment of the Ground and Structural Stability Based on PS-InSAR*

The PS-InSAR results consisted of the displacements in the LOS and the vertical direction. As the test site was represented by an urban area, the PS-InSAR technique was optimal for the assessment of the building displacements and consequently of the cultural heritage condition. Five intervals were defined for the representation of the results, namely, vertical displacement velocity values higher than +3.50 mm/year were identified as the uplifting process, the values between +1.50 ÷ +3.50 mm/year were attributed to moderate uplifting, the values between +1.50 ÷ −1.50 mm/year were considered to indicate stability, the values from the −3.50 ÷ −1.50 mm/year interval were attributed to moderate subsidence, and the values lower than −3.50 mm/year were associated with subsidence.

It can be observed that subsidence in Figure 19 is concentrated mainly in the south of Alba Iulia city, whereas in Figure 20 it is located in the east part. The explanation resides in the fact that the data acquired from descending orbits is more sensitive to the deformation that occurs in the west, northwest, and southwest directions, and the data from ascending orbits is more sensitive to deformation along the east, northeast, and southeast directions, as mentioned in [65].

The PS-InSAR results revealed that subsidence was encountered mainly in the districts that were developed on the former marshland areas and that are currently affected by high groundwater levels that flood the foundations of the buildings when heavy rains occur. The use of SAR data acquired from different Sentinel-1 orbits (i.e., descending and ascending) enabled the identification of both areas of concern. Figure 21 illustrates a detailed view of the vertical displacements of the buildings located on the eastern side of the fortress and the evolution of Alba Iulia city starting in the 18th century. The displacement map (Figure 21a) contains the PSs generated from the processing of the Sentinel-1A/B datasets corresponding to this particular area of interest and with a vertical displacement velocity value lower than −1.50 mm/year in order to better emphasize the subsidence areas. The urban planning of Alba Iulia city (Figure 21b) depicts different stages of development, namely, the areas elevated in the 18th century, 19th century, between 1900 and 1945, and after 1945.

**Figure 19.** The vertical displacement map obtained from Sentinel-1A data (© background data source: OSM).

**Figure 20.** The vertical displacement map obtained from Sentinel-1B data (© background data source: OSM).

**Figure 21.** The subsidence areas within the district affected by high groundwater levels: (**Left**) the vertical displacement map obtained from Sentinel-1A/B (© background data source: OSM) and (**Right**) urban development evolution of Alba Iulia city after 1700 (© [44]). As previously mentioned, PS-InSAR enables the assessment of the structural preservation state also at the individual building level.

An analysis of the vertical displacements within Alba Iulia Fortress revealed that overall the area is stable. Some of the PSs indicated the presence of subsidence and require further investigation using either a larger Sentinel-1 time series (with a shorter revisiting time for the acquisition of the images and also by adding new data in order to observe whether the trend remains unchanged) or multitemporal SAR imagery acquired by other satellite missions (e.g., TerraSAR-X, COSMOSkyMed) over the same time frame.

Figure 22 illustrates the results of the PS-InSAR processing for an emblematic building of Alba Iulia Fortress, namely Coronation Cathedral.

**Figure 22.** Coronation Cathedral vertical displacement velocity map obtained from Sentinel-1A/B (color symbols as above figures; © background data source: Google 2020).

The values on the map represent the vertical displacement velocities derived from the processing of the two Sentinel-1 datasets. Although the values were mainly within the stability interval (as defined above), the results indicate a slight subsidence trend (on the northern and southern sides, to the west) as well as a slight uplift trend (on the northern side, to the east) for the perimeter buildings. The previous ground examinations revealed that the cathedral presented deformations caused by subsidence.

It could be argued that the results are not completely reliable and their interpretation is difficult to some extent. However, the results are valuable and accurate, as the displacements were measured at the individual building level. The PS-InSAR technique signaled the potential susceptible areas that require further analysis using satellite data or ground investigations (e.g., terrestrial measurements, visual inspections, field surveys). From this perspective, PS-InSAR is excellent for the early detection of the potential threats to the cultural heritage sites. In addition, the results' level of detail is adequate also for archaeological research. From the cultural heritage management perspective, the early detection of building instability (e.g., inclination, displacement, deformation) is extremely important. For example, an inclination of more than 2 degrees suggests an ongoing deterioration process. In addition, horizontal displacements of 1–2 cm signify structural instability or subsidence. These types of processes cannot be observed during a visual inspection, as they are imperceptible to the human eye. Hence, PS-InSAR provides essential information for cultural heritage management, as timely specific intervention actions could be undertaken to safeguard the buildings. It is also important to find the underlying cause of these

processes. For example, if the inclination is caused by vibrations, then the cultural heritage building requires reparation work using filling materials. If the inclination is caused by displacements, then major intervention work must be performed at the foundation of the affected building. In these cases, a detailed analysis regarding the building integrity is also performed locally by construction experts. As the displacements of cultural heritage are critical and require high restoration costs, early identification leads to diminished risks and lower intervention costs.

#### *3.4. Relevance of the Results for Cultural Heritage Management*

The results obtained within the three cases studies strengthen cultural heritage management by supporting the prioritization of intervention actions. In addition, the results contribute to the elaboration of the site intervention documentation and to the assessment of cultural heritage's structural integrity. In Alba Iulia, the cultural heritage management plan is elaborated both at the level of the territorial administrative unit and individualized for each area. Hence, each cultural heritage site has a distinct management plan. The satellite-derived results are pertinent to the revision of the management plans.

Accordingly, the users of these results are the authorities in charge of the conservation and protection of the cultural heritage sites. Alba Iulia City Hall, in collaboration with the National Museum of Unification Alba Iulia, manages the local cultural heritage sites. A new department will be established within City Hall with the aim of enabling an improved management of cultural heritage by using an integrated approach. At national level, the National Heritage Institute and the Ministry of Culture represent potential users of the results.

A diagram showing the overall concept of the research and the relevance and usefulness of the results in the context of cultural heritage management is presented in Figure 23. Summing up, the results of the Alba Iulia study consist of an analysis of the urban expansion that consequently might lead to increased pollution and more demanding archaeological interventions within the areas with changed land use, an investigation regarding the ascending temperature trend that represents a key element of climate change that directly affects the conservation of the cultural heritage, and an evaluation of the ground and structural integrity, which have a clear value for preservation measures. Without exception, the results of the study can be exploited for improved individual management plans for cultural heritage sites.

**Figure 23.** Relevance of the satellite-derived urban development multifaceted analysis in the context of cultural heritage management.

#### **4. Discussion and Conclusions**

The main objective of the present research was to demonstrate and promote the use of satellite remote sensing as a convenient approach for the generation of sustainable results in the process of cultural heritage stewardship, preservation, and management.

The study contributes to a better understanding of the extent and impact of the primary and secondary factors that affect Alba Iulia's cultural heritage sites. Urban expansion is associated with increased residential areas and commercial and industrial development. Consequently, these factors generate an increased use of transport infrastructure, air pollution, temperature rise, and altered environmental or biological factors such as humidity and dust. Additionally, climate change also contributes to local temperature rise and increased hazard risks. Moreover, other human activities such as illegal construction threaten the integrity of cultural heritage. In general, all the related urban development factors are recognized as having a major impact on heritage sites [74]. In this context, sustainable urban development should be strongly linked to the protection of cultural heritage [75], and thus information regarding the current state of conservation and the challenges faced by urban heritage is essential [76]. The urban analysis performed in this study based on multitemporal satellite imagery generated significant results that can be further exploited by the responsible authorities (i.e., Alba Iulia City Hall) to adapt the individual management plans for cultural heritage sites.

Specifically, the study enabled a detailed analysis of the temporal and spatial evolution of the Alba Iulia territorial administrative unit's urban growth. Besides an accurate quantification of the magnitude of urban expansion, the study identified the specific time frame in which the most significant changes occurred. Furthermore, the study provided evidence that upholds the analysis and monitoring of the vicinities of the investigated cultural heritage site, as these areas usually offer essential information that should be taken into consideration by cultural heritage site managers when revising protection and conservation strategies. Therefore, the accurate assessment of the urbanization trend obtained based on the LULC analysis expands the knowledge related to the risk factors that can irreparably destroy cultural heritage. The upward trend in built-up land can contribute to an upward trend in pollution caused either by intensified vehicle traffic or newly developed industrial areas. Over time, pollution represented a key factor in the degradation of the surfaces of historical buildings and monuments. Pollutants emitted into the atmosphere can have a major and sometimes irreversible impact on the objectives of cultural heritage, through phenomena of corrosion caused by chemicals or soiling caused by particles [48]. The early identification and anticipation of this threat contribute to the development of strategies and the adoption of techniques for safeguarding cultural heritage. In the past, the third gate of the fortress was affected by pollution due to various atmospheric agents that produced the physical disintegration of the cultural monument, for which restoration work was required. Thus, the early identification of these threats contributes to risk mitigation and implicitly lowers intervention costs.

The urban heat islands that were delineated inside the study area represent useful information for cultural heritage authorities. In this way, the intervention measures that can be undertaken (e.g., establishment of new vegetation areas, extension of the existing ones) would be more effective when using the information related to the location of the heat islands.

The remote sensing analysis also exposed the potential vulnerable areas within the study area. Currently, PS-InSAR represents the most advantageous tool for the accurate and rapid identification and assessment of ground and structural displacements. PS-InSAR results flag the buildings that might require further on-ground investigations and, in this manner, strengthen the preventive measures within the cultural heritage management process.

Overall, the results obtained in the three case studies contribute to a better characterization of the Alba Iulia territorial administrative unit that encloses the Alba Iulia Fortress cultural heritage site. Remote sensing enables the analysis and monitoring of cultural

heritage sites over impressive time spans. In addition, due to the recent Earth Observation data policies that allow the free use of data, studies have become more complex and the integration of different approaches has enriched the variety of the results and enhanced their join interpretation. The outcomes of the study strengthen archaeological research and integrate a wider spectrum of information that is fundamental for advanced cultural heritage management decisions and actions.

Considering the impact of air pollution on cultural heritage conservation, future work will focus on the quantification and monitoring of air pollution levels based on satellite data. The recently launched Sentinel-5P satellite mission is dedicated to air quality monitoring and enables the identification of trace gases and pollutants [77]. Using an interdisciplinary approach, the study will target the identification of the pollutants and particles with a high impact on cultural heritage, as well as the improvement of conservation and risk mitigation solutions.

Furthermore, future work involves the use of very high-resolution optical remote sensing data for the identification and monitoring of the urban changes at a very detailed scale. Such an analysis would enable the detection of individual buildings that are constructed within the boundaries of the protected cultural heritage site. For example, the early identification of illegal constructions based on satellite imagery can prevent the deterioration of the cultural heritage site by adopting counteractive measures. In addition, in some cases, ground inspections cannot be performed due to various reasons; accordingly, satellite monitoring would represent a feasible source of accurate information. Not only can illegal constructions affect cultural heritage sites, but also the buildings that receive formal approval for construction within the protected area of the archaeological sites. For example, changes to the skyline can have a major impact on the historic landscape, consequently affecting the value of the cultural heritage site. Therefore, the continued monitoring of cultural heritage sites using very high-resolution satellite imagery would provide essential information for management authorities to change the General Urban Plan. For more accurate results, processing very high-resolution satellite data will be performed based on modern classification algorithms such as Support Vector Machine (SVM) or object-based image analysis (OBIA).

Likewise, future work involves the application of PS-InSAR based on very highresolution SAR imagery (e.g., TerraSAR-X, COSMOSkyMed). The purpose of this approach targets multiple aspects, such as the generation of more PSs for an improved examination of cultural heritage buildings, and the cross-validation of the displacement trends obtained from the exploitation of the Sentinel-1 data. In addition, in order to continue and complement the analysis of the multitemporal SAR series, the processing of newly acquired Sentinel-1 data is foreseen. A special emphasis will be placed on the validation of the SAR interferometric processing results based on in situ observations (e.g., geodetic measurements performed in reference test sites, integration of data collected by the Global Navigation Satellite Systems (GNSS) permanent stations).

Lastly, the study demonstrated in all respects the potential of Earth Observation data to provide meaningful information for cultural heritage management, especially for the assessment of the conservation state, the early detection and impact estimation of different threats/factors, and the planning of intervention actions. In addition, multitemporal and multisensor remote sensing analysis indirectly contributes to target 11.4, "Strengthen efforts to protect and safeguard the world's cultural and natural heritage" [78], of the United Nations Sustainable Development Goals (UN SDGs).

**Author Contributions:** Conceptualization, C.M., I.D.N., C.E.M., A.M.L. and A.L.D.; methodology, C.E.M., A.M.L. and A.L.D.; writing—original draft preparation, C.M., I.D.N., C.E.M., A.M.L., A.L.D., G.T.R. and I.C.I.; writing—review and editing, C.M., I.D.N. and A.B.; supervision, A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by a grant from the Romanian Ministry of Research and Innovation, CCCDI—UEFISCDI, project number PN-III-P1-1.2-PCCDI-2017-0413/contract number 50PCCDI/2018, within PNCDI III.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available due to ongoing researches in this field.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Photogrammetric Measurement of Erosion at the Sabbath Point Beothuk Site in Central Newfoundland, Canada**

**James Williamson 1,\* and Ionut Cristi Nicu 2,3,4**


Received: 4 August 2020; Accepted: 11 September 2020; Published: 14 September 2020

**Abstract:** Erosion at archaeological sites in Central Newfoundland, Canada is a major concern, which is compounded by the fact that there has been a dearth of archaeological research in this region. While more than 70 house pits are known, very few excavations have examined whole features in the Exploits River Valley (ERV), and the archaeology of many has not been examined yet. The aim of this study is to examine the rate of erosion at the Sabbath Point house pit, a recently recorded archaeological site, located on the bank of Red Indian Lake (RIL), and to describe a low-cost methodology for analysing site level bank changes. This site is particularly important, as it represents an example of a late Beothuk residential feature about lifeways practiced in this region. The surveys employed here were carried out using image-based modelling. GRASS GIS was used to measure the diachronic difference between bank edges. The Digital Elevation Models (DEMs) were then compared, and the differences were measured using a transect based method. The erosion measurement has shown that Sabbath Point is in danger of being completely eroded. This shows that a salvage excavation program covering the entire feature is necessary within the next few years, as the feature itself will begin to erode.

**Keywords:** erosion; Beothuk; cultural heritage; GRASS; photogrammetry; UAV; Newfoundland; GIS

#### **1. Introduction**

The erosion of cultural heritage sites has been analysed using geospatial technology since the 1990s, although studies have primarily focused on the use of satellite imagery [1] and spatial statistical methods [2]. These studies tend to apply either analytics focused on forecasting erosion [3] or measuring ongoing erosion from prior data [4–6]. These efforts are directed towards different natural disasters like landslides [3,7], gully erosion [8,9], rockfalls [10], debris flow [11,12], flooding [13], coastal erosion [14], and anthropic causes including: infrastructure works [15], urban sprawl [16], intensive agriculture [17], etc. Cultural heritage represents an intersection between history and society [18], and archaeological cultural heritage has proven to be an important nexus of identity and social cohesion as well as a valuable economic resource via tourism [19]. Therefore, assessing the current state of threatened cultural heritage is crucial, given its importance.

The stochastic nature of erosion has also been widely discussed, as has the increasing severity of weather patterns [20]. Studies have shown that archaeological sites within the arctic and sub-arctic are under greater threat from anthropic effects due to global warming. Whether satellite imagery or geophysical prospection are used [21], the results are promising, and the methods are continually being improved upon.

The use of unmanned aerial vehicle (UAV) photogrammetry in quantifying the dynamics, morphology and erosion of landforms has become increasingly common, as it offers high-resolution local imagery of archaeological sites. UAVs have been used to monitor river morphology, bank erosion [22] and to create data for use in river management and flood analysis [23]. The use of UAV photogrammetry in archaeology began properly around 2010, when a proliferation of archaeological papers appeared which analysed the usefulness of UAVs for understanding sites [24] and to document their erosion [25]; however, this technique is still in its infancy, and quality standards have not yet been formally agreed upon. Image-based modelling has been used to record archaeological features and excavations [26], and to record artefacts [27]. Verhoeven [28] began the widespread discussion of the use of aerial image-based matching in 2011, and the technique has been widely tested against traditional aerial photography [29]. The use of both aerial imagery and terrestrial photography for image-based modelling has also been tested in Greece, where a combination of aerial and terrestrial photography was used to prepare 3D models of the features [30]. The combination of terrestrial and aerial photography has become more useful as image-based modelling algorithms have also improved [31]. This allows 3D models to be prepared using different acquisition methods as needed. Early attempts to use photogrammetry for measuring erosion confirmed that non-metric cameras can be used to prepare 3D models to examine changes in landscapes [32].

Recent investigations in the use of photogrammetry for measuring erosion have shown that an external error of 3 cm is likely when compared to other techniques; however, the method has also been shown to be effective in analysing subtle landscape changes, even in slow moving rivers, as the internal error is often much lower [33]. There have been several archaeological studies focusing on the use of photogrammetry for measuring the erosion of cultural heritage and the use of UAV imagery to analyse erosion [25,34–37], but these have focused on map based comparisons and the analysis of photos and orthophotos. These have been conducted at different scales, from studies on footprints [36] to studies of monumental sites [25,37]. This paper introduces a methodology for studying bank erosion using a DTM analysis, focusing on the indigenous archaeology of the Boreal forest, an area which has traditionally been ignored, as archaeological sites are often obscured by vegetation and thus may remain un-recorded. Archaeological techniques, including erosion studies, must be altered to work in these conditions [38]. This study is novel, as it applies a photogrammetric method using DTMs for studying a bank edge at a cultural heritage site. As similar photogrammetric methods have been used in environmental studies [32,39], this technique can be implemented for studying erosion at archaeological sites. The use of consumer grade cameras for measuring erosion has also been tested, with the result that these can be effective for analysis [39].

Canada's cultural heritage, especially the archaeological sites located on the coastal areas in the Arctic, is threatened by thermokarst activity [40] and coastal erosion [40–42]. Archaeological sites from the coastal area of Newfoundland are vulnerable to increased erosion resulting from future sea-level rise and increasingly intense weather caused by global warming [43]. As studies on erosion at cultural heritage sites located on the shores of inland man-made reservoirs are scarce [4–6,44,45], this study is necessary to assess the danger to archaeological features at Sabbath Point, Newfoundland, Canada. The study involved the recording and measurement of the bank at the Sabbath Point site using image based modelling, referencing the models, and then measuring the bank by outlining the edge and also by using an elevation profile method in GRASS GIS [46]. The use of image-based modelling, a photogrammetric technique based around matching images to create 3D models [47], has become a common technique for archaeological recording, but not for measuring erosion.

The current study will focus on the rate of erosion on the bank beside the Sabbath Point house pit, which is gradually eroding into RIL, central Newfoundland, Canada [48]. This feature was part of a larger historic indigenous use of this area by the Beothuk, who have since become extinct [49]. This study was carried out as part of a larger project "Beothuk Settlement in the Exploits River Valley" (ERV). This project focuses on surveying the archaeological landscape of the ERV through UAV imagery and geospatial statistics, to examine the social characteristics of these features.

#### *Regional Settings and Archaeological Background*

Sabbath Point is a historical indigenous residential feature on the bank of RIL, in Central Newfoundland [50] (Figure 1), north of the mouth of the ERV (Figure 2), and is among the later sites inhabited by the Beothuk, who were the indigenous people of this region. This site is believed to have been inhabited in the 18th century, as it has a hexagonal outline, which fits with Marshall's house pit typology [49].

**Figure 1.** (**A**) Geographical location of the study area in Newfoundland context; (**B**) Location of Sabbath Point Site in regional context; (**C**) Detail on the house-pit location on the shore of Red Indian Lake.

**Figure 2.** The location of all the archaeological sites in the ERV.

The Sabbath Point house pit has seen three separate excavations (the areas excavated are visible in Figure 3), which have each focused on the interior of the house pit. These took place in July 2018 [50], September 2018 [51], and June 2019 [52]. As they focused on the interior of the feature, they did not examine the exterior architecture, aside from the September 2018 excavation, which showed that they were partially stone based. The area around the site contains important indigenous cultural heritage (CH) sites and is within a settled landscape. Historical evidence suggests that this area was almost exclusively inhabited by the Beothuk between the 16th and 19th centuries [53]. The relationships between different indigenous groups in Newfoundland have been antagonistic [54] but recent discussions have called this into question [55] and suggest that a more nuanced understanding may be necessary. The wider cultural landscape has also been negatively affected by anthropic changes to the landscape. Erosion measurement has been used here to show how erosion is proceeding at the Sabbath Point house pit, and that erosion on RIL is a danger to other regional features.

The ERV and RIL in central Newfoundland were cultural hot spots throughout human history in the area and have seen 5000 years of human habitation [56]. The site immediately to the north of Sabbath Point, Indian Point, was occupied by the Maritime Archaic Peoples, followed by the Little Passage peoples (the ancestral Beothuk) and then the Beothuk themselves. Within the ERV, barring the Dorset Paleo-Inuit, there are sites associated with every indigenous group to live in Newfoundland [57–59]. The area was culturally central to the Beothuk [53]. The shore of RIL also held the funerary hut of both Nonosabasut and Demasduit, in 1820, although these features have since disappeared, and are likely to have eroded into RIL [49]. The indigenous archaeology of this region is eroding into RIL at an accelerated rate, and many features are likely to have been lost.

**Figure 3.** View of the areas excavated, showing the house pit itself, and the bank edge from 07/09/2018 (excavations took place on the interior part of the house pit and through a single wall which is not in immediate danger of eroding).

The underlying geology of the area is within the Dunnage zone of the Newfoundland Appalachians. Volcanic and sedimentary rocks within Dunnage zone are remnants of continental and intra-oceanic areas and are composed of back-arcs and ophiolites deposited within the Iapetus Ocean during the Cambro-Ordovician period. The Red Indian Line separates rocks of the Buchans and Robert Arm groups of the Notre Dame Subzone from the Victoria Lake Supergroup of the Exploits Subzone. The Red Indian Line is located along the south shore of RIL and extends on a north-northeast direction along Joes Lake and eastern edge of Dawes Pond and continues northeast into Notre Dame Bay [60]. The climate of Newfoundland is typically northern Atlantic, with short summers and long, mild winters. The average temperature in central Newfoundland is about 17 ◦C in the summer and −6 ◦C during the winter. Mean annual precipitation ranges from 700–900 mm/yr and the mean annual snowfall ranges between 275–325 cm [61].

The erosion has been caused by the raising and lowering of the water-level for the Grand Falls-Windsor power plant downstream of this site, and this has impacted the cultural heritage of the region. Since 1980, Beothuk house-pits at different areas within Central Newfoundland have also eroded into the river [58,62,63]. At Red Indian Falls (RIF) 5, 23 km to the East of this site, a house pit has eroded into the Exploits River [63], at RIF 1 the features are disappearing because they are regularly flooded [64], and at Boom Island several hearth features from the ancestral Beothuk are almost permanently underwater [58,65]. Erosion can destroy features, and submergence under a fast-flowing river will likely remove subtle features than can allow surface archaeological analyses. There have been no prior studies of erosion at archaeological sites in the interior of the island. While the Coastal Archaeological Resources Risk Assessment (CARRA) Project [66] focused on coastal archaeological resources, very little thought was given to the state of inland sites. The reason for this is that interior archaeological sites on the island of Newfoundland have largely been ignored in archaeological literature, as they do not fit the current coastal model of life on the island prior to contact [67].

A method for studying the erosion at this site, since its discovery in 2016 [68], was necessary to provide an approximate baseline for how erosion should be expected to continue at other sites [59]. Image based modelling was used to prepare measurements of this feature, as it is a scalable method for measurement which can be applied to different areas easily [69]. These models were processed using Agisoft Metashape, which has recently been shown to have the most accurate photogrammetric algorithm [31], when compared with open source tools, and has been applied to archaeology at several different scales and in varied situations [24,27,70]. Agisoft Metashape is both more user friendly, more accurate, and has a more efficient workflow than its nearest rival, MicMac ENSG [71].

#### **2. Materials and Methods**

The measurement of erosion at Sabbath Point has been carried out using a photogrammetric technique, image-based modelling, and analysed using GRASS, an open source GIS [46]. A comprehensive flow chart methodology is outlined in Figure 4.

**Figure 4.** Methodological workflow used in this study (in blue and green at the top is the raw data, the red shows the photogrammetric processing, and the blue in the middle shows the raw rasters produced in Agisoft; finally, the pink shows the steps involving GRASS GIS).

Image-based modelling was used to measure the bank edge as field measurements would have provided less precise models of the bank edge, and the extra time spent standing on the edge of the bank may have damaged it. The data was gathered on the 07/09/2018, 07/11/2018, 26/06/2019, and 11/10/2019, during each fieldwork visit to the region. The photos were processed in Agisoft Metashape, using points to tie the photos together, then to build a dense point cloud through interpolation, and finally to create a mesh [53]. A final model, which did not cover the edge of the bank, was acquired on the 24/06/2019 (and was used to reference the other models against the grid set up using the total station (Nikon 300 Series Total Station). The meshes were then referenced to the local grid using this model. A Digital Terrain Model (DTM), and an Orthophoto were processed from these models. These are then imported into GRASS GIS, where they are processed to produce a terrain visualisation (aspect, or the direction of slope beginning at North, with 0 degrees). The aspect is then used to digitise the edge of the bank. Then a line along the bank edges was drawn to create a guide for transects to be drawn to measure the edges along their length. The measurements are then compared in a spreadsheet.

#### *2.1. Photograph Acquisition*

Two different photogrammetric acquisition methods were used to prepare the 3D models. The methods of acquisition used were low altitude UAV (using a Phantom 3 Pro) and commercial handheld photography (using a WB36F). These tools, DJI Phantom UAVs and consumer grade cameras, have both been used separately for analysing erosion [5,39]. The photographic acquisition was carried out whenever fieldwork in the area was taking place, to maintain a record of the position of the bank during this time. The handheld camera was used during two of the 3D modelling attempts, it was not possible to prepare 3D models using a UAV and a handheld digital camera was used instead on the 07/11/2018 and 26/06/2018. The aim of the acquisition scheme was to maximise coverage, and to ensure that each part of the bank was covered by at least 9 photos. During the 07/09/2018 acquisition 182 phots were taken, on 07/11/2018 177 photos were taken, on the 26/06/2019 acquisition 432 photos were taken, and on the 11/10/2019 acquisition 665 photos were taken. The variation in the number of photos was due to the size of the area covered by each model, and the coverage of the house-pit. The 24/06/2019 acquisition included 163 photos and was used to reference the models. This did not cover the bank and focused on the interior of the house-pit.

The acquisition flight with the UAV was performed manually, as autonomous flight paths rely upon the UAVs internal GPS, and can often waver, making planned flights unsafe in areas under tree canopy. When using the camera, photos were taken from 1–2 m from the bank edge, separated by around 50 cm, or a single step. The camera-based acquisition method involved taking a photo pointing at the bank, or the house-pit, and another at an oblique angle, to increase the possibility of photo matching. As handheld photogrammetry is not particularly suitable for large horizontal features [64] a focus on the bank was necessary, as it provided a vertical dimension, that made it easier to prepare a model. The photographic acquisition method with the UAV involved taking photos 2–4 m above the feature flown in a rough grid, with photos taken directly facing the feature and oblique angles. The useful geometric data captured in oblique photos is often uneven throughout the photo but allows for several photos to be matched more easily, as they will include geometric data from a wider area.

#### *2.2. Image-Based Modelling Method*

The photographs were processed using Agisoft Metashape [72]. Agisoft has the most accurate image-based modelling algorithm [31] when compared with other photogrammetric programs and is the industry standard in archaeology because of its user friendliness, the fact that data remains on the computer of its user, and is cheaper than the nearest competitor. MicMac ENSG, the nearest open source competitor, has recently been in transition between two different algorithms, and was not stable enough to be used here. Image-based modelling works by comparing photographs of the same object, which share geometric features, to produce a set of points which tie together each photograph [64,72]. These are then added to by interpolating points from the photographs around this, and these are then

used to prepare the whole 3D point cloud. This point cloud is then sampled to a fifth of the dense cloud density for processing to create a mesh. A mesh, or triangles joining every point to form a surface are produced. The mesh was used to create the Digital Elevation Model (DEM), and an Orthophoto Mosaic. The 3D models were then compared against each other after referencing the models against the model prepared on the 24th of June 2019. As the site is largely under thick tree canopy, Real Time Kinematic (RTK) points in certain areas of the site had high levels of error [53]. The 24 June 2018 model was aligned using the corners of the June 2019 excavation, measured using the total station. It had an error of 3 cm from the grid. This model was then used to pick out points in common with each different model, and these were evaluated according to their error within Agisoft. The 07/09/2018 model had an error of 2.8 cm, the 07/11/2018 model had an error of 1.5 cm, the 26/06/2019 model had an error of 1.9 cm, and the 11/10/2019 model had an error of 1.2 cm. The greatest error present was 3.5 cm on a control point in the 07/09/2018 model. Each of the models used for measurement had 6 or more control points, spread throughout the model. As the area around the bank is changing, points from near the edge of the bank, could not be used. The bank is subsiding, as it is being undermined. The subsidence of the bank meant that preparing a linear measurement using the total station with a high level of detail (which may have required 30 measurements), could have caused further damage to the edge of the bank. The DEMs and Orthophotos are then exported from Agisoft for import into GRASS GIS.

#### *2.3. Geographic Information System (GIS) Processing*

The DEMs and Orthophotos were then imported to GRASS GIS which was used to outline the edge of the bank by comparing different visualisations of the 3D models and drawing the line of the bank over these. GRASS GIS was used because it has powerful and varied tools for raster analyses [72]. The raster processing involved testing the usefulness of different visualisations of the model, which included using comparisons of slope and DEMs, hillshade (or relief), and composite images. The comparison of different visualisation has been shown to be an effective tool in the past and has been central to resolving feature edges in other studies [39].

Finally, the Aspect visualisation was used (prepared using the r.slope.aspect module), as this highlighted the edge of the bank more clearly. The edge was digitised as the aspect changed between sloping North West to sloping South East. The edge was recognised by a solid change in colour which corresponded to the difference in direction, and the edge of the bank was highlighted along this. In Figure 5, the four aspect views of the site can be seen, and the fact that the change in the aspect marks where the model is used. The aspects show how a change in the colour gradient shows where the bank should be drawn on. This was overlaid on the DEM, to show a terrain visualisation that could be used to draw on the edges using the GRASS GIS digitisation tool, over the first major change in direction of the bank.

Similar 2D methods are widely used to measure erosion [4,6,14,41]. 2D digitisations are commonly used to analyse 3D models [34,73], and a reliance on satellite imagery and orthophotos can be seen throughout the analysis of erosion on banks. There are limitations to 2D methods, and 3D models offer the possibility of volumetric analysis through voxel analysis in GRASS GIS. This method may allow more detailed results to be rendered, as dynamic processes involved in erosion should show that the bank has expanded. However, this would reduce the ability of researchers to use historic data, as historic maps tend to have either limited or non-existent 3D information. At best this data are often present as contours or relief maps. This is also unnecessary in the case of the present study, as the data required is 2D, to prepare an analysis of the features to show the distance that the bank is retreating.

A line following the broad trend of these models was then drawn along the bank, following the trend of the other models, to provide a common orientation for the measurements. To prepare measurements, transects were drawn at 0.1 m. The distance between the bank edge at 07/09/2018 and the later bank edges was measured along the transects, to show how the bank edge has retreated. The layering of the transect and the measurement lines, when compared with the bank edge can be seen in the diagram below to create a visualisation of how the data are being gathered (Figure 6). This shows how different layers can be used to prepare the measurements. Finally, the GRASS GIS data was imported into the QGIS GRASS Module, where it was used to prepare the graphics used here.

**Figure 5.** Visualisation of the aspect of the bank on (**A**) 07/09/2018, (**B**) 07/11/2018, (**C**) 26/06/2019, and (**D**) 11/10/2019.

**Figure 6.** Components of the measurement method showing (**A**) the line drawn between the end points of the measurable bank and the transect at 25 cm distance; (**B**) different bank edges measured; (**C**) merged layers (background: orthophoto of the site from 07/09/2018).

#### **3. Results**

The changes in the bank can be seen in Figures 5 and 7, with each line representing the bank at the given dates, with the image below showing the bank at different areas. The background shows the orthophoto of the site based on the 07/09/2018 imagery captured using a Phantom 3 Pro. To the immediate south of the bank, the house-pit wall can be seen. This shows how this site is in danger. Much of the architectural remains of the houses and zooarchaeological evidence in surrounding middens are placed outside the house [26,29,50]. The exterior archaeological features are in particular danger.

**Figure 7.** Orthophotos of the bank on (**A**) 07/09/2018, (**B**) 07/11/2018, (**C**) 26/06/2019, and (**D**) 11/10/2019.

#### *3.1. Bank Edge Erosion*

The lines that show the bank as it has progressed throughout this period. As can be seen in Figure 8, the bank edge has retreated at different speeds between 07/09/2018 and 11/10/2019 and that the rate of erosion is increasing. Between the first and last data gathered the bank has moved back by 72 cm at the widest point. The shortest distances are below the level of error. In one area, the distance between 07/11/2018 and 26/062018 shows that the bank has subsided and slide forward. This is also an artefact of the analysis, because a 2D analysis, as is commonly used in erosion studies, does not show volumetric changes. This study shows how the bank is retreating, and thus, needed only to show a 2D measurement. This suggests that there are dynamic processes occurring around the bank, which is being undermined by the annual water-level fluctuation between Summer and August. The effect of the snow on the site is less than the effect of soil falling out during the late summer and early autumn. The difference between these stages can be seen in Figure 8. Other changes have taken place in the bank as well, as it is being undercut by rising water levels. Holes have opened in the bank, which were partially overgrown by moss the 11/10/2019 measurement. These holes can be seen in Figure 9. This shows the changes were more complex than bank regression and shows that the feature itself may be undermined soon.

**Figure 8.** The graph is showing the distance eroded (m) from the 07/09/2018 bank and the future bank measurements.

**Figure 9.** Holes that began appearing in the bank on 26/06/2019.

#### *3.2. Diachronic Analysis and Spatial Averages*

The bank edge has eroded an average of 32 cm in across the measurements, between 07/09/2018 and 11/10/2019. Because the bank edge has retreated an average of 10 cm in two months, 17 cm in 9 months, and 30 cm in 13 months with the greatest difference being 70 cm, it seems as though erosion is accelerating. However, the average difference per month is greatest in the first two (07/09/2018 and 07/11/2018), measurements. When the area between 50 cm and 160 cm from the left is counted, the average change in the bank is greatest monthly change is between the 3rd and 4th models (26/06/2019 and 11/10/2019). Between the 50 cm and 160 cm (from south west to north east) measurements along the transect line, the monthly change between each model is 4 cm between the 07/09/2018 and 07/11/2018 models, 3 cm between the 07/11/2018 and 26/06/2019 models, and 9 cm between the 26/06/2019 and 11/10/2019 modes. This shows how the rate of erosion is spatially varied throughout the transect area. This can also be seen in Figure 5. In the areas in which 2D erosion has slowed down, holes have opened underneath the bank edge, as can be seen in Figure 9, and this suggests that erosion is continuing, and that the process is more complex.

#### **4. Discussion**

#### *4.1. Regional Implications*

The deterioration of the bank at Sabbath Point draws parallels with other Beothuk residential features. The feature is not only facing erosion at the slope beside the water line but is also being undercut by the edge of the bank. The raised walls, which have not been excavated on this side of the site (the excavations are outlined in Figure 3), could show information about the architecture [31], and midden heaps, normally outside the features [56,57], may have also been lost.

A map of the unlabelled archaeological features in the ERV can be seen in Figure 2. These have been left unlabelled to offer some protection to the sites. Important Beothuk features on the bank of RIL including the burial hut of Nonosawbasut, and Demasduit are already likely to have been lost, as RIL water-level oscillates by at least 5 m between June and September (measured from the furthest tree stump visible on the beach) [28]. This relates to a serious issue in studying the geospatial aspects of each Beothuk site, in that the original bank of the river from the Contact period has now been lost. This is also a direct danger that relates the feature to other sites within the ERV, where other house pits have been eroded, as at RIF 5 [52], where erosion has destroyed a house pit, and at RIF 1, where the features are regularly submerged. Other features such as the Boom Island hearths have been submerged for many years [45] and are now only visible from the air. This means that the risk of losing features is not simply limited to RIL, but also occurs throughout the ERV. The erosion throughout the year is variable, and between September and November, most of the erosion occurs. This shows that specific mechanisms taking place in early autumn drive accelerated erosion. The erosion is also variable throughout the measured area. This suggests a more nuanced model of erosion is necessary for archaeological sites in the ERV. Erosion measurement has been carried out at several different sites using UAVs. Erosion at different sites must therefore be analysed at different sites in the future with this consideration, as it could inform site maintenance plans. Two example studies have been chosen for comparison, both of which focus on reservoirs, where much of the damage is anthropic [4–6].

A primary difficulty in drawing comparisons between these discussions is the fact that most studies into erosion take place at large scales, of at least 100 m, while this example has required a grid of 1 m for visualization and measurements at 10 cm. The difference between this study and that at the Kuibyshev dam is two orders of magnitude [6] as the study area in this study is 3 m long, while that of the Kuibyshev dam at the Beganchik site is more than 100 m. Archaeological measurements and geospatial discussions are often scalable [58], but, the implications of erosion are different at separate scales. At sites where erosion is considered at an intersite, or almost regional level, measurements can be used for statistical techniques to forecast erosion due to the stochastic nature of this process [3]. The effect of the scalar difference is that measuring erosion at small sites may produce outliers when compared with large scale measurements.

#### *4.2. Comparison with the Kuibyshev and Stânca-Costes,ti Reservoirs*

The man-made Kuibyshev Reservoir in European Russia draws a useful comparison, as it is also a case of anthropic erosion [6]. The site of Beganchik, represents much of the Palaeolithic evidence for this region, and has been surveyed using UAVs. Erosion at this site has approached a speed of 2–3 m per year since it was surveyed in 2017. Other features within the area of the reservoir are also being submerged. This site shows how anthropic features may gradually eliminate sites but is also disappearing much faster than other regions. The difference in the rate of erosion is also noticeable; however, this may be due to the fact that they are each subject to different mechanical changes, as Beganchik is now on a cliff, while Sabbath Point house pit is now on a lower terrace and has the protection of trees.

A comparison with the Ripiceni—*Holm* (an Eneolithic site in Romania) on the bank of the Stânca-Costes,ti reservoir may also help discuss the dangers man-made dams pose for heritage at the reservoir [4]. The bank edge at Ripiceni—*Holm* has seen a change of 96 cm per year. A direct comparison with erosion at Sabbath Point shows that the Ripiceni—*Holm* site is eroding much faster. As the site is on a plain, the soil in the region is more mobile, which means that the features are in more danger. As the forests of RIL help protect it from erosion events, the vegetation should be protected. The study on the Ripiceni—*Holm* site shows that erosion is causing more damage on the left bank of the river being fed by the reservoir. This has serious implications for sites in the ERV and suggests that future research must focus on forecasting erosion in this area, as the level of erosion may vary throughout the region.

#### *4.3. Vegetation and Bank Deterioration*

Both the Kuibyshev and Stânca-Costes,ti reservoir have archaeological sites on open, unforested land and have faster rates of erosion [4,6] than Sabbath Point. Forests stabilise soil, and tree roots tend to protect sites from erosion [23,59]. The issue with preserving the tree cover features is that many features cannot be properly surveyed without removing the tree cover, and that fast growing alders can destroy or damage archaeological features [47]. Examples of this difficulty in surveying features in the ERV include the sites RIF 5 and RIF 2, where previous surveys lead to recording the existence of features, but were unable to analyse their morphology, leading to changes in their interpretation after brush cutting [47,60]. Brush cutting has since become a regular part of surveying house pits. Sustainable management procedures suggested at Arctic sites have included the fact that vegetation helps protect features from erosion [61] but the rapid growth of alder thickets has been seen as a possible danger for features in Central Newfoundland [47]. This creates a paradoxical situation, and the necessity to manage our brush cutting activities, while preserving vegetation to protect sites. This will however require further research into methods for sustainably managing sites to prepare policies for these sites.

#### **5. Conclusions**

Measuring erosion using a combination of photogrammetry and GRASS GIS has been successful. As the method is based on cost-effective tools, consumer grade cameras, UAVs, and common GIS-based analyses and visualisations, the method could be easily used in future analyses of erosion rates of eroding archaeological sites. The Beothuk house pit at Sabbath Point could help to answer questions about the nature of hexagonal house pits, their interior architecture, and the resources that the Beothuk relied on, but has only been partially excavated until today. Erosion within the study area has seen a maximum of 70 cm, with a measured average of 30 cm in 13 months, however, this does not include the undercutting of the bank seen in Figure 9. The house pit will begin to erode within the next two years, with deterioration in the north wall, however exterior evidence, which has been present at other Beothuk sites. Tree cover in this area has helped to maintain a lower rate of erosion than that seen at either Stânca-Costes,ti or the Kuibyshev reservoirs. Tree cover is also an issue in this area, as it makes

photogrammetry more difficult while also protecting the edge of the bank. Research into erosion using sequential satellite maps and historic maps focusing on statistical methods for erosion analysis is a priority for pursuing a necessary heritage protection scheme.

**Author Contributions:** Conceptualization, J.W.; methodology, J.W.; software, J.W.; validation, J.W.; formal analysis, J.W. and I.C.N.; investigation, J.W.; resources, J.W.; data curation, J.W.; writing—original draft preparation, J.W.; writing—review and editing, J.W. and I.C.N.; visualization, J.W. and I.C.N.; supervision, I.C.N.; project administration, J.W.; funding acquisition, J.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** For J.W.: this research was funded by the Provincial Archaeology Office of Newfoundland and Labrador, the ISER Research Grant Program and the Smallwood Research Grant Program. For I.C.N., part of this work is performed according to the Russian Government Program of Competitive Growth of Kazan Federal University. The APC was funded by the Department of Archaeology, Memorial University of Newfoundland, and the Chair of Northern Indigenous Community Archaeology, Memorial University.

**Acknowledgments:** For J.W.: I would like to acknowledge the help given to me by my supervisors Lisa Rankin and Peter Whitridge in this project, and thanks to Laurie McLean, Don Holly, and Chris Wolff for including me in the excavations at Sabbath Point. I would also like to thank Don Pelley for help during fieldwork in the Exploits River Valley. The constructive comments of two anonymous reviewers are kindly acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **The Pavilions at the Alhambra's Court of the Lions: Graphic Analysis of Muqarnas**

#### **Antonio Gámiz-Gordo 1,\*, Ignacio Ferrer-Pérez-Blanco <sup>2</sup> and Juan Francisco Reinoso-Gordo 3,\***


Received: 16 July 2020; Accepted: 11 August 2020; Published: 13 August 2020

**Abstract:** This research documents and graphically analyzes the pavilions muqarnas at the Court of the Lions in the Alhambra in Granada, a World Heritage Site. In order to cast some light on the understanding and preservation of these 14th century architectural elements, after a brief report of historical data on catastrophes and restorations, a novel methodology for the case study based on three complementary graphic analyses is presented here: First, there is a review of outstanding images ranging from the 17th to the 20th centuries; subsequently, new CAD (computer-aided design) drawings from pavilions muqarnas testing the theoretic principles from their geometric grouping are accomplished for the first time; and finally, a 3D laser scanner is used to understand the precise present-day state from the point cloud obtained. Comparing drawings allows us to assess the muqarnas relevance while proving, for the first time, that the muqarnas of both pavilions have distinct configurations and different amounts of pieces. Besides, this process reveals geometric deformations existing in the original Nasrid muqarnas compositions, identifying small pieces hitherto unknown, plus additional deformations resulting from adjustments after important threats that both pavilions and their muqarnas overcame for centuries, despite their fragile construction.

**Keywords:** muqarnas; Alhambra; graphic analysis; drawings; 3D laser scanner; historical images; cultural heritage; UNESCO; Spain

#### **1. Introduction**

#### *1.1. Background. Brief Historical Data on Catastrophes and Restorations in the Court of the Lions*

The study target in this research is located at the Alhambra Royal House in Granada (southern Spain), a Medieval urban complex -Alhambra, Generalife and Albayzín- included by UNESCO in the World Heritage List taken into account several reasons, particularly the outstanding criterion IV [1]:

"The Alhambra and Generalife bear exceptional testimony to Muslim Spain the 13th and 15th centuries. They form a remarkable example of the palatine residences of medieval Islam, neither destroyed nor changed by the vicissitudes of time, as with the examples in Maghreb. The attributes contained in the inscribed property justify their exceptional position in the Islamic architectural tradition of the Early Middle Ages, and they express the authenticity in a reliable way. Since its conception as a palatine city, its architecture began from a proportional system, following the principles of area compartmentalization, no exteriorization and the typical acclimatized design of the Islamic culture. Together with this, it comes to fruition in a decorative program based upon geometry, epigraphy and vegetable decoration that attain its most characteristic expression in Mocárabe [muqarnas] vaults".

No other place in the Alhambra is as interesting as the Lions Palace and its famous courtyard, a masterpiece of Nasrid and Islam art in Spain, built during the reign of Muhammad V (1354–1359 and 1362–1391) [2]. The groupings of muqarnas pieces used in this environment are undoubtedly one of the most outstanding features of its architectural identity. The Court of the Lions has a rectangular plan measuring 28.50 × 15.70 m and formed by two perpendicular axes. Its elegant perimeter gallery rests on 124 white marble columns with ringed shafts that are grouped in subtle compositional rhythms (Figure 1).

**Figure 1.** Owen Jones; Jules Goury. 1842. "Plan of the Royal Arabian Palace in the ancient fortress of the Alhambra" (red squares mark the Pavilions in the Court of Lions; top of the image facing North) [Plans, elevations, sections, and details of the Alhambra Vol. 1, plate 3].

In the middle of the smaller sides of the patio stand out two almost square-shaped temples which are considered authentic masterpieces of the Nasrid art (Figure 2a,b). Above their columns and abacuses there are brick pillars covered with fine plaster decoration, which support the beams on which eaves and roofs rest. This lintel-based layout freed the arcs from a structural or load-bearing function, and facilitated the inclusion of a beautiful openwork ornamentation that allows the passage of light and the incorporation of sophisticated compositions of muqarnas, the target of this study.

The pavilions's inner floor is made of marble and covered with stunning hemispherical wooden domes, which are one of the finest examples of 14th century carpentry works in Granada. Plaster muqarnas pendentives, which are also a target of this study, were designed as a transition from the square plan to the domes circular plan.

**Figure 2.** Pictures from the Court of the Lions [taken by the authors]: (**a**) View to the western Pavilion (2004); (**b**) view to the eastern Pavilion (2018).

The Nasrid dynasty carefully chose locations for their fortresses and foundations, taking into account the important seismic activity in its kingdom and particularly in the citadel of Granada. The Alhambra underground, called "Formación Alhambra" by current Geology experts is quite suitable to reduce seismic risk. Both site and geological disposition have been crucial factors for the Palace's conservation thus far [3,4]. A historical document [5] reported the flood of the river Darro on 5 March 1600 that affected the base of the Alhambra hill, leading to significant landslides in what is nowadays known as 'Tajo de San Pedro', an area which keeps a part of the monument at geological risk.

The risk of ruin threatened the Court of the Lions many times, as indicated by the historical data outlined below. According to the most prominent Alhambra curator architect, Leopoldo Torres Balbás, the Nasrid decadence had a turning point on 27 June 1431, when a major earthquake collapsed some walls and towers, a few days after the entrance of John II of Castile II in the Vega de Granada starting the Battle of La Higueruela [6].

Repairs have been incessant in the Court of the Lions since both town and landmark passed into Christian hands in 1492 [7]. In documents from 1541 and 1542, mention was made of the replacement of plasterwork and the placement of braces to control the movements of the columns. In 1590, a fire and explosion occurred in a powder mill near the river Darro, on the north slope of the Alhambra, severely affecting the Palace of the Lions, blowing up glasses, doors, and windows. In the Hall of the Mocarabes, next to the western pavilion, its muqarnas vault was destroyed; and other emblematic architectural elements of the monument were also greatly affected, such as the nearby Comares Tower [8].

In 1626, 16 braces were placed in the Court of Lions due to major collapses; and between 1688 and 1691 the ruin was complete. For this reason, the pavilions' roofs were modified between 1691 and 1694, so that they lost their original aspect. The west pavilion nowadays keeps the roof built in the late 17th century, which replaced a previous one with a steeper slope, hipped and with glazed tile trestles [9]. In 1729, some repairs took place due to the arrival of Philip V of Spain, and new repairs occurred in 1744 and 1757. New reports dating from 1784 informed about an imminent ruin, aggravated by a lightning strike occurred on November 5, 1787.

In 1812 and 1820 some light works on the roofs were performed, and also at the beginning of the 19th century a garden was created in the Court and the irrigation water induced earth movements, so it was removed around 1844–1846 [10]. As described later, in 1858, an unfortunate restoration of the eastern pavilion took place; in 1889 it was rebuilt; and in 1910 the subsoil was repaired [11]. Towards 1934, Torres Balbás restored it again in the form it still has nowadays. In addition, in 1966, the architect Francisco Prieto-Moreno consolidated and restored the western pavilion.

Finally it should be remembered that along the history, there have been many natural hazards that have affected the archaeological sites or monuments, for example, of the Seven Wonders of the Ancient World, only one of them—the Great Pyramid of Giza—has survived until our days. Regarding the other six, three have been damaged by earthquakes. The Colossus of Rhodes collapsed around 227 BC, the same occurred to the Pharos Lighthouse in Alexandria in the 14th century AD. Finally, the Mausoleum in Halikarnassos, destroyed by floods and earthquakes and rebuilt several times, eventually disappeared in the 15th century [12,13].

It is surprising how the Court of the Lions has been preserved until our days with such a fragile structure. As indicated above, throughout the history of the Alhambra, movements caused by earthquakes, explosions, or water leaks in the subsoil have been aggravated by neglect or lack of maintenance due to shortage of economic resources, by the effect of weather elements, or as a consequence of inadequate repairs. All these factors have caused on many occasions the detachment of plaster ornaments and the collapse of the columns.

#### *1.2. General Objectives and Methodology*

The pavilions at the Court of the Lions and their fragile muqarnas have barely been studied and, until now, no digital surveys have been made to reflect with precision their formal complexity. For this reason, the present research aims to document and graphically analyze the muqarnas, to facilitate their study and improve their future conservation.

In order to tackle the last objective, we propose a novel methodology for the case study based on an architectural graphic analysis including three types of drawings that complement each other. After reviewing the historical data, some important images from the pavilions extending from the 17th to 20th century have been compared in order to document and illustrate their transformations, as an initial reference for the CAD (computer-aided design) drawings. Subsequently, a digital model is produced to allow the comparison between the reality and the theoretical geometric principles that rule the assembly of muqarnas. Finally, a 3D laser scan provided new images from the point cloud, thus allowing a more precise understanding of their current state. A later comparison between drawings unveils new contributions about both the configuration and current state of the muqarnas, remarkable 14th-century architectural pieces.

#### **2. Materials and Methods. Graphic Analysis of the Pavilions Muqarnas**

#### *2.1. The Pavilions and their Muqarnas Historical Images*

In the first known views of the Court of the Lions, drawn and engraved by Louis Meunier around 1668 [14], the pavilions show a rudimentary but expressive aspect. Their roofs were quite steep, although the eaves representation is strange and the columns grouping patterns is inaccurate, casting doubt on the reliability of the details. Nevertheless, the pavilions seem completely symmetrical (Figure 3).

**Figure 3.** Louis Meunier, h. 1668. "El patio de los Leones in Grenada" [Spanish National Library, invent/19565].

This apparent symmetry is supported by a valuable drawing by the outstanding architect Juan de Villanueva toward 1766–1767, published in "Las Antiguedades Árabes de España" (1787) and edited by the Royal Academy of Fine Arts of San Fernando in Madrid [15]. There is a longitudinal section towards the Hall of the Abencerrajes (Figure 4) where the shading interior of the pavilions are drawn for the first time, but the roofs no longer had the hipped layout depicted by Meunier. The representation of the columns pattern is accurate and the arches ornamentation is schematically drawn, as well as the low-detailed pendentives of interior muqarnas, barely showing their groupings, but with no differences between both pavilions.

**Figure 4.** Juan de Villanueva, 1766–1767. Longitudinal section of the Court of the Lions. [Royal Academy of Fine Arts of San Fernando, Madrid, MA-0540].

Many other draftsmen accomplished outstanding views from the Court of the Lions in the late 18th century and the first half of the 19th century, similarly to what occurred with drawings from other Arab landmarks such as the Mosque of Cordoba [16]. Henry Swinburne published a drawing of the paved courtyard and the east pavilion in his book "Travels through Spain in the years 1775 and 1776 ... " (1779). The accurate drawings by William Gell around 1801, hosted in the British Museum, show the aspect of the courtyard before the garden was created.

It is also worth noting the posthumous work of the Irish architect James Cavanah Murphy, "The Arabian Antiquities of Spain" (c. 1813–1815); and the detailed images by Alexandre Laborde included in his book "Voyage Pittoresque et Historique de l'Espagne" (1808–1820). It should be mentioned the beautiful views by John Frederick Lewis included in "Sketches and drawings of the Alhambra" (1835); de Joseph-Philibert Girault de Prangey in "Souvenirs de Grenade et de l'Alhambra" (1836); or by Nicolas-Marie-Joseph Chapuy in "Le Moyen Age Monumental et Archeologique" (nº 251, h. 1841–1846), separating from other authors that drew simplified muqarnas and did not understand their formal genesis.

The architects Owen Jones and Jules Goury deserve to be paid special attention, as they performed a scientific approach analyzing the ornaments of the Alhambra and researched the geometric construction laws of the muqarnas in order to draw them rigorously. Their splendid longitudinal section from the Court of the Lions pointing to the Hall of the Two Sisters (in the opposite direction to the section by Villanueva from 1766–1767) was published in 1838 (Figure 5), and included in their monumental book "Plans, Elevations, Sections and Details of the Alhambra" (1842).

Still nowadays no CAD plans have exceeded the exquisite level of detail of Owen Jones and Jules Goury's drawings. The pavilions were symmetrically drawn, with hipped roofs different from those drawn by Villanueva; the muqarnas are accurate and very detailed. In addition, Jones and Goury yielded many other drawings of muqarnas on cornices, arches, pendentives or domes, proving a deep knowledge and ability to represent them [17].

**Figure 5.** Jones and Goury (dib.); Kennion (grab.). Longitudinal section of the Court of the Lions, published by Owen Jones, 1838 [private collection].

In 1859 and 1862, Nicomedes de Mendívil y Quadra produced splendid drawings from a pavilion, both preserved at the Royal Academy of Fine Arts of San Fernando in Madrid. They were published as isolated print plates in "Calcografía Nacional" belonging to the outstanding "Los Monumentos Arquitectónicos de España" (1856–1882).

The drawings fill two plates, the first plate representing the plan, details from the muqarnas and the half dome, and the other one representing the section [15]; the pavilion depiction is very reliable and was carefully shaded (Figure 6). The author does not indicate whether if it is the east or the west pavilion the one that is drew on the plate. The plan only contains two out of four arcades, probably because he assumed they were symmetric, although this is not true as we will explain later. The plate published by Nicomedes de Mendívil y Quadra neither finished the other two arcades.

**Figure 6.** Nicomedes de Mendívil y Quadra. Pavilion at the Court of Lions: (**a**) Interior details, 1859; (**b**) section, 1862 [Royal Academy of Fine Arts of San Fernando, Madrid. MA-0210; MA-0211].

Many photographs were taken from the Court of the Lions shortly after the appearance of this technique, all of great documentary interest [18]. A recent exhibition has brought together many of those photographs by various authors [19]: Pablo Marés (h. 1852), Felix Alexander Oppenheim (1852), Alphonse de Launay (h. 1854), Charles Clifford (1854), John Gregory Crace (1855), Jules Palanpin-Dufresne (1856), Eduardo García Guerra (h. 1856), Joaquín Pedrosa (1857), Alexis Gaudín (1858), Jakes August Lorent (1858), Gustave de Beaucorps (1858), and others.

These photographs show the braces placed to prevent the columns from collapse, a threat aggravated by the irrigation of the garden created at the beginning of the 19th century and eliminated around 1844–1846 because the water would cause serious problems in the foundations [20]. The disturbing state of conservation in the Court of the Lions led to a delirious proposal by architect Salvador Amador in 1846, which proposed the complete demolition and reconstruction of the Arab palace, which fortunately was never undertaken [21]. An expressive photograph by Gustave de Beaucorps taken in 1858 shows the east pavilion fully propped up to reinforce its foundation and fix its tilt (Figure 7a).

The aforementioned photographs also show the aspect of the east pavilion prior to the radical intervention of Rafael Contreras, Alhambra curator from 1858 to December 1859 [22], who removed the 17th century roof and the upper plaster walls in order to build a new dome covered with glazed ceramic (Figure 7b). The main problem of this intervention, and others accomplished by Contreras, is that it tried to beautify the monument without a scientific basis and without any rigor, so that original elements were replaced by unrealistic ones, something that fully clashes with the current principles of restoration in architectural heritage.

The new external dome built by Contreras caused rainwater drainage problems that damaged the interiors. After renovating the structure, Torres Balbás commissioned a new tile which was quite controversial due to the marked slope roof [23], removing the upper plaster walls by Contreras that Torres Balbás did not rebuilt. From then until now, two different roof solutions cover the eastern and western pavilions.

**Figure 7.** Eastern pavilion: (**a**) Gustave de Beaucorps, 1858 [Legacy from Ortiz Echagüe, University of Navarra]; (**b**) Purger & Co., c. 1902, after the restoration of Rafael Contreras in 1859 [private collection].

#### *2.2. New CAD Drawings Following Theoretical Principles on Muqarnas Grouping*

Muqarnas are small prismatic pieces grouped on complex patterns or sequences shaping cornices, capitals, arches, pendentives, or vaults which were a symbol of the Alhambra Nasrid architectural identity. It is necessary to carefully draw the plans of muqarnas to understand their complex geometry, their proportions and compositions. The aforementioned plan does not require a scale number to obtain the valuable information referred to: Pieces identification and their grouping.

This research has accomplished the first CAD drawings of plans of the pavilions of the Lions and all their muqarnas pieces using the software ArchiCAD. To achieve this goal, essential references concerning the rules about muqarnas layout have been used: Diego López de Arenas' treatise published in Seville in 1633 [24], Fray Andrés de San Miguel's manuscript from the first half of the 17th century, published in 1969 [25], and a recent reinterpretation of this text by Enrique Nuere [26,27]. References to Eastern muqarnas geometry can be seen in Dold-Samplonius (2015) [28] and Garofalo (2010) [29]. Both papers refer to a treatise by the mathematician Al-Kashi from 1425, which is a fundamental reference for studying muqarnas in the Eastern world.

A review of Owen Jones and Jules Goury's muqarnas drawings has been undertaken, as they propose an elementary grammar on their grouping [17] starting from three flat Figures A, B and C, two triangles and a rectangle, that make up the base for seven key muqarnas pieces, A1, A2, A3, B, C1, C2 and C3, which allow various muqarnas compositions (Figure 8a).

The pavilion plan drawing by Nicomedes de Mendivil includes muqarnas elementary pieces, but the muqarnas grouping does not match the real ones at the pavilion. Mendivil draws six pieces identified from A to F, but he did not sort them according their flat figure shape, nor did he indicate the linkage between the sides from different pieces (Figure 8b).

It should be noted that other authors have proposed a different number of basic pieces or consider that some of them are divided in others [30]. In any case, the three flat figures proposed by Jones and Goury do not resolve the great diversity of muqarnas groupings existing in the Alhambra, where other auxiliary geometries, such as small triangles or squares, frequently appear to enrich many compositions.

**Figure 8.** Key muqarnas pieces: (**a**) Jones and Goury (dib.) 1834 [Plans, Elevations, Sections, and Details of the Alhambra, plate 10, 1842]; (**b**) Nicomedes de Mendívil y Quadra (dib.) 1862; Esteban Buxó (eng., 1864) [Los monumentos arquitectónicos de España, no. 27, 1865].

From the aforementioned precedents and a careful observation of reality, this research identified the existing pieces in both pavilions to consequently draw their plans and basic volumes. Different colors have been assigned to each type of piece to facilitate its location in the subsequent drawings of groupings (Figure 9).

**Figure 9.** Pieces library from the pavilions at the Court of the Lions [authors' own drawing].

During the graphic analysis process, it was detected that the arcades' muqarnas are smaller than those of the pendentives and friezes placed in a higher position and at a greater distance from the observer. It has also been found that in higher ceilings, such as that of the Hall of the Kings, the pieces sizes are even larger.

Finally, it has been observed that the western pavilion displays 2258 pieces, with 1032 located in its pendentives and 1226 in its arcades; on the other hand, the eastern Pavilion has 2222 pieces, with 948 in its pendentives and 1274 in its arcades (Figure 10).

**Figure 10.** Muqarnas plans: (**a**) Western pavilion; (**b**) eastern pavilion [author's own drawing].

#### *2.3. Images Rendered from 3D Scanners*

Thanks to the valuable collaboration from Patronato de la Alhambra y Generalife, the scanning of the muqarnas groupings from both pavilions could accurately assess their current state. The point cloud was captured using a BLK360 laser scanner with a precision of 6 mm when the object is within a distance of 10 m, and 8 mm for 20 m. Three point densities are allowed by the BLK360: Low, medium and high that capture points separated 20 mm, 10 mm and 5 mm from each other when the object is 10 m away. This research worked with medium density and obtained data from the 360◦ environment. All scans were registered in the same reference system (Figure 11). To facilitate the 3D modelling, HDR pictures were taken setting the scanner to HDR mode.

Recap Pro software from Autodesk was employed to obtain the registration in both automatic and manual forms; the manual method would be used only if the automatic option failed. Three variables define the quality register for each scan: (a) Common volume between the current scan (not registered yet) and the scan/s previously registered, (b) the balance between the number of points scanned through the 3 orthogonal spatial directions, the optimal scene is when the points equally expand the three spatial directions, and (c) percentage of points with a fitting error smaller than 6 mm.

Data analysis has covered various phases. Volumes of interest have been selected using AutoDesk Recap. Plans, elevations and sections have been obtained from the point cloud using Cloudcompare software. Then, key points on flat projections have been identified and line drawings have been made on the muqarnas assembly contours and the internal main axes. A 3d model could have been obtained by photogrammetry from a series of photos with targets, which coordinates were known, but the accuracy would be equivalent, as shown in [31].

(**a**) (**b**)

**Figure 11.** Views from the 3D laser scanner [authors' own drawing]: (**a**) Western pavilion; (**b**) eastern pavilion.

#### **3. Results. Comparing Drawings**

#### *3.1. Comparing Di*ff*erent Pendentives*

After drawing the muqarnas geometry, it was observed that their composition in the pavilions pendentives is different. It is easier to identify and visualize the number of pieces of each type and their locations when using different colors for each type.

For example, the naked eye may detect the different locations a trapezoidal muqarnas (marked in brown color) occupies on each pavilion (Figure 12). This muqarnas was identified in Jones and Goury's drawings and frequently occur in other Alhambra rooms.

**Figure 12.** Pendentives composed of muqarnas [authors' own drawing]: (**a**) Western pavilion; (**b**) eastern pavilion.

#### *3.2. Comparing Di*ff*erent Arches*

Arcade muqarnas were made up by partial groupings separated by strips that layout the whole set. These strips are similar to those ones called "medinas" in the text by López de Arenas, but they were not drawn. In each arcade, it is a straightforward observable fact to identify both inner and outer muqarnas rows.

On the other hand, it should be considered that each pavilion includes three arches between columns, one central and two sides identical to each other. This local symmetry in each arcade is one of the few features common to all of them.

The 1226 muqarnas in the western pavilion arcades are located as follow: Northern 309, southern 309, western 307, eastern 301; while in the arcades in eastern pavilion there are 1274: Norhern 317, southern 317, western 301, eastern 339. Thus, considering both pavilions, there are six different types of arches (Figure 13).

**Figure 13.** Different arches of muqarnas [own drawing]: West (western pavilion); east (western pavilion); north and south (western pavilion); west (eastern pavilion); east (eastern pavilion); north and south (eastern pavilion).

#### *3.3. Comparing Di*ff*erent Floor Plans*

Differences between muqarnas floor plans from both pavilions are not so easily perceived due to their geometric complexity. Neither the models derived from the point cloud unveil the differences. For this reason, a layout displaying a color for each unique muqarnas group is provided in order to indicate the differences between pavilions (Figure 14).

Although both have different pendentives and arches, the plan shows an axis of symmetry that matches with the longitudinal axis of the courtyard. It should be noted that there are some groupings duplicated in both pavilions. According to the plan provided in this research (Figure 14) the inner central arches (red) and the northern and southern ones (green) are the same in both pavilions.

Groupings in the western pavilion are simpler than the eastern ones, the differences are in small muqarnas groups with similar contour condition. In this case, every side arch in all arcades is identical (yellow). Central arches in western arcades (red) and eastern (dark blue) are different; although they are very similar to those of the eastern pavilion and just differ in two groups composed of four muqarnas replaced by a flat piece.

The eastern pavilion is more diverse with three types of side arches (yellow, orange and light green in Figure 14). The green arches are only in one arcade and the side arches are duplicated in the northern and southern arcades.

**Figure 14.** Different groups of muqarnas [own drawing]: (**a**) Western pavilion; (**b**) eastern pavilion.

#### *3.4. Comparing Pendentives Elevations*

Comparative analysis between different pendentive elevations is of great interest to assess their accuracy level and their similarity to real works (Figure 15). Jones and Goury drawn identical pavilions more accurately than the previous drawing by Juan de Villanueva, as they drawn the muqarnas rows and a "shell" piece in the center, which is supported by another piece which was not featured, so that the shell piece was enlarged. Nevertheless, the lower level from muqarnas profile is represented with skillful precision, providing high thoroughness and amazing level of detail, particularly taking into account the drawing scale and the fact that it encompasses the whole Court of the Lions.

The comparison between Nicomedes de Mendivil's drawing and the CAD production shows the high accuracy achieved by Nicomedes, and proves that the pavilion produced by Nicomedes was the western one, a fact hitherto unknown, as the artist never indicated it. However, one outstanding finding, observed when comparing the 3D scans, is the significant geometric deformation of elements analyzed in the following section.

**Figure 15.** Inner muqarnas pendentives at the pavilion according to Jones and Goury, Mendivil, and authors' own production: (**a**) Western pavilion; (**b**) eastern pavilion.

#### **4. Discussion. Geometric Deformations**

#### *4.1. Original Layout Deformations*

In a recent paper on the Hall of the Bark in the Alhambra, the difficult geometric problem posed when designing a muqarnas pendentive was explained, as it has to be adapted to the circular shape from an upper dome starting from a square shape [32]. Nasrid artisans did not follow in this case the theoretical process described by Jones and Goury or later authors, but instead deformed and adjusted the muqarnas close to the circle, creating with great mastery and intuitively pieces shaped ad-hoc, barely noticeable to the naked eye.

In order to assess this, a new plan drawing has been provided attempting to follow the geometric rules. It proven impossible to configure the upper levels without using special pieces, because in some areas the pieces overlap each other and other areas become empty. This made necessary the creation of small deformed pieces (darker color in Figure 16). As in the pendentives of the Hall of the Bark, the deformed muqarnas pieces were composed symmetrically in both pavilions.

**Figure 16.** Deformed muqarnas in western Pavilion pendentive [own drawing]: (**a**) Plan; (**b**) axonometry.

#### *4.2. Deformations by Construction*

The point cloud from 3D laser scanner allows derive the pendentives contours from both pavilions, so that we can draw plans and sections in order to take accurate measures (Figure 17). Then we find out dimensional differences or deformations the ideal design (Figure 10) and the current conserved construction.

Analyzing both plans dimensional schemes, it was observed that they are not perfect squares, as expected, but instead they are stretched towards the courtyard center around 10–13 cm; although their angles are close to 90 degrees. Nevertheless, their diameters, i.e., the inner dome basis, have important deformations, close to 13 cm in the eastern pavilion, and 5 cm in the western one. Important angular deformations were also observed in sections that are obviously visible and close to 4 degrees in some of the cases (Figure 18).

**Figure 17.** Plan deformations [authors' own drawing]: (**a**) Western pavilion; (**b**) eastern pavilion.

**Figure 18.** Sections deformations [author's own drawing]: (**a**) Western pavilion; (**b**) eastern pavilion.

Some of these deformations were already described by architect Francisco Prieto-Moreno after the consolidation and restoration of the west pavilion in 1966 [33]. This revealed that the interior was deformed and turned towards the courtyard, probably because the nave pushed the pavilion and it was not possible to restore it to its initial location. For this reason, it was decided to advance the columns bases until achieving their verticality and, subsequently, to eliminate the braces.

A similar operation had to be carried out in previous consolidation works over the centuries, which in many cases were provisionally performed with braces. As plaster is a fragile material, there would be fissures, fractures, and even loss of parts. Thus, the current deformed state is a consequence of successive repairs in which fissures were filled or even new parts were inserted, consolidating their distortions. In this way, a material that could seem ephemeral has facilitated easy repair and conservation over time [34]. These types of repairs allowed other Nasrid buildings in Granada to survive, such is the case of the Corral del Carbón, currently showing important deformations [35].

#### **5. Conclusions**

Drawing is an essential tool for the analysis and sustainable conservation of architectural heritage. Without a deep knowledge about the reality to be preserved or restored, it is not possible to perform suitable actions for its protection. Throughout history, documentation techniques of the architectural heritage have been used with different accuracy levels depending on the needs, the means available, and the characteristics of the portrayed reality. A monument graphically documented in a suitable manner increases its probability to survive over time. In the worst scenario, when, unfortunately, the architectural heritage has been destroyed, its memory can be maintained if there are images that facilitate its virtual recreation or material reconstruction.

The aim of this research was to graphically document and analyze the muqarnas of the pavilions at the Court of the Lions taking into account their historical background. To achieve this goal, after providing a brief account on threats, catastrophes, and restorations over the centuries, an innovative methodology was carried out based on three complementary graphic analyses.

Firstly, this research reviews the main known historical images from the 17th to 20th centuries that clearly illustrate all the transformations undergone by the monument. In this sense, the drawings made with manual techniques in the 19th century are particularly remarkable. Jones and Goury were the first authors to accurately understand, identify, and draw pieces of the Alhambra muqarnas, although they only represented the ones at the western pavilion and assumed that both pavilions were symmetrical. The drawing by Nicomedes de Mendivil is of great quality, although he did not indicate whether the referred pavilion was the western or the eastern one, perhaps because he assumed that both were identical and locally symmetrical, and therefore he only drew two arcades not taking into account that the others are different. Nevertheless, this study concludes that the pavilion drawn is the western one.

This research produced the 3D CAD models or CAD plans for the pavilions muqarnas for the first time, illustrating with precision the formal complexity of muqarnas groupings, considering that valuable drawings from 19th century only included the western pavilion. After analyzing the pavilions in detail, and taking into account the theoretical geometric principles from muqarnas groupings, it has been possible to identify for the first time all their pieces: The western pavilion displays 2258 pieces (1032 on pendentives and 1226 on arcades), while the eastern pavilion has 2222 (948 on pendentives and 1274 on arcades). It has been observed that the pieces located at a higher height have a larger size than those located at a shorter distance from the observer. Schemes provided explain the subtle differences and symmetries between the various pendentives and arches in both pavilions, otherwise too complex to be observed, by means of photographs or 3D laser scannings.

Significant deformations, hitherto unknown, have been detected based on the laser scanner 3D point cloud results. It has been observed that the original design involved deformed pieces that Nasrid artists from the 14th century virtuously built to solve a difficult geometric problem: Adapting the pendentives upper rows to the dome circular base.

Besides, other significant deformations had been reflected in plan and sections schemes, resulting from columns tilts probably due to earthquakes, explosions, lightning, differential seatings or water leaks in the subsoil, aggravated by poor maintenance or by inadequate restorations. It is utterly remarkable how some pavilions supported on such slender columns and muqarnas built on fragile materials have survived over the centuries.

Finally, the graphic analysis provided by this research will hopefully contribute to a sustainable conservation approach towards both pavilions and their muqarnas. Similarly, the article is expected to contribute to the consolidation of the muqarnas significance as outstanding identity symbols of this beautiful Medieval palace in the Alhambra of Granada (Figure 19).

**Figure 19.** Photographs showing the current state in the Court of the Lions [authors' own production]: (**a**) Western pavilion view from interior of eastern pavilion; (**b)** eastern pavilion view from interior of western pavilion.

**Author Contributions:** I.F.-P.-B., A.G.-G., and J.F.R.-G. took part in the entire researching process. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The present research has been made possible thanks to the collaboration of the specialized group of research and dissemination of the Council Patronato de la Alhambra y Generalife as well as the dedicated personnel who surveil the monument day after day. The authors thank the support from the research group Expregrafica. Lugar Arquitectura y Dibujo (PAIDI-HUM-976), Programa de Doctorado en Arquitectura, Instituto Universitario de Arquitectura y Ciencias de la Construcción from the University of Seville; The Survey Modeling Laboratory (SMLAB) from the University of Granada; and the research group Ingeniería Cartográfica (PAIDI-TEP-164) from the University of Jaen. We thank Eugenia Fernández-Beltrán for her work in translating this paper, as well as her dedication.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Coastal Erosion A**ff**ecting Cultural Heritage in Svalbard. A Case Study in Hiorthhamn (Adventfjorden)—An Abandoned Mining Settlement**

**Ionut Cristi Nicu 1,\*, Knut Stalsberg 2, Lena Rubensdotter 3,4, Vibeke Vandrup Martens <sup>5</sup> and Anne-Cathrine Flyen <sup>5</sup>**


Received: 24 February 2020; Accepted: 13 March 2020; Published: 16 March 2020

**Abstract:** Hiorthhamn is an abandoned Norwegian coal mining settlement with a loading dock and a lot of industrial infrastructure left in the coastal zone. In this study, changes in the position of 1.3 km of the Hiorthhamn shoreline, which affect cultural heritage, is described for a time-period spanning 92 years (1927–2019). The shoreline positions were established based on a map (1927), orthophotos (2009) and a topographic survey with differential Global Positioning System (GPS) (summer 2019). Detailed geomorphological and surface sediment mapping was conducted to form a framework for understanding shoreline-landscape interaction. The shoreline was divided into three sectors to calculate the erosion/stability/accretion rates by using the DSAS (Digital Shoreline Analysis System) extension of ArcGIS. The DSAS analysis showed very high erosion in Sector 1, while Sectors 2 and 3 showed moderate accretion and moderate erosion, respectively. Sector 1 is geologically composed of easily erodible sorted beach sediments and protected remains from the mining industry such as wrecks of heavy machines, loading carts, wagons and rusty tracks that are directly exposed to coastal erosion. The all-sector average shoreline erosion rate (EPR parameter) for the 92 years period was −0.21 m/year. The high shoreline erosion rates in Sector 1, together with the high potential damage to cultural heritage, supports the urgent need of continued coastal monitoring and sustainable management of cultural heritage in Hiorthhamn.

**Keywords:** coastal erosion; shoreline; monitoring; geomorphological mapping; cultural heritage; Svalbard; DSAS; high Arctic

#### **1. Introduction**

Coastal areas of seas and oceans are among the most dynamic landforms and are permanently under threat of change from natural processes and anthropogenic pressure. A significant proportion of the world's population, at present and throughout history, lives in coastal areas (between 15% and 40%), which are under the direct effects of climatic change-related [1,2] processes such as sea-level rise [3] and changes in the intensity and frequency of storms [4]. Over the last decades, coastal areas have been subjected to increased rates of erosion [5], and this will most likely cause significant economic

losses in the future. Arctic coastal areas can experience erosion rates similar to, or higher than, those in temperate regions due to the added influence of thawing permafrost and extreme temperatures.

Coastal areas located at high latitudes are even more affected by the changes in the environment (e.g., air temperatures, major river discharges and open water season length have increased, and storm tracks and intensities are changing) [6]. The Arctic coastal zone is defined as the region both seaward and landward of the coastlines of the Arctic shelf areas, including all archipelagos and islands [7]. Arctic coastal areas can experience erosion rates similar to or higher than those in temperate regions due to the added influence of thawing permafrost and extreme temperatures, even though the erosional processes are usually still limited to a few months per year [6]. Arctic areas have become important hot spots for studying the effects of a changing climate [8], which is felt earlier there than elsewhere on Earth [9]. The present ongoing research focuses on modelling the velocity of glaciers and ice caps [10], land cover and ice-wedge polygon mapping [11,12], the surface morphology of fans [13], retrogressive thaw slumps triggering [14,15], coastal erosion [16], human impact [17], etc. One of the most exposed Arctic areas is Svalbard, which is experiencing amplified climate change when compared to the global average [18,19]. Svalbard's coastal area is under high pressure from natural [20] and in some area's anthropogenic changes [21,22].

Studying the physical dynamics of coastal areas is a challenging task with a huge potential to be applied in future coastal management plans and sustainable development of coastal areas. With the technological development, tracking coastal erosion and the associated coastal hazards became easier; new tools are being developed and implemented, along with their integration in GIS (Geographic Information System). Various GIS-based tools have been developed over the years in order to realise different analyses of coastal areas; they have been used in different studies, as follows: XBeach modelling [23], Ntool [24], CERA (Coastal Erosion Risk Assessment) [25], CHW (Coastal Hazard Wheel) [26,27], CESM (Coastal Erosion Susceptibility Model) [28], CVI (Coastal Vulnerability Index) [29], DSAS (Digital Shoreline Analysis System) [30–32].

Aside from the economic activities and the local population located on the coast that is under threat from climatic changes, another significant asset is often neglected when studies are made in the coastal areas—cultural heritage. Globally speaking, coastal areas have been used since Prehistory for human settlement, due to their abundance in natural resources needed for survival and development [33]. Coastal cultural heritage represents an important part of cultural resources in the coastal areas [34]; regardless of the location of cultural heritage, on the coasts of the inland shorelines of large man-made reservoirs [35,36], open seas and ocean coastal areas [37–39], there is a high risk of erosion [40].

The world's cultural heritage is under direct threat from climate change [41]. Cultural heritage located in Arctic areas is even in greater danger [42,43], both from natural causes (such as erosion or burial) [44] and the increasing pressure of Arctic tourism [22,45], which has grown especially in Svalbard [46]. Significant efforts are being made to minimise the impact of visitors on the environment and to move the tourist industry towards a sustainable approach [47]. Over the years, the coastal areas of Svalbard have been subject to intense research (Table 1, Figure 1).


**Table 1.** List of studies dealing with coastal erosion in Svalbard (numbers in column 1 correspond to the red numbers in Figure 1.


#### **Table 1.** *Cont.*

\* DSAS—Digital Shoreline Analysis System.

Studies that approach the coastal processes affecting cultural heritage sites are of high significance; as cultural heritage in the Arctic is in great danger [61–63]. The present study aims to assess and quantify the main changes along an approximately 1.3 km stretch of shoreline located in Hiorthhamn—on the eastern shore of Adventfjorden, Svalbard, where cultural heritage is under threat from coastal erosion. The DSAS tool [64] is implemented to evaluate the erosion and accretion rate of the coastline over time. The quantification of change is complemented with a detailed geological and geomorphological mapping of the affected and adjacent area that will bring new insights in understanding how different substrata and geomorphological processes interact and are affected by coastal processes in a changing Arctic climate. Along with the assessment of coastal erosion, field-based geological mapping was employed, in order to understand the associated geological processes and to map them. The area under interest was chosen within a research project financed by the Research Council of Norway—CULTCOAST [65].

**Figure 1.** Distribution of studies dealing with coastal erosion in Svalbard (red numbers correspond to numbers from Table 1, column 1) (base map from Norwegian Polar Institute) [66].

#### *The CULTCOAST Project*

CULTural Heritage Sites in COASTal Areas. Monitor, Manage and Preserve Sites and Landscapes under Climate Change and Development Pressure—CULTCOAST is a project financed by the Research Council of Norway (RCN) for the period 2019–2023 under the programme MILJØFORSK [65]. The primary objective of this project is to monitor, manage and preserve cultural heritage sites, environments and landscapes under climate change and development pressure. The project brings together an international multidisciplinary team specialised in the evaluation, assessment, protection and conservation of heritage, from research and education institutions in Norway (Norwegian Institute for Cultural Heritage Research—NIKU; University of South-Eastern Norway—USN; Geological Survey of Norway—NGU; Foundation for Industrial and Technical Research—SINTEF), UK (University of St Andrews), and Russia (Kazan Federal University).

#### **2. Study Area and Cultural Heritage Background**

#### *2.1. Study Area*

The study area is located on Spitsbergen, which is the largest island of the Svalbard archipelago (Figure 2a), governed by Norway as it was established by the Spitsbergen Treaty from 9 February 1920. The Svalbard archipelago is situated between the North Pole and Northern Norway [49]. Svalbard thus forms a terrestrial node separating the Barents Sea, the Greenland Sea and the Arctic Ocean. This position gives a complex climatic situation where different oceanic currents influence both long-term climate and short-term meteorology of the islands. The focus of this study is the shoreline of Hiorthhamn—on the northeastern shore of Adventfjorden (Figure 2b), approximately 3 km north of Longyearbyen (Figure 2c), which is the major administrative, tourist and scientific centre in Svalbard, also named "European Gateway" to the Arctic [21].

**Figure 2.** (**a**). Geographical location of the study area in the context of Svalbard; (**b**) Local context; (**c**) Location of the protected cultural heritage (base orthophoto from Norwegian Polar Institute) [66].

The Hiorthhamn fjord side is step-shaped with south-west-oriented slopes running down from the top of the Hiorthfjellet mountain to the Advent fjord shore. The bedrock setting of the study area is dominated by the subhorizontal layering of sedimentary rocks of Early Permian to Eocene age, including some coal seams that formed the base for the mining activities [67]. Sediment of varying thickness is draped over the bedrock slopes and steps, primarily consisting of weathering material together with different types of slope-process deposits and some patches with glacial till. The highest coastline after deglaciation was approximately 62 m above the present sea level [68] and the lowermost parts of slope form a low angle beach landscape.

The climate is a Polar tundra climate, according to the Köppen climate system. The mean annual air temperature in central Spitsbergen is −5 ◦C, the mean annual precipitation is 190 mm (concentrated mainly during the summer season), and dominant winds come from southwest or west. Svalbard is an area with continuous permafrost, meaning all but the surface is permanently below 0 degrees Celsius. A permafrost landscape has a so-called active layer, meaning the surface of the soils, often ca 0.5–1 m thick, which thaws in the summer and refreezes in the winter. The arctic climate influences the landscape through high activity of frost weathering together with surface cryogenic processes, such as solifluction, which are associated with permafrost and the active layer. The cryogenic surface processes are all also effected by gravity, resulting in a net downslope movement even on only gently sloping ground. Other important processes are rock fall, snow avalanches and debris flows, that also bring material from higher to lower elevations [69]. These climatically linked processes have in central Svalbard resulted in a "smoothing" of the landscape, with few sharp topographic features and a rounded topography. At the start of the 20th Century, the records have shown a strong warming trend (starting with the 1920s), where within a period of 5 years the mean annual air temperature changed from −9 ◦C to −4 ◦C. This warming is generally interpreted as representing the end of the Little Ice Age (LIA) [12].

Politically, Svalbard is governed by a representative of the federal government, Sysselmannen (the Governor) of Svalbard. Sysselmannen must ensure that activities on Svalbard are in line with the Norwegian national safety and security goals, as well as being the administrator of all cultural heritage and natural values of the archipelago [70]. In the Norwegian national heritage database, Askeladden—Riksantikvaren [71], there are 2025 listed cultural heritage sites in Svalbard; 342 built heritage sites and 1683 archaeological sites. An overwhelming majority of these sites are coastal, thus being exposed to geohazards, development and wear and tear from tourism.

#### *2.2. Cultural Heritage Background*

The Arctic environment is ideal for long-term preservation of archaeological remains, both for artefacts and environmental proxies [43]; however, studying cultural heritage in the Arctic is quite challenging due to a changing warmer climate [63]. The discovery of Spitsbergen was made by a Dutch crew in 1596 led by Willem Barentsz; this led to a flow of Russians and Europeans craving for glory, money and fame [49]. There were different activities that brought explorers to Svalbard; whaling (starting from the 17th Century), hunting and scientific exploration (starting with the 18th Century), industrialisation (19th and 20th Centuries) [72].

All the "traces" that were left behind by these activities came in time to be listed cultural heritage. Nowadays, Svalbard's cultural heritage is under the protection of the Svalbard Environmental Protection Act, which states that all remains from before 1946 are automatically protected as listed cultural heritage sites. Moreover, all traces of human graves, including crosses and other grave markers, as well as bones and bone fragments found on or below the surface of the ground are automatically protected regardless of their age. Same rules apply to skeletal remains at slaughter sites for walruses and whales and in connection with self-shooting traps for polar bears; the protected cultural heritage sites have a 100-m buffer protected area around the sites. The buffer area is granted the same importance as the site itself [73].

From 1955, numerous archaeological expeditions took place in Svalbard; they were organised by the Soviet Union, the Netherlands, Denmark, Sweden, Finland, Poland and Norway [74]. Early explorers of Svalbard have shown a great interest in cultural heritage, like the Dutch whaler Cornelius Gisbert Zoorgdrager, who in his book "Bloyende Opkomst der Aloude en Hedendaapsche Groenlansche Visschery" from 1720, mentions all the places on the coast of Svalbard where he has seen remains of buildings and blubber ovens during his whaling expeditions around 1700. The first real archaeological survey was undertaken by Christian Keller in 1967, as an experimental project; and the first scientific excavation took place in 1958 on Midthuken in Bellsund by the Finnish scholar Tegengren. It is worth to mention the fact that many Svalbard cultural heritage sites have been threatened by the effects of erosion since their existence; only in 1989 their destruction begin to be acknowledged, resulting in a recommendation that research should be focused on threatened cultural heritage [75], along with the need for geological mapping of Svalbard coasts [49].

Starting from 1978, the Russian archaeologist V.F. Starkov was leading extensive excavations in Svalbard of the Russian hunting stations; he was followed by Norwegian, Dutch and Polish archaeologists [76,77]. Russian settlements represent the most numerous group of historical sites on Spitsbergen, being spread all over the territory of the archipelago [78]. Exploration of Svalbard's natural resources (oil, petroleum, coal) left behind significant traces; coal was found and used as early as 1610 by whalers in Svalbard on their ships. Permanent coal-mining communities were established in Longyearbyen and Hiorthhamn (mining active 1906–present) (Figure 3a–c), Grumant (1910–1962), Ny-Ålesund (1916–1962), Bjørnøya (1916–1925), Svea (1917–2016) (Figure 3d–f), Barentsburg (1920–present) and Pyramiden (1927–1998) [79].

**Figure 3.** (**a**) Former coal-mining infrastructure in Longyearbyen; (**b,c**) Former coal-mining infrastructure in Hiorthhamn; (**d–f**) Remains of the coal-mining exploitation in Svea; (**g–i**) Current images from the Russian coal-mining city of Barentsburg.

Extensive whale and walrus hunting took place in Svalbard starting from 1611. The hunting activities started with the English and Dutch whalers, followed by Russian and Norwegian. This lasted about three hundred years and lead to killing of thousands of Greenland right whales and thousands of Atlantic walruses [80]. The remains of the hunting stations, slaughter sites, blubber ovens and the whalers' cemeteries are all listed heritage sites.

Cultural heritage in Svalbard also includes remains like plane wrecks from WWII, and even these young sites experience the same threats, exposure to geo-hazards, tourism and development. The historic remains at Svalbard are considered international cultural heritage [81] and the Governor of Svalbard states that conserving cultural heritage is amongst the most important environmental goals. Thus, monitoring and maintenance of cultural heritage sites is considered an important task [82].

#### *2.3. Hiorthhamn Mining Settlement*

In 1917, a Norwegian company, A/S De Norske Kulfelter Spitsbergen, started constructing the mining facilities of what was to become Hiorthhamn city, but coal mining itself did not start up until 1919. Buildings and machines were moved from nearby Advent City, an English coal mining city that was closed down in 1909. The mine at Hiorthhamn was in use only for a few years, from 1919 until 1921 and again from 1938 until 1940. Nevertheless, an entire mining town was built up at the foot of the mountain Hiorthfjellet (Figure 4) [83]. At the start of the 20th century, the business activities in Svalbard was characterised by great optimism. However, many of them failed due to lack of knowledge.

**Figure 4.** Detail of Hiorthhamn mining settlement from 1936, with emphasis on the three sectors analysed in this study (image courtesy: Norwegian Polar Institute) [84].

The central cable car station in Hiorthhamn was built in 1939 (Figure 3c) and is listed in the National heritage database Askeladden under the ID 93040-6. It consists of a large wooden loading construction. The inside machinery is in good shape, with wheels and gears on cast foundations. A wooden staircase leads to the top of the construction. Remains of two railway tracks come out from the southern side of the building directly into the sea. From the side towards the mountains, rests of cable car wires are hanging [70].

Today the former mining settlement is the second largest gathering of listed buildings at Svalbard (next to Ny–Ålesund) and is considered an important cultural heritage site. The listed historic site counts a total of 13 standing buildings and remains from mining infrastructure, such as cable car poles, a cable car station and railway tracks. The mine itself is situated high up in the mountain, 580 meters above sea level, and the coal was transported on an aerial cableway from the mine to the seaside. Except a few buildings up in the mountains, the whole site is situated close to the coastline, vulnerable to coastal processes [85].

#### **3. Materials and Methods**

Coastline dynamics were evaluated on the basis of modern GIS and RS techniques. In this study, we used the remote sensing data available from this area: a topographic map from 1927 (scale 1:50,000), aerial images from 2009, and field surveys from August 2019 to assess the coastal change over a period of 92 years. The aerial images from 2009 were obtained by the Web Map Services (WMS) developed by the Norwegian Polar Institute (NPI) [69]. All the data was integrated into a GIS using the WGS\_1984\_UTM\_Zone\_33N reference system. The 1927 topographic map was georeferenced in ArcGIS ver. 10.4 (ESRI, Redlands, CA, USA) by adding six control points by using first order polynomial, auto-adjust transformation. The total RMS error was 4.5 m, which is close to the global errors mentioned in Bourriquen et al. [55].

The geomorphology and surface sediments of the study area was mapped through a combination of remote sensing and four days of field validation in beginning and end of August 2019. Field data was collected using Getac F110 rugged field tablets with integrated GPS. The base for the geomorphological map is the 2009 aerial images available together with elevation data (DEM and contour lines) and the map was made in ArcGIS ver. 10.6 software.

Change rates between shorelines extracted from old maps, aerial images and field surveys were computed using the DSAS v.5 extension of ArcGIS [65]. Different parameters were calculated: Shoreline Change Envelope (SCE, expressed in m), Net Shoreline Movement (NSM, in m), End Point Rate (EPR, m/yr), and Linear Regression Rate (LRR, m/year). SCE describes the variability of each transect considering the maximum spatial recorded displacement of shoreline, but without considering the time span. NSM deals only with the dates of two shorelines, reporting the distance between the oldest (1927) and newer (2019) shorelines for each transect, whereas this movement may not be the maximum shoreline displacement recorded. EPR is calculated in m/yr by dividing the distance of shoreline movement by the time elapsed between the oldest and the most recent shorelines. A total of 146 transects at 10 m alongshore intervals were analysed using the EPR parameter of DSAS to obtain annual mean rates of shoreline change. LRR is based on the overall minimum of the squared distance to the known shoreline using all available data to find the best-fit line and is being recognised as a very useful tool for computing long-term rates of shoreline change. The confidence interval in DSAS was at 95%. The shoreline was divided into three sectors in order to get a better picture of the final geomorphological and DSAS analyses.

Field surveys were made in August 2019 using a GPS system comprised of an Altus APS–NR2 antenna and a Carlson controller. Two initial points were surveyed with the GPS system (±0.5 m accuracy), then the survey was conducted with a Trimble S5 Series Motorized total station, along with a Trimble TSC3 controller (Figure 5a). In total, a number of four fixed points were established and recorded for present and future monitoring of the shoreline (Figure 5b), to gather data to support the geomorphological mapping (Figure 5c) as well as for monitoring of surface soil movements over time (Figure 5d).

**Figure 5.** (**a**) Total station and controller on one of the fixed points; (**b**) Coastline survey; (**c**) Geological mapping of the upper area of Hiorthhamn; (**d**) Survey of surface soil movements.

#### **4. Results and Discussion**

#### *4.1. Geomorphological and Surface Sediment Mapping*

Results from the geomorphological mapping are shown in Figure 6. In addition to the result from the present study, it includes some sub-surface interpretations by [86]. The map coastline is based on aerial photos of 2009, whereas the present coastline based on GPS field surveying (2019) is shown as a purple line.

The gently sloping plain down to the shoreline in the northern part of the study area (Sector 1) is mainly covered by peat and grass vegetation, and the low wave-cut section (Figure 7a,e) clearly shows how today's organic material overlies gravelly beach sediments. A series of low subdued gravelly ridges oriented obliquely from the present-day beach and inward is only visible in the aerial photos and is interpreted as relict beach ridges formed during earlier beach migration towards south east. This development requires a net sediment supply to Sector 1 and long-shore sediment drift from northwest towards southeast, which is confirmed by the ongoing build-up of a small spit bar shifting the small tidal channel outlet towards southeast between 2009 and 2019 (south-eastern end of Sector 1 and Figure 7a).

The wave-erosion shoreline scarp (Figure 7f) which cuts the old subdued beach ridges at an angle, together with observed in-land displacement of over-wash deposits parallel to the beach, indicate a change in shoreline environment towards increased erosion. This change may be the result of a shift in dominating wind-, and hence, wave direction during especially erosive storm events. The detailed land survey data of the present shoreline (Figure 6) supports this interpretation.

**Figure 6.** Map showing geomorphology and surface sedimentology in the study area. Simplified genetic sediment classifications of Sector 3 are based on the study of Lønne and Nemec [86].

**Figure 7.** (**a**) Overview of primarily the active fan with boxes marking the position of interpreted coastal sectors. Note that Sectors 1 and 3 extend outside the photo. Blue arrow indicates the long-shore sediment transport direction forming a small spit-bar in the eastern end of Sector 1; (**b**) Photo towards north west, from the coastal escarpment of the relict alluvial fan in Sector 3, towards sectors 1 and 2. Dotted line delineates the post 2009 debris flow advance of the coastline in Sector 2; (**c**) Photo from the middle of Sector 2 towards east showing shoreline erosional edge (left of person) and the recent debris flow deposit (dotted line); (**d**) Sector 3. Fresh small slump from the shoreline escarpment in the relict alluvial fan; (**e**) The peat covered relict beach ridge plain in Sector 1; (**f**) The active erosion edge of the shoreline in Sector 1 at high tide (rifle for scale). Note the coarse sorted beach sediment in the section and the thin strip of beach separating the shoreline scarp from the water line.

Sector 2 of the study area includes the distal part of an active alluvial fan system which is subjected to both erosional and depositional processes. A braided fluvial system is transporting coarse and poorly sorted sediments towards the shoreline and episodic debris flows or debris torrents transport larger masses, which are later remodelled by the more continuous fluvial processes. According to [86] detailed descriptions of sedimentary processes and environment, fluvial channels migrate laterally over the fan throughout the summer due to melting of a snowpack that deflects the water flow. Dynamic shifts in location of the fluvial transport and debris flows down to the swash zone today (Figure 7b,c) result in an irregular shoreline with significant yearly and seasonal changes (Figure 5b).

The south-eastern Sector 3 of the study area forms a long erosional escarpment through a raised abandoned alluvial fan–delta [86]. Wave action forms the steep escarpment by undercutting the raised delta sediments. Slumping of primarily sandy sediments in the over-steepened escarpment moves sediments down to the beach (Figure 7d) where they are moved into the sea and further eastwards by wave action and long-shore transport. The shoreline displacement between 2009 and 2019 (Figure 6) clearly highlights the shoreline retreat resulting from the undercutting and erosion.

Solifluction is a ubiquitous process transporting sediments downward from the slopes above both Sectors 1 and 2. The same net downward movement of surface material is the result of active layer detachments forming scars and slump (Figure 6). Although these slope processes have no direct impact on the shoreline sedimentary environment today, they contribute to the total down slope sediment flux, and enhanced thawing of permafrost due to a warmer future climate may increase this effect, subsequently resulting in increased net shoreline input.

Lønne and Nemec [86] described the onshore area of Sectors 2 and 3 in detail and characterized the combined feature as a Holocene fan delta and that the depositional system may hypothetically extend under water up to 400 m offshore. The postglacial lowering of the relative sea level and a simultaneous decline in glaciofluvial activity enhanced fluvial erosion of the delta fan and formed the younger, and still active, incised braided system feeding Sector 2 with fluvial-, and debris-flow sediments. Apart from some snow avalanche deposition in the uppermost part, there is today no sediment input from higher up to the raised delta fan of Sector 3. The shoreline in Sector 3 has no significant terrestrial sediment input and future shoreline changes are more dependent on marine-controlled processes and erosional potential. Since the escarpment today is dependent on the internal permafrost in the raised delta-sediments, future increase in soil temperature and active layer thickness may lead to increased speed of undercutting and subsequent shore–line retreat.

#### *4.2. Shoreline Displacement*

Using the DSAS analysis of the coastline, the following parameters were computed EPR (Figure 8a), SCE (Figure 8b). EPR shows the highest rates between the oldest (1927) and the most recent (2019) shoreline in the northwestern part of our area of interest. The shoreline from 2009 was not included in the analysis and is shown in Figure 8a only for orientation purpose.

The final rates were reclassified into five classes by using the Natural Breaks (Jenks) method: very high (−0.76 to −0.57 m/year), high (−0.56 to −0.34 m/year), moderate (−0.33 to −0.20 m/year), low (−0.19 to 0.06 m/year) and very low (0.07 to 0.45 m/year). The average erosion rate for the entire coastline is −0.21 m/year. The highest EPR values are located in the north-western part of the coastline (Sector 1), where the most important part of the cultural heritage (the loading dock) is located, making it very vulnerable to coastal erosion. Positive values are recorded in the middle part of the interest area (Sector 2), where sediments are being deposited from the alluvial fan. The southeastern part of the coastline (Sector 3) can be framed into the average erosion rate of −0.21 m/year. However, Sector 3 has no cultural heritage assets. SCE (m) values represent the greatest distance among all the shorelines that intersect a given transect and they are visible in Figure 8b. The highest interval (48.22–70.84 m) is located within Sector 1 and the lowest values (2.82–17.95 m) in Sector 2.

NSM, expressed in m is visible in Figure 9a. Sector 1 is characterised by very high and high erosion rates, while Sector 2 and Sector 3 by moderate and high accretion and moderate and high erosion, respectively.

**Figure 8.** (**a**) End Point Rate and (**b**) Shoreline Change Envelope of the coast in Hiorthhamn area as calculated by DSAS (1927, 2009 and 2019 shoreline position used).

The LRR parameter is detailed in Figure 9b and highlights the areas with erosion, stability and accretion, divided into five classes, as follows: very high erosion = −0.73 to −0.52; high erosion = −0.51 to −0.33; moderate erosion = −0.32 to −0.15; stability = −0.14 to 0.07; accretion = 0.06 to 0.45. The average LRR is −0.16 m/year. There was a total number of 112 transects, out of which 76 (67.86%) are erosional and 36 accretional (32.14%). Transect number 102 has the highest value for an erosional transect and is located in Sector 1, in the same area as the loading dock is located. Whilst the transect with the highest accretion value is number 93, located in the upper part of Sector 2. A significant percentage of the shoreline (48%) has transects with statistically significant erosion. This highlights the vulnerability of the coastal cultural heritage.

A study that focused on approximately 62,000 km of Arctic coasts [7] reported an average erosion rate of −0.5 m/year. The highest erosion rate in the Arctic was registered on Alaskan coast (Drew Point) with a value of −8.4 m/year, while the lowest erosion rate was reported in Svalbard. Other reported erosion rates in Svalbard range from −0.5 to −4.5 m/year for Longyearbyen [21], −0.26 m/year for Isbjørnhamna [48], −0.34 m/year and −0.47 m/year for ice-poor cliffs and ice-rich cliffs, respectively [32]. With a value of EPR of −0.21 m/year, our study area can be framed in around the average erosion rate previously reported in Svalbard. A recent study [51] reported EPR values of −0.19 m/year, which has increased from −0.07 m/year for period 1936–2007. This means that large portions of Svalbard coasts, strictly referring to those surveyed in the past, have a stable erosion rate. This is due to the fact that in many cases high erosion rates compensate with high accretion rates. The complex interaction between different sediment transporting processes and shoreline processes is well exemplified in this study and indicates that understanding of the upland systems is important for a predictive interpretation of the shore-line changes into the future.

**Figure 9.** (**a**) Coastal stability derived from Net Shoreline Movement analysis (in m) (very high erosion = NSM −70.84 to −53.33; high erosion = NSM −53.32 to −30.41; moderate erosion = NSM −30.40 to −9.01; stability = NSM −9.00 to 5; moderate accretion = NSM 5.01 to 24.12; high accretion = NSM 24.13 to 41.96; (**b**) LRR parameter expressed in m/year.

#### *4.3. Cultural Heritage*

As shown by Hollesen et al. [43], the cold wet climate of the Arctic is the perfect environment in which cultural heritage can be preserved; following the climate variables and tendencies nowadays, the "friendly" environment is slowly disappearing. This leads to a systematic degradation of cultural heritage assets in the Arctic, which is the case of the present study, where coastal cultural heritage is in great danger of being destroyed.

The most vulnerable element of protected cultural heritage is the cable car station or loading dock (Figure 10a,b). Other remains that are washed away by the sea are visible in Figure 10c from 2011; those remains are no longer present in Figure 10d from 2019. This shows the power that sea exerts on the shoreline, along with longshore currents. To this is added the increase in global sea level, which is caused by an increase in temperature of the oceans and significant mass loss of the Svalbard glaciers [87]. Other cultural heritage remains exposed are the exploitation rail tracks (Figure 10e), along with wrecks of heavy machines, loading carts, small loading wagons, metal barrels and other industrial waste that is scattered in Sector 1, but very densely seen in Sector 2.

**Figure 10.** Details over Sector 1; (**a**,**b**) Details highlighting the erosion of what was in the past the loading dock; (**c**) Figure from 2011 highlighting remains from the exploitation equipment on the beach (courtesy of A.-C. Flyen); (**d**) The only remained equipment on the beach in 2019; (**e**) Exploitation rail tracks affected by coastal erosion.

#### **5. Conclusions**

Historical evolution and temporal morphodynamics of shoreline position is of great importance in evaluating the spatial dynamics of the coastal system behaviour and of cultural heritage assets located in coastal areas. DSAS analysis of shoreline changes in years 1927–2019 shows both erosion and accretion for different sectors of the coastline. Sector 1, which is located in the north-west part of the investigated coastal area and has a length of 460 m is characterised as having high and very high erosion; Sector 2 has a length of 230 m has a moderate accretion, and Sector 3 with a length of 610 m showed moderate erosion, respectively. Potential changes in wind and wave direction and possibly intensity since 1927 have interacted with sediment availability and processes transporting material down to the beach zone. The type of shoreline sediments exposed to erosion also affects the resulting rate of erosion with already well-sorted and fine-grained beach sediments being more susceptible to abrasion than more course, poorly sorted debris flows and braided river sediments. Hence, the same change in wave energy might have a different result in shoreline displacement at different positions along the shore—causing a varied erosional hazard to pre-existing cultural heritage. If the active-layer thickness increases in a warming climate, this will contribute significantly to the potential erosion rates of sedimentary coast-lines. The fact that almost half of the shoreline has statistically high erosion transects, highlights the danger for the coastal cultural heritage and the need from local authorities of a future sustainable management plan.

**Author Contributions:** Conceptualisation, I.C.N., L.R.; methodology, I.C.N., K.S., L.R.; software, I.C.N.; validation, I.C.N.; formal analysis, I.C.N.; data curation, I.C.N., K.S., L.R.; writing—original draft preparation, I.C.N.; writing—review and editing, I.C.N., K.S., L.R., A.-C.F., V.V.M.; visualization, I.C.N.; supervision, I.C.N.; project administration, V.V.M.; funding acquisition, V.V.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Research Council of Norway, CULTCOAST—Monitor, Manage and Preserve Sites and Landscapes under Climate Change and Development Pressure (MILJØFORSK), project number 294314.

**Acknowledgments:** The authors are grateful to the Norwegian Polar Institute (NPI) personnel for the assistance in the logistic fieldwork. Two anonymous reviewers are kindly acknowledged for their fruitful comments.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **GIS-Based and Statistical Approaches in Archaeological Predictive Modelling (NE Romania)**

#### **Ionut Cristi Nicu 1,\*, Alin Mihu-Pintilie <sup>2</sup> and James Williamson <sup>3</sup>**


Received: 7 October 2019; Accepted: 25 October 2019; Published: 27 October 2019

**Abstract:** Archaeological predictive modelling (APM) is an important method for archaeological research and cultural heritage management. This study tests the viability of a new statistical method for APM. Frequency ratio (FR) is widely used in the field of geosciences but has not been applied in APM. This study tests FR in a catchment from the north-eastern part of Romania to predict the possible location(s) of Eneolithic sites. In order to do that, three factors were used: soils, heat load index and slope position classification. Eighty percent of the sites were used to build the model, while the remaining 20% were used to externally test the model's performance. The final APM was made with the help of GIS software and classified into four susceptibility classes: very high, high, medium and low. The success rate curve and the prediction rate curve reported values of the area under curve (AUC) of 0.72, and 0.75 respectively. The Kvamme's Gain value for the model has a value of 0.56. Therefore, the final APM is reliable, so FR is a viable technique for APM. The final map can be successfully used in archaeological research, cultural heritage management and protection, preventive archaeology and sustainable development.

**Keywords:** cultural heritage; frequency ratio; AUC; predictive modelling; GIS; Kvamme's Gain; north-eastern Romania

#### **1. Research Aims**

The present work aims to present the results obtained using new statistical method (FR—frequency ratio) for archaeological predictive modelling (APM) in a catchment from north-eastern (NE) Romania. Despite the use of APM globally, only one study was undertaken in Romania, focusing on post-Roman sites [1]. This study aims to fill the gap and to propose and test a new statistical method—frequency ratio applied in APM. This study also aims to prepare the first APM for Eneolithic sites in the NE part of Romania. The final APM will be of great help for cultural heritage management, preventive archaeology, planners, stakeholders, etc. In combination with landslide/gully erosion susceptibility maps, this could represent a powerful tool for cultural heritage mitigation, vulnerability assessment and improving the sustainability of cultural heritage resources.

#### **2. Introduction**

The use of GIS and statistical modelling to depict the possible locations of archaeological sites has faced an upward trend over the last decades [2,3]. This has not been matched by an improvement in the statistical methods applied until recently, when a proliferation of different methods can be seen. Currently, APM represents a powerful tool for preventive archaeology [4], cultural heritage management [5,6] and improving national-scale archaeological inventories [7,8]. APMs have been used successfully in different geographical areas of the globe including Africa [9,10], Europe [1,11–13], Asia [14,15] and America [16–18]. Different intuitive (qualitative) and quantitative statistical methods have been used to identify the locations of archaeological settlements. The most commonly used method for building qualitative predictive models is the analytic hierarchy process [9]; the main issues of this method are related to the subjective judgement of the expert(s), and that intuitively finished models can easily be influenced by prior biases in the archaeological record.

Quantitative methods can limit bias, which helps to reduce the effects of this problem and can improve the replicability of the workflow for use in other locations. The most commonly used geospatial formulae for predictive modelling are maximum entropy (MaxEnt) [19], logistic regression [20], evidential reasoning [14], fuzzy logic [21], weight of evidence [22] and multivariate statistics [23].

Statistical modelling of the landscape to find suitable areas for prehistoric settlement choices has been well used by geographers and archaeologists [14,24–26]. This method can effectively be conceptualised as a form of archaeological prospection. The most commonly tested era is the Neolithic, likely due to the pivotal nature of this period, and the fact that landforms were like modern features. Different areas (e.g., Italy [4], Greece [12], Scotland [27]) have distinct landscape features which may have been decisive for prehistoric people when choosing where to place a settlement. For example, in Italy [4], the favoured features were alluvial plains, terraces, caves (lithic industry), and low slope areas, while in Greece [12], the area around Thessaly represented an important connection between the islands of the Aegean Sea with the northern and southern parts of the country, being at the same time characterised by permanent inhabitation. In mainland Scotland [27], on the other hand, the locations of chambered cairns, timber halls and megalithic stone tombs influenced the way Neolithic people chose to place their settlements, as the megalithic stone tombs were perceived as land tenure and symbolize territorial control.

While this method has been widely applied in research archaeology, issues have been reported in CRM use in the Netherlands [28]. These, however, are extremely specific and may stem from the lack of periodisation or landscape specialisation. The lack of widespread use of this model in a European context suggests that testing this methodology using different methods, and for different periods, is extremely important and a desideratum.

#### **3. Archaeological Background and Study Area**

Cucuteni culture is considered to be the last great Eneolithic civilisation of Old Europe. It is part of the Cucuteni–Ariusd–Trypillia Cultural Complex, covering an area of approximately 350,000 km2, comprising Romania, Ukraine and the Republic of Moldova [29]. The north-eastern part of Romania is well known for its high density of Eneolithic sites; many studies have approached the mobility of Eneolithic settlements and their vulnerability [30] towards natural hazards [31,32] and anthropogenic impact [33]. The discovery of the Cucuteni Culture in the first half of the 19th century marked the beginning of archaeological research in north-eastern Romania. Many archaeological and interdisciplinary studies aimed to better understand how prehistoric people lived [25,34] and moved across the landscape in search of new resources [35–37]. As shown by [38], Romania has an incomplete registry/inventory of archaeological sites. This means that local authorities and cultural heritage planners are incapable of mitigating threats to archaeological sites [39], including looting [40]. The method proposed in this study can be used as a starting point to produce APMs for the entire Moldavian Plain and Plateau and to be used along with the natural hazards susceptibility maps in order to establish the present state of the cultural heritage.

The study area is in the North-eastern part of Romania (Figure 1a), at the contact of Moldavian Plain with Suceava Plateau; it has an area of approximately 340 km2. In geomorphological literature, the area is known as Dealul Mare–Harlau Coast (Coasta Dealul Mare–Harlau) [41] (Figure 1b). Previous studies have shown the tremendous archaeological potential of this area [25,29–39]. Within this area

we created a database of 100 Eneolithic sites, using the online databases of the Institute of Cultural Memory (CIMEC), National Archaeological Registry (RAN) and National Heritage Institute (INP), as well as numerous field trips.

**Figure 1.** (**a**) Geographical location of the study area in Romania; (**b**) Location of Eneolithic sites, divided into training sites and testing sites.

#### **4. Materials and Methods**

Within the study area, there are a total number of 100 Eneolithic sites. Out of these, 80 sites (80%) were randomly selected for model training, while the other 20 sites (20%) were used to test the performance of the model (Figure 1b). This was done to provide an external dataset for testing the model, and due to the small area, number of sites, and the fact that an internal random selection was not to be performed, this high ratio of excluded to included sites was acceptable.

In this study, three factors were used to determine the most likely areas to host Eneolithic sites: soil type, heat load index (HLI) and slope position classification (SPI). Originally, several other possible

factors were chosen and tested in combination with each other, including slope, aspect, landforms and distance to rivers; these lowered the predictive strength of the model, and so were abandoned.

Soil types were extracted from The Soil Map of Romania (scale 1:200,000) (Figure 2a); the classification was updated after [42]. The Chernozem soils seem to have been the most important because they are very fertile. As a result of the importance of farming, these areas would have been extremely important. Soil represents a significant factor used in APM [20].

In order to prepare the topographically derived factors for the model (HLI and SPI), a Lidar derived digital elevation model (DEM) with a resolution of 1 × 1 m was used. HLI is a novel method in the geosciences, and this is the first time it is being used in APM; it was calculated with the help of Geomorphometry & Gradient Metrics Toolbox of ArcGIS [43]. By using this factor, we have included slope and aspect as factors, because these are components in the HLI calculation. HLI takes these into consideration by "folding" the aspect so that the highest values are south-west, and the lowest values are north-east [44]. The method also considers the steepness of the slope. In this study, HLI (Figure 2b) was computed and reclassified into five classes 0.097–0.59, 0.59–0.64, 0.64–0.67, 0.67–0.71, 0.71–1.04.

The last factor used to compute the final APM was SPI [44]. This extension computes TPI (Topographic Position Index) grids from elevation grids, providing a simple method to classify the landscape into slope position and landform category using the TPI values [25,45]. SPI was computed with the help of Topography Tools Toolbox of ArcGIS (Figure 2c); prehistoric people's preference for midslope and flat areas is once again highlighted following this classification. Placing the settlements on midslopes was the first criteria, this being enough high to protect them from flooding [25] and other natural hazards of which they were aware [46].

**Figure 2.** The factors used to compute the archaeological predictive modelling; (**a**) Soil type; (**b**) HLI (Heat Load Index); (**c**) TPI (Topographic Position Index).

FR represents a quantitative statistical method used in various branches of geosciences; it offers better prediction and validation results when compared to other statistical methods. For example, it is commonly used in landslide susceptibility mapping [32,46–48], gully erosion susceptibility mapping [49], groundwater assessment and contamination [50,51], flood susceptibility [52], deforestation susceptibility [53], etc. Despite its broad applicability, it has not been used in APM. FR has several advantages over other methods, as it allows all factors to be examined while prioritising those with the greatest applicability. This study tests the view that there is a casual relationship between Eneolithic settlement choice and the features of each location, which can be described through the metrics used in this study. In order to determine the frequency ratio for each class factor, the ratio between Eneolithic sites occurrence and non-occurrence was calculated (Table 1, column 7).

The Eneolithic sites used in this study are representative for Eneolithic culture (chronological framework: ca. 5000–3500 BCE) from the NE part of Romania [54]. Figure 3a shows the site Cucuteni (Cetăt,uia), the eponymous site of the Cucuteni culture, known as the last great Eneolithic civilisation of Old Europe, part of Cucuteni–Trypillia Cultural Complex [55]. It represents one of the most important archaeological discovery in the 19th century by the German archaeologist H. Schmidt [56]. Figure 3b,c illustrate the sites from Giurges,ti (Dealul Mănăstirii/Chetrosu) and Costes,ti (Cier/La S, coală), respectively. They represent sites iconic for the Cucuteni culture, and each year archaeological investigations take place [57].

**Figure 3.** Representative Neolithic sites from the study area; (**a**) Cucuteni (Cetăt,uia) site, the eponymous site of the Cucuteni culture; (**b**) the site of Giurges,ti (Dealul Mănăstirii/Chetrosu); (**c**) the site of Costes,ti (Cier/La S, coală).

Out of 100 Eneolithic sites within the study area, 13 sites belong to the Precucuteni period (ca. 5000–4600 BCE), 35 sites belong to Cucuteni A period (ca. 4600–4100 BCE), 8 sites are framed to the Cucuteni A-B period (ca. 4100–3850 BCE), 24 sites belong to the Cucuteni B period (ca. 3850–3500 BCE), and 20 sites belong to the same culture, but with an unknown phase (hereafter named Cucuteni unknown) [54–57]. For each site, representing a point, a buffer area of 100 m was made to obtain an average surface of 3 ha for each site, as this is the average surface of a Cucuteni site [29]. This was done to prevent the site from being described as a single point within the GIS, which is a common fiction in geospatial studies [58].


**Table 1.** Frequency ratio values for the conditioning factors.

The weight of each factor was calculated, and by summarising the weights we obtained the APM using the following equation:

$$\text{APM} = \sum \text{FR}\_i \tag{1}$$

where FR*<sup>i</sup>* is the FR of each factor type and FR is the area where Eneolithic sites occurred.

The FR method selected soils as being the most important conditioning factor, followed by SPI and HLI. The final susceptibility maps were reclassified into four susceptibility classes low, medium, high and very high. The resulting model must be tested and validated. In addition to the use of the classic procedure to validate the APM [59] (Equation (2)), we used the area under curvature (AUC). Based on the issues identified by [60] which argues the use of topographical features in APM, we chose to address this issue by using in our modelling just three factors.

$$\text{Kvamme's Gain} = 1 - \text{(Area Percentage/Percentage of Sites)} \tag{2}$$

#### **5. Results and Discussion**

The final APM made with the help of the FR method is shown in Figure 4a; the classes of the APM are low probability, medium, high and very high probability. The final raster (1 × 1 m/pixel) was reclassified by using the Natural Breaks (Jenks) method. Calculating the statistics of the APM (Table 2), low probability area covers 0.9% of the total area, while medium, high and very high probability areas cover 17.1%, 46% and 36% of the study area, respectively. High and very high probabilities are related to upper, toe and flat areas, as highlighted by [25]. Medium and low probability areas are mostly associated with very steep slopes, usually affected by geomorphological processes (e.g., landslides). A significant number of sites (33) are located in areas with medium probability, fact that was also highlighted by [25]. It was shown that following the SPI analysis, a significant number of Eneolithic sites were located on middle slopes, as Eneolithic people were aware of the impact of floods [25] and landslides [46]. This is highlighted in this study by the ratio obtained in Table 1 for midslope, with a value of 0.18; higher values have sites located at the toe slope and upper slope with values of 0.21 and 0.25, respectively.

The results of this study complete the conclusions obtained by [25]. Further research is needed in order to frame the Cucuteni unknown sites into one of the four periods (Precucuteni, Cucuteni A, Cucuteni A-B, and Cucuteni B). This represents one of our future goals in our endeavour of fully understand the mobility of Eneolithic populations.

When it comes to the APM validation, reasonably good results were obtained. The success rate (Figure 4b) of the APM produced with FR has been assessed by applying the receiver operating curve (ROC) with the area under curve (AUC), having a value of 0.72. The prediction rate curve has a value of 0.75, which is considered as good. Other studies [61] have reported AUC values of 0.65, therefore our results are pertinent, taking into consideration the resolution of our data used in the study (1 × 1 m/pixel).

A particularly widely applied measure of predictive model quality is Kvamme's Gain, which gives a ratio of the precision of a model against its accuracy [62]. The goal of this is to show that the area described is small enough when compared with its accuracy. The result of this is usually between 1 and −1, with a positive result suggesting that the model works, while a negative result would suggest that the model is predicting areas without features. The Kvamme's Gain in this case is 0.56, which suggests that this has been a success.

$$\text{Gain} = 1 - (\text{36/23}) = 1 - 1.56 = 0.56$$


**Table 2.** Statistical results following the use of frequency modelling (FR) method for archaeological predictive modelling (APM).

The method proposed in this study can be useful in making APMs for the whole area of Moldavian Plain; this represents an area with a high density of Eneolithic sites, much of them being still undiscovered and in this way cannot be protected or included in the cultural heritage database. Based on this model, this problem can be solved in an easy and sustainable way. This can be the base for future space-time APMs. Discovering more Eneolithic sites could bring us close to understand their mobility. Using APMs together with natural hazards susceptibility maps will have a significant impact in reducing the effect of hazards on cultural heritage, risk reduction, land use planning, hazard mitigation and for developing sustainable methods and practices for cultural heritage preservation.

The results of this study shows the fact that FR could be successfully used and applied to identify areas likely to host archaeological sites, in this case Eneolithic sites. Following the good validation results, both using ROC and a well-known tool to validate APMs, the Kvamme's Gain, the method proposed can stand alone for future use. Some previous studies reported Gain values of 0.26 [16], 0.70 [15], 0.15 [63] and 0.92 [9], which is placing our method proposed and results obtained as being pertinent. Used in combination with other statistical methods, like AHP [9], fuzzy logic [21] and MaxEnt [19], better predictive models can be produced in the future.

**Figure 4.** (**a**) The final APM generated with the help of FR; (**b**) Success and prediction rate curves with associated area under curve (AUC) values.

Another significant use of the APM is in the preparation of future infrastructure projects (building roads and motorways). As highlighted by [39], the road network in this area has increased significantly from 0.67 km/km2 in 1894 to 2.64 km/km<sup>2</sup> in 2012. The future A8 motorway [64], which is projected to intersect the southern part of the study area, will be one of the biggest infrastructure projects from the north-eastern part of the country (Figure 5). Currently, the projected path of the motorway does not intersect any known Eneolithic sites; however, it overlaps with areas identified as having a high and very probability of "hosting" archaeological sites.

**Figure 5.** The A8 Motorway overlapping the APM.

#### **6. Conclusions**

Over recent years, APM has become a powerful tool in searching new locations of archaeological interest. There are very few studies (only one) employing statistical modelling of archaeological site locations in Romania. This study aims to fill that gap and to propose a new method—FR to APM. To do this, a set of three factors were used, i.e., soils, HLI and SPI. The statistical modelling procedure involved the selection of 80% of the sites to build the model and 20% of the sites were used to test the model's performance. The final APM was divided into four probability classes: low, medium, high and very high. With a success rate of 72%, a prediction rate of 75% and a value of 0.56 for Kvamme's Gain, FR proves to be a reliable method in determining areas with a high archaeological potential. Final maps can be used in preventive archaeology in studies that deal with the degradation of archaeological sites from natural hazards, cultural heritage management, understanding prehistoric people mobility and behaviour and finding sustainable development measures for cultural heritage.

**Author Contributions:** Conceptualisation, I.C.N. and A.M.-P.; methodology, I.C.N., A.M.-P. and J.W.; software, I.C.N. and A.M.-P.; formal analysis, I.C.N., A.M.-P. and J.W.; investigation, I.C.N., A.M.-P. and J.W.; resources, I.C.N. and A.M.-P.; data curation, I.C.N. and A.M.-P.; writing—original draft preparation, I.C.N., A.M.-P. and J.W.; writing—review and editing, I.C.N., A.M.-P. and J.W.; visualisation, I.C.N., A.M.-P. and J.W.; supervision, I.C.N. and A.M.-P.

**Funding:** This research received no external funding. The APC was funded by NIKU.

**Acknowledgments:** The authors would like to thank to the employees of the Romanian Waters Agency Bucharest, Prut-Bîrlad Water Administration Ia¸si, who kindly provided the LiDAR data used in this study. Two anonymous reviewers and the editor are kindly acknowledged for their fruitful comments.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Sustainability* Editorial Office E-mail: sustainability@mdpi.com www.mdpi.com/journal/sustainability

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18