Using UAV Survey, High-Density LiDAR Data and Automated Relief Analysis for Habitation Practices Characterization during the Late Bronze Age in NE Romania

Abstract

The characterization of prehistoric human behavior in terms of habitation practices using GIS cartography methods is an important aspect of any modern geoarchaeological approach. Furthermore, using unmanned aerial vehicle (UAV) surveys to identify archaeological sites with temporal resolution during the spring agro-technical works and automated mapping of the geomorphological features based on LiDAR-derived DEM can provide valuable information about the human–landscape relationships and lead to accurate archaeological and cartographic products. In this study, we applied a GIS-based landform classification method to relief characterization of 362 Late Bronze Age (LBA) settlements belonging to Noua Culture (NC) (cal. 1500/1450-1100 BCE) located in the Jijia catchment (NE Romania). For this purpose, we used an adapted version of Topographic Position Index (TPI) methodology, abbreviated DEV, which consists of: (1) application of standard deviation of TPI for the mean elevation (DEV) around each analyzed LBA site (1000 m buffer zone); (2) classification of the archaeological site’s location using six slope position classes (first method), or ten morphological classes by combining the parameters from two small-DEV and large-DEV neighborhood sizes (second method). The results indicate that the populations belonging to Noua Culture preferred to place their settlements on hilltops but close to the steep slope and on the small hills/local ridges in large valleys. From a geoarchaeological perspective, the outcomes indicate a close connection between occupied landform patterns and habitation practices during the Late Bronze Age and contribute to archaeological predictive modelling in the Jijia catchment (NE Romania).

1. Introduction

GIS-based cartography is an essential tool for any modern geoarchaeological approach [1,2] and contributes to archaeological predictive modelling [3] for better protection and management of local or regional cultural heritage [4,5]. Therefore, the GIS-based relief analysis applications well-documented by [6,7] (e.g., combinations of geomorphometric parameters, supervised and unsupervised classification, probabilistic clustering algorithms, double ternary diagram classification and object-oriented image analysis) have proven to be very effective not only in the Earth science disciplines such as geo-pedology [8,9,10,11], terrestrial geomorphology [11,12,13,14,15], seafloor mapping [16,17,18], hydrology [19,20], climatology [21], landscape mapping [22] and landscape ecology [23]; but also in geoarchaeological investigations [24,25]. This statement is supported by the fact that landform characteristics are factors which certainly influenced prehistoric human behavior in terms of habitation practices [25]. In this context, there are many geoarchaeological studies which have focused on the human–environment relationship, testifying to the significant influence of landform features on settlement patterns during the different prehistoric and early historic cultural phases [24,25,26,27,28,29], among others.

Therefore, even if the automated [24] or semi-automated [6] GIS delineation of small-scale relief features based on Light Detection and Ranging (LiDAR)-derived Digital Elevation Models (DEM’s), such as the Topographic Position Index (TPI)-derived methodology abbreviated in this study as DEV after [26], can be considered just an alternative to developing innovative digital geoarchaeological maps [4,30], this method can also solve many other archaeological issues related to human–landscape interactions, especially in areas with high densities of cultural heritage spread across heterogeneous landscapes [24,25,26]. This is the case for the plateau–plain transition zone of the Jijia catchment (NE Romania), which is characterized by an impressive number of archaeological settlements (>2000 sites) dated throughout the most representative prehistoric periods for eastern Europe (e.g., Neolithic, Chalcolithic, Bronze Age and Iron Age) [25,31,32].

In this framework, due to the fact that each prehistoric cultural phase identified in northeastern Romania left particular habitation traces in the landscape, some of them well-documented, such as Neolithic and Chalcolithic cultures [3,25,31,32], we chose to analyze the spatial patterns and geomorphological characteristics of the settlements belonging to the Late Bronze Age (LBA, chronological framework: cal. 1500/1450-1100 BCE) [33]. This selection was mainly motivated by the relatively low amount of research into LBA sites, and secondly, by the need to complete the existing archaeological database (i.e., National Archaeological Record of Romania) which so far has been lagging for the LBA, in comparison with other prehistoric periods relevant to the region. In this regard, in the present study, we provide the first TPI-based landform classification using the LiDAR-derived DEM, performed on the archaeological sites of Noua Culture (NC), the most representative communities that inhabited the Jijia catchment (NE Romania) during the Bronze Age (BA). The main purposes of this study are (a) to emphasize the role of geomorphological conditions and GIS-based detailed landform mapping as a key tool for investigating the possible influences of landscape features on the spatial evolution of NC settlements and (b) to characterize the preservation status of LBA sites in the current landscape configuration of the Jijia catchment.

2. Study Area

2.1. Geography of the Jijia Catchment

Jijia River drains a territory of 5757 km² located in the northeastern part of Romania (watershed centroid: 47°30′N/27°00′E) and is the main tributary (right-bank) of Prut River—a natural border between Romania and the Republic of Moldova [34,35]. The watershed overlaps the lower area of the Moldavian Plateau, also known as the Jijia Plain [25]. The elevation range between 18 m a.s.l. and 584 m a.s.l. (average elevation—150 m a.s.l.) (Figure 1a) induces a relief energy ranging from 0.5–50 m/km² (Jijia floodplain) to 50–235 m/km² (average relief energy—70.1 m/km²) (Figure 1b). The highest elevation values (>400 m a.s.l.) indicate the contact area between plateau–plain transition zone in the western (Moldavian Plain–Central Moldavian Plateau) and southern (Moldavian Plain–Suceava Plateau) flanks. The average slope is 5.71°, where the highest declivity (>35°) corresponds to the front of the cuestas, frequently affected by landslides, and the active riverbanks consumed by erosion [25,34,36] (Figure 1c). The climate is temperate continental with a mean annual temperature of 8–10 °C and an average annual precipitation ranging from 460 mm (<150 m a.s.l.) to 670 mm (>500 m a.s.l.) [36].

The general morpho-structure is a monocline with Miocene-Pleistocene dipping strata from north-west to south-east, from Suceava Plateau to the middle Prut Valley [36]. The geological structure consists of a succession of thin layers of limestone and sandstone (2–30 m thick; Lower and Medium Sarmatian deposits) and sands with clays layers (200–300 m thick) and thin layers (2–5 m thick) of limestones and andesitic cinerites (Upper Sarmatian deposits) [37], covered by a loess layer (Pleistocene deposits) with thicknesses between 1.0–2.5 m (most frequent) and 15–30 m (e.g., fluvial terraces, reverse cuesta slopes) [38]. The valleys corresponding to major watercourses (e.g., Jijia, Sitna, Miletin and Bahlui) are characterized by recent alluvial deposits (Holocene period), of which the middle and lower sector of the Jijia floodplain is most developed in the landscape [34,39] (Figure 1d).

The large-scale landforms which dominate the heterogenous landscape of the Jijia catchment are well-contoured valleys separated by interfluves that are most often of the cuesta type. From a geoarchaeological point of view, the cuesta landforms are the most representative, in terms of prehistoric habitation practices, because they produce two different types of slope–site relationships [25]: (1) the first type consists of the cuesta dip slopes characterized by low roughness—preferred for agro-pastoral practices; and (2) the second type is represented by cuesta scarp slopes, generally affected by deep stream incisions at the base, diffuse and well-defined gully erosion along the slopes and large landslides—preferred for settlement locations due to the dominance in the relief and for the high visibility [25,40]. Generally, this typical morpho-structure along with the open slopes and local ridges in the major floodplains are the main small-scale landforms used by prehistoric populations for placing the settlements. However, more details on the geological, geomorphological, hydrological and climatological aspects can be found in other various studies related to the geoarchaeological context of the studied area [31,37,40,41,42].

2.2. Archaeological Context: Late Bronze Age—Noua Culture

From a cultural point of view, the end of the Bronze Age (BA) in Romania is characterized by the emergence of two new important cultural complexes, namely, Zimnicea-Plovdiv and Noua–Sabatinovka–Coslogeni. The latter was documented, so far, in eastern Romania and Transylvania, and Transcarpathian Ukraine and Republic of Moldova, reaching the middle and upper Dniester [43,44,45]. The main characteristics of these human groups is represented by the presence of the so-called ashmounds (grey spots, visible on aerial photographs and on site, with diameters of 20–30 m and small elevations) that sparked interest amongst specialists since the end of the 19th and beginning of the 20th centuries [46,47] (Figure 2).

Until recently, ashmounds were mostly considered linked to possible remains of burnt dwellings or hearths [43,44,48,49,50] (Figure 2a) and sacred/cultic areas of the settlements [51,52,53,54] (Figure 2b). Lately, due to the interdisciplinary studies conducted on such structures from Republic of Moldova [55], a more plausible explanation was issued, specifically the one of household pits, whose organic content led to physicochemical changes in the soil, producing the ash-like color, representing evidence of the type of work practiced, namely, cattle shepherding. The latter is well proven by the high number of animal osteological remains (mostly belonging to cattle or sheep), some of which were processed and converted in different types of tools, present in all the settlements specific to LBA (Figure 2c). Besides the presence of ashmounds and numerous animal bones, the most important cultural markers for Noua Culture (NC) communities are the double-handled kantharoi (present especially in funerary contexts); the curved knives made of flint (krummesser); crenated scapulae, used for processing animal skin; and tupik sickles, made out of large animals’ mandibles and utilized in agriculture [48].

3. Data and Methods

3.1. Archaeological Data: LBA Sites Inventory

In order to identify the relationship between the LBA settlements and the geomorphological features of the occupied landscape, a geo-referenced database has been generated, firstly in Google Earth Pro 7.1.5., and exported in ArcGIS 10.3 as shapefile data (Table S1). The cartographic and aerial image products used for validation of sites location consisted of old maps; military topographic maps (1:25,000 scale); orthophotoplans existing on portals, such as Inis Viewer and Atlas Explorer; and LiDAR-derived DEM [56] and aerial images collected during the field work. Additionally, the LBA sites inventory was compiled by consulting the archaeological monographs [43,44,57,58] and archaeological repertories existing for the counties Iași and Botoșani [59,60,61,62], whose territories overlap, partially, the Jijia catchment. In addition to cartographic documentation, due to the fact that NC settlements are characterized, in most cases, by the presence of the so-called ashmounds, we were able to identify 70 new sites during the field surveys using UAV technology [63,64,65] (Figure 3). To this end, we used a drone (Phantom 4 Pro v.2) in order to obtain aerial photographs. Two techniques were used: oblique photography and vertical photography, using the mission planner available in the DJI Pilot application. For the vertical one, we used a 70% overlap between each photo. The flight altitude used was between 70 and 100 m. The photos were later imported in Agisoft Metashape in order to obtained large ortorectified images [64].

Overall, regardless of the inventory method used, the site mapping was done only for the discoveries that consisted, among others, of remains from dwellings of certain settlements. Therefore, we have compiled a database consisting of more than 400 certain Bronze Age sites (BA-chronological framework: ca. 2900/2800–1100 BCE), from which 362 belong to the LBA period (out of which 195 present ashmounds) of the Noua Culture (ca. 1500/1450–1100 BCE) (Figure 4a). It should be noted that, besides the LBA settlements, we also identified and mapped the funerary contexts, the hoards and other types of discoveries (isolated or uncertain) belonging to the same chronological interval, but these will be analyzed in the future. In this framework, the last step of the archaeological database construction was to generate a 1000 m buffer zone for each selected LBA site used for automated relief analysis and landform classification purposes (Figure 4b; Table S1).

3.2. Elevation Data: LiDAR-Derived DEM

The elevation data were obtained from National Administration “Romanian Waters”—Prut-Bîrlad Water Administration (NARW-PBWA) which, through the implementation of SMIS-CSNR number 17945: Works to reduce the flood risk in the Prut–Bîrlad River Basin (PBRB) (2013) [56], managed to scan the entire northeastern territory of Romania using high-density airborne LiDAR technology, a popular remote sensing method used for measuring the exact length of an object on the Earth’s surface [25,34,35]. Broadly, the LiDAR method uses light in the form of a pulsed laser to measure ranges (variable distances) of the Earth. These light pulses—combined with other data recorded by the airborne system—generate precise, three-dimensional information about surface characteristics [25]. Therefore, the elevation data used in this work consisted of more than 700 raster files generated based on raw ground point elevation data collected at spatial density between 4 point/m² (out of built-up areas) and 16 point/m² (built-up areas).

The LiDAR-derived DEM used for the relief analysis within the buffer zones (1000 m radius) around each LBA site identified in the Jijia catchment was achieved by spatially processing the raster files, and following these steps: (1) the raster’s were generated in grid formats through Inverse Distance Weighting (IDW) interpolation [66,67], with 1 m cell sizes; (2) the resulting small-scale DEMs were filtered using flow direction, sink and fill tools, to reduce the errors generated by merging the .tiff files [68,69]; (3) the slope raster was generated using Spatial Analyst Tools via ArcGIS 10.3; (4) the delineation of landform units was performed using an adapted TPI-landform classification tool from Relief Analysis Toolbox for ArcGIS abbreviated by [26] with DEV [8,24,25].

3.3. Automated Relief Analysis: TPI-Based DEV

The GIS algorithm used for landform classification of LBA site locations was based on the Topographic Position Index (TPI) tool developed by Weiss, A. D. [70] and implemented as an ESRI ArcView 3.x. extension by Jenness J. [71]. The TPI tool calculates the difference between elevations at the central point $z_0$ (Equation (1)), which in our case was the central point of each LBA site, and the average elevation $\overline{z}$ (Equation (2)) around it within a known radius R [25]: R = 100 m, R = 300 m, R = 600 m, R = 1200 m and R = 2000 m around the LBA sites selected for this study. After this process, the positive values of TPI indicate that the central point is located higher ($z_0$>$\overline{ z}$) than its average surroundings and the negative values of TPI ($z_0$< $\overline{z }$) indicate the opposite.

(1)$TPI=z_0-\overline{z}$

(2)$\overline{z}={\frac{1}{n}}_R{{\displaystyle \sum }}_{i\in R}{z}_i$

Based on the initial TPI algorithm, De Reu J. and his collaborators [72] developed a new computational equation (abbreviated DEV) which uses the TPI calculation and the standard deviation (SD) of the surrounding elevation of $z_0$ (Equation (3)). According to the authors [72], DEV improves the results because it measures the TPI as a fraction of local relief normalized to local surface roughness (Equation (4)) [25,26]. However, even if both TPI and DEV tools are frequently used in the landscape archaeology studies [24,25,26,72], we considered it most appropriate to use DEV instead of TPI due to the higher potential accuracy of landform classification in the study area, in the catchment area of Jijia River.

(3)$DEV=\frac{z_0-\overline{z}}{SD}$

(4)$DEV=\sqrt{\frac{1}{n_R-1}{{\displaystyle \sum }}_{i=1}{\left(z_i-\overline{z}\right)}^2}$

The next step was to divide the landscape into the six discrete slope position classes described in Table 1 (first method) or the ten morphological classes described in Table 2 (second method), both methods using the GIS-based methodology adapted to TPI-derived DEV [70,71]. Therefore, the TPI-based DEV data processing involves classification of the topographic surface into a complex landscape feature by combining the parameters from two small-DEV (e.g., R = 300 m) and large-DEV (e.g., R = 1200 m) neighborhood sizes by using the TPI-landforms class tool from the Relief Analysis Toolbox [8]. As this method has been applied before in a previous study for another prehistoric interval, for the Chalcolithic cultures Precucuteni and Cucuteni [25], we already know that the TPI-based DEV results highlight all landform types [70,71] occupied by the prehistoric settlements in the heterogenous landscape of northeastern Romania. The landform classification accuracy has been verified based on two methods: first, by visual interpretation of morphological limits using the UAV imagery for large-scale landforms (e.g., valleys and hills) and the LiDAR database for small-scale morphological features (e.g., deeply incised streams and small hills in plains); secondly, by comparing the TPI-based DEV results with the specific geomorphological features of LBA site locations identified and characterized during the successive field surveys.

4. Results

4.1. Landform Classification Accuracy and Optimal Neighborhood Sizes Combination

The descriptive statistics of LBA archaeological sites placement classified into six slope position classes, using the first method (see Table 1) for four candidate radii (R = 100 m; R = 300 m; R = 600 m; R = 1200 m), are shown in Table 3 and Figure 5. The descriptive statistics of LBA archaeological sites placement, classified into ten landform classes (second method, see Table 2) for the four combined versions of small TPI-based DEV and large TPI-based DEV neighborhood sizes (DEV 100 and DEV 300 m; DEV 300 and DEV 1200 m; DEV 300 and DEV 2000 m; DEV 600 and DEV 2000 m), are shown in Table 4 and Figure 6.

The TPI-based DEV landform classification accuracy has been assessed using the visual interpretation of UAV aerial imagery collected during the field work and by comparing the automated relief analysis results of LBA sites’ locations with the specific morphological features of more than 200 prehistoric site locations provided by previous geomorphological and archaeological studies [31,32,37]. Furthermore, based on a previous study [25], in which we applied a similar automated relief classification within the heterogenous landscape of NE Romania (plateau–plain transition zone of Moldavian Plain), we were able to use, in this case also, the same TPI-based DEV thresholds and various neighborhood sizes.

The optimal neighborhood size combinations identified in the previous study [25], which are also applicable in this work, indicate that for the landform classification in six slope position classes (first method), the 300 m radius value around z₀ are the most appropriate for the rest of our approach (Figure 5). We support this statement due to the fact that R = 300 m has a low density of patches and discriminates the various features with less fragmentation, in comparison with R = 100 m, and also does not generalize the relief features such as R = 600 and R = 1200 m (Table 3). For the second relief analysis method in which we classified the landscape into ten morphological classes, the results of small (R = 300 m) and large (R = 1200 m) combined neighborhood sizes highlight the most accurate, dominant landform types occupied by LBA settlements (Figure 6). Other neighborhood size combinations emphasize the U-shaped valleys and discriminate the deeply incised streams (e.g., small-R = 100 m and large-R = 600 m) or upper slopes and headwaters classes (e.g., small-R = 300 m and large-R = 2000 m), or as in cases of small-R = 600 m and large- R = 2000 m, reduce the classes of landform features from ten landform classes to just two dominant relief features, hilltops/high ridges and deeply incised streams (Table 4). Therefore, the geoarchaeological interpretation of LBA site placement per landform classes was achieved based only on this optimal neighborhood size combination.

4.2. Slope Position Classification of LBA Sites’ Locations

According to the slope position classification resulting from R = 300 m of TPI-based DEV and LiDAR-derived slope combination, 38.12% of LBA settlements were placed on the concave landforms, 115 sites in valleys (Sp1) and 23 sites on lower slopes (Sp2); 7.45% of LBA sites were located on the flat areas (Sp3) (27 sites); and over 54% of the LBA settlements were placed on the convex landforms (44 sites on the middle and upper slopes—Sp4/Sp5; 153 sites on the ridges—Sp6). The preference manifested by the LBA communities for placing their settlements on the top of cuesta (42.2% sites located on the ridges—Sp6) can be interpreted, firstly, as a necessity to provide defense for at least one or two sides of the settlement, as in the case of Cucuteni Culture [25]; and secondly, as the necessity to have a wide perspective on the valleys. In this context, the settlements which occupied the concave landforms, such as valleys (Sp1) and lower slopes (Sp2) (>140 LBA sites), can be interpreted as agro-pastoral locations during the vegetation season. This affirmation is supported by the fact that the LBA sites which have in their structures the so-called ashmounds (195 sites) occupied a relatively lower number of top/summits/ridges landforms, compared with the ones that do not present the mentioned structures (167 sites). However, the locations of LBA sites on the tops of the hills but near the steep slopes and/or inside of the large valleys indicate human behavior in close connection with agro-pastoral activities.

4.3. Landform Classification of LBA Sites’ Locations

According to the TPI-based DEV landform classification using small-DEV 300 m and large-DEV 1200 m combined neighborhood sizes, 47.8% of LBA sites are located on convex landforms: 77 on hilltops, high ridges (L10); 4 sites on middle slope ridges, small hills on plains (L9); 72 on local ridges/hills in valleys (L8); and 20 on upper slopes (L7). Next, 1.38% of the LBA settlements were located on the flat or gentle slope areas (<5° plains, flat areas (L5)—2 sites; >5° open slopes (L6)—3 sites). The remaining 50.82% of LBA sites were located on the concave landforms: 69 sites in U-shaped valleys (L4); 20 on upland drainages/headwaters (L3); 4 sites on middle slope drainages/shallow valleys (L2); and 91 on deeply incised streams (L1). The most represented landform classes are: deeply incised streams (L1)—25.14%; hill tops and high ridges (L10)—21.27%; local ridges/hills in valleys (L8)—19.88%; and U-shaped valleys (L4)—19.06%. The least represented landform classes are the flat or gentle slope areas, most likely due to the wetlands which occupied the floodplains of the Jijia River and its tributaries (e.g., Bahlui, Miletin, Bașeu), areas not suitable for habitation but very important for pastoral activities. The different slope position classification outcomes, presented in the previous subchapter, based on the TPI-based DEV landform classification using small-DEV 300 m and large-DEV 1200 m, reveal that half of the LBA sites located on the ridges (Sp6) are actually located on the local ridges/hills in the valley (L8). These results are also highlighted by the high number of the LBA sites with ashmounds (>50 sites), which occur on the local ridges/hills in valley (L8). Therefore, the second method of landform selection for settlements’ placements during the LBA was finding the locations in the vicinity of natural channels, such as stream meanders or steep banks inside the major floodplains, or on small hills.

5. Discussion

5.1. Characterization of Habitation Practices during the Late Bronze Age

In this study, an alternative solution for automated relief analysis based on slope positions and landforms patterns of the location of Noua Culture settlements, in the landscape of the Jijia catchment (NE Romania), is presented [24,25,70,71,72]. The outcomes related to the relationship between the archaeological site locations and the geomorphological features indicate four specific landforms classes preferred by the prehistoric communities throughout the entire Late Bronze Age period: L10—hilltops, high ridges, L8—local ridges/hills in valleys, L4—U-shaped valleys and L1—deeply incised streams. Therefore, the settlements located on the hilltops and ridges, and the settlements located on the local ridges and hills in valleys, indicate the necessity to provide defense for at least one or two sides of the LBA sites and to gain better views of the valleys. These habitation practices characteristics, also identified at other prehistoric communities that occupied the Jijia River basin earlier [31,32], indicate that this was the first criterion used in selecting permanent site locations based on local topography [25]. However, what is specific to Noua Culture and different than other prehistoric cultures in the region is the high number of settlements located inside of U-shaped valleys of the Jijia River and its major tributaries (Bahlui, Miletin and Sitna rivers). For example, during the Late Bronze Age, more than 24% of the sites were identified in the floodplain areas of Jijia River, compared with Precucuteni–Cucuteni Culture which occupied the study area between cal. 5000–3500 BCE, where only 8% of the settlements were located within floodplain areas [25]. These significant differences in habitation characteristics can be explained, firstly, by climatic conditions (the climate changed during the mid-Holocene—a drier period towards the end of Bronze Age period which caused a probable reduction in flood events [73,74,75,76,77]), and secondly, due to the necessity of LBA communities for wide river plains areas that could support the work practiced, namely, cattle shepherding [24,43,44,45,46,47,48,49,50,51,52,53,54,55] (Table 5).

Another particularity of the LBA sites is the presence of the ashmounds (see Figure 3), which were observed in 195 sites out of the 362 sites investigated in this study. In this study, we wanted to find out if there were significant differences regarding landform selection between the two different types of LBA sites (see Figure 3). To investigate this issue, we compared the results of the slope position (Figure 7a,b) and landform classification (Figure 7c,d) between the LBA settlements with ashmounds and the ones without.

According to the slope position classification results, there are not significant differences between the two categories, even if the flat areas (Sp3) and lower slopes (Sp2) classes are relatively more representative for LBA sites without visible ashmounds (Table 6). According to the landform classification based on combination of small-DEV and large-DEV, several differences can be observed in U-shaped valleys (L4) and local ridges/hills in valley (L8) classes (Table 7). There are 15% more LBA sites with ashmounds located on the local ridges/hills in valleys (L8) than LBA sites without these structures; and there are 9% more sites without ashmounds located in U-shaped valleys (L4) landforms then LBA sites with ashmounds. The most plausible explanation for these outcomes is the ashmound’s structure and how they appeared in the LBA settlements as a result of cattle shepherding activities [76]. Therefore, the reason why there are more LBA sites with ashmounds located on the local ridges/hills in valley (L8) than in U-shaped valleys (L4) is the fact that the wide and open valleys are good for grazing but unsuitable for habitation due to weather exposure and floods events [25].

Although, usually, the individuals that inhabited the Jijia catchment during the end of the Bronze age are considered to have been mere shepherds, concerned only with performing the mentioned activity, the totality of the preferences manifested by them in selecting the locations for their settlements reveals different behavior. Thus, the appetite for proximity to water sources finds its explanation in the good visibility and control gained over the river found in its vicinity, especially since rivers do not represent only a source of water supply but also a means of communication/transportation. Additionally, this characteristic allows good artefactual mobility, a fact proven by the presence of a large number of LBA sites in the immediate proximity of Jijia River, the most important watercourse in the entire workspace. Bearing in mind the risks of flood events, the LBA human groups chose to be located on the highest landforms present in the Jijia Valley, and the explanation for this lies in the very high visibility they gained, thereby controlling the entire middle sector of the river.

Last but not least, we have to take into consideration some of the archaeological excavations performed on settlements with no ashmounds visible on the surface, located in areas that were not subject to agricultural processes [44,58]. In some cases, the invasive research conducted at these sites has highlighted the presence of the ash-like soil layer at a depth of 20–30 cm. The presence of the LBA cultural layer so close to the topsoil could represent evidence of the level of destruction [78,79] and an appeal to the necessity of preservation and conservation of the archaeological remains belonging to the end of the Bronze Age.

5.2. Automated Relief Analysis: Applicability and Limitations for Archaeological Studies

In this study, an alternative GIS-based solution for analyzing the spatial patterns and geomorphological characteristics of the settlements belonging to Late Bronze Age (LBA, chronological framework: cal. 1500/1450-1100 BC) was presented. To this end, we used an adapted version of Topographic Position Index (TPI) methodology [70,71], abbreviated DEV after [72], which consists of: application of the standard deviation of TPI for the mean elevation (DEV) around each analyzed LBA site (1000 m buffer zone); classification of the archaeological sites’ locations using six slope position classes (first method) and ten morphological classes by combining the parameters from two small-DEV and large-DEV neighborhood sizes (second method). The outcomes produced new and valuable information regarding the 362 archaeological sites belonging to the Noua Culture, which flourished during the Late Bronze Age (cal. 1500/1450-1100 BCE) in the heterogenous landscape of the Jijia catchment (NE Romania). However, even if the automated relief analysis using TPI-based DEV methodology could not replace geomorphological expert opinions, the results bring new insights regarding remote sensing applied in cultural heritage assessments and archaeological predictive modelling.

According to [72], the TPI, or adapted DEV, is an example of a topographic metric that became institutionalized by its integration into widely used tools, such as the ESRI products [72]. This resulted not only in the popularity of this metric parameter, but also in the uncritical application by users without solid geomorphological backgrounds [72], such as archeologists [25]. In this context, there are many examples of when the GIS-based landform classification was easily integrated into various geo-archaeological studies due to its large applicability and low-cost performance [24,25]. Overall, most studies have focused on these automated modeling techniques, especially when the investigation of large archaeological databases was needed for various purposes (e.g., cultural heritage management, habitation practices characterization during different cultural periods, impact assessments of natural and anthropogenic hazards on the archaeological sites and so on) [24,25,26,27]; but there are also examples of when the TPI algorithm was applied in small-scale situations (intra-site) when the digital elevation models with high resolution (e.g., LiDAR-derived DEM’s) were available [2,3,25]. Therefore, one of the technical limitations of the GIS-based methodology proposed in this study is the low resolution of DEM.

6. Conclusions

The automated relief investigation using the LiDAR-derived DEM and TPI-based DEV landform classification has produced new and valuable information regarding the 362 archaeological sites (195 sites with ashmounds) belonging to the Noua Culture, which flourished during the Late Bronze Age (cal. 1500/1450-1100 BCE) in the heterogenous landscape of the Jijia catchment (NE Romania). Therefore, the main habitation practices characteristics derived from relief analysis of the LBA sites’ locations are:

• According to the slope position classification (first method) resulting from R = 300 m of TPI-based DEV and LiDAR-derived slope combination, over 38% of LBA settlements were placed on the concave landforms (Sp1—valleys and Sp2—lower slope); 7.45% of LBA sites were located on the Sp3—flat areas (≤5°); and over 54% of LBA settlements were placed on the convex landforms (Sp4—middle slope, Sp5—upper slope and Sp6—ridges, tops of the hills).

• According to the TPI-based DEV landform classification using small-DEV 300 m and large-DEV 1200 m combined neighborhood sizes, over 47% of LBA sites were located on convex landforms (L10—hilltops, high ridges, L9—middle slope ridges, small hills in the plains, L8—local ridges/hills in valleys and L7—upper slopes), 1.38% of the LBA settlements were located on the flat or gentle slope areas (L5—plains/flat areas and L6—open slopes) and 50.82% of LBA sites were located on the concave landforms (L4—U-shaped valleys, L3—upland drainages/headwaters, L2—middle slope drainages/shallow valleys and L1—deeply incised streams).

• The very high-density of NC sites located just on the four specific landforms, of which two are convex landforms (L10—hilltops, high ridges and L8—local ridges/hills in valleys) and two are concave landforms (L4—U-shaped valleys and L1—deeply incised streams), indicates habitation practices based on agro-pastoral activities manly induced by the suitability of the local topography. Additionally, the landform patterns highlight a specific eco-cultural niche just for the LBA communities, different from those of the other prehistoric cultures (e.g., Precucuteni-Cucuteni Culture) which flourished in the same workspace.

From the perspective of the GIS-based methodology applied in this study, the automated relief analysis using TPI-based DEV combined with LiDAR-derived DEM and other geomorphological variables (e.g., terrain slope) can be integrated very easily into various geo-archaeological studies (e.g., paleo-environmental reconstructions, eco-cultural niche modelling and archaeological predictive modelling) due to its great applicability. Furthermore, there are many improvements that the GIS-based techniques used in this work bring to the conventional geo-archaeological surveys: the fast and low-cost performance of relief analysis at both small and large scales, the ability to divide the landscape into ten landform classes and replicate the analysis for various archaeological contexts and also the ability to describe prehistoric human behavior based on certain geographical datasets.

Author Contributions

Conceptualization, A.M.-P. and C.B.; methodology, A.M.-P.; software, A.M.-P., C.B. and C.C.S.; validation, A.M.-P. and C.B.; formal analysis, A.M.-P. and C.B.; investigation, A.M.-P. and C.B..; resources, A.M.-P., C.B. and C.C.S.; data curation, A.M.-P. and C.B.; writing—original draft preparation, A.M.-P. and C.B.; writing—review and editing, A.M.-P.; visualization, A.M.-P., C.B. and C.C.S.; supervision, A.M.-P. and C.B.; project administration, A.M.-P. and C.B.; funding acquisition, A.M.-P., C.B. and C.C.S. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

LiDAR-derived DEM data were provided by NARW-PBWA based on the institutional collaboration protocol between PBWA and UAIC. All information was processed in the Geoarchaeology Laboratory of Institute of Interdisciplinary Research, Science Department (Arheoinvest Centre), University “Alexandru Ioan Cuza” of Iaşi (UAIC).

Conflicts of Interest

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:

GIS	Geographic Information System
TPI	Topographic Position Index
SD	Standard Deviation
DEV	Standard Deviation of Topographic Position Index
LiDAR	Light Detection and Ranging
DEM	Digital Elevation Model
IDW	Inverse Distance Weighting
R	Radius
$\overline{z}$	Average elevation pixels for various candidate radii
$z_0$	Central pixel elevation for various candidate radii
BA	Bronze Age
LBA	Late Bronze Age
NC	NOUA Culture
NARW-PBWA	Romanian Waters–Prut-Bîrlad Water Administration
PBRB	Prut–Bîrlad River Basin
UAV	Unmanned Aerial Vehicle
Sp1	Slope position class: Valley
Sp2	Slope position class: Lower slope
Sp3	Slope position class: Flat area
Sp4	Slope position class: Middle slope
Sp5	Slope position class: Upper slope
Sp6	Slope position class: Ridge
L1	Landform class: Deeply incised streams
L2	Landform class: Middle slope drainage, shallow valley
L3	Landform class: Upland drainage, headwaters
L4	Landform class: U-shaped valleys
L5	Landform class: Plains, flat areas (<5°)
L6	Landform class: Open slopes (>5°)
L7	Landform class: Upper slopes
L8	Landform class: Local ridges/hills in valley
L9	Landform class: Middle slope ridges, small hills in plains
L10	Landform class: Hill tops, high ridges

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

Figures and Tables

Figure 1. Geographic location of Jijia river basin in northeastern Romania with (a) elevation, (b) relief energy, (c) slope and (d) geological sketch maps.

Figure 2. Examples of archaeological remains belonging to Noua Culture (LBA) in the Jijia catchment (NE Romania): (a) inside of an LBA site from the Jijia lowland—an excavated ashmound which was visible on aerial photographs and on site; (b) a tomb belonging to Noua Culture discovered during the excavation of a Chalcolithic settlement from Jijia-Siret upland; (c) LBA tools and ceramic fragments collected from various locations in the Jijia basin: c1, c2 and c3—crenated scapulae used as tools for processing animal skin; c4—a ceramic fragment of a cooking pot; c5—a stone-axe.

Figure 3. Aerial images (left—during the vegetation season; right—during the spring agro-technical works) with unexcavated LBA sites belonging to Noua Culture in the Jijia catchment: (a); Site 115—Dumești, Holmului Hill; (b) Site 125—Erbiceni, Spinoasei Hill/Valea Lungă; (c) Site 132—Fântânele, Cimitirul Ortodox; (d) Site 21—Aroneanu, Șapte Oameni; (e) Site 23—Dorobanț, La Chișcă; (f) Site 31—Boureni, Popa Mort Hill; (g) Site 36—Valea Oilor, Mădârjești Pond; (h) Site 37-Valea Oilor, North of the village; (i) Site 66—Ceplenița, Ion Clacă Hill; (j) Site 127—Spinoasa, Drumul Poștei Hill; the so-called ashmounds (circular gray spots visible on the field) with temporal resolution during the spring agro-technical works (right image of each example).

Figure 4. Noua Culture (NC) sites distribution in the Jijia River basin (see Table S1) used for habitation practices characteristics during the Late Bronze Age: (a) 362 LBA settlements, of which 167 sites are without ashmounds and 195 sites are with ashmounds (see Figure 3); (b) 1000 m buffer zone around each LBA sites used for automated relief analysis and landform classification.

Figure 5. Slope position classification based on TPI-derived DEV of the LBA archaeological sites (Noua Culture) in the Jijia catchment, with six morphological classes (see Table 1); the relief analysis was provided for four TPI-derived DEV neighborhood sizes: (a) 100 m, (b) 300 m, (c) 600 m and (d) 1200 m.

Figure 6. Landform classification based on TPI-derived DEV of the LBA archaeological sites (Noua Culture) in the Jijia catchment, with ten landform types (see Table 2); the relief analysis was provided for the combined of two TPI-derived DEV neighborhood sizes (a) 100 and 600 m, (b) 300 and 1200 m, (c) 300 and 2000 m and (d) 600 and 2000 m.

Figure 7. Automated relief analysis of LBA sites in the Jijia catchment: slope position classification based on TPI-derived DEV with six morphological classes (see Table 1) of the (a) LBA archaeological sites without ashmounds (167 sites) and (b) LBA sites with ashmounds (195 sites); landform classification based on TPI-derived DEV with ten landform types (see Table 2) of the (c) LBA archaeological sites without ashmounds (167 sites) and (d) LBA sites with ashmounds (195 sites).

Table 1

Description and abbreviation of slope position classes used for habitation practices’ characterization during the LBA in the Jijia catchment, obtained based on TPI-based DEV calculated for various neighborhood sizes and terrain slope (first method).

Slope Position ClassesDescription 1	DEV Threshold	Landform ClassesAbbreviation
Ridge	TPI > 1 SD	Sp6
Upper slope	0.5 SD < TPI ≤ 1 SD	Sp5
Middle slope	−0.5 SD < TPI ≤ 0.5 SD	Sp4
Flat area	−0.5 SD < TPI ≤ 0.5 SD	Sp3
Lower slope	−1 SD < TPI ≤ −0.5 SD	Sp2
Valley	TPI ≤ −1 SD	Sp1

¹ Slope position classes adapted after [70,71] for specific landscape characteristics of the study area.

Table 2

Descriptions and abbreviations of landform classes used for habitation practices characterization during the LBA in the Jijia catchment, obtained based on combined small-DEV and large-DEV neighborhood sizes (second method).

Landform Classes Description 1	Small-DEVNeighborhood Size	Large-DEVNeighborhood Size	Landform ClassesAbbreviation
Hill tops, high ridges	Z0 > SD	Z0 > SD	L10
Middle slope ridges, small hills in plains	Z0 > SD	0 ≤ Z0 ≤ SD	L9
Local ridges/hills in valley	Z0 > SD	Z0 < −SD	L8
Upper slopes	−SD ≤ Z0 ≤ SD	Z0 > SD	L7
Open slopes (>5°)	−SD ≤ Z0 ≤ SD	0 ≤ Z0 ≤ SD	L6
Plains, flat areas (<5°)	−SD ≤ Z0 ≤ SD	−SD ≤ Z0 < 0	L5
U-shaped valleys	−SD ≤ Z0 ≤ SD	Z0 < −SD	L4
Upland drainage, headwaters	Z0 < −SD	Z0 > SD	L3
Middle slope drainage, shallow valley	Z0 < −SD	0 ≤ Z0 ≤ SD	L2
Deeply incised streams	Z0 < −SD	−SD ≤ Z0 < 0	L1

¹ Landform classes description after [70,71] for specific landscape characteristics of the study area.

Table 3

Number of LBA archaeological sites occurring over six slope position classes in the Jijia catchment (see Figure 5).

1 Slope PositionClasses Description	DEV Threshold	LandformCode	2 R = 100 m	2 R = 300 m	2 R = 600 m	2 R = 1200 m
Sp6: Ridge	TPI > 1 SD	Sp6	86	153	126	107
Sp5: Upper slope	0.5 SD < TPI ≤ 1 SD	Sp5	42	22	7	4
Sp4: Middle slope	−0.5 SD < TPI ≤ 0.5 SD	Sp4	102	22	16	7
Sp3: Flat area	−0.5 SD < TPI ≤ 0.5 SD	Sp3	66	27	7	5
Sp2: Lower slope	−1 SD < TPI ≤ −0.5 SD	Sp2	30	23	7	6
Sp1: Valley	TPI ≤ −1 SD	Sp1	36	115	199	233

¹ Slope position classes adapted after [70,71] for specific landscape characteristics of the study area; ² R—radius value around z_0.

Table 4

Number of LBA archaeological sites occurring over ten specific landform classes in the Jijia catchment (see Figure 6).

1 Landform ClassesDescription	Small-DEVNeighborhood Size	Large-DEVNeighborhood Size	Combined 2 Small-DEV and 3 Large-DEV
			100 m and 600 m	300 m and 1200 m	300 m and 2000 m	600 m and 2000 m
L10: Hill tops, high ridges	Z0 > SD	Z0 > SD	50	77	81	85
L9: Middle slope ridges, small hills in plains	Z0 > SD	0 ≤ Z0 ≤ SD	3	4	5	5
L8: Local ridges/hills in valley	Z0 > SD	Z0 < −SD	33	72	67	36
L7: Upper slopes	−SD ≤ Z0 ≤ SD	Z0 > SD	71	20	12	8
L6: Open slopes (>5°)	−SD ≤ Z0 ≤ SD	0 ≤ Z0 ≤ SD	22	2	5	2
L5: Plains, flat areas (<5°)	−SD ≤ Z0 ≤ SD	−SD ≤ Z0 < 0	11	3	9	3
L4: U-shaped valleys	−SD ≤ Z0 ≤ SD	Z0 < −SD	136	69	68	24
L3: Upland drainage, headwaters	Z0 < −SD	Z0 > SD	5	20	14	24
L2: Middle slope drainage, shallow valley	Z0 < −SD	0 ≤ Z0 ≤ SD	1	4	3	3
L1: Deeply incised streams	Z0 < −SD	−SD ≤ Z0 < 0	30	91	98	172

¹ Landform classes description after [70,71] for specific landscape characteristics of the study area; ² R—small radius value around z₀ (100, 300 and 600 m); ³ R—large radius value around z₀ (600 m, 1200 and 2000 m).

Table 5

Differences in habitation practices characteristics between Noua Culture (Late Bronze Age) and Precucuteni–Cucuteni Culture (Chalcolithic) [25] highlighted by the archaeological sites occurring over ten specific landform classes in the Jijia catchment (see Figure 6); the main differences are underlined in bold.

1 Landform ClassesDescription	2 Small-DEVNeighborhood Size	3 Large-DEVNeighborhood Size	Relative Frequency (%)
			Noua Culture(Late Bronze Age)	Precucuteni–Cucuteni Culture (Eneolithic)
L10: Hill tops, high ridges	Z0 > SD	Z0 > SD	21.27	45.86
L9: Middle slope ridges, small hills in plains	Z0 > SD	0 ≤ Z0 ≤ SD	1.10	2.51
L8: Local ridges/hills in valley	Z0 > SD	Z0 < −SD	19.89	10.05
L7: Upper slopes	−D ≤ Z0 ≤ SD	Z0 > SD	5.52	4.92
L6: Open slopes (>5°)	−SD ≤ Z0 ≤ SD	0 ≤ Z0 ≤ SD	0.55	0.73
L5: Plains, flat areas (<5°)	−SD ≤ Z0 ≤ SD	−SD ≤ Z0 < 0	0.83	0.73
L4: U-shaped valleys	−SD ≤ Z0 ≤ SD	Z0 < −SD	19.06	7.84
L3: Upland drainage, headwaters	Z0 < −SD	Z0 > SD	5.52	2.62
L2: Middle slope drainage, shallow valley	Z0 < −SD	0 ≤ Z0 ≤ SD	1.10	0.63
L1: Deeply incised streams	Z0 < −SD	−SD ≤ Z0 < 0	25.14	24.09

¹ Landform classes’ descriptions after [70,71] for specific landscape characteristics of the study area; ² R = 300 m (small radius value around z₀); ³ R = 1000 m (large radius value around z₀).

Table 6

Descriptive statistics of LBA sites (see Figure 7a) vs. LBA sites with ashmounds (see Figure 7b) occurring over six slope position classes in the Jijia catchment.

1 Slope PositionClasses Description	2 DEV Threshold	Landform Code	Number of LBA Sites	Number of LBA Sites with Ashmounds
Sp6: Ridge	TPI > 1 SD	Sp6	60 (35.93%)	93 (47.69%)
Sp5: Upper slope	0.5 SD < TPI ≤ 1 SD	Sp5	10 (5.99%)	12 (6.15%)
Sp4: Middle slope	−0.5 SD < TPI ≤ 0.5 SD	Sp4	9 (5.39%)	13 (6.67%)
Sp3: Flat area	−0.5 SD < TPI ≤ 0.5 SD	Sp3	20 (11.98%)	7 (3.59%)
Sp2: Lower slope	−1 SD < TPI ≤ −0.5 SD	Sp2	15 (8.98%)	8 (4.1%)
Sp1: Valley	TPI ≤ −1 SD	Sp1	53 (31.74%)	62 (31.79%)

¹ Slope position classes adapted after [70,71] for specific landscape characteristics of the study area; ² R = 300 m (radius value around z₀).

Table 7

Descriptive statistics of LBA sites (see Figure 7c) vs. LBA sites with ashmounds (see Figure 7d) occurring over ten landform classes in the Jijia catchment.

1 Landform Classes Description	2 Small-DEVNeighborhood Size	3 Large-DEVNeighborhood Size	Number of LBA Sites	Number of LBA Sites with Ashmounds
L10: Hill tops, high ridges	Z0 > SD	Z0 > SD	41 (24.55%)	40 (20.51%)
L9: Middle slope ridges, small hills in plains	Z0 > SD	0 ≤ Z0 ≤ SD	2 (1.2%)	3 (1.54%)
L8: Local ridges/hills in valley	Z0 > SD	Z0 < −SD	17 (10.18%)	50 (25.65%)
L7: Upper slopes	−SD ≤ Z0 ≤ SD	Z0 > SD	6 (3.59%)	6 (3.08%)
L6: Open slopes (>5°)	−SD ≤ Z0 ≤ SD	0 ≤ Z0 ≤ SD	3 (1.8%)	2 (1.03%)
L5: Plains, flat areas (<5°)	−SD ≤ Z0 ≤ SD	−SD ≤ Z0 < 0	5 (2.99%)	4 (2.05%)
L4: U-shaped valleys	−SD ≤ Z0 ≤ SD	Z0 < −SD	40 (23.95%)	28 (14.36%)
L3: Upland drainage, headwaters	Z0 < −SD	Z0 > SD	8 (4.79%)	6 (3.08%)
L2: Middle slope drainage, shallow valley	Z0 < −SD	0 ≤ Z0 ≤ SD	1 (0.6%)	2 (1.03%)
L1: Deeply incised streams	Z0 < −SD	−SD ≤ Z0 < 0	44 (26.35%)	54 (27.69%)

¹ Landform classes description after [70,71] for specific landscape characteristics of the study area; ² R = 300 m (small radius value around z₀); ³ R = 1000 m (large radius value around z₀).

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/rs14102466/s1. Table S1: Noua Culture sites’ distribution in the Jijia River basin used in this work for landform classification of settlements’ locations and habitation practices characterization for the Late Bronze Age.

References

1. Butzer, K.W. Challenges for a cross-disciplinary geoarchaeology: The intersection between environmental history and geomorphology. Geomorphology; 2008; 101, pp. 402-411. [DOI: https://dx.doi.org/10.1016/j.geomorph.2008.07.007]

2. Verhagen, P.; Drăguț, L. Object-based landform delineation and classification from DEMs for archaeological predictive mapping. J. Archaeol. Sci.; 2012; 39, pp. 698-703. [DOI: https://dx.doi.org/10.1016/j.jas.2011.11.001]

3. Nicu, I.C.; Mihu-Pintilie, A.; Williamson, J. GIS-Based and Statistical Approaches in Archaeological Predictive Modelling (NE Romania). Sustainability; 2019; 11, 5969. [DOI: https://dx.doi.org/10.3390/su11215969]

4. Biscione, M.; Danese, M.; Masini, N. A framework for cultural heritage management and research: The Cancellara case study. J. Maps; 2018; 14, pp. 576-582. [DOI: https://dx.doi.org/10.1080/17445647.2018.1517699]

5. Gioia, D.; Bavusi, M.; Di Leo, P.; Giammatteo, T.; Schiattarella, T. Geoarchaeology and geomorphology of the Metaponto area, Ionian coastal belt, Italy. J. Maps; 2020; 16, pp. 117-125. [DOI: https://dx.doi.org/10.1080/17445647.2019.1701575]

6. Piloyan, A.; Konečný, M. Semi-Automated Classification of Landform Elements in Armenia Based on SRTM DEM using K-Means Unsupervised Classification. Quaest. Geogr.; 2017; 36, pp. 93-103. [DOI: https://dx.doi.org/10.1515/quageo-2017-0007]

7. Simensen, T.; Halvorsen, R.; Erikstad, L. Methods for landscape characterization and mapping: A systematic review. Land Use Policy; 2018; 75, pp. 557-569. [DOI: https://dx.doi.org/10.1016/j.landusepol.2018.04.022]

8. Deumlich, D.; Schmidt, R.; Sommer, M. A multiscale soil-landform relationship in the glacial-drift area based on digital terrain analysis and soil attributes. J. Plant. Nutr. Soil Sci.; 2010; 173, pp. 843-851. [DOI: https://dx.doi.org/10.1002/jpln.200900094]

9. Illés, G.; Kovács, G.; Heil, B. Comparing and evaluating digital soil mapping methods in a Hungarian forest reserve. Can. J. Soil Sci.; 2011; 91, pp. 615-626. [DOI: https://dx.doi.org/10.4141/cjss2010-007]

10. Mora-Vallejo, A.; Claessens, L.; Stoorvogel, J.; Heuvelink, G.B.M. Small scale digital soil mapping in southeastern Kenya. Catena; 2008; 76, pp. 44-53. [DOI: https://dx.doi.org/10.1016/j.catena.2008.09.008]

11. Pracilio, G.; Smettem, K.; Bennett, D.; Harper, R.; Adams, M. Site assessment of a woody crop where a shallow hardpan soil layer constrained plant growth. Plant Soil; 2006; 288, pp. 113-125. [DOI: https://dx.doi.org/10.1007/s11104-006-9098-z]

12. Dębniak, K.; Mège, D.; Gurgurewicz, J. Geomorphology of Ius Chasma, Valles Marineris, Mars. J. Maps; 2017; 13, pp. 260-269. [DOI: https://dx.doi.org/10.1080/17445647.2017.1296790]

13. Drăguț, L.; Blaschke, T. Automated classification of landform elements using object-based image analysis. Geomorphology; 2006; 81, pp. 330-344. [DOI: https://dx.doi.org/10.1016/j.geomorph.2006.04.013]

14. McGarigal, K.; Tagil, S.; Cushman, S. Surface metrics: An alternative to patch metrics for the quantification of landscape structure. Landsc. Ecol.; 2009; 24, pp. 433-450. [DOI: https://dx.doi.org/10.1007/s10980-009-9327-y]

15. Tagil, S.; Jenness, J. GIS-based automated landform classification and topographic, landcover and geologic attributes of landforms around the Yazoren Polje, Turkey. J. App. Sci.; 2008; 8, pp. 910-921. [DOI: https://dx.doi.org/10.3923/jas.2008.910.921]

16. Wilson, M.F.J.; O’Connell, B.; Brown, C.; Guinan, J.C.; Grehan, A.J. Multiscale terrain analysis of multibeam bathymetry data for habitat mapping on the continental slope. Mar. Geod.; 2007; 30, pp. 3-35. [DOI: https://dx.doi.org/10.1080/01490410701295962]

17. Wright, D.J.; Heyman, W.D. Introduction to the special issue: Marine and coastal GIS for geomorphology, habitat mapping, and marine reserves. Mar. Geod.; 2008; 31, pp. 223-230. [DOI: https://dx.doi.org/10.1080/01490410802466306]

18. Zieger, S.; Stieglitz, T.; Kininmonth, S. Mapping reef features from multibeam sonar data using multiscale morphometric analysis. Mar. Geod.; 2009; 264, pp. 209-217. [DOI: https://dx.doi.org/10.1016/j.margeo.2009.06.002]

19. Lesschen, J.P.; Kok, K.; Verburg, P.H.; Cammeraat, L.H. Identification of vulnerable areas for gully erosion under different scenarios of land abandonment in southeast Spain. Catena; 2007; 71, pp. 110-121. [DOI: https://dx.doi.org/10.1016/j.catena.2006.05.014]

20. Liu, H.; Bu, R.; Liu, J.; Leng, W.; Hu, Y.; Yang, L.; Liu, H. Predicting the wetland distributions under climate warming in the Great Xing’an Mountains, northeastern China. Ecol. Res.; 2011; 26, pp. 605-613. [DOI: https://dx.doi.org/10.1007/s11284-011-0819-2]

21. Bunn, A.; Hughes, M.; Salzer, M. Topographically modified tree-ring chronologies as a potential means to improve paleoclimate inference. Clim. Change; 2011; 105, pp. 627-634. [DOI: https://dx.doi.org/10.1007/s10584-010-0005-5]

22. Guitet, S.; Cornu, J.-F.; Brunaux, O.; Betbeder, J.; Carozza, J.-M.; Richard-Hansen, C. Landform and landscape mapping, French Guiana (South America). J. Maps; 2013; 9, pp. 325-335. [DOI: https://dx.doi.org/10.1080/17445647.2013.785371]

23. Fei, S.; Schibig, J.; Vance, M. Spatial habitat modeling of American chestnut at Mammoth Cave National Park. For. Ecol. Manag.; 2007; 252, pp. 201-207. [DOI: https://dx.doi.org/10.1016/j.foreco.2007.06.036]

24. Argyriou, A.V.; Teeuw, R.M.; Sarris, A. GIS-based landform classification of Bronze Age archaeological sites on Crete Island. PLoS ONE; 2017; 12, e0170727. [DOI: https://dx.doi.org/10.1371/journal.pone.0170727] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28222134]

25. Mihu-Pintilie, A.; Nicu, I.C. GIS-based Landform Classification of Eneolithic Archaeological Sites in the Plateau-plain Transition Zone (NE Romania): Habitation Practices vs. Flood Hazard Perception. Remote Sens.; 2019; 11, 915. [DOI: https://dx.doi.org/10.3390/rs11080915]

26. De Reu, J.; Bourgeois, J.; De Smedt, P.; Zwertvaegher, A.; Antrop, M.; Bats, M.; De Maeyer, P.; Finke, P.; Van Meirvenne, M.; Verniers, J. et al. Measuring the relative topographic position of archaeological sites in the landscape, a case study on the Bronze Age barrows in northwest Belgium. J. Archaeol. Sci.; 2011; 38, pp. 3435-3446. [DOI: https://dx.doi.org/10.1016/j.jas.2011.08.005]

27. Noviello, M.; Cafarelli, B.; Calculli, C.; Sarris, A.; Mairota, P. Investigating the distribution of archaeological sites: Multiparametric vs probability models and potentials for remote sensing data. Appl. Geogr.; 2018; 95, pp. 34-44. [DOI: https://dx.doi.org/10.1016/j.apgeog.2018.04.005]

28. Schmaltz, E.; Märker, M.; Rosner, H.-J.; Kandel, A.W. The integration of landscape processes in archaeological site prediction in the Mugello basin (Tuscany/Italy). Proceedings of the 42nd Annual Conference on Computer Applications and Quantitative Methods in Archaeology CAA; Giligny, F.; Djindjian, F.; Costa, L.; Moscati, P.; Robert, S. Paris 1 Panthéon-Sorbonne University: Paris, France, 2014; pp. 451-458.

29. Verhagen, P.; Whitley, T.G. Integrating Archaeological Theory and Predictive Modeling: A Live Report from the Scene. J. Archaeol. Method Theory; 2012; 19, pp. 49-100. [DOI: https://dx.doi.org/10.1007/s10816-011-9102-7]

30. Gioia, D.; Bavusi, M.; Di Leo, P.; Giammatteo, T.; Schiattarella, M. A geoarchaeological study of the Metaponto coastal belt, southern Italy, based on geomorphological mapping and GIS-supported classification of landforms. Geogr. Fis. E Din. Quat.; 2016; 39, pp. 137-148. [DOI: https://dx.doi.org/10.4461/GFDQ2016.39.13]

31. Asăndulesei, A. Inside a Cucuteni Settlement: Remote Sensing Techniques for Documenting an Unexplored Eneolithic Site from Northeastern Romania. Remote Sens.; 2017; 9, 41. [DOI: https://dx.doi.org/10.3390/rs9010041]

32. Brigan, R.; Weller, O. Neo-Eneolithic settlement pattern and salt exploitation in Romanian Moldavia. J. Archaeol. Sci. Rep.; 2018; 17, pp. 68-78. [DOI: https://dx.doi.org/10.1016/j.jasrep.2017.10.032]

33. Dietrich, L. Die Mittlere und Späte Bronzezeit und die Ältere Eisenzeit in Südostsiebenbürgen Aufgrund der Siedlung von Rotbav; European Verlag Dr. Rudolf Habelt GMBH: Bonn, Germany, 2014.

34. Huţanu, E.; Mihu-Pintilie, A.; Urzica, A.; Paveluc, L.E.; Stoleriu, C.C.; Grozavu, A. Using 1D HEC-RAS Modeling and LiDAR Data to Improve Flood Hazard Maps Accuracy: A Case Study from Jijia Floodplain (NE Romania). Water; 2020; 12, 1624. [DOI: https://dx.doi.org/10.3390/w12061624]

35. Stoleriu, C.C.; Urzica, A.; Mihu-Pintilie, A. Improving flood risk map accuracy using high-density LiDAR data and the HEC-RAS river analysis system: A case study from north-eastern Romania. J. Flood Risk Manag.; 2020; 13, e12572. [DOI: https://dx.doi.org/10.1111/jfr3.12572]

36. Bacauanu, V. Câmpia Moldovei. Studiu Geomorfologic; Editura Academiei Romane: Bucuresti, Romania, 1968; pp. 1-222.

37. Niculiță, M.; Mărgărint, M.C.; Santangelo, M. Archaeological evidence for Holocene landslide activity in the eastern Carpathian lowland. Quat. Int.; 2016; 415, pp. 175-189. [DOI: https://dx.doi.org/10.1016/j.quaint.2015.12.048]

38. Haase, D.; Fink, J.; Haase, G.; Ruske, R.; Pécsi, M.; Richter, H.; Altermann, M.; Jäger, K.-D. Loess in Europe—its spatial distribution based on a European loess map, scale 1:250,000. Quat. Sci. Rev.; 2007; 26, pp. 1301-1312. [DOI: https://dx.doi.org/10.1016/j.quascirev.2007.02.003]

39. Romanescu, G.; Cîmpianu, C.I.; Mihu-Pintilie, A.; Stoleriu, C.C. Historic flood events in NE Romania (post-1990). J. Maps; 2017; 13, pp. 787-798. [DOI: https://dx.doi.org/10.1080/17445647.2017.1383944]

40. Mărgărint, M.C.; Niculiță, M. Landslide Type and Pattern in Moldavian Plateau, NE Romania. Landform Dynamics and Evolution in Romania; Radoane, M.; Vespremeanu-Stroe, A. Springer Geography, Springer: Cham, Switzerland, 2017; pp. 271-304. [DOI: https://dx.doi.org/10.1007/978-3-319-32589-7_12]

41. Nicu, I.C. Is overgrazing really influencing soil erosion?. Water; 2018; 10, 1077. [DOI: https://dx.doi.org/10.3390/w10081077]

42. Nicu, I.C. Natural risk assessment and mitigation of cultural heritage sites in North-eastern Romania (Valea Oii river basin). Area; 2019; 51, pp. 142-154. [DOI: https://dx.doi.org/10.1111/area.12433]

43. Florescu, A.C. Contribuții la cunoașterea culturii Noua. Arheol. Mold.; 1964; II–III, pp. 143-216.

44. Dascălu, L. Bronzul Mijlociu și Târziu în Câmpia Moldovei; Editura Trinitas: Iași, Romania, 2007.

45. Sava, E. Așezări din Perioada târzie a Epocii Bronzului în Spațiul Pruto-Nistrean (Noua-Sabatinovka); Biblioteca Tyragetia: Chișinău, Moldova, 2014; Volume XXVI.

46. Zaretskyi, I.A. Za-metka o drevnostyakh Khar’kovskoy gub. Bogodukhovskogo uyezda, slobody Likha-chevki. Khar’kovskiy Sb.; 1888; 2, 229.

47. Gorodčov, V.A. Dnev-nik arkheologicheskikh issledovaniy v Zen’kovskom uyezde Poltavskoy guber- nii v 1906 g. Issledovaniye Bel’skogo gorodishcha. Tr. 14 Arkheologi-Cheskogo S”Yezda V Chernigove; 1911; 3, pp. 93-161.

48. Petrescu-Dîmbovița, M. Contribuții la problema sfîrșitului epocii bronzului și începutului epocii fierului în Moldova. Stud. Și Cercet. De Istor. Veche Și Arheol.; 1953; IV, pp. 443-486.

49. Motzoi-Chicideanu, I. Cenușar. Enciclopedia Arheologiei și Istoriei Vechi a României, I; Preda, C. Editura Enciclopedică: București, Romania, 1994.

50. Gerškovič, J.P. Fenomen zol’nikov belogrodovskogo tipa. Ross. Arheol.; 2004; 4, pp. 104-113.

51. Comșa, E. L’évolution des types d’habitation du territoire de la Roumanie (depuis l’énéolithique jusqu’à la fin de l’ âge du bronze). Dritter Internationaler Thrakologischer Kongress zu Ehren W. Tomascheks, 2–6 Juni 1980 Wien; Peschew, A.; Popov, D.; Jordanov, K.; Fol, A. Swjat: Sofia, Bulgaria, 1984; pp. 121-137.

52. Berezanskaja, S.S. Komarówskaja kul’tura. Arheol. Ukrajn’skoj RSR; 1985; I, pp. 428-437.

53. Arnăut, T. Spații sacre și Practici Funerare din Mileniul I a.Chr. în Arealul Carpato-Balcanic; Casa Editorial-Poligrafică Bons Offices: Chișinău, Moldova, 2014.

54. Pieniąžek, M. Kultische Landschaften in der Steppe. Zu den Anfängen sakraler Architektur im Nordpontikum. Prähistorische Z.; 2011; 86, pp. 8-30. [DOI: https://dx.doi.org/10.1515/pz.2011.002]

55. Sava, E.; Kaiser, E. Poselenie s „Zolnicami” u Acela Odaia-Miciurin, Respublica Moldova (Arheologhicesne i Estestvennonaucinie issledovaniia) / Die Siedlung mit „Aschenhügeln” Beim Dorf Odaia-Miciurin, Republik Moldova (Archäologische und Naturwissenschaftliche Untersuchungen); Biblioteca Tyragetia: Chișinău, Moldova, 2011; Volume XIX.

56. SMIS-CSNR 17945 (Water Administration Prut—Bîrlad, Romania) Works for Reducing the Flood Risk in Prut—Bîrlad Basin. Available online: http://www.romair.ro (accessed on 28 March 2022).

57. Diaconu, V. Cultura Noua în Regiunea Vestică a Moldovei. Ph.D. Thesis; Institute of Archaeology, Romanian Academy: Iași, Romania, 2014.

58. Florescu, A.C. Repertoriul culturii Noua-Coslogeni din România. Cultură și Civilizație la Dunărea de Jos; Museum of Dunarea de Jos: Călăraşi, Romania, 1991; Volume IX.

59. Chirica, V.; Tanasachi, M. Repertoriul Arheologic al Județului Iași, I; Institutul de Istorie și Arheologie ‘A. D. Xenopol’: Iași, Romania, 1984.

60. Chirica, V.; Tanasachi, M. Repertoriul Arheologic al Județului Iași, II; Institutul de Istorie și Arheologie ‘A. D. Xenopol’: Iași, Romania, 1985.

61. Păunescu, A.; Șadurschi, P.; Chirica, V. Repertoriul Arheologic al Județului Botoșani; CIMEC: București, Romania, 1976.

62. Șovan, O.L. Repertoriul Arheologic al Județului Botoșani; 2nd ed. CIMEC: Botoșani, Romania, 2016.

63. Themistocleous, K. The Use of UAVs for Cultural Heritage and Archaeology. Remote Sensing for Archaeology and Cultural Landscapes; Chini, M.; Ehlers, M.; Lakshmi, V.; Mueller, N.; Refice, A.; Rocca, F.; Skidmore, A.; Vadrevu, K. Springer: Cham, Switzerland, 2020; pp. 241-269. [DOI: https://dx.doi.org/10.1007/978-3-030-10979-0_14]

64. Muzirafuti, A.; Cascio, M.; Lanza, S. UAV Photogrammetry-based Mapping the Pocket Beach of Isola Bella, Taormina (Northeastern Sicily). Proceedings of the 2021 IEEE International Workshop on Metrology for the Sea (MetroSea 2021); Reggio Calabria, Italy, 4–6 October 2021; [DOI: https://dx.doi.org/10.1109/MetroSea52177.2021.9611617]

65. Manfreda, S.; McCabe, M.F.; Miller, P.E.; Lucas, R.; Madrigal, V.P.; Mallinis, G.; Dor, E.B.; Helman, D.; Estes, L.; Ciraolo, G. et al. On the Use of Unmanned Aerial Systems for Environmental Monitoring. Remote Sens.; 2018; 10, 641. [DOI: https://dx.doi.org/10.3390/rs10040641]

66. Zimmerman, D.; Pavlik, C.; Ruggles, A.; Armstrong, M.P. An experimental comparison of ordinary and universal Kriging and Inverse Distance Weighting. Math. Geol.; 1999; 31, pp. 375-390. [DOI: https://dx.doi.org/10.1023/A:1007586507433]

67. Lu, G.Y.; Wong, D.W. An adaptive inverse-distance weighting spatial interpolation technique. Comput. Geosci.; 2008; 34, pp. 1044-1055. [DOI: https://dx.doi.org/10.1016/j.cageo.2007.07.010]

68. Zhou, T.; Popescu, S.; Malambo, L.; Zhao, K.; Krause, K. From LiDARWaveforms to Hyper Point Clouds: A Novel Data Product to Characterize Vegetation Structure. Remote Sens.; 2018; 10, 1949. [DOI: https://dx.doi.org/10.3390/rs10121949]

69. Doneus, M. Openness as visualization technique for interpretative mapping of airborne LiDAR derived digital terrain models. Remote Sens.; 2013; 5, pp. 6427-6442. [DOI: https://dx.doi.org/10.3390/rs5126427]

70. Weiss, A.D. Topographic Position and Landforms Analysis. Proceedings of the ESRI User Conference; San Diego, CA, USA, 9–13 July 2001; pp. 227-245. Available online: http://www.jennessent.com/downloads/TPIposter-TNC_18x22.pdf (accessed on 28 March 2022).

71. Jenness, J. Topographic Position Index (tpi_jen.avx) Extension for ArcView 3.x, v. 1.2. Jenness Enterprises. Available online: http://www.jennessent.com/arcview/tpi.htm (accessed on 28 March 2022).

72. De Reu, J.; Bourgeois, J.; Bats, M.; Zwertvaegher, A.; Gelorini, V.; De Smedt, P.; Chu, W.; Antrop, M.; De Maeyer, P.; Finke, P. et al. Application of the topographic position index to heterogeneous landscapes. Geomorphology; 2013; 186, pp. 39-49. [DOI: https://dx.doi.org/10.1016/j.geomorph.2012.12.015]

73. Watanabe, S.; Hajima, T.; Sudo, K.; Nagashima, T.; Takemura, T.; Okajima, H.; Nozawa, T.; Kawase, H.; Abe, M.; Yokohata, M. et al. MIROC-ESM 2010: Model description and basic results of CMIP5-20c3m experiments. Geosci. Model Dev.; 2011; 4, pp. 845-872. [DOI: https://dx.doi.org/10.5194/gmd-4-845-2011]

74. Hijmans, R.J.; Cameron, S.E.; Parra, J.L.; Jones, P.G.; Jarvis, A. Very High Resolution Interpolated Climate Surfaces for Global Land Areas. Int. J. Climatol.; 2005; 25, pp. 1965-1978. [DOI: https://dx.doi.org/10.1002/joc.1276]

75. Ivanova, M.; De Cupere, B.; Ethier, J.; Marinova, E. Pioneer farming in Southeast Europe during the early sixth millennium BC: Climate-related adaptations in the exploitation of plants and animals. PLoS ONE; 2018; 13, e0197225. [DOI: https://dx.doi.org/10.1371/journal.pone.0197225]

76. Fick, S.E.; Hijmans, R.J. WorldClim 2: New 1km spatial resolution climate surfaces for global land areas. Int. J. Climatol.; 2017; 37, pp. 4302-4315. [DOI: https://dx.doi.org/10.1002/joc.5086]

77. Gent, P.R.; Danabasoglu, G.; Donner, L.J.; Holland, M.M.; Hunke, E.C.; Jayne, S.R.; Lawrence, D.M.; Neale, R.B.; Rasch, F.J.; Vertenstein, M. et al. The Community Climate System Model Version 4. J. Climate; 2011; 24, pp. 4973-4991. [DOI: https://dx.doi.org/10.1175/2011JCLI4083.1]

78. Serjeantson, D. Intensification of animal husbandry in the Late Bronze Age? The contribution of sheep and pigs. The Earlier Iron Age in Britain and the Near Continent; Haselgrove, C.; Pope, R. Oxbow Books: Oxford, UK, 2017; pp. 80-93. [DOI: https://dx.doi.org/10.2307/j.ctvh1dwqj.7]

79. Brașoveanu, C. Perioada târzie a Epocii Bronzului în bazinul Jijiei (România). Habitat și materialitate. Ph.D. Thesis; Alexandru Ioan Cuza University of Iași: Iași, Romania, 2021.