Integrating GIS-MCDA and Machine Learning for Geoarchaeological Habitation Suitability Modeling in the Karst Mountain Environment of Biokovo Nature Park

1. Introduction

Understanding the geographical factors that determine habitation suitability remains one of the central challenges of landscape archaeology, particularly in mountain environments where pronounced spatial heterogeneity strongly influences patterns of human occupation [1,2,3,4]. Contemporary mountain archaeology no longer considers mountains as marginal areas but as dynamic socio-environmental systems shaped by the long-term interaction between natural processes and human activities [1,3,5,6,7].

Consequently, archaeological research has shifted from documenting individual sites towards investigating the spatial organization of cultural landscapes and the processes that structured them through time [7,8,9].

This perspective is especially relevant in Mediterranean karst mountain regions, where rugged topography, limited surface water availability and pronounced environmental gradients create highly variable conditions for human occupation and land use [2,4,10,11]. In such environments, habitation patterns result from the interaction of multiple geographical characteristics rather than a single environmental factor, making their reconstruction dependent on the integration of archaeological evidence and spatial analysis [1,12,13].

Habitation suitability modelling has become an important analytical approach for investigating the relationship between archaeological landscape features and landscape characteristics [14]. Based on the premise that the spatial distribution of archaeological landscape features reflects the interaction between human behavior and environmental conditions rather than random processes, predictive modelling integrates archaeological and geographical data within a GIS environment to evaluate habitation suitability and test hypotheses concerning the factors that influenced the spatial organization of human occupation [14,15,16,17,18]. Beyond identifying areas of archaeological potential, these models provide a quantitative framework for evaluating long-term human–environment relationships and interpreting the processes that shaped cultural landscapes [18,19,20].

Biokovo Mountain provides a suitable case study for applying this approach. Its complex karst geomorphology, rugged relief, pronounced environmental gradients and the close relationship between the Adriatic coast and the mountain interior have created highly heterogeneous conditions for long-term human occupation and land use. Archaeological research has documented numerous anthropogenic landscape features, including burial mounds, communication routes, dry-stone structures, seasonal pastoral settlements, sacred sites and cave localities with evidence of human activity, indicating recurrent use of the mountain from prehistory to the modern period [5,21]. Despite this rich archaeological record, the geographical factors associated with habitation suitability across Biokovo have never been systematically evaluated, nor has an integrated spatial model been developed to examine their role in high-mountain habitation patterns.

To address this research gap, this study presents the first GIS-based habitation suitability model for Biokovo Mountain integrating archaeological landscape features with geographical variables. Unlike approaches based solely on the spatial distribution of archaeological sites, the proposed model combines multiple categories of archaeological landscape features representing different forms of past human activity with the delineation of suitable and unsuitable areas derived from landscape characteristics. This modelling approach is used to test the hypothesis that the spatial distribution of archaeological landscape features reflects long-term relationships between human activity and geographical conditions, to evaluate the relative importance of individual geographical variables in shaping habitation suitability, and to provide a new interpretative framework for understanding the organization of the high-mountain cultural landscape.

2. Materials and Methods

2.1. Study Area

Biokovo Mountain is located along the eastern Adriatic coast within the Dinaric mountain system and represents one of the most prominent karst massifs in southern Croatia (Figure 1). Extending in a northwest–southeast direction parallel to the coastline, the massif rises abruptly from sea level to its highest peak, Sveti Jure (1762 m a. s. l.), over a relatively short horizontal distance. This pronounced altitudinal gradient creates distinct environmental contrasts between the coastal zone, transitional slopes and the high-mountain plateau [22].

The mountain is composed predominantly of Mesozoic limestones and dolomites, while its present-day morphology is the result of prolonged tectonic activity and karstification [22,23]. The landscape is characterized by steep slopes, rocky ridges and numerous surface and subterranean karst landforms, including karren, dolines, karst ridges, pits and caves [22,23]. The summit plateau exhibits well-developed polygonal karst, with densely distributed dolines separated by sharp interdoline ridges and pyramidal peaks, producing exceptional microtopographic variability [22]. Water resources are scarce, consisting primarily of springs, ponds, wells and natural rock pools [24]. Biokovo is characterized by pronounced climatic variability resulting from the interaction of Mediterranean and continental influences [25]. Temperature, precipitation and vegetation change markedly with altitude, affecting pasture productivity and the seasonal suitability of different parts of the mountain for human occupation [25,26].

2.2. Methodological Workflow of Habitation Suitability Modelling

The methodological workflow consisted of seven main stages (Figure 2). First, input spatial datasets were compiled and used to derive 21 morphometric, hydrogeomorphological, climatic, and archaeological criteria. All criterion layers were harmonized to a common 5 m raster grid, and reference data were mapped and divided into training and independent validation subsets. Training data were then used for criterion diagnostics, response-curve definition, sensitivity analysis, and spatial cross-validation. Following the analysis, two GIS-MCDA (Equal Weight and AHP) and two machine-learning (Random Forest (RF) and Extreme Gradient Boosting (XGBoost)) models were subsequently developed. Finally, all models were independently validated and compared in terms of predictive performance and spatial robustness.

2.2.1. Input Data and Geoarchaeological Criterion Generation

The predictive modelling incorporated 21 geoarchaeological criteria, divided into four groups: morphometric, hydrogeomorphological, climatic, and archaeological. All criteria were generated using a digital elevation model (DEM) derived from airborne LiDAR data provided by the State Geodetic Administration of Croatia [27]. Due to computational constraints, the DEM was resampled to a common 5 m analysis grid prior to criteria generation.

Morphometric criteria

Elevation (ELV) represents the vertical position of each raster cell and ranges from approximately sea level to 1762 m a.s.l. It was included because altitude influences accessibility, climate, vegetation, snow persistence, and the seasonal suitability of mountain habitation areas [3,4,28].

Slope (SLO) expresses terrain gradient in degrees, ranging from nearly flat surfaces to slopes exceeding 80°. It directly affects accessibility, construction potential, erosion, and the availability of usable ground, with gentler terrain generally considered more favourable for habitation [29,30,31].

Aspect (ASP) describes slope orientation and influences solar radiation, wind exposure, soil moisture, vegetation, and local microclimate. It was included to account for differences in environmental conditions associated with terrain orientation [28,29,32].

Plan curvature (PLAN) describes horizontal terrain curvature and distinguishes convergent, divergent, and approximately planar surfaces. It affects runoff concentration and drainage dispersion and therefore helps identify hollows, ridges, and locally planar terrain potentially relevant to habitation [29].

Profile curvature (PROF) describes terrain curvature along the direction of maximum slope. It influences flow acceleration, erosion, sediment accumulation, and slope stability, distinguishing locally convex, concave, and linear terrain positions [29].

Sky View Factor (SVF) measures local topographic openness as the proportion of visible sky above each raster cell. It captures differences in enclosure that can influence radiative cooling, cold-air accumulation, humidity, and exposure [33].

Terrain Ruggedness Index (TRI) quantifies local elevation variability, with higher values representing rougher and more dissected terrain [34,35]. Lower ruggedness generally corresponds to smoother surfaces that are easier to access and more suitable for construction and recurrent human use [30,31].

Topographic Position Index (TPI) describes the elevation of a location relative to its surroundings. Positive values represent locally elevated positions, negative values indicate depressions, and values close to zero represent terrain with an elevation similar to its surroundings, allowing differentiation of ridges, depressions, margins, and relatively neutral terrain positions [34].

Hydrogeomorphological criteria

Drainage density (DD) represents the local concentration of drainage features and serves as a proxy for runoff concentration, terrain dissection, erosion potential, and moisture conditions [36]. In karst terrain, it should be interpreted cautiously because surface drainage is often discontinuous [37].

Convergence Index (CI) describes the tendency of surface flow to converge or diverge [38]. It differentiates terrain where runoff and sediment may accumulate from better-drained or more divergent positions and therefore provides information on moisture, erosion, and local surface stability.

Length-Slope Factor (LSF) combines slope steepness and contributing flow length as a proxy for potential erosion intensity [39]. High values represent terrain with greater runoff-related erosion potential and generally less favourable conditions for stable surfaces suitable for habitation.

Stream Power Index (SPI) estimates the potential erosive power of concentrated runoff based on slope and accumulated flow [40]. High values identify areas potentially affected by stronger runoff or gullying and are interpreted as a proxy for runoff energy rather than permanent stream activity in the Biokovo karst.

Topographic Wetness Index (TWI) combines contributing area and slope to estimate potential moisture accumulation [41]. It was included because local moisture conditions may influence water availability, vegetation productivity, and the suitability of locations for habitation or pastoral use.

Archaeological criteria

Distance to ponds (lokve) (DL) represents proximity to mapped ponds, which constitute important potential water sources in the karst mountain environment. Access to such water sources may have influenced habitation, pastoral activity, and seasonal occupation [5,21,42,43].

Distance to mountain paths/roads (DR) measures the distance to the mapped mountain path and road network and was used as an indicator of landscape accessibility. Proximity to communication routes may have influenced habitation, pastoral infrastructure, movement, and the spatial distribution of archaeological remains [44].

Least-cost distance to the coast (LCDC) represents terrain-mediated accessibility from the coast. Slope was transformed into a walking-cost surface, with walking speeds constrained between 0.5 and 6.0 km h⁻¹ and accumulated least-cost travel time was calculated from the coastline [44,45].

Density of speleological objects (DSO) represents the spatial concentration of known caves and other speleological features [46]. Cave-rich areas provide potential shelter, water-related resources, or other archaeologically relevant activity spaces, although the pattern can also reflect differences in survey intensity.

Climatic criteria

Heat Load Index (HLI) was derived from slope and aspect to represent relative topographic heat exposure [47]. It captures differences in thermal conditions that influence snow persistence, soil moisture, vegetation, and the seasonal suitability of mountain areas for human occupation.

General Wind Shelter Index (GWSI) represents topographic shelter from wind independent of a single dominant direction [48]. Sheltered locations may provide more favourable conditions for habitation, livestock management, and outdoor activity by reducing wind exposure and associated thermal stress.

Bora Wind Shelter Index (BWSI) represents topographic protection from bora, a strong and cold regional wind [48,49]. Higher values identify terrain better protected by surrounding relief, whereas lower values indicate more exposed slopes and ridges.

Jugo Wind Shelter Index (JWSI) represents topographic shelter from jugo, a regional southerly to southeasterly wind [48,50]. It captures local differences in exposure that may affect humidity, precipitation, thermal comfort, and the suitability of locations for habitation or pastoral activity.

2.2.2. Reference Data Design

The reference dataset comprised 804 manually delineated reference polygons: 447 locations interpreted as habitation-related TRUE cases (e.g. remains of stone-built houses, shepherds’ dwellings and seasonal pastoral settlements) and 357 surveyed FALSE cases without evidence of habitation (Figure 3). The FALSE class intentionally included not only clearly unsuitable terrain (e.g. steep cliffs, gullies and deep dolines) but also surveyed locations that appeared environmentally plausible yet lacked known habitation evidence. This design reduced the risk that the models would learn only an artificial contrast between archaeological locations and inherently unsuitable for habitation. Criterion values were extracted for each polygon as exact area-weighted raster-cell means, producing a complete 21-variable modelling table.

Before any criterion screening or model development, the reference dataset was divided into spatially stratified training and independent-validation subsets using a 70:30 split. A regular 2 km spatial grid was used to distribute both classes throughout the study area while retaining the target class proportions. The training subset contained 563 polygons (313 TRUE and 250 FALSE), and the independent validation subset contained 241 polygons (134 TRUE and 107 FALSE). The validation subset was then sealed and was not used for criterion screening, response standardization, AHP weighting, machine-learning tuning, model fitting, or threshold selection.

2.2.3. Training-Only Criterion Diagnostics, Response Curves and Sensitivity Analysis

Criterion relevance and potential redundancy were assessed exclusively using the 563 polygons of the training reference subset. Individual criteria were evaluated in terms of their ability to distinguish habitation-related TRUE locations from FALSE locations using descriptive statistics, univariate discrimination measures, and a fixed five-fold spatial cross-validation procedure. Potential redundancy among predictors was examined using rank correlations, with aspect treated as a circular variable. Correlation was used as a diagnostic rather than as an automatic criterion-exclusion rule.

Spatially cross-validated response curves were subsequently derived for each criterion to characterise the form of its relationship with habitation suitability. These empirical relationships, together with the observed training-data distributions and the geoarchaeological interpretation of individual variables, were used to define transparent criterion-specific GIS-MCDA membership functions. All 21 criteria were standardised to a common 0–1 suitability scale using linear, optimum-based, circular, or nonlinear functions depending on the observed response and the expected underlying environmental process.

2.2.4. GIS-MCDA Modelling

Two different GIS-MCDA models were constructed from the standardised geoarchaeological criteria. The Equal-Weight model (EW GIS-MCDA) assigned an identical weight to each of the 21 criteria:

$$w_i={\displaystyle \frac{1}{21}}=0.047619$$

Accordingly, habitation suitability at location $x$ was calculated as:

$$H{S}_{EW}\left(x\right)={\displaystyle \frac{1}{21}}\sum _{i=1}^{21}{S}_i\left(x\right)$$

where $S_i\left(x\right)$ denotes the standardised suitability value of criterion $i$ at location $x$.

The second GIS-MCDA model was based on the Analytic Hierarchy Process (AHP), one of the most used methods for criteria weighting in GIS-MCDA [51,52,53,54]. Pairwise judgments were defined independently of the training AUC values and without access to the independent validation data. Comparisons were first conducted among the four thematic groups and subsequently among the criteria within each group. The final group weights were 41.55% for morphometric criteria, 29.26% for archaeological criteria, 18.49% for hydrogeomorphological criteria, and 10.70% for climate criteria. The five largest individual criterion weights were assigned to DR (13.32%), SLO (12.84%), TRI (8.50%), DSO (7.69%), and LSF (7.67%). All pairwise comparison matrices satisfied the conventional consistency requirement, with a maximum consistency ratio of 0.0263.

The final AHP suitability surface was calculated as:

$$H{S}_{AHP}\left(x\right)=\sum _{i=1}^{21}{w}_i{S}_i\left(x\right)$$

with the criterion weights constrained to sum to one:

$$\sum _{i=1}^{21}{w}_i=1$$

2.2.5. Predictive Machine-Learning Modelling

RF and XGBoost were implemented as complementary nonlinear predictive models using the same 21 criteria employed in the GIS-MCDA analysis. Unlike the GIS-MCDA models, the machine-learning models were trained on the original harmonized criterion values rather than the standardized 0–1 suitability score. Because aspect is a circular variable, ASP was represented by its sine and cosine components:

$$ASP\_sin = sin\left(ASP\right), ASP\_cos = cos\left(ASP\right)$$

This transformation produced 22 numerical model features while preserving the original set of 21 conceptual criteria.

Both algorithms were tuned exclusively within the fixed five-fold spatial cross-validation design. For RF, 24 predefined combinations of $mtry$ and minimum node size were evaluated, with 1000 trees grown per forest. The selected RF configuration was:

$$mtry = 7, min.node.size = 3, n\_trees = 1000$$

For XGBoost, 24 predefined combinations of learning rate ($\eta$), maximum tree depth, minimum child weight, and number of boosting rounds were evaluated. The remaining sampling and regularization parameters were fixed as follows:

$$subsample = colsample\_bytree = 0.80, \gamma = 0, \lambda = 1, \alpha = 0$$

The selected XGBoost configuration was:

$$\eta = 0.03, max\_depth = 3, min\_child\_weight = 1, n\_rounds = 250$$

Criterion importance was evaluated by permutation within the held-out spatial folds. For each original criterion, deterioration in predictive performance after permutation was quantified as the increase in Brier score and the reduction in AUC. The sine and cosine components of ASP were permuted jointly so that aspect remained a single conceptual predictor. Finally, RF and XGBoost were refitted using all 563 training polygons and their locked hyperparameters and then applied to the harmonized raster stack. The resulting probability surfaces were restricted to the same common analysis mask as the GIS-MCDA models and checked for consistent raster geometry, valid-cell count, and the expected probability range:

$$0 \le p\left(x\right) \le 1$$

2.2.6. Independent Validation and Spatial Robustness Assessment

Only after all four habitation suitability models had been finalized was the independent 30% validation dataset used for their accuracy evaluation. Mean model scores were extracted for the 241 validation polygons using exact polygon–cell overlap. Model discrimination was assessed using receiver operating characteristic (ROC) curves and area under the curve (AUC) values with 95% DeLong confidence intervals. The AUC direction was specified a priori so that higher scores represented the TRUE class, and scores were not reversed after inspection of the validation results.

Because all models were evaluated using the same validation polygons, pairwise differences in AUC were tested using paired DeLong tests. Holm correction was applied across the six pairwise model comparisons. RF and XGBoost were additionally evaluated as probabilistic models using the Brier score, log-loss, calibration intercept, and calibration slope. These calibration metrics were not applied to the GIS-MCDA outputs because their 0–1 values represent suitability indices rather than calibrated event probabilities. No classification threshold was derived from the validation data.

A supplementary spatial robustness analysis examined whether model performance was sensitive to the distance between the independent validation polygons and the training data. The minimum polygon-to-polygon Euclidean distance between each validation polygon and the 563 training polygons was calculated in EPSG:3765. Two distance-restricted validation subsets were defined a priori: polygons located ≥250 m and ≥500 m from the nearest training polygon. ROC analysis was performed only when at least ten observations were available in each response class. Consequently, the ≥500 m subset was retained as a sample-size diagnostic but was not used for AUC estimation when too few TRUE cases remained.

3. Results

3.1. Criterion Relevance, Redundancy and Response Structure

Training-only spatial cross-validation showed clear differences in the standalone discriminatory power of the 21 criteria (Table 1). Overall, terrain morphology, accessibility, and climate criteria provided the strongest individual separation between TRUE and FALSE reference polygons, whereas several hydrogeomorphological and resource-proximity criteria showed comparatively weak standalone discrimination. However, low univariate performance was not considered sufficient justification for criterion exclusion, as predictors with limited individual discriminatory power may still contribute useful information through nonlinear effects and interactions within multivariate models [55,56].

Redundancy analysis identified several strongly correlated predictor pairs, particularly among terrain-structure and exposure-related variables, indicating partial overlap in the environmental information represented by individual criteria. Correlation was likewise treated as a diagnostic rather than an automatic exclusion criterion, because related predictors may contribute differently to GIS-MCDA standardisation and nonlinear machine-learning models. Accordingly, all 21 criteria were retained for subsequent habitation suitability modelling.

3.2. Accuracy of Habitation Suitability Models

Four habitation suitability models were produced for Biokovo Nature Park: two GIS-MCDA models based on Equal Weight (EW GIS-MCDA) and the Analytic Hierarchy Process (AHP GIS-MCDA), and two machine-learning models based on RF and XGBoost. Although all four models identified broadly similar areas of favorable and unfavorable terrain, the machine-learning outputs showed a more spatially contrasted pattern, with high-suitability areas more sharply separated from the surrounding landscape (Figure 4).

All four models showed strong discrimination between habitation-related TRUE locations and FALSE locations (Figure 5; Table 2). The two machine-learning models achieved the highest predictive accuracy, with XGBoost reaching an AUC of 0.9886 and RF an AUC of 0.9876. Their confidence intervals largely overlapped, and the difference between them was not statistically significant, indicating essentially equivalent discriminatory performance.

The GIS-MCDA models also performed well, although their accuracy was lower than that of the machine-learning approaches. The AHP model achieved an AUC of 0.9447, compared with 0.9152 for the Equal-Weight model. Pairwise statistical comparisons confirmed that AHP performed significantly better than Equal Weight, while both RF and XGBoost significantly outperformed the two GIS-MCDA models. These results indicate that structured expert weighting improved predictive performance relative to simple equal weighting, while nonlinear machine-learning algorithms provided a further increase in discrimination.

Probability-based performance measures (Brier score and log-loss), which are applicable only to the RF and XGBoost outputs, also slightly favored XGBoost (Table 2). However, the difference between the two algorithms was small, reinforcing the conclusion that both models provided similarly strong predictive performance.

Supplementary distance-based validation was used to examine whether the high model accuracy could be explained primarily by validation locations situated close to the training data. Restricting the analysis to validation polygons located at least 250 m from the nearest training polygon did not reduce model discrimination, with the relative ranking of the four models remaining unchanged. This provides additional evidence that the high validation accuracy was not driven solely by spatial proximity between training and validation samples. At distances of at least 500 m, however, too few TRUE cases remained for reliable ROC analysis. The independent 241-polygon validation dataset therefore remains the primary measure of model accuracy, while the distance-restricted analysis provides supplementary evidence of spatial robustness.

3.3. Final Habitation Suitability Model of Biokovo Nature Park

Based on the independent validation results, the XGBoost model was selected as the final spatial representation of habitation suitability in Biokovo Nature Park (Figure 6). The continuous model output was classified into five suitability classes using Jenks natural breaks classification: very low, low, medium, high, and very high. Of the approximately 193.17 km² modelled area, 120.77 km² (62.52%) was classified as very low suitability, 24.92 km² (12.90%) as low suitability, 15.17 km² (7.85%) as medium suitability, 13.27 km² (6.87%) as high suitability, and 19.03 km² (9.85%) as very high suitability. Thus, very low and low suitability together accounted for approximately 75.42% of the modelled area, indicating that most of Biokovo represents environmentally demanding terrain with limited potential for habitation. In contrast, high and very high suitability cover only 16.72% of the study area.

The final model indicates that favourable habitation environments on Biokovo are relatively scarce and occur as a discontinuous network of spatially restricted patches embedded within a predominantly unsuitable mountain karst landscape. The pronounced dominance of very low suitability reflects the generally harsh geomorphological conditions of the massif, where steep slopes, high terrain ruggedness, exposed ridges, and strongly dissected karst terrain constrain the availability of locations suitable for habitation. Very low and low suitability predominated particularly along the steep southwestern exposed ridges and highly dissected slopes. In contrast, high and very high suitability occurred in smaller, spatially fragmented patches, mainly within the interior karst plateau and on locally gentler terrain.

4. Discussion

4.1. Spatial Patterns and Environmental Controls of Habitation Suitability

The results demonstrate that habitation suitability on Biokovo cannot be explained by a single geographical characteristic but rather reflects the combined influence of terrain morphology, accessibility and local climatic conditions. Although individual criteria differed considerably in their discriminatory capacity, slope (SLO), distance to mountain paths and roads (DR), and Heat Load Index (HLI) consistently emerged among the most influential predictors in both machine-learning models. Their consistent importance in both independently tuned algorithms suggests that these relationships are not model-specific.

From an archaeological perspective, this combination is particularly meaningful in a high-mountain karst environment. Slope affects not only the physical suitability of terrain for habitation but also movement, construction, livestock management and the availability of usable surfaces [7]. Accessibility, represented by distance to the mapped network of mountain paths and roads, further emphasises the importance of connectivity within the mountain landscape. Heat load, in turn, reflects local differences in thermal conditions that may influence snow persistence, soil moisture, vegetation and the seasonal suitability of particular locations. Favourable areas were therefore characterised not by a single resource or environmental advantage, but by the convergence of several conditions conducive to repeated or seasonal human activity.

The spatial organisation of these favourable areas further indicates that habitation potential should not be understood as a continuous altitudinal or topographic zone, but rather as a network of localised favourable areas embedded within the karst massif. Suitable locations generally occur as spatially restricted patches separated by considerably more rugged and restrictive terrain. In mountain landscapes, such spatial fragmentation may reflect the organisation of different activities around dispersed but potentially complementary areas rather than within large and continuously occupied settlement zones [58].

In several parts of the mountain, areas of high and very high suitability correspond to landscapes characterised by intensive anthropogenic modification, particularly dense networks of dry-stone structures associated with pastoral land use. Houses, livestock enclosures, walls, paths and other dry-stone constructions represent the material expression of the organisation of mountain space and of the long-term investment required to maintain pastoral activities in a demanding karst environment. Their concentration within these areas suggests repeated use of favourable parts of the mountain for habitation and economic activities.

This spatial organisation is expressed in different forms across the interior karst plateau (Figure 7). In the highest part of the massif around Sveti Jure (Figure 7A), favourable terrain forms a fine-grained mosaic associated with locations such as Studenci, Baškovića staje, Lokva and Radov dolac. In contrast, the Lemišni doci–Plužine–Podglogovik area (Figure 7B) contains a more extensive and spatially coherent concentration of favourable terrain, while Kupušnjak near Vošac (Figure 7C) represents a more elongated and spatially constrained favourable zone within otherwise rugged terrain. Suitable environments on Biokovo therefore occur at different spatial scales and in different topographic settings.

Podglogovik, located within this more extensive favourable zone (Figure 7B), provides a clear example of how modelled suitability corresponds with the archaeological use of mountain space. The earliest documented human activity within this relatively restricted area is represented by a prehistoric necropolis. The same area also contains a pastoral settlement and its associated infrastructure, water sources, sacred structures and communication routes of different periods. The landscape was also repeatedly modified and reinterpreted: the church of St Elijah was constructed on a prehistoric burial mound, while individual tumuli were later reused for observation positions and trenches during the First World War [5]. Podglogovik is therefore significant not only for the density of archaeological remains, but also for the repeated organisation and reuse of the same spatial unit for different purposes.

A comparable pattern can be observed at Lokva, within the fine-grained mosaic of favourable terrain identified in the wider Sveti Jure area (Figure 7A). The prehistoric tumulus situated at approximately 1460 m a.s.l. contained a burial attributed to the Cetina cultural group and radiocarbon dated to 2410–2198 BC, demonstrating that the highest parts of Biokovo were already incorporated into patterns of human mobility and landscape use during the third millennium BC [5,21]. The same tumulus was reused several millennia later for two burials dating to the fourteenth century. The wider landscape around Lokva contains a pond, a clay source, remains of dry-stone structures and enclosed dolines, as well as several speleological objects, including Dusa Ice Pit nearby and Stara ledenica along the communication route between Vošac and Lokva. At Lokva, the repeated use of this broader spatial unit coincides with access to pasture, water, raw materials and communication routes.

Importantly, both the Podglogovik necropolis and the Lokva tumulus are located within areas of very high suitability for human habitation and landscape use. This correspondence should not be interpreted as evidence that burial locations themselves represent habitation sites. Rather, their placement within favourable zones suggests that funerary places formed part of wider landscapes that were accessible and repeatedly used for different economic and social activities. The occurrence of archaeological features from different periods within these same favourable spatial units reflects the long-term importance of particular parts of the mountain, although the character and intensity of their use changed through time.

Figure 8 shows how the modelled patterns correspond with habitation-related features at Lemišni doci, Podglogovik and Lokva, where dwellings, agricultural or pastoral surfaces, pathways and water-related features occur within broader karst micro-landscapes. Their correspondence with high- and very-high-suitability zones indicates that the model identifies broader environmental settings rather than individual settlement locations.

Accessibility was particularly important, with distance to mountain paths and roads among the strongest predictors in the models. This result requires chronological caution because the mapped communication network cannot automatically be considered contemporaneous with all archaeological features included in the analysis. Nevertheless, mountain routes are strongly constrained by topography and frequently follow passes, gentler slopes and naturally accessible corridors [59]. Later paths may therefore preserve or approximate older patterns of movement, or reuse earlier communication corridors, even when direct chronological continuity cannot be demonstrated. Historical evidence from Biokovo also illustrates the importance of such corridors: in the Sveti Jure area, two routes ascending from the coast, one via Lokva and the other via Vošac, converged with paths from Župa, Krstatice and Vrgorac, while an old route documented cartographically in 1968 continued towards Mucića and Crna ledenica ice pits [60].

In the high karst zones of Biokovo, permanent surface water is scarce, and pastoral communities developed multiple strategies for securing water, including the use of ponds, wells, stone troughs and natural depressions. The relatively weak standalone discriminatory capacity of distance to ponds should therefore not be interpreted as evidence that water was unimportant. Rather, it indicates that proximity to mapped ponds alone was insufficient to determine habitation suitability, which depended on the interaction of several geographical characteristics. Historical and ethnographic evidence also documents the use of snow and ice from natural ice pits to supplement conventional water sources during periods of scarcity [60,61,62,63].

Pastoral settlements associated with communities from Župa, Krstatice, Zagvozd, Brela, Bast, Tučepi and Podgora were distributed across the high-mountain landscape, some situated close to speleological objects such as Dusa Ice Pit, Crna ledenica, Mucića ledenica, Jarova rupa and Ledenica pored barake. Pastures, pastoral settlements, communication routes and different water resources thus formed interconnected components of the high-mountain pastoral landscape.

The chronological depth of this evidence requires a distinction between continuity of activity and persistence of spatial significance. The presence of prehistoric burial mounds, medieval burials, pastoral settlements, sacred structures, historical routes and more recent modifications within the same favourable areas does not demonstrate uninterrupted occupation or the continuation of identical practices. Instead, it suggests that particular spatial units were repeatedly selected and reinterpreted by communities with different social and economic requirements [63]. This is particularly evident at Podglogovik and Lokva, where evidence of different forms of activity accumulated over several millennia. Viewed from a longue durée perspective, this recurrent use reflects the interaction between relatively persistent geographical conditions and changing human practices over time [64,65].

The different modelling approaches provide complementary views of the identified spatial pattern. GIS-MCDA provides a transparent representation of the contribution assigned to individual landscape characteristics, whereas Random Forest and XGBoost capture nonlinear relationships and interactions between predictors. Although the machine-learning models achieved significantly higher discriminatory performance than the GIS-MCDA approaches, different modelling approaches converged on broadly comparable spatial patterns and identified several of the same influential predictors. This convergence suggests that the resulting suitability pattern reflects meaningful landscape structure rather than the behaviour of a single modelling technique. Taken together with the archaeological evidence, the results indicate that geographical conditions influenced human choices, while repeated human activity progressively transformed favourable areas into a structured cultural landscape.

5. Conclusions

This study demonstrates the potential of integrating GIS-MCDA and machine-learning approaches for modelling habitation suitability in complex karst mountain landscapes. Based on the results obtained for Biokovo Nature Park, the following conclusions can be drawn:

• The integration of 21 morphometric, hydrogeomorphological, climatic and archaeological criteria with a reference dataset of 804 polygons (447 habitation-related TRUE cases and 357 surveyed FALSE cases) provided a structured framework for identifying environmental conditions associated with past mountain habitation. The separation of training and independent validation datasets ensured that model performance was evaluated on data not used during model development.
• All four developed models showed strong independent discriminatory performance. The AHP GIS-MCDA model (AUC = 0.9447) significantly outperformed the Equal-Weight model (AUC = 0.9152), indicating that structured weighting improved the predictive capability of the GIS-MCDA approach. Machine-learning models achieved substantially higher discrimination, with RF and XGBoost reaching AUC values of 0.9876 and 0.9886, respectively, while their performances were statistically comparable.
• GIS-MCDA provides a transparent representation of criterion standardisation and weighting, whereas RF and XGBoost are better able to capture nonlinear relationships and interactions among environmental variables. Their combined use therefore brings together the interpretability of GIS-MCDA and the higher predictive performance achieved by the machine-learning models.
• The final XGBoost model indicates that favourable habitation environments occupy a relatively limited part of Biokovo. High- and very-high-suitability zones account for only 16.72% of the modelled area and occur mainly as spatially discontinuous patches within the interior karst plateau. Their distribution around dolines, gentler surfaces and other locally favourable micro-landscapes supports an interpretation of Biokovo as a selectively used mountain environment rather than a uniformly occupied landscape. Archaeological evidence from Podglogovik and Lokva further indicates the repeated use of some of these favourable areas in widely separated chronological periods, without implying continuous occupation.
• The resulting suitability model should not be interpreted as a deterministic reconstruction of settlement locations, but as a spatial framework for identifying areas with increased potential for recurrent human use. It can therefore support targeted archaeological prospection, particularly within predicted high- and very-high-suitability areas that have not yet been systematically surveyed. Prospective field validation will be essential for testing its archaeological usefulness and assessing its transferability to other Mediterranean and Dinaric karst mountain environments.