GIS-based Landform Classification of Eneoli thic Archaeological Sites in the Plateau-plain Transition Zone (NE Romania): Habitation Practices vs. Flood Hazard Perception

The landforms of the Earth’s surface ranging from large-scale features to local topography are factors that influence human behavior in terms of habitation practices. The ability to extract geomorphological settings using geoinformatic techniques is an important aspect of any environmental analysis and archaeological landscape approach. Morphological data derived from DEMs with high accuracies (e.g., LiDAR data), can provide valuable information related to landscape modelling and landform classification processes. This study applies the first landform classification and flood hazard vulnerability of 730 Eneolithic (ca. 5000–3500 BCE) settlement locations within the plateau-plain transition zone of NE Romania. The classification was done using the SD (standard deviation) of TPI (Topographic Position Index) for the mean elevation (DEV) around each archaeological site, and HEC-RAS flood hazard pattern generated for 0.1% (1000 year) discharge insurance. The results indicate that prehistoric communities preferred to place their settlements for defensive purposes on hilltops, or in the close proximity of a steep slope. Based on flood hazard pattern, 8.2% out of the total sites had been placed in highly vulnerable areas. The results indicate an eco-cultural niche connected with habitation practices and flood hazard perception during the Eneolithic period in the plateau-plain transition zone of NE Romania and contribute to archaeological predictive modelling.


Introduction
The ability to describe the geomorphological setting based on GIS-landform classification is an important aspect of any environmental analysis or landscape modelling effort [1,2]. Landforms are defined as specific geomorphic features on the Earth's surface, ranging from large-scale features (e.g., plains, mountain ranges) to small-scale features (e.g., individual hills, valleys), as well as their component landforms, such as hilltops, valley bottoms, exposed ridges, flat plains, and upper or lower slopes [1][2][3]. In recent decades, the development of GIS software and free access to datasets has attracted the interest of researchers in implementing new computer algorithms to approach the morphometric attributes and topography of the Earth's surface [2,4]. Also, there has been an upward trend in the use of GIS-based analysis to classify landforms over various scientific fields such as geo-pedology [5][6][7][8], geomorphology and seafloor mapping [9][10][11][12][13][14][15], hydrology [16][17][18], climatology [19], landscape mapping and ecology [20,21], and archaeology [4,[22][23][24]. The morphological features provide useful information for archaeological studies because the data can be interlinked with settlement distribution [25], and evolution for reconstructing the paleo-landscapes [26,27]. In this context, one of the first approaches to find application in archaeological studies was the analysis of convexity, concavity, and flatness of the topographic surface [28][29][30][31]. Despite this, at present only a few studies with applications in the field of archaeology based on GIS-landform classifications have been developed [4,24,28,29]. In recent years, there have been intensive applications of LiDAR data in the analysis of floods in a GIS environment [32,33], yet no applications towards prehistoric flood hazards.
Understanding the connections between the small-scale features, large-scale landforms, flood hazard perception, and the types of archaeological settlement is an important method applied in the study of the prehistoric peoples because the landscape can reveal insights into settlement distribution and dynamics over time [4,27]. This paper provides the first landform classification of 730 Eneolithic sites, using the TPI (Topographic Position Index) [3,5,28,55], and the SD (standard deviation) of the mean elevation, abbreviated as DEV by [29], around archaeological sites [3,4,29], which can classify the landscape in terms of slope position and landform categories and morphological classes based on the geomorphology [1,4,56,57]. The results can provide insights into factors favoring human habitation during the Eneolithic period in the plateau-plain transition zone of NE Romania and contribute to archaeological predictive modelling at regional-scale based on small-scale morphological features and flood hazard patterns [45,49,[58][59][60][61].

Regional Setting
The study area (8789 km 2 ) is located in the north-eastern part Romania, between the Siret floodplain in the west and the Prut floodplain in the north and east, where the Prut River is a natural border between Romania and the Republic of Moldova to the east and northeast, and Romania and Ukraine to the north [45,46] (Figure 1a). The southern limit is a structural one, represented by the contact area between the Moldavian Plain and the Central Moldavian Plateau [42]. The Moldavian Plain, known as the Jijia Plain, represents 88.65% (7880 km 2 ) of the study area, and the Suceava Hills (Suceava Plateau) on the western flank, represent the remaining 11.35% (909 km 2 ) (Figure 1b). The elevation range between 20 m a.s.l. and 591 m a.s.l., with an average of 163.5 m a.s.l. (Figure 1c). The relief energy does not exceed 150-200 m/km 2 , with an average of 63.4 m/km 2 , where the highest values indicate the contact area between plateau-plain transition zone and the smallest values indicate the floodplain of the Siret, Jijia and Prut rivers (Figure 1d). The slope average ranges between 0° and 38.88°, with an average value of 4.74° (Figure 1e). The highest declivity corresponds to the fronts of the cuestas, generally with an eastern or north-eastern slope aspect, and the gentle slope with the backs of the cuestas with southern or south-western aspects (Figure 1f,g).
The general morpho-structural setting consists of a monocline, dominated by cuesta landforms and deeply incised valleys, where the strata are gently dipping from northwest to southeast (Miocene-Pleistocene deposits) [42]. In the Suceava Hills (western flank) and the Moldavian Plain, the lithology is characterized by successions of clays with sands (200-300 m thick), and thin layers (2-30 m thick) of limestone and sandstone (Lower and Medium Sarmatian deposits) [42]. Over these, a loess layer lies over the entire study area, with thicknesses frequently less than 2 m [62], but which can reach a thickness of 15-30 m on the cuesta slopes and on fluvial terraces [48]. The gravel deposits occur in the alluvial plains of Siret and Prut rivers [45,46]. The landscape dominated by the cuesta landforms produces two different types of slope: (i) cuesta dip slopes characterized by a low roughness; (ii) cuesta scarp slopes, generally affected by deep stream incision at the base, diffuse and well-defined gully erosion along the slopes, and landslides [43,44,48]. Generally, this typical morpho-structure along with the headwaters and local ridges in the valley of the Baseu, Jijia, and Bahlui rivers, are the main small-scale landforms used by prehistoric populations for the placement of settlements in this region [35][36][37][38]47,49]. Regarding flood hazards in NE Romania, the temperate continental climate with heavy rains creates favorable conditions for extreme floods on the main watercourses, especially along the Siret and Prut rivers. This phenomenon also occurs in the secondary rivers, but with a lower frequency due to recently created ponds. In recent decades, there have been many historical flood events both at a regional and local level which exceeded the 0.1% (1000 year) discharge insurance [45,46]. For the chronological framework of this study, three weather intervals based on flood events reconstruction have been identified by [63]: ca. 5300 BCE, ca. 4000 BCE, and ca. 3500-3000 BCE. This evidence indicates a hydrological activity more or less similar to the nowadays with 600-640 mm annual rainfall, but between these wet intervals, the climate was probably drier than in the present [48]. However, the flood hazard has always been present near to the main watercourses in the study area.

Inventory of Archaeological Sites
To identify the relationship between settlements placement and morphological features, a geo-referenced database was created using Esri ArcGIS 10.3, based on field surveys ( Figure 2), along with relevant archaeological documentation and registries. It should be noted that the analysis was based only on certain settlements, where archaeological documentation was well-grounded in reports, radiocarbon-based chronology, scientific articles, and geo-archaeological maps [34][35][36][37][38][39][40][41]

Elevation Data
The elevation model based on high-density airborne LiDAR (Light Detection and Ranging) data used to analyses the area within the buffer zones (1000 m radius) around each archaeological site was achieved by spatially processing in the ArcGIS software of 730 .tiff files. The rasters were generated in grid formats through Inverse Distance Weighting (IDW) interpolation, with cell sizes of 1 m used for large-scale analysis, and 25 m used for small-scale representations of the study area [64][65][66] (Figure 3). The resulting small-scale DEMs were filtered using flow direction, sink and fill tools, to reduce the errors generated by merging the .tiff files [64,67]. In the large-scale maps, the 25 m DEM resolution was used, where vegetation, buildings, and other artificial structures were filtered before processing the archaeological data. The slope pattern was generated using Spatial Analyst Tools (ArcGIS), and delineation of landform units based on TPI and DEV was performed using Relief Analysis Toolbox for ArcGIS [3,5,55] (Figure 4).

Flood Hazard Data
To generate the flood hazard pattern in the area of interest, the HEC-RAS v 5.0.1 software (Hydrologic Engineering Centers-River Analysis System), which is an auxiliary module for ArcGIS 10.2, was used [45,68]. The hydrological risk assessment process comprised three steps: (i) the pre-processing step, involving the generation of the thematic layers (thalweg, banks, flow paths, and cross-section vector) based on DEM with a resolution of 0.5m/pixel; (ii) the processing step, involving the export of the thematic layers to the HEC-RAS software and introducing the parameters required to run the flood simulation (Manning roughness coefficient; Flows with an insurances of 0.1%); and (iii) the post-processing step, involving the export of the HEC-RAS result into the ArcGIS software and the generation of the flood extent in order to obtain flood bands with an insurance of 0.1 % (1000 years). The large-scale flood hazard assessment was chosen instead of frequent floods with lower insurance (e.g., 1% or 100 year) because the topographic surface has been anthropogenically modified over time, and this fact disturbs the flood pattern at small-scale.

TPI and DEV
TPI is based on the algorithm developed by Weiss A.D. [3] and implemented as an extension for ESRI ArcView 3.x. by Jenness J. [55]. This tool calculates the difference between elevations at the central point 0 (Equation 1) and the average elevation (Equation 2) around it within a known radius R [3,28,29,55,69]. Positive or negative values of TPI indicate that the central point is located higher ( 0 > ) or lower ( 0 < ), respectively, than its average surroundings. In this equation, the range of TPI depends not only on distinguishes between 0 and , but also with respect to R, because a large R-value generally reveals large-scale landform units (major valleys, mountains, hills), while smaller values highlight small-scale features (stream valleys, headwaters, local ridges) [4,5,29] (Figure 4a).
In addition to the basic algorithm, DEV measures the 0 using TPI and the standard deviation SD of the elevation (Equation (3)) [29,55]. DEV improves the results because it measures the topographic position as a fraction of local relief normalized to local surface roughness (Equation (4)) [28,29]. As with TPI results, positive values of DEV ( 0 > ) indicate that the central point is situated higher than its neighbourhood, and negative values ( 0 < ) indicate that it is situated lower [3,69] ( Figure 4b).
TPI and DEV are two complementary methods frequently used in archaeological landscape research [70]. However, we consider it most appropriate to use DEV instead of TPI due to the higher potential accuracy of landform classification and the ability to identify the topographic preferences of archaeological settlements in a heterogeneous landscape [29] (Figure 4c).

Validation of Landform Classification Accuracy for Various Neighbourhood Sizes
The summaries of archaeological site placement classified into six slope position classes using Criterion (1) for five candidate radii are shown in Table 1 and Figure 5. The summaries of archaeological site placement classified into 10 landform classes for the four combined versions of small-TPI and large-TPI neighbourhood sizes according to Criterion (2) are shown in Table 2 and Figure 6. were extracted using a 1000 m buffer zone around an Eneolithic site selected randomly from the study area: (a) calculation of TPI rasters for 100 m, 300 m, 600m, 1000 m, and 2000 m thresholds using the algorithm developed by [3] and [55]; (b) calculation of standardized TPI for each threshold rasters based on SD and Mean after the ArcGIS algorithm described by [75]; (c) generate the DEV models for each threshold rasters based on standardized TPI after [29]; (d) classification of landscape features into six slope position classes using the DEV and slope for each threshold rasters (Method 1) after [3];   The accuracy of the results has been verified using the visual interpretation of aerial imagery and by comparing the landform classification generated by TPI with the specific morphological features of over 100 settlement locations provided by previous geomorphological and archaeological surveys in the study area [39,[47][48][49][50][51][52][53][54]. Also, according to [28,29], which applied the same methodology of landscape classification within a heterogenous landscape in Belgium, in this case, the same DEV thresholds and various neighbourhood sizes were used; the results obtained by [28,29] are applicable for our study area, being a plateau-plain transition zone.   Figure 4d). Open slopes (>5°) −SD ≤ Z0 ≤ SD 0 ≤ Z0 ≤ SD  (Figure 6). The statistical analysis and geoarchaeological interpretation of archaeological site placement per landform classes were achieved based on these two results.

Classification of Archaeological Site Placement Based on Slope Position
According to the slope position classification generated using SD of TPI (DEV) and R = 300 m, over 65% of settlements were placed on the convex landforms: 470 sites on ridge, summit, or hilltops; 37 sites on upper slopes; 58 sites on middle slopes (>5°); and 28 sites (≤5°) on flat areas. Thirty percent of the settlements were placed on concave features: lower slope, foot slope, 32 sites on toe slopes; and 240 sites in valleys (Figure 7). In the Precucuteni period (ca. 5000-4600 BCE), 48.3% of settlements were located on the top of the hills, 13.3% on upper and middle slopes and the remaining 38.4% of sites preferred flat areas or lower landforms like foot slopes and valleys ( Figure  7a). During the Cucuteni period (ca. 4600-3500 BCE), regardless of the cultural phases (Cucuteni A, A-B, or B), the location of settlements seems to follow the same slope position pattern: an average of 54.4% on ridges, summits, and tops 4.2% on upper slopes; 6.7% on middle slopes (>5°); 3.4% on flat areas (≤5°); 3.7% on lower slopes (foot slope, toe slope); and 27.5% in valleys (Figure 7b-e). Overall, the preference of the Eneolithic communities for placing their settlements on the top of cuesta near the steep slope is determined by the necessity to provide defence for at least one or two sides of the settlement. This is the first criterion used by prehistoric communities in the selection of habitation locations based on local topography (Figure 7).

Classification of Archaeological Site Placement Based on Landform Units
Based on TPI-landform classification using 300 m and 1000 m combined neighbourhood sizes, 59.5% of sites are located on convex landforms: 348 sites on hilltops, high ridges; 24 sites on middle slope ridges, small hills in the plains; 96 sites on local ridges/hills in the valley; and 47 sites on upper slopes. 1.7% of the sites are located in the flat areas or on the gentle slope surfaces: plains, flat areas (<5°)-7 sites; open slopes (>5°)-7 sites (Figure 8). The remaining of 38.8% of sites overlie concave landforms: 211 sites on deeply incised streams; 6 sites on middle slope drainages, small hills in the plains; 25 sites on upland drainages, headwaters; and 94 sites in U-shaped valleys. The high ridges and hills remained the main landform classes used by prehistoric communities for the location of settlements (average 39%), regardless of the cultural period (Figure 8a-e). The least represented classes are the flat or gentle slope areas (average 0.7%), most likely due to the wetlands, which occupied the flood plain of main rivers and are not suitable for habitation but important for hunting and fishing. Of the concave landform classes, the deeply incised streams (average 24.22%) and U-shaped valleys (average 12.12%) are the most represented, confirming that the second method of defending settlements was to be located in the vicinity of natural channels like stream meanders or steep banks. For the same reason, in 10.9% of cases, the sites are located on the local ridges or small hills in the valley (Figure 8).

Habitation Practices During the Eneolithic Period
This work has quantified the landform variations of Precucuteni (PC) and Cucuteni (CA, CA-B, CB, and CU) settlement locations in the landscape between Siret and Prut rivers (Moldavian Plain, NE Romania). A general trend is observed throughout the entire Eneolithic period, when the prehistoric communities preferred to place the settlements on the hilltops or in the culmination area.  In this way, the front of cuestas provides a good defence of the settlement, offers a wider perspective on the landscape, and a very good (inter)visibility [47][48][49][50][51][52][53][54]. The sites located in the valleys have a lower relative frequency, and it is possible to correlate their placement with seasonal mobility for the practice of agriculture [35][36][37]75].
The increasing number of sites located on the top of the hills since the end of the Precucuteni period is connected with demographic growth and with a slight climate change which occurred at the beginning of the Cucuteni A period. Compared to the previous Eneolithic phase, the dry period during the Cucuteni caused a decline of subsistence agriculture, forcing human communities to migrate towards the northern part of the Moldavian Plain in search of new fertile lands, pastures and resources (see Cucuteni A-B settlement distribution) [35,75]. Even though in Cucuteni B some of the communities return to the old territories, the slope position pattern remains the same. Due to the heterogeneous landscape which characterizes the area between Siret and Prut rivers, the habitat preferences and the cultural practice, especially those related to subsistence agriculture, did not change significantly during the Eneolithic. From this point of view, the results of this approach highlight the eco-cultural niche occupied by settlement in the Precucuteni and Cucuteni periods in the study area.
Regarding the accuracy of landform classification into discrete slope position classes using the DEV, and by combining the parameters from two neighbourhood sizes using large-TPI and small-TPI, the results indicate a high agreement with the archaeological surveys and landscape description achieved within the study area [34][35][36][37][38]. Even this approach could not replace the classic geomorphological expert opinion, it does bring new insights in remote sensing applied in cultural heritage assessment, preventive archaeology and archaeological predictive modelling.

Flood Hazard Perception During the Eneolithic Period
The placement of flooding areas likely indicates that the territories occupied by wetlands were used by prehistoric communities for hunting and fishing, but they are also emphasized as being a very inappropriate habitation place due to floods. According to HEC-RAS flood hazard pattern provided using 1,000-year discharge insurance [45,46,49], only 8.2% of the total sites were placed in vulnerable areas (Figure 9). During the Precucuteni and Cucuteni A phases, the settlements potentially affected by high flood events do not exceed 6.5% (Figure 9a,b), in Cucuteni A-B, this value increases to 9.9% (Figure 9c), and during Cucuteni B, the number of vulnerable settlements reaches 12.3% (Figure 9d,e). According to slope position classification generated using SD of TPI (DEV), the most potentially affected settlements are those that were located on the low areas like: valleys (37 sites), lower/middle slopes (12 sites), and flat areas with ≤ 5º slope (13 sites) ( Table 3). According to the TPI-landform classification by combining two neighbourhood sizes, the most potentially affected settlements by floods were built on the concave landforms: deeply incised streams (37 sites), U-shaped valleys (25), and plains with ≤ 5° slope (3 sites) ( Table 4). These sites may correspond to temporary settlements used for a particular activity (e.g., fishing, hunting, clay exploitation, flint processing), and could indicate the seasonal mobility of prehistoric communities generated by the annual hydrological regimes [34][35][36][37][38]41,75].
Overall, the habitation practices deduced from settlement landform patterns attest that prehistoric communities had a high awareness of flood hazard, especially during the Precucuteni and Cucuteni A periods (ca. 5000-4100 BCE). The relative increase of the number of sites placed in the vulnerable areas in the second half of the analysed period (ca. 4100-3500 BCE) can be explained by the migration of inhabited territory to the north and north-east (Cucuteni A-B phase; ca. 4100-3850 BCE) near to the Prut floodplain, the river with the highest hydrological activity in the lowland region. Also, the Cucuteni B phase (ca. 3850-3500 BCE) corresponded with the increase of the prehistoric population, and implicitly, with the territorial expansion of habitation places during the drier periods [48,63]. However, the low-presence of Eneolithic sites in floodplain of the main rivers indicate a relative cyclicality of hydrological events associated with overflow (e.g., seasonal flooding, flash floods), and also highlights the areas affected by excess of humidity (e.g., wetlands or riparian zone), which are unfit for permanent habitation.

Conclusions
The geoarchaeological investigation of the heterogeneous landscapes of the Moldavian Plain (NE Romania) using GIS landform classification and flood hazard assessment has produced valuable information regarding the distribution of 730 Precucuteni and Cucuteni settlements during the Eneolithic period. The habitation practices and flood hazard perception results based on DEV (SD of TPI) and HEC-RAS modelling technique applied in this approach are: • According to slope position classification based on DEV 300 m, over 65% of settlements were placed on the convex landforms (e.g., ridge, summit, hill top), <5% of total settlements on flat areas with a slope ≤ 5°, and 30% of settlements were placed on concave features (e.g., valleys).

•
According to the TPI-landform classification by combining two neighbourhood sizes, in this study DEV 300 and DEV 1000, 59.5% of sites are located on positive landforms (e.g., hill tops, high ridges, small hills in plains, local ridges/hills in valley), 1.7 % sites are on the flat areas or on the gentle slope surfaces (< 5°), and 38.8% sites overlap on the negative landforms (e.g., U-shaped valleys, headwaters, shallow valley, deeply incised streams).

•
According to flood hazard pattern generated for an extent with 0.1% insurance (1000 years), 8.2% of sites are located in vulnerable areas which indicate a high flood hazard perception during the Eneolithic period.

•
The high-density settlements built on specific landforms (e.g., ridge, top of cuestas) indicate a habitation practice during the Eneolithic based on local topography and highlight a specific eco-cultural niche for the prehistoric communities in the plateau-plain transition zone of NE Romania.
Regarding the methodology applied in this approach, the GIS landform classification based on TPI and DEV combined with other morphological variables (e.g., slope) can be integrated very easily into future paleo-environmental, archaeological predictive modelling, and cultural heritage management studies. Furthermore, the difference between conventional archaeological surveys and the GIS techniques used in this work is made by rapid, low-cost and the ability to perform the analysis both at small and large scale.