Probabilistic Models for Predicting Archaeological Site Locations in Marj Bisri (Southern Lebanon): Comparing Frequency Ratio and Shannon’s Entropy

3.1 Archaeological and geospatial data

The establishment of the archaeological predictive model for the study area is based on the data provided by the archaeological survey conducted by The Polish Centre for Mediterranean Archaeology at the University of Warsaw between 2004 and 2005 (Jakubiak 2011; Jakubiak & Neska 2005; 2007). The data is publicly available in the Environmental and Social Impact Assessment of the Greater Beirut Water Augmentation Project report (Dar Al-Handasah 2014). Seventy-five sites were recorded during these campaigns spanning from the Palaeolithic through to the present day, 67 of which were incorporated into the present study. Seventeen sites fall within the area of the dam reservoir, 13 sites within a 100 m distance and 25 further sites within a 1 km distance. The remaining 12 sites lie within a distance ranging between 1 and 6 km from the project dam reservoir boundaries (Table 1). The Bisri valley and its surroundings offer a suitable environment for human settlement based on the recorded archaeological sites. The archaeological evidence points to an extensive occupation during the Roman period with 17 sites. One of the major landmarks during the period are the drums of four grey granite columns, which seem to have belonged to a tetrastyle prostyle temple located at the crossing point of the two rivers, the Nahr el-Barouk and the Nahr ‘Aray which form the Nahr el-Awwali (Boulanger 1955; Khalil 2009; Tallon 1967). According to Aliquot (2009), the dimension of the portico (14.70 m) is equivalent to that of the largest rural temples of Lebanon. The existence of a rural settlement connected to the temple is an assumption also considered by Aliquot (2009). The temple is believed to have been dedicated to Eshmoun, the chief deity of Sidon as argued by Khalil (2009; 2015). Evidence of a road network connecting the area of Marj Bisri to the Bekaa is reported by Tallon (1967) and further evidence is provided by Khalil (2015). Since the surveyors did not specify the surface area of each recorded sites. The spatial extent for each site was determined by a bounding square (115.35 × 115.35 m) corresponding to 16 pixels, each matching the spatial resolution of the ALOS Prism 30 (AW3D30) (JAXA, n.d.) This method is particularly effective for overlaying archaeological site surfaces onto geo-environmental raster data. It enables the extraction of values from the centres of 16 raster cells in overlapping areas of archaeological sites using the “extract value to table” function in ArcGIS Pro (Nicu, Mihu-Pintilie & Williamson 2019). The total surface area of the 67 archaeological sites recorded within the study area amounts to 1,032 pixels (476,175 m²), instead of 1,072 pixels (510,459 m²), due to the overlapping spatial extents of some sites given their proximity (Table 1).

The creation of the archaeological susceptibility map necessitates various geospatial data layers. To achieve this, we utilized the ALOS Prism 30 (AW3D30) DSM, which offers a horizontal resolution of approximately 30 meters. This data is derived from images captured by the Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM) on the Advanced Land Observing Satellite (ALOS) between 2006 and 2011. From this DEM, we extracted and classified the following geomorphometric parameters: slope, aspect, elevation, and difference from mean elevation (DIFF). The drainage network was created from the Digital Elevation Model using flow accumulation analysis with a stream threshold based on the mean accumulation value (Tang 2000). The soil map was digitized from the 1:50,000 soil map of Lebanon (Darwish et al. 2006) and georeferenced using the provided control points. Subsequently, it was converted from stereographic coordinates to UTM 36 N.

3.2 Model variables

Six geo-environmental variables were selected based on theoretical and empirical studies to be implemented in the process of generating the archaeological predictive models (Balla et al. 2014; Danese et al. 2014; Graves 2011; Heilen et al. 2013; Mink et al. 2009; Nsanziyera et al. 2018). The topographic variables cost distance to streams, slope, difference from mean elevation (DIFF), aspect, and elevation, were all created from the ALOS Prism 30 (AW3D30) (JAXA, n.d.). Proximity to water sources is one of the most crucial factors in establishing settlements. To assess this, a cost surface for traversing to streams was created using the Tobler hiking function (Tobler 1993) and the Path Distance Spatial Analyst tool, following the workflow by Tripcevich (2009). The author uses the following formula: 0.00016666 * exp (3.5 * abs(s + 0.05)) where s is the mathematical slope to convert speed to time and compute the vertical factor cost. A vertical factor table is then used as input in the ArcGIS Pro Path distance tool to generate an accumulative cost surface based on the vertical relative moving angle (VRMA), which is the angle of slope from the From cell to the To cell. According to this model, the minimum time required to traverse one meter occurs at a downslope gradient of 0.049957916 (–2.86 degrees). This approach measures the cumulative time cost of traversing each raster cell in the study area to reach the nearest available water source. The cost distance-to-streams map was reclassified into 10 different decimal hours’ classes using the equal interval method (Figure 2A). Slope is also a commonly used factor in establishing APM given that settlements are mostly established on areas with gentle slopes. The slope was rescaled logarithmically using the Rescale by Function tool in ArcGIS Pro Spatial Analyst. It was then reclassified into 10 different classes using the Natural Breaks (Jenks) method (Figure 2B). Difference from mean elevation (DIFF) was also implemented in this process. This index measures the difference in elevation between each cell in the DEM and the mean elevation within a specified neighbourhood around that cell. Cells with positive values indicate areas higher than their average neighbourhood such as ridges, while negative values indicate lower areas than the average neighbourhood such as valley floors (Wilson and Gallant 2000; De Reu et al. 2011). The DIFF was computed using the workflow determined by De Reu et al. (2011) using a circular neighbourhood of a 300 m radius. Six distinct morphological classes were identified and classified according to the methodology adapted from Weiss (2001), which is based on the standard deviation of the DIFF values: 1- Valley, 2- Lower slope, 3- Flat slope, 4- Middle slope, 5- Upper slope, and 6- Ridge. This classification followed the workflow outlined by De Reu et al. (2011; Figure 2C). Aspect is also considered as a very influential variable considering that sun exposure and wind direction are crucial factors determining site location (Figure 2D).

Figure 2

Elevation or altitude above sea level was also integrated into the computation process and reclassified using equal intervals (Figure 2E). Soil types are also deemed relevant to settlement choice and thus were included in the modelling process. The soil map digitized from 1:50,000 maps shows 13 different types or combinations of soil types within the study area (Figure 2F; Darwish et al. 2006).

3.3 Investigating Archaeological site Distributions in Relation to geo-Environmental Categories

The initial step involves evaluating the spatial distributions of archaeological sites in terms of their relation to selected geo-environmental variables. The purpose of this is to determine whether there are notable statistical differences between the spatial distributions of archaeological sites observed and the background environment. In other terms, we seek to examine, for instance, if sites tend to be closer to water sources or on gentler slopes than expected by chance (Kvamme 1990).

In this assessment, we use the Kolmogorov-Smirnov goodness-of-fit test, originally developed by Kvamme (1990) for archaeological contexts, and subsequently applied in various case studies (Nsanziyera et al. 2018; Wheatley 1995). A comparison is made between the expected surface area of archaeological sites, treated as a statistical population, which is calculated based on the percentage of surface area for each class within each geo-environmental variable on one hand, and on the other hand, the observed surface area occupied by archaeological sites within each class, treated as the sample (Kvamme 1990). For this purpose, the cumulative distribution of the sample and the statistical population is plotted, and their difference (D) is calculated. A critical value (d) is then compared to the maximum value of D (Dmax). The critical value d of a two-tailed test using a significance level of 0.05 is calculated using the equation 1.36/√n (Thomas 1986).

Two hypotheses can be formulated in this context:

The H0 hypothesis is rejected if Dmax exceeds the critical value d, indicating a statistically significant difference between the sample and population distributions. In Table 2, a statistically significant difference is observed for all geo-environmental variables where Dmax exceeds the critical value d of 0.022. This implies that archaeological sites within the analysis area are likely located in relation to geo-environmental variables. However, as Wheatley (1995) states, the result of this test should be interpreted with caution, as it reflects association rather than causation, and the lack of repeated trials may affect the reliability of the findings.

3.4 Methods

4.1 Validation

To validate the predictive model, the simple Gain statistics after Kvamme (1988) were employed, and the results were expressed using the following formula:

where G is the gain mode, PS is the percentage of the area characterized by the highest probability of hosting archaeological sites, and GS is the percentage of observed archaeological sites within this area. The gain values extend from 0 to 1. A value closer to 1 indicates higher accuracy in predicting archaeological sites. The calculated gain for the FR model is 0.64 indicating a good performance with 24.71% of the observed sites located in the area characterized by the highest potential and covering 8.98% of the study area’s total surface. The SE model also shows a similar performance with a gain value of 0.64. Since Kvamme’s Gain is typically calculated using the same dataset that was used to build and train the model primarily it is therefore considered as an internal validation metric susceptible to overfitting and does not inherently test the model’s ability to generalize to new, unseen data.

External validation is therefore essential to ensure that the model can generalize effectively to unseen data. To achieve this, k-Fold Cross-Validation was implemented, a technique particularly valuable when dealing with a small dataset, as it enables the model to be trained and tested multiple times on different data subsets, optimizing the use of available data (Kohavi 1995; Walker et al. 2023; Wang, Shi & Oguchi 2023). In this case, the dataset was divided into five equal-sized folds. The model was trained on four folds while being tested on the remaining fold, with this process repeated five times. The overall model performance was then averaged across the five test folds. By using k-Fold Cross-Validation, a robust estimate of the model’s performance is obtained, as each site is used for both training and testing in different iterations. In the present case, the 67 archaeological sites were divided into approximately five equal folds by assigning random values to the attribute table in ArcGIS Pro, which were then used to allocate each site to a particular fold. Five APMs were then produced using four folds and tested using the remaining fold and the Kvamme’s Gain for each fold of cross validation was calculated for both methods. The average Kvamme’s Gain is than calculated and compared to Kvamme’s Gain calculated using all sites. If the two value are similar, the model generalizes well. The results in Table 5 indicate that average Kvamme’s Gain values for both methods are not statistically significant as indicated by the results of the T-test were the p-values (0.202 and 0.414) for respectively the FR and Shannon methods are greater than 0.05.

A Monte Carlo simulation using random site distributions was also used to validate Kvamme’s Gain. The aim was to assess how well the Archaeological Predictive Model performs compared to random chance by simulating 100 random distributions of archaeological sites across the study area, matching the actual number of sites in the dataset. The objective was to determine whether the Kvamme’s Gain statistics observed for both APMs are significantly higher than what would be expected from random site distributions. For each iteration of the Monte Carlo simulation, Kvamme’s Gain was calculated. The “observed” Kvamme’s Gain from both APMs was then compared to the distribution of Kvamme’s Gain values generated from the random simulations. As illustrated in Figures 6 and 7, the observed Kvamme’s Gain for both APMs is much higher than the distribution of the Monte Carlo simulated gains, which are mostly clustered around 0. This confirms that the APM models are valid and perform significantly better than random chance.

Figure 6

Figure 7