Between Land and Sea: An Airborne LiDAR Field Survey to Detect Ancient Sites in the Chekka Region/Lebanon Using Spatial Analyses

The interdisciplinary project “Between Land and Sea” combines geological, geomorphological and paleoenvironmental approaches to identify archaeological remains of the Chekka region (Lebanon). In order to record the topography of this area, the first ever scientific airborne LiDAR data acquisition in Lebanon was conducted in autumn 2018. This work describes not only the acquisition and processing of the LiDAR data, but also the attempt to derive possible archaeological sites from the generated elevation model based on methods for spatial analysis. Using an “inverted mound” (iMound) algorithm, areas of possible settlement structures could be identified, which were classified regarding their probability of a possible ancient site using a deductive predictive model. A preliminary validation of some of the detected favoured areas using high-resolution aerial images has shown that the methods applied can provide hints to previously undiscovered sites. It was possible to identify probable ancient wall remains at several detected locations. In addition, least-cost path analyses were performed to reconstruct possible trade and transport routes from the Lebanon Mountains to the Mediterranean coast. The combination of the results of the iMound detection and classification as well as the calculated path system could point to the strategic location of the modern village of Kfar Hazir as a kind of traffic junction. Moreover, reconstructed main transport routes provide indications of heavily frequented roads and may form the basis for further investigations. To validate the results, upcoming field surveys will be realized on site.


Introduction
Due to the civil war in Lebanon (1975Lebanon ( -1990 archaeological activities came to a halt and only in the last decade was intensive scientific research in this country resumed. While some regions are better known, others are still waiting to be investigated. Belonging to the latter is the region between Ras Chekka and Enfeh and its hinterland up to the Mount Lebanon massif (ca. 300 km²). Although it is known from ancient texts that people living in this area during the Bronze and Iron Ages were in contact with the important civilizations of their time (Kopetzky, 2010;Cohen-Weinberger & Goren, 2004), hardly any physical evidence of these periods was detected in this part of Lebanon. Recent publications concentrate on the Tells of Fadous-Kfarabida and Batroun south of the study area (Badreshany & Genz, 2009;Genz, 2016;Höflmayer et al., 2014;Genz et al., 2018;Genz, 2010), on Enfeh north of it (Panayot-Haroun, 2015;Semaan & Salama, 2019) and on Roman temples in the hinterland (Fares, 2010).
The project "Between Land and Sea: The Chekka region in Lebanon", funded by the Austrian Science Fund, intends to investigate the ecological and economical basis of an ancient cultural landscape, which for millennia connected the timber growing areas in the mountains to the shores of the Mediterranean. The Bay of Chekka is a natural harbor and thus an ideal mooring place for ships. Here, two ancient settlements, Tell el-Heri and especially Tell Mirhan, are evidence for this Mediterranean trade (Kopetzky et al., 2019). However, aside from some Roman and medieval remains, no earlier archaeological features are known so far.
This area is a typical Mediterranean landscape with a complex topography covered with macchia and olive groves, which makes a traditional survey on foot nearly impossible. In addition to the invisibility of the ground, modern constructions of buildings and infrastructure have changed this ancient landscape severely at such a rapid speed that conventional methods cannot keep up with and it is feared that archaeological sites will be lost. Thus, it was decided to support the regional field survey with LiDAR technology to detect endangered ancient structures more quickly.
The use of airborne LiDAR (Light Detection and Ranging) in archaeology has become more and more important in the last few decades. Especially in highly vegetated areas, this method can provide detailed surface information with an option to see "through" dense vegetation. Consequently, it has been used for archaeological prospection in forested regions around the world like in Belize (Chase et al., 2011), in Cambodia (Evans et al., 2013), in Slovenia (Štular et al., 2012), in the United States (Davis et al., 2019) and in the Kingdom of Tonga (Freeland et al., 2016).
Spatial analyses that go beyond a visual inspection and interpretation of a LiDAR derived model are rarely used in archaeological prospection, some of those focus on (automated) feature detection (e.g. Freeland et al., 2016;Masini et al., 2018;Davis et al., 2019;Verhagen & Drăguţ, 2012;Cerrillo-Cuenca, 2017). Others concentrate on visibility analyses and calculation of viewsheds in order to recognize spatial patterns (e.g. Vinci & Bernardini, 2017;Doneus & Kühteiber, 2013). One of the most common spatial analyses of digital elevation models in archaeology is the reconstruction of road or path networks with least-cost paths (LCP) (e.g. Diwan & Doumit, 2017;Verhagen & Jeneson, 2012;Seifried & Gardner, 2019;Palmisano, 2017;Howey, 2007) in order to reconstruct the network between known settlements and between settlements and their hinterland. The calculation of LCPs is based on the assumption that people optimize frequently used paths over a long period of time. Because of the variety of new GIS tools and fast developing computing capacities, LCP analyses have been frequently used since the beginning of this millennium (Herzog, 2013). A broad overview of this method including multiple variations and possible problems and errors is described in Herzog (2014).
The presented work deals with the first scientific airborne LiDAR data acquisition in the Lebanon as well as the processing of the data and their evaluation regarding possible archaeological sites in the Chekka region. The main objective of this study was the detection of favoured areas for possible ancient settlements using an inverted Mound (iMound) algorithm and their classification with a deductive predictive model (Verhagen & Whitley, 2012;McCoy & Ladefoged, 2009;Danese et al., 2014) based on several spatial derivations of the LiDAR model. As a result, locations with a high probability of archaeologically interesting findings were identified, which can be taken into consideration in the planning of future regional surveys. As a second objective these very local areas were put in a larger context performing LCP analyses in order to get an idea of possible trade and transport routes of timber (especially cedar) from their source areas (mountains) to the coast, as they were the main commodity in ancient times.

Study Area
The study area is part of the Lebanese Levantine coast at the eastern part of the Mediterranean Sea, and is located within the triangle of Byblos in the south, Tripoli in the north and Baalbek in the east (Fig. 1). It forms the hinterland of the modern city of Chekka, a predominantly industrial town with a large cement factory and many quarrying areas close by. The region extends from the coast around Chekka to the highest peaks of the Lebanon Mountains (more than 3000 m) over an altitudinal range of nearly 1800 m, including different landscape units. Beginning from the West, the narrow coastal plain merges into badlands, which have been deeply incised into a large alluvial plain. This plain was formed by the directly linked river Abu Ali, which drains the upper part of the catchment. Two other major rivers (Nahr el Jaouz and Nahr el Ouadi) with their tributaries drain the southern parts of this region and have formed a complex hilly landscape with steep slopes and deeply incised valleys.
The region around Chekka is dominated by Mesozoic limestone and tertiary sediment deposits. While the sedimentary deposits with their fine to coarse grain sizes are dominated by fluvial and gravitational (landslides) geomorphic processes (badlands), most of the limestone parts show strongly pronounced karst forms. These limestone layers include embedded phosphate nodules and chert bands, which are known as the Chekka Formation (Walley, 1997). This material is now used on a large scale for cement production (Kopetzky et al., 2019). Beside the quarrying activity, the whole area has been massively transformed and changed by anthropogenic impacts, as e.g. terracing for olive groves, road constructions and buildings.
Today's climate is characterized by cool, wet winters and hot, dry summers. Mean annual rainfall usually ranges between 700 and 1000 mm with a high variance in time and space and with the highest values in the mountain region. Since the beginning of the Holocene, drier and cooler as well as wetter and warmer phases have alternated in the eastern Mediterranean (Psomiadis et al., 2018). These changing conditions have had a strong impact on the vegetation (e.g. the tree line) as well as on the geomorphic dynamics and thus on human activity (Broodbank, 2013) within the Levant.
Without anthropogenic disturbances, the natural vegetation would be composed of macchia with deciduous species (e.g. oaks) at lower altitudes evergreen species with increasing altitudes. The Lebanon cedar (cedrus libani) would occur under the current climatic conditions from elevations between 800 and 2100 m (Messinger et al., 2015). Due to overexploitation in the past, it can only be found in twelve reserves in Lebanon today (Hajar et al., 2010). It can be assumed, however, that the lower "cedar line" was subject to fluctuations due to the changing climatic conditions during the Holocene.

Data Acquisition
The ALS data acquisition was extremely challenging, as there was no previous experience in this region for such a kind of investigation. The main problem were the restrictions that apply to the Lebanese airspace. Thus, airborne data could only be collected by the Lebanese military using UH-1D helicopters (Fig. 2).
A Riegl VP-1 Helicopter Pod system was used, consisting of a VUXSYS LR laser scanner, including an inertial mounting unit (IMU), a global navigation satellite system antenna (GNSS; the system used is able to acquire GPS and GLONAS satellite information) and two Sony Alpha 6000 cameras (riegl.com, see Tab. 1). The system was connected via cable to a toughbook inside the helicopter, so that the scanner could be powered (by on-board electronics), operated and the data could be streamed to the PC, where they were stored in real time. Since there was no applicable mount for the VP-1 POD available for this type of helicopter, a prototype for such a mount had to be developed in Eichstätt, which could then be adapted to an exhibition helicopter of the German Museum Munich. After transporting the mount to Lebanon, it had to be approved by a security commission of the Lebanese Air Force in Beirut. After successful approval, a flight plan was drawn up in cooperation with the military leadership, which divided the area into single flight sections. In the end, the area could be surveyed in a total of 14 flying hours on three consecutive days (divided into 7 flights of 2 hours each) along parallel flight strips (strip orientation north-south, 400 m distance between strips) and along profile flights through the main valleys. An eighth flight was scheduled but could not be carried out due to a defect of the helicopter.
Besides the airborne part of the data acquisition, the processing of the trajectory (see 3.2.1) also required the recording of differential GNSS raw data (dGNSS) with a temporal resolution of 1 Hz by a ground team within the study area. These data were used in a postprocessing procedure to correct the raw information of the GNSS antenna on the helicopter in order to eliminate errors, which are caused by e.g. tropospheric/ionospheric refraction. Due to the large differences in elevation within the study area and as specified by the flight plan, a total of 3 locations (see Fig. 1) were chosen as dGNSS reference stations. Two systems were deployed redundantly on each location (Stonex S9III and Stonex S9i).
The large helicopter used turned out to be very suitable for data acquisition and the pilots benefited from a radar system able to determine the flight altitude above ground, unlike commercial helicopters. This enabled the constant checking of the pre-calculated ideal flight parameters. Nevertheless, due to the extreme differences in altitude (valleys and hills) in certain areas, it was not always possible to adjust the flight height quickly to the topography, which resulted in varying point resolutions of the point cloud. An attempt to compensate for this was made by sufficient overlaps between the flight strips and by additional flight routes through the valleys. Despite the overall successful flight mission, problems with the GPS satellite availability occurred during single flights. Dips in the GPS signal occurred and were  concentrated in a relatively limited area, probably the result of the deliberate manipulation of the GPS signal on the ground by military units. This problem ultimately affected individual narrow sections of three flight strips. In certain cases, the resulting errors could be eliminated by the postprocessing of the trajectory and a final strip adjustment, but one flight strip had to be partly removed from the data set.

Raw Data Processing
The processing of the raw data followed a three-step workflow, including the calculation of a precise trajectory (flight line), the combination of the trajectory with the raw scan data and a final strip adjustment. The trajectory was derived by the software package PosPac MMS (Applanix), using the information acquired during the flight by the GNSS antenna and the IMU, both mounted on the VP-1. This unprocessed trajectory information from the scanner system was combined afterwards with the data of the dGNSS ground stations, which were operated during the flights within the single flight areas (see Fig. 1), resulting in a processed trajectory. Then this highly accurate reconstruction of the flight path was transferred to the software package Riegl Riprocess.
Within Riprocess, the raw scan data were linked to the processed trajectory and point clouds were derived for every single flight strip of the mission. These strips as well as the processed trajectory information were then exported and a final strip adjustment was performed by using the approach of Glira et al. (2015) as implemented in the point cloud processing software OPALS (Pfeifer et al., 2014). In order to serve as reference data for the adjustment, the most reliable strips were selected based on quality criteria and adjusted in advance. Consequently, the rest of the strips could be adjusted using more flexible spline trajectory correlation models (Glira et al., 2016). These were necessary due to the significantly varying trajectory quality also within single strips. As a final result about 2.7x10 9 points within 180 flight strips were exported as pointclouds (LAS files) and served as data base for all following processing steps.

Filtering and Classification
After the postprocessing of the raw data, the resulting pointclouds had to be filtered regarding positive and negative outliers and classified, especially concerning ground points. For the following steps the GIS-software SAGA (System for Automated Geoscientific Analyses) (Conrad et al., 2015) with the extension LIS Pro 3D of Laserdata (www.laserdata. at) was used (Petrini-Monteferri et al., 2009) in combination with the programming language Python and the statistics software R.
The workflow for the processing of the ALS data in order to generate a Digital Terrain Model (DTM) was created based on the morphological classification in Hilger (2017) and is shown in Figure 3. At first, the raw pointcloud was divided into 50 x 50 m tiles. This small size of 2500 m² per part was chosen in order to reduce the computation time for the following steps. The tiling of the pointcloud was followed by a removal of the outliers. All points with errors in z-direction (low and air points) were filtered using the SAGA LIS tool "Remove Isolated Points (PC)". The next and most important step was the ground classification itself, which is based on the work of Hilger (2017). Within each cell of the tiled pointcloud (50 x 50 m), the lowest point was declared as a ground point (seed point). These seed points were meshed in a triangulated irregular network (TIN), which worked as a preliminary ground surface. Next, the cell size of the tiles was reduced step by step until it reached a size of 10 x 10 m. With every step, new lowest points were detected and integrated into the TIN. After that procedure, a previously unclassified point was only declared as ground, if the angle between the TIN plane and the line connecting the point and the closest valid ground point is less than 40° and the point is located within a distance of 1.5 m towards the TIN mesh. After new points were verified as ground, the TIN mesh was updated and the threshold parameters were checked again for all unclassified points until finally no further points could be assigned to the class "ground" (Ma et al., 2018). After the ground classification, the classified pointcloud was carefully checked and corrected manually. The last step included the creation of a DTM using only ground points. To exclude even low vegetation from the model, the z-value of only the lowest point per cell was used to determine the cell value. The created DTM with a cell size of 1 m² is the base for all following evaluation steps.
Between Land and Sea: An Airborne LiDAR Field Survey to Detect Ancient Sites in the Chekka Region ...

Data Evaluation by Spatial Analysis
The derived DTM was analysed using spatial calculations in order to detect possible ancient settlements and other traces of anthropogenic structures. These analyses were used in a final step to detect favoured areas for possible archaeological sites. The basic workflow of the spatial analyses is shown in figure 4. First of all, areas had to be extracted from the input DTM that could be considered as potential sites for settlement structures (as described below, the term iMound is used for these areas). Then those areas were classified regarding their potential as a possible archaeological site. In the end LCP analyses intend to reconstruct possible trade and transport routes from the Lebanese Mountains towards the coast in order to link the coastal sites with their hinterland.

iMound Detection
Regarding the settlements, the main focus of the presented study was on tells and settlements on hills, as a preselection had to be made due to the size of the study area. Many of the known inland sites in the region are either situated at elevated locations (temples of Hardine and Ain Akrine) or on partly artificial tells (Tell el-Heri and Tell Arqa). Consequently, the first task was to extract hill-shaped landforms from the DTM, which were analysed afterwards whether they would have been a suitable place to settle for ancient people or not.
To accomplish this an inverted Mound (iMound) algorithm was used, which was originally developed by Freeland et al. (2016) in order to detect burial mounds in the Kingdom of Tonga. Figure 5 illustrates the detection algorithm of the

254
Jakob Rom et al.
mounds step by step. After a low-pass filter was applied to the DTM in order to remove noise from the model, the DTM was flipped: Hills became pits and vice versa. Now the sinks were filled using the algorithm of Wang and Liu (2006), known from hydrologic modelling toolboxes. The next step was to subtract this filled model from the original one so that only the filled sinks (which are in fact inverted mounds) remain. By using thresholds for minimum mound height (1 m) and maximum plain area (8 ha) the valid iMounds were extracted (Fig. 5) and the DTM was flipped back. All the following investigations were then conducted on these hill-shaped landforms.

iMound Classification
The detected iMounds are primarily landforms and not archaeological features. In order to establish whether those areas may have been suitable locations for ancient settlers or not, they were subjected to a classification with a deductive predictive model. Several parameters will be introduced in the following paragraphs based on which a score value was calculated (see section 3.3.2.5) estimating the probability of an ancient settlement at this location.

Visual Inspection
The beginning of the analyses was a visual examination that resulted in the designation of a first visual parameter (a) and a second parameter (r). The iMounds were manually checked in a multi-stage procedure and evaluated with respect to their location and possible recognizable structures. The first inspection included a visual check of all iMounds using satellite images. For incorrect or completely unsuitable hills, both visual parameters for score calculation (a and r) were set to 0. These include, for example, modern anthropogenic hills in quarrying areas or areas on the narrow ridges of geological structures. For all remaining iMounds the parameter a was set to 1. The second check was carried out only for these valid areas (a = 1). Each iMound was scanned for possible ancient (wall) structures based on different visualization methods of the DTM. These included Analytical Hillshades (simple relief shading) by one light source (cf. Kokalj et al., 2013), primarily using the SAGA tool "Combined Shading", which creates the impression of diffuse light scattering in combination with the slope. The height of the light source above the horizon was set between 15° and 35° and the result was displayed in a grey scale with a linear histogram stretch (Fig. 6A). The slope gradient in degree was also displayed using a linear color scheme (yellow to red). Due to the very high values, both anthropogenic walls (including agricultural terraces, building remains) as well as natural steep structures could be identified (Fig. 6B). In addition, Local Relief Models (Fig. 6E) were calculated using the improved methodology from Hesse (2010) and by varying the used filter radii between 6 and 15 m to be able to detect small-scale features on the flat to medium relief iMound areas. The Sky View Factor (cf. Kokalj et al., 2011;Zakšek et al., 2011) was also applied and was mostly calculated using 16 sectors and radii between 6 and 20 m (Fig. 6F). Finally, positive and negative openness (cf. Doneus, 2013) were also derived in order to visualize the individual iMounds ( Fig. 6C and 6D). Both the number of sectors and the used radii were in accordance with the Sky View Factor.
Only in the case of recognizable structures using the described visualization methods or a particularly suitable location of the iMound, the parameter r was set to 1; in the case of no recognizable structures, the corresponding iMound remains at r = 0. In the end, three types of iMounds can be distinguished based on the visual inspection: Incorrect or unsuitable iMounds (a = 0, r = 0), suitable areas without any indications of archaeological structures (a = 1, r = 0) and suitable areas with features indicating past settlements or a very favourable spatial location (a = 1, r = 1).

Availability of Fresh Water
The availability of fresh water was crucial for ancient settlements, which is why this factor was probably one of the most important in the choice of a settlement location. This assumption was included in the calculation of the score, as it is a vital component in other predictive models as well (cf. Brandt et al. 1992). The current river system was derived indirectly from the DTM using a flow accumulation approach based on the relief parameters, whereby a grid cell must have a potential hydrological catchment area of at least 10,000 m² in order to be classified as a river of sufficient size (cf. Wilson & Gallant, 2000;Quinn et al., 1991). The derived river network served as a basis for further calculations. Two restrictions must be taken into account: (i) the river network has changed more or less strongly since ancient times, although at least in the deeply incised valleys of the larger rivers, especially at higher altitudes in the study area, no drastic changes in flow path are expected. And (ii) surely not all derived flow paths indicate perennial rivers. Hydrological data from the American University of Beirut show that there are only three perennial river sections (Abu Ali, Nahr el Ouadi and Nahr el Jaouz) throughout the year under the present day climatic and hydrological conditions, which surely have changed during the Holocene. However, the basis and accuracy of this dataset remains unclear and it was decided to use the before mentioned minimum hydrological catchment size (10,000 m²) for a probable perennial river.
To get an idea of which iMounds are in a good position regarding the availability of fresh water, LCPs of each iMound (starting from the centre of the polygon) to the nearest calculated river or known spring were computed. A map of presently active springs was provided by Professor Joanna Doummar from the Department of Geology at the American University of Beirut. If the LCPs are calculated based on the DTM with a cell size of 1 m², the computational steps will require a high computing power with long processing times. Another problem with small cell sizes is that the LiDAR data set includes modern roads, which can be easily recognized with a 1 m grid. Since the modern roads usually have more favourable inclination angles than their immediate surroundings, the LCP algorithms prefer to follow the modern roads automatically. As a consequence, the DTM was converted to a cell size of 10 x 10 m in order to prevent this bias, by using a mean value upscaling method. This 10 m grid does no longer contain any relicts of modern roads or other small-scale features and the LCPs could be calculated purely based on the relief parameters.
The parameter D describes the accumulated costs of each iMound to reach the next river or spring using the 6 th degree polynomial function by Herzog (2010) (see chapter 3.3.3.1) for calculating the cost grid. Olaya (2004) proposed to exclude high slope values from the computation of LCPs because very steep slopes can be dangerous for hikers. As the presented calculations offer the possibility to cross hillsides, this threshold was set to a high value (45°).

Area and Circularity
Unrealistic large iMounds probably covering more than one hilltop were avoided by setting the maximum iMound size to 8 ha, the size of the largest Bronze and Iron Age sites known in this region. Nevertheless, the detection of a mound with a large area is a better indication for a valid ancient settlement than a small area, because settlements had to find enough space on the hilltops. However, the shape of an iMound does not necessarily represent the extent of a possible tell. The latter may well extend beyond the boundaries of the digitized mound. In addition, contrary to the previous parameter, no ranking order can be established. This means for example that an 8 ha iMound is not necessarily more suitable than a 5 ha one. Only the fact that too small areas are more likely to indicate natural structures than tells can be considered in the score calculation. To achieve this, the area parameter (A) for iMounds larger than 0.3 ha was set to 1. For smaller iMounds, the value 0.5 was assumed for this parameter.
Tells are often nearly circular, in particular manmade tells (cf. Menze et al., 2006). This applies also to the outlines of several known tells close to the study area, for example Tell el-Heri, Tell Arqa and Tell Ardeh. For these reasons, high circularity of an iMound is considered a significant indicator of a likely settlement location. The circularity (parameter C) was calculated after Skau (1951): where A is the area and p the perimeter of the iMound polyg ones: Parameter s represents the slope of a grid cell in percent. A similar calculation of the cost grid was suggested for example by Bell and Lock (2000).
where A is the area and p the perimeter of the iMound polygon. This formula calculates values close to 1 for more or less round features whereas irregularly shaped areas receive values close to 0.

Visibility
Visibility probably played a major role in the choice of a location for ancient settlements. Therefore, viewshed analyses are part of many predictive models (cf. Danese et al., 2014). If one could overlook large parts of the surrounding landscape when standing on the hill, it was possible to recognize approaching threats at an early stage. It was also possible to control and monitor trade routes when visibility was wide. The radius of the viewsheds computed for the iMounds was set to 2100 m because the field experiments of Fábrega-Álvarez and Parcero-Oubiña (2019) indicate that it is not possible to detect motion of a man-sized object in a vegetated environment at a larger distance.
To be able to perform visibility calculations in an acceptable computing time, the DTM with a cell size of 10 x 10 m was used and the viewshed for each iMound was calculated only from the highest topographic point plus an estimated eye height of 1.60 m. The resulting binary grid was buffered using the established radius of 2100 m detection range (Fig.  7) and the visible area within this buffer radius was summed up for each iMound, resulting in parameter V.

Score Calculation
All parameters of the spatial analysis were combined in a predictive model to calculate a score for each iMound. This score indicates the probability of archaeological features to be expected in the areas according to the variables computed here. The value of the score (S) can range from 1 (very likely) to 0 (very unlikely) and is composed of the parameters described above, whose values can also vary from 1 to 0. The parameters D and V were normalized and adjusted so that the value 1 was assigned to the best and the value 0 to the worst fitting iMound. The normalized parameters are indicated by the subscript "no" in the formula (D no , V no ). Exceptions are, as described above, the area (A) and also the circularity (C), the latter directly resulting from the formula in chapter 3.3.2.3.

S=[(D no +A+C+V no +r)*a]/5
If indications of structures were registered during the visual inspection (r = 1), this can indicate possible archaeological sites. The parameter a in this formula ensures that a score value was only assigned to valid iMounds.

LCP Analysis
The iMound classification in order to detect potential settlement areas on hills was followed by simple LCP analyses to reconstruct former trade routes and paths from mountainous regions to the coast. Three main destinations on the Mediterranean coast were selected. These include the well-known tells of Mirhan and el-Heri as well as a spot further south representing the ancient sites of modern Batroun (cf. Markoe, 2000). Batroun itself is located outside the area, which was surveyed by LiDAR and is therefore not included in the DTM, but the distance from the selected destination to the modern city center is only 4 km and the path is quite flat as it follows the plain of Nahr el Jaouz. As starting points for the LCP calculations, the midpoints of only those iMounds were used, whose score value exceeds 0.5. This ensures that only areas with a high probability for ancient settlements were included in this calculation. The computation of the LCPs was again based on an upscaled DTM with a cell size of 10 x 10 m (c.f. 3.3.2.2).
Four different slope-dependent cost functions were applied for calculating the LCPs, these cost functions are described within the next chapter. As all of these depend on the slope of the DTM only, the resulting path networks likely concentrate on the incised and low-slope flood plains of the main river sections. Even though these plain areas seem well suited for paths, narrow passages and river meander prevent easy progress and require frequent crossing of the rivers. As this leads to considerable problems, especially in the case of large rivers, it was decided to adjust the calculation of LCPs there. For this purpose, the main river sections, according to the data of the American University of Beirut (Abu Ali, Nahr el Ouadi, Nahr el Jaouz), were buffered (Buffer radius: 25 m) and provided with very high costs, so that crossing the river was made considerably more difficult, but not impossible.

LCP Methods
Four different methods were used to calculate the main cost grids. In order to get near to an anisotropic (direction dependent) calculation of the cost surface, the aspect grid was used to set the direction of maximum cost in the SAGA tool "Accumulated Cost" (Olaya, 2004). Once the costs for each of the destinations were accumulated, the final LCPs could be created.

LCP based only on the slope
The first method is the simplest. The slope gradient is directly included in the calculation of the cost grid. In principle, flat routes are preferred to steeper ones:

Cost(s)=s
Parameter s represents the slope of a grid cell in percent. A similar calculation of the cost grid was suggested for example by Bell and Lock (2000).

LCP based on Tobler (1993)
The Tobler hiking function is the most common method for calculating LCPs in archaeology (Herzog, 2010). The formula was developed on an empirical basis and was e.g. published in Tobler (1993). The formula of the speed of a walker as a function of the slope (V(s)) is: where s represents the mathematical slope. To calculate a cost grid for LCP applications, the reciprocal of this formula must be used. It describes how long a hiker needs for a certain distance depending on the slope and is expressed in time per distance (e.g. seconds per meter). (2007) The third method of constructing the cost grid is based on the considerations of Llobera and Sluckin (2007), specified by Herzog (2013):

LCP based on Llobera and Sluckin
2010). The formula was developed on an empirical basis and was e.g. published in Tobler (1993). The formula of the speed of a walker as a function of the slope (V(s)) is: where s represents the mathematical slope. To calculate a cost grid for LCP applications, the reciprocal of this formula must be used. It describes how long a hiker needs for a certain distance depending on the slope and is expressed in time per distance (e.g. seconds per meter).
2) LCP based on Llobera and Sluckin (2007) The third method of constructing the cost grid is based on the considerations of Llobera and Sluckin (2007), specified by Herzog (2013): (2010) Again, s represents the mathematical slope. The formula is a mathematical adaptat where s represents the slope gradient and š a critical inclination value. This formula was developed primarily with regard to wheeled vehicles, which were most likely used in ancient times in the region around Chekka. For the use of carts š must not be set too high. Herzog (2013) recommends values between 8 and 16%, in this work a threshold value of 10% was used.

LCP based on Herzog (2010)
The last method is a 6th degree polynomial function developed by Herzog (2010): ones: ( ) = Parameter s represents the slope of a grid cell in percent. A similar calculation of the cost grid was suggested for example by Bell and Lock (2000).

1) LCP based on Tobler (1993)
The Tobler hiking function is the most common method for calculating LCPs in archaeology (Herzog, 2010). The formula was developed on an empirical basis and was e.g. published in Tobler (1993). The formula of the speed of a walker as a function of the slope (V(s)) is: where s represents the mathematical slope. To calculate a cost grid for LCP applications, the reciprocal of this formula must be used. It describes how long a hiker needs for a certain distance depending on the slope and is expressed in time per distance (e.g. seconds per meter).
2) LCP based on Llobera and Sluckin (2007) The third method of constructing the cost grid is based on the considerations of Llobera and Sluckin (2007), specified by Herzog (2013): (2010) Again, s represents the mathematical slope. The formula is a mathematical adaptat Again, s represents the mathematical slope. The formula is a mathematical adaptation to physiological data from Minetti et al. (2002), which were obtained at slope gradient values between -0.45 and +0.45, and is thus based on the energy that a pedestrian must expend.

Validation of the LCPs
The best LCP method was selected by comparing the LCPs with the recent main traffic routes derived from OpenStreetMap (OSM) data (www.openstreetmap.org). Trails, residential roads, and similar pathways were removed from the data set. It is quite likely that the locations of the main traffic roads today are more or less comparable to the main routes in ancient times. To be able to make a statement about the quality of the LCP methods, a buffer of radius 50 m was created around the LCPs. If OSM roads were detected within these buffer zones, the corresponding path sections of the LCPs were marked. The ratio of total length and corresponding path sections can provide information about the quality of the different LCP methods (Tab. 2). A hit rate of 45% for the Llobera and Sluckin method means that almost half of the road sections within a radius of 50 m actually contain modern main roads as recorded in the OSM data set.
Based on this validation, the best fitting method to calculate LCPs in the considered area around Chekka turned out to be the method after Llobera and Sluckin (2007). Thus, it was used for the reconstruction of the ancient paths. In order to provide a comparable alternative to the Llobera and Sluckin paths, the method based on Herzog (2010) was additionally used for the reconstruction.

DTM
The calculated DTM with a cell size of 1 m² effectively covers an area of about 290 km² and extends from the Mediterranean coast to altitudes of 1780 m in the southeast of the study area (Fig. 8). The arithmetic mean of the density of all recorded points is 9.977 pts/m², but in central areas and regions well covered by the flight strips, significantly higher density values were achieved. Due to the GPS problems during the recording of the data (see chapter 3.1), only minor inaccuracies in the strip adjustment of the affected flight strip occurred.
The classification of the pointclouds led to a reliable detection of the ground points in large parts of the investigation area. For an automated approach of filtering and classification, the overall result is satisfying. Even in densely forested areas, it was possible to detect the ground beneath the leaf canopy using the LiDAR method. Only a few of the 1 m² cells did not contain any ground points resulting in small gaps in the DTM. These originate from buildings or appear  within very dense vegetation. The resulting gaps were thus filled by using a spline-interpolation. Figure 9 shows the differences between the unfiltered Digital Surface Model, the DTM (only ground points), and the final gap-filled DTM. In many publications concerning archaeological prospection LiDAR derived DTMs are used because of the numerous advantages of this method. Archaeological features are often made visible quickly and easily by a simple hillshading. Examples can be found in Chase et al. (2011) for Central America, where the DTM clearly identifies Mayan structures, as well as in Evans et al. (2013) for the Khmer temples of Angkor in Cambodia. In the work presented here, however, no archaeologically relevant features can be made visible without great effort or just using simple visualization methods. One reason is that the region in the hinterland of Chekka has been strongly and almost completely reworked in modern times, whereas, for example, Chase et al. (2011) found almost unchanged landscapes in the rainforest of Belize since Mayan times. The second reason is that a very long time has passed since the Bronze Age, when settlements can be assumed to have been established due to the intensive cedar wood trade in this area. The remains of settlements have either been destroyed by modern construction, or have decayed to a point where they only have a minimal influence on the topography. Because of the elapsed time, the strong modern overprinting and the partly strong relief, a direct detection of settlement remains from ancient times was not feasible.

iMound Detection and Classification
590 areas in the study area were classified as iMounds. Each iMound is clearly identified by an individual identification number (ID). By visual inspection of the respective areas 314 mounds could be assigned to level 1 (a = 0, r = 0), 157 to level 2 (a = 1, r = 0) and 119 to level 3 (a = 1, r = 1). Many of the mounds are concentrated in the areas east of the coastal plain. This strongly fragmented landscape shows a diversified geomorphological setting including many natural hill structures, which were identified by the iMound algorithm.
For iMounds classified to levels 2 and 3, the probability of an ancient settlement at this location was estimated by a score (see 3.3.2.5). A high score reflects suitable parameter values for an ancient settlement. The locations of the iMounds as well as their individual score values are presented in figure 10. The calculated score is theoretically distributed between 0 and 1, the achieved values vary between 0.258 (iMound ID 89) and 0.869 (iMound ID 243). In addition, known and presently active springs are mapped in figure 10 as well. Table 3 shows the iMounds with the 10 highest scores and their parameters for calculating the score value. It is remarkable that seven of the ten highest scores are located in the alluvial plain of the Abu Ali River, mainly near the modern towns of Kfar Hazir and Amioun. The reasons for these high scores in this region are, on the one hand, low slope values and, on the other hand, short and relatively direct paths to rivers. In addition, the visibility in plains shows high values, because on hills -in an otherwise flat area -it is possible to overlook large parts of the surroundings. Due to the large number of iMounds with a high score, probably only a part of them indicate tell locations in this landscape.
As it is very difficult to detect archaeologically interesting structures only by visualizing the DTM, the iMound detection and classification was supposed to point to areas of possible ancient settlements. However, the calculated areas must not be considered as complete. It is very likely that villages were also formed between the detected hilltops, which are not covered by this methodology. However, the iMounds as an automated calculation result provide largearea indications of favourable settlement locations, especially tells or settlements on hilltop structures, without the need to search the DTM manually with great time expenditure.

LCP Analysis
For all iMounds with a score > 0.5 (n = 164) LCPs were calculated to one of the three destinations Tell Mirhan, Tell el-Heri or Batroun. The algorithm automatically selects the best destination for each path. Figure 10 shows the results of these analyses. Based on the validation (see 3.3.3.2), the paths were calculated using the method according to Llobera & Sluckin (blue) and Herzog (red). The intensity of the color indicates how many paths overlap at this point. In this way, possible main routes can be determined. Numerous coinciding LCPs, computed by using both methods, indicate that these apparently very suitable paths were also heavily used in past times. When looking at the reconstructed paths in figure 10, two routes are particularly noticeable: First, the southern route from the high mountains in the southeast to the destination at Batroun near the coast. Especially in the last third, many paths of both methods merge beside the valley of Nahr el Jaouz, and even cross the river at the same location. Shortly before arriving at the destination, 32 Llobera & Sluckin paths and 38 Herzog paths follow the same trail. It should be stated again that high costs were attributed to traverse the river Nahr el Jaouz due to the river morphology. Nevertheless, the possibility of ancient routes along sections of the floodplain remains.
The second remarkable route with large overlaps of both methods leads from the alluvial plain of the Abu Ali River through the badlands towards Tell Mirhan. The reason for this large number of overlaps of the LCPs is that many paths starting in the eastern study area join at Kfar Hazir in order to cross the terrain step from the alluvial plain through the badland area towards the coastal plain. As a result, 67 Herzog paths and 79 Llobera & Sluckin paths follow a similar route towards Tell Mirhan. The fact that nearly all LCPs of those iMounds located east of the plain meet at Kfar Hazir is very remarkable. It implies that today the easiest way from the plain towards the coast without any real alternative is through Kfar Hazir. In fact, the modern main road is located exactly at this point. It has to be considered that the landscape, especially the geomorphic very active badlands, has changed a lot during the last millennia, but that even in ancient times the best way out of the mountains towards Tell Mirhan was probably through this narrow passage. The village of Kfar Hazir is therefore of great importance when considering the LCP results. With the help of the calculation of digital relief parameters from the DTM, it was also possible to derive least-cost Voronoi polygons for the three destinations. This means that for each pixel of the elevation model, the destination that can be reached most easily on the basis of the underlying LCP formula (in this case the formula according to Llobera and Sluckin, see 3.3.3.1) was calculated. The calculation itself is similar to the derivation of hydrological catchment areas from a DTM. However, it is not based on the elevation model itself, but on the cost grid from the LCP calculation. Figure 10 clearly shows that especially the least-cost Voronoi polygon for Tell Mirhan covers large parts of the study area, though the size of the least-cost Voronoi polygons of Tell Mirhan and Batroun can only be determined when considering additional destinations and extending the study area. This observation supports the thesis that Tell Mirhan was of great importance as a harbor town in ancient periods (Kopetzky et al., 2019). In contrast, the Voronoi polygon for Tell el-Heri is relatively small and only covers areas up to a maximum of 9 km inland. For all areas further east, it seems to be easier to move to the coast towards Tell Mirhan or Batroun. This could be an indication that el-Heri was not a large port city but served other purposes. The LCP methods used are based on the slope-dependent cost functions combined with high costs in the buffers around the main rivers. A more detailed analysis including other topographic components (vegetation cover, soil properties, smaller river systems) or even social and cultural components (visibility) could provide further information for the accurate reconstruction of ancient road systems (Herzog, 2013).

Examples of Favoured Areas
The iMound classification has led to the identification of many potentially archaeologically interesting areas. Using highresolution RGB aerial images (ground resolution ca. 0.05 m), which were taken by the two Sony Alpha 6000 cameras during the LiDAR flight, these sites have been examined visually. However, due to the intensive modern construction (buildings, infrastructure, terracing) of some areas, it must be assumed that some sites are not shown in the aerial images. Therefore, on-site visits to verify the results are crucial and will be conducted in the course of the project.
Nevertheless, in some places the aerial images complement the results of the methods presented. An example is the iMound ID 508 about 2 km northwest of the city of Douma. With a score of 0.645, it achieves an above-average value and is moreover not situated in any plain but in complex topographical relief, which further strengthens the score. Already during the visual inspection of the DTM at this location, linear structures were detected and a spring has been mapped in the immediate vicinity (Fig. 10). The visual observation of aerial images of iMound 508 confirms the evidence. In the elevation model as well as in the aerial image the linear structures of ancient wall remains are clearly visible (Fig. 11). The obvious building foundations indicate a probable ancient settlement on this site, which had been unknown to the authors until this point. A survey on site is necessary for further information on this settlement as well as on the nearby iMound at Douma (ID 500), which reaches a top-ten score. Unfortunately, this area is completely covered with modern constructions, thus aerial images cannot provide additional information. Analogous to this procedure, other structures could be identified on aerial images. Among others, wall remains are clearly visible in iMounds 256, 344 and 554. All of these iMounds achieve very high score values. With iMound ID 365 near the villages of Bahbouch and Bziza, structures can be identified on aerial images (Fig. 11). Whether the traces of an ancient settlement still slightly influence the terrain surface or whether these formations are of different origin cannot be determined by remote sensing methods.
Due to the reasons mentioned above not all iMounds could be verified by aerial images. Among these numerous unconfirmed sites, the area around the town of Kfar Hazir will be discussed in more detail. Here, at the transition from the Abu Ali flood plain to the heavily fragmented areas of the badlands, the iMounds with the highest score values can be found. In addition, all LCPs running in this direction roughly coincide in this area. Based on today's relief, Kfar Hazir seems to be the optimal place to overcome the terrain step from the high-altitude flood plain towards the coastal plain. A settlement at this point would have had an outstanding strategic importance in ancient times, as the roads could be monitored and large parts of the alluvial plain could be overseen. As a consequence, there is much to suggest that modern Kfar Hazir can look back on a long settlement history. However, the strong modern imprinting of the landscape has the consequence that hardly any old structures can be seen in aerial images. More reliable information about Kfar Hazir as a potential archaeological site can only be obtained by validation on site.

Conclusion
The DTM created from the acquired LiDAR data of the Chekka region provides a good database for archaeological investigations. The spatial resolution as well as the solid filtering and classification of the pointcloud allow further evaluations. In addition to the spatial analyses presented here, areas of potential archaeological sites could be delimited without time-consuming examination of all the DTM visualizations.
By identifying iMounds and estimating the probability of an ancient settlement at this location, new possible site locations have been detected automatically. As the examples of iMounds 256, 344, 365 and 508 show, some of the detected areas can be verified as archaeologically interesting spots by aerial images. For other calculated areas a verification on site is necessary.
In order to put the locally restricted iMounds into a larger context, connections between numerous favoured areas and destinations on the Mediterranean coast were established by calculating LCPs. Thus, possible main trade and transport routes from the Lebanon Mountains to the coast could be reconstructed. It is noticeable that many of these routes roughly coincide in some places due to the topography. This is most obvious in the village of Kfar Hazir close to a terrain step. Such probable past traffic junctions are of archaeological importance because strategic settlements were often located nearby.
The reconstructed paths, as well as the derived river systems lead to the question of the transport of timber as an important commodity from the higher altitudes to the coast. The calculated DTM shows parts of the large river courses of the Abu Ali, Nahr el Ouadi and Nahr el Jaouz. All three river channels are characterized by narrow river valleys with hardly any straight sections. The problem of log jams is very serious in large parts of the river sections, which would have made transport via the river systems very labour-intensive. The reconstructed route system could have been used for the timber transport. A presumed main road from Kfar Hazir to the coast leads directly to Tell Mirhan. These reasons speak for a preferred transport of timber via a land route, which is also in line with the results of Semaan (2015). However, a combination of fluvial and terrestrial transport is also conceivable, with a higher probability of water transport in the upper mountain regions where the Lebanese cedar grows.
It has been shown that LiDAR data acquisition and subsequent spatial analysis can help to answer archaeological questions. In politically unstable countries like Lebanon, however, such surveys are difficult to conduct. A long and intensive preparation period and negotiations with political, military and scientific authorities must be expected. Shortterm changes in the data collection process and other unexpected problems may occur. Nevertheless, the recording and evaluation of a high-resolution DTM leads to valuable results and makes a first evaluation of relatively large areas possible, which could be shown within this study.