Probabilistic Modelling using Monte Carlo Simulation for Incorporating Uncertainty in Least Cost Path

The movement of past peoples in the landscape has been studied extensively through the use of Least Cost Path (LCP) analysis. Although methodological issues of applying LCP analysis in Archaeology have frequently been discussed, the effect of vertical error in the DEM on LCP results has not been fully assessed. This research proposes the use of Monte Carlo simulation as a method for incorporating and propagating the effects of vertical error on LCP results. Random error fields representing the vertical error of the DEM are calculated and incorporated into the documented and reproducible LCP modelling process using the R package leastcostpath. By incorporating vertical error into the LCP modelling process the accuracy of the LCP results can be understood probabilistically, with the likelihood of obtaining an LCP result quantified. Furthermore, the effect of incorporating vertical error on the LCP results can be expressed through the use of probabilistic LCPs, allowing for a graphical representation of the uncertainty in the LCP calculation, as well as identifying the most probable location of the „true‟ least cost path. The method of understanding LCP results probabilistically is applied to a Roman road case study, finding that the accuracy of the LCP from south-to-north without incorporating vertical error is not representative of the LCP population with vertical error accounted for. In contrast, the accuracy of the LCP without incorporating vertical error from north-to-south is representative of the LCP population. The use of probabilistic LCPs suggests that the location of the Roman road in the southern section of the study area was selected to minimise the difficulty of moving up and down slope, irrespective of the direction of movement. However, the identification of two corridors of similar likelihood of containing the „true‟ location of the LCP in the northern section when modelling movement south-to-north suggests that the input data and parameters used in the LCP analysis are unable to discern which corridor contains the most probable „true‟ location of the LCP. Therefore, this research suggests that different input data and parameters are used and tested.

Uncertainty in LCP analysis has been addressed in numerous studies. Conolly and Lake (2006), Branting (2012), Harris (2000) and Herzog (2014) discuss how the accuracy of the LCP can be affected by using DEMs of different resolutions (also see Becker et al., (2017) for site catchment analysis); Lock and Pouncett (2010) assess how different computational methods for calculating slope can affect the LCP results; and finally Herzog and Posluschny (2011) focus on how the number of neighbouring cells used in the LCP calculation can affect the LCP results. Despite this, there has been little archaeological discussion on the vertical error of DEMs and its implications for the stability of LCP results (Ehlschlaeger and Shortridge, 1996 for non-archaeological example; Herzog and Yépez, 2015). In many analyses the DEM is used as an error-free model of the topography (Wechsler and Kroll, 2006), despite the existence of vertical error (Carlisle, 2005;Fisher and Tate, 2006;Heuvelink, 1998;Holmes et al., 2000;Veregin, 1997;Wechsler and Kroll, 2006). By ignoring the vertical error of the DEM, the stability of the computed LCP is not statistically evaluated. More specifically, the computed LCP is not assessed as to whether it is representative of the LCP population with random error accounted for, or whether the LCP result is a simple random sample that is statistically unlikely and therefore not representative of the LCP population. Due to this, archaeological interpretations based on a single LCP realisation may not be statistically stable, and so inferences such as the reconstruction of the layout of ancient paths (i.e. predictive (Parcero-Oubiña et al., 2019); (e.g. Parcero-Oubiña et al., 2019;Verbrugghe et al., 2017;Verhagen and Jeneson, 2012)) or understanding the factors that may have influenced the location of known routes (i.e. postdictive (Parcero-Oubiña et al., 2019); e.g. (Fonte et al., 2017;Güimil-Fariña and Parcero-Oubiña, 2015;Herzog, 2017)) may be susceptible to change if the LCP modelling process is run with a slightly modified, but equally likely, realisation of the true elevation surface. This paper introduces probabilistic modelling of LCPs via the use of Monte Carlo simulation. The application of probabilistic LCPs will be demonstrated through a case study focusing on understanding the factors that may have influenced the location of a Roman road in Cumbria, England.

Monte Carlo simulation for Uncertainty Propagation
The objective of uncertainty propagation is to assess how error in the model input data propagate to the model outputs (Heuvelink, 1998;Zhang and Goodchild, 2003, p. 94). For clarity, the term "error" used here refers to the unknown difference between reality and a representation of reality, while "uncertainty" is the acknowledgement and modelling of error through the use of probability distributions (Heuvelink et al., 2007;Wechsler and Kroll, 2006). The vertical error (or accuracy) of the DEM refers to how close the DEM"s elevation values represent "reality" (Carlisle, 2005;Wechsler and Kroll, 2006). The vertical error is often summarised by a single global statistic such as the root mean square error (RMSE; Carlisle, 2005;Hunter and Goodchild, 1997;Wechsler and Kroll, 2006), which quantifies the average deviation between ground-based observations and DEM values at a set of control points. Although the RMSE statistic has been criticised for not adequately representing vertical error as it does not account for spatial variance, error distributions are generally unavailable and therefore cannot be incorporated into the analysis (Aguilar et al., 2010;Canters et al., 2002;Fisher, 1998;Heuvelink, 2006;Holmes et al., 2000;Kyriakidis et al., 1999;Zhang and Goodchild, 2003, p. 94). The two main approaches for handling vertical error in the absence of higher accuracy elevation data or ground truth field survey data are of interest here: (1) ignore it; (2) model the vertical error distribution and propagate the error through Monte Carlo simulation (Carlisle, 2005;Gesch, 2018;Hunter and Goodchild, 1997;Oksanen and Sarjakoski, 2005). The use of Monte Carlo simulation can also be applied to DEMs with high spatial resolution, as high resolution does not guarantee greater certainty in the DEM accuracy (Temme et al., 2009). However, Monte Carlo simulation is computationally demanding (Heuvelink, 2006(Heuvelink, , 1998, and often necessities the use of lower resolution DEMs (Fereshtehpour and Karamouz, 2018).
Whilst the majority of archaeological LCP studies mention the DEM used, the vertical error is often ignored and its effect on the stability of the LCP not assessed (e.g. Diwan and Doumit, 2017;Moutsiou and Agapiou, 2019;Rademaker et al., 2012;Rosenswig and Martínez Tuñón, 2020). Herzog and Posluschny (2011) and Herzog and Yépez (2015) recommend assessing the impact of the DEM errors on the stability of LCP results by adding random error, with Herzog and Yépez (2014) adding random values from a uniform distribution to the DEM. However, the study is limited to comparing only two LCPs. The low number of computed LCPs may reflect the high computational requirement of Monte Carlo simulation (Hengl et al., 2010;Heuvelink, 2006), particularly in archaeological cases when multiple factors such as different cost functions or barriers are incorporated into the LCP modelling process (e.g. Herzog and Yépez, 2014). Nonetheless, there remains the need for incorporating vertical error into the LCP modelling process and understanding how this uncertainty affects the stability of the LCP results through the use of probability theory. In order to model the vertical error of the DEM and to propagate the error when calculating an LCP, this research proposes the use of Monte Carlo Simulation (Brouwer Burg et al., 2016;Zhang and Goodchild, 2003, p. 94). Monte Carlo simulation has been used extensively for the modelling of DEM uncertainty in viewshed analysis (e.g. Fisher, 1992Fisher, , 1991Nackaerts et al., 1999), flooding (e.g. Cooper et al., 2015;Gesch et al., 2020;Gesch, 2018;Hawker et al., 2018;Leon et al., 2014), and slope failure prediction (e.g. Davis and Keller, 1997;Holmes et al., 2000;Zhou et al., 2003), however as far as the author is aware, Monte Carlo simulation has not been applied to least cost path analysis in Archaeology. This approach is called the "probabilistic" method and aims to encapsulate the vertical error uncertainty through the incorporation of random error fields that matches the error distribution characteristics of the DEM (Gesch, 2018;Hunter and Goodchild, 1997). Each random error field realisation is added to the DEM (Figure 1: B) which is then used to compute the cost surface (Figure 1: C) and finally the LCP (Figure 1: D). Therefore, each LCP can be viewed as a deterministic realisation from an equally likely realisation of the true elevation surface. From this, the stability of the LCP can be evaluated probabilistically, with the expected value of the LCP population mean being approximated by the sample means of the simulated LCPs incorporating vertical error (Figure 1: E).
Furthermore, the LCP results from the Monte Carlo simulation can be visually communicated probabilistically, that is the likelihood of an LCP crossing a raster cell, by calculating the number of times the least cost paths  The simplest random error field representing DEM vertical error can be defined as random values drawn from a normal distribution of mean = 0 and standard deviation = RMSE of the DEM (Fisher, 1998;Fisher and Tate, 2006;Wechsler and Kroll, 2006). Multiple realisations of the random error field can be generated, resulting in the DEM incorporating slightly different realisations of error (Fisher, 1998). However this ignores the spatial autocorrelation of vertical error (Fisher, 1991;Hunter and Goodchild, 1997;Wechsler and Kroll, 2006;Zhang and Goodchild, 2003, p. 94), with overestimated elevation at a particular location also being overestimated at its neighbouring locations (Temme et al., 2009). The lack of spatial autocorrelation in the random error field can be overcome by applying a mean-low-pass 3x3 filter over the surface (Wechsler and Kroll, 2006). This method only requires a global measure of DEM vertical error such as RMSE and therefore makes it suitable for many

The R package leastcostpath
In order to make the incorporation of vertical error accessible to other researchers, the method has been implemented into the R package leastcostpath (Lewis, 2020). The package is available at https://cran.rproject.org/package=leastcostpath and includes the function add_dem_error(). The add_dem_error() function creates a random error field based on the supplied RMSE and adds it to the DEM (i.e. it automates steps A and B in Figure 1).
The R package leastcostpath is built on the classes and functions in the R package gdistance ( van Etten, 2017). According to van Etten (2017), the results from gdistance are comparable to other software such as ArcGIS and GRASS. gdistance converts the raster map representing cost into a graph ( van Etten, 2017). The centre cells of the raster map are represented as nodes in the graph, with the weight of the edges corresponding to the permeability between adjacent nodes ( van Etten, 2017). The use of permeability is in contrast to other GIS software, where calculations are done using cost, friction, or resistance values ( van Etten, 2017). The permeability graph can be represented as a matrix (van Etten, 2017). Matrices allow for efficient storage using sparse matrices ( van Etten, 2017). Sparse matrices only record the non-zero values and are therefore suitable in LCP analysis as the connection of cells in a raster map is often limited to adjacent cells only, with the vast majority of cells not connected, and thus having a conductance value of zero (van Etten, 2017) (Figure 2). The R package leastcostpath utilises these conductivity matrix surfaces to represent permeability through a landscape. In particular, the matrix method allows for the calculation of both isotropic permeability, whereby permeability across a surface is neither benefited nor hindered by directionality, and anisotropic cost, where the permeability is direction-dependent (Conolly and Lake, 2006;Wheatley and Gillings, 2002). For example, the anisotropic property of slope can be calculated and used with cost functions such as Tobler"s Hiking function (Tobler, 1993) when assessing the difficulty of moving up and down slope. From this, the least cost path can be calculated using Dijkstra"s algorithm (Dijkstra, 1959;Herzog, 2014).

Issue with Vector to Raster conversion of Least Cost Path
The resultant Least Cost Path, represented as a Vector Line of zero width, can be rasterised via the create_lcp_density() function in the R package leastcostpath. That is, the cells of the raster that intersect with the LCP are assigned a value of one, with everywhere else assigned a value of zero. If multiple LCPs are supplied the values are summed; resulting in a Raster that indicates how many times an LCP crosses a cell.
However, the conversion of the LCP to Raster cells is not straightforward. As the LCP is calculated by connecting nodes in the graph based on the supplied adjacent cells in the Raster, increasing the number of adjacent cells to greater than eight neighbouring cells allows the LCP to "jump" Raster cells. Whilst the effects of this "jumping" on LCP results has been discussed when incorporating barriers to movement such as rivers (e.g. Herzog, 2014Herzog, , 2013van Leusen, 2002), the impact of converting the LCP to Raster cells has not been explored.
Although increasing the number of adjacent cells used within the LCP modelling process, of which the R package leastcostpath allows 48, reduces the elongation error (the difference in length between the LCP and the optimal straight line route) (Goodchild, 1977;Harris, 2000;Herzog and Posluschny, 2011;Huber and Church, 1985), it also increases the distance and possibility for the LCP to "jump" Raster cells. This issue becomes a concern when rasterising LCPs as the nodes in the graph chosen when calculating the LCP represent the vertices in the Vector Line, with the coordinates of each vertex being at the centre of the corresponding Raster cell (Figure 3: A and B). Therefore, when the LCP "jumps" Raster cells, the trajectory of the LCP between the centres of two adjacent-but-not-touching Raster cells will cross nearby Raster cells, and will be included when rasterising the Vector Line (Figure 3: C). This results in the rasterised LCP having a width greater than one cell where the LCP has "jumped" Raster cells, and so can lead to inaccurate assessments of how many times an LCP crosses a Raster cell. The impact is assumed to be negligible, however further investigation in how to overcome the issue is outside the scope of this research. The impact of incorporating vertical error on the LCP results will be assessed through the use of probabilistic LCPs. By considering LCP results as expressions of confidence (i.e. probabilistically), the probability of obtaining an LCP within a certain maximum distance from the known route of the Roman road can be statistically assessed. Through this, greater confidence can be given that the probabilistic LCPs generated reflect the most probable location of the "true" LCP, and therefore interpretations on whether the Roman road was located to minimise difficulty when moving through the landscape are more stable and less susceptible to change when different realisations of vertical error are introduced into the DEM. Furthermore, the probability of an LCP crossing a cell will be visualised, providing a graphical representation of the uncertainty in the LCP calculation due to vertical error, as well as visually identifying the most probable location of the "true" LCP.

Background
During the late first century AD and early second, forts north of the Tyne-Solway isthmus were abandoned (Breeze and Dobson, 1985). The forts of Brougham and Ambleside were subsequently built and occupied (Breeze, 1988) (Figure 4). the High Street Roman road (Figure 5), likely built during the late first or early second century A.D (Breeze and Dobson, 1985, p. 6;Hindle, 1998;Shotter, 2004, p. 53), is thought to have connected these two forts and provided connection to the main-western route to Carlisle (Codrington, 1919, pp. 155-157;Cool, 2007, p. 55;Shepherd, 2004, p. 29).

Data and methods
The data and code to reproduce the analysis in this case study are available at: https://doi.org/10.5281/zenodo.4074956.

High Street Roman road
In order to assess how closely the computed LCPs align with the Roman road, the known route of the High Street Roman road is required. Despite the location of the Roman road being available from other sources (e.g. Bishop, 2014;Talbert and Bagnall, 2000), the extant route of the Roman road was recorded by Whitehead and Elsworth (2008) during an archaeological evaluation of the area. This represents the most accurate polygonal representation of the Roman road available, whilst also only including that which can be identified in the field (see Collingwood, 1930, p. 118;Nicholson, 1861, p. 7 for speculation on northern and southern sections). The location of the Roman road was retrieved from the Historic England Scheduled Monuments via historicengland.org.uk/listing/the-list/data-downloads/ (List Entry -1003275).

Computing the Cost Surface and Least Cost Paths
In order to represent the real-world topography and to calculate slope within the study area, the OS Terrain 50 DEM (Ordnance Survey, 2020) was used. Although the spatial resolution of the OS Terrain 50 DEM is greater than the commonly used DEMs such as the SRTM (Federal Geographic Data Committee, 1998) (50m compared to approximately 30m), the RMSE is less (4m compared to 9.73m). Therefore, the OS Terrain 50 elevation values are more accurate and consistent with the true elevation of the study area, and so the impact of the vertical error on the DEM values is reduced.
In order to incorporate vertical error into the LCP results, 1,000 realisations of random error fields were created and added to the DEM (steps A and B in Figure 1). Although Heuvelink (1998)

noted that 100 Monte
Carlo realisations should be enough, the calculation of 1,000 realisations ensures that the LCP results are stable and the estimated population distribution is reliable (Heuvelink, 1998, p. 45). As simulated DEMs have been shown to have unrealistic geomorphology due to small sinks that in reality do not exist (Temme et al., 2009), the DEMs with random error were modified using a sink removal algorithm as implemented in the R package DMMF (Choi, 2018;Wang and Liu, 2006). The R package leastcostpath was used to compute 1,000 conductivity matrix surfaces, with each surface created from a different realisation of the random error field added to the DEM. Nearly all archaeological LCP studies are based on slope as it is critical for all pedestrian travel (Herzog, 2014(Herzog, , 2010Howey, 2007). Slope is commonly calculated by recording the slope value in the direction of steepest descent (Herzog, 2014). However, this does not reflect the anisotropic properties of slope, with the slope value being dependent on the direction of descent (Conolly and Lake, 2006;Herzog, 2014).
Effective slope, which takes into account the direction of descent, was computed in leastcostpath by calculating the difference in elevation between centre cells and their sixteen neighbouring cells, and then correcting for distance (Herzog, 2014;van Etten, 2017;Yu et al., 2003). By calculating the effective slope rather than slope in the direction of steepest descent, the LCPs from origin to destination and the destination to origin may differ (Herzog, 2014). Tobler"s Hiking cost function (Tobler, 1993) which calculates the difficulty of moving between cells in terms of time taken (Güimil-Fariña and Parcero-Oubiña, 2015;Herzog, 2010;Tobler, 1993), was applied to the slope values. Finally, to incorporate the effect of effective slope on the calculated LCP, the 1,000 slope-based conductivity matrix surfaces were used to calculate 1,000 LCPs from north-to-south and 1,000 from south-to-north. The calculation of LCPs from both north-to-south and south-to-north reflect that the intended direction of movement when building the High Street Roman road is not known. However, by calculating the LCPs in both directions, the LCP results may offer evidence to suggest that a certain direction of movement was in mind.

Assessing the stability of Least Cost Path results
In order to assess whether the computed LCPs without vertical error are representative of the LCP population that takes into account vertical error of the DEM, the maximum straight-line distance deviation was identified by calculating the maximum distance between the LCP and the Roman road polygonal vector. This is in contrast to the commonly used buffer method , whereby the proportion of the LCP within a buffer distance from the true linear feature (i.e. the known road) is calculated (e.g. Güimil-Fariña and Parcero-Oubiña, 2015;Lynch and Parcero-Oubiña, 2017;Zohar and Erickson-Gini, 2020). Instead, the maximum straight-line distance can be thought of as representing the worst-case accuracy of the LCP, as well overcoming the subjectivity in the buffer distance used.
The accuracy of the north-to-south and south-to-north LPCs without incorporating vertical error and the two LCP population distributions of the accuracy of the LCPs incorporating vertical error from north-tosouth and south-to-north were statistically assessed using a Monte Carlo hypothesis test. Monte Carlo Hypothesis testing works by ranking a summary statistic (in this case the maximum straight-line distance between the LCP and the Roman road polygonal vector) describing the observed data against a summary statistic from random samples (Hope, 1968). From this, the probability of obtaining an LCP accuracy value equal to less than the accuracy of the LCP from north-to-south and south-to-north without incorporating vertical error was calculated. Furthermore, probabilistic LCPs, whereby the probability of an LCP crossing a cell is calculated, were used to visualise the uncertainty on the most probable location of the "true" LCPs.

Least Cost Paths without incorporating vertical error
The LCPs calculated when not incorporating vertical error represents the paths that minimises the difficulty of moving up and down slope to travel from north-to-south (Figure 6: A) and south-to-north (Figure 6: B). It is apparent that the computed LCP from the south-to-north more closely aligns with the known location of the High Street Roman road, with the maximum distance of the LCP from south-to-north being 607.68m less than the maximum distance of the LCP from north-to-south (Table 1).

Figure 6. Least Cost Paths based on DEM without incorporating vertical error (A) from North to South and (B)
from South to North

Statistical assessment of the stability of Least Cost Paths without incorporating vertical error
The use of Monte Carlo hypothesis testing to compare the accuracy of the 2,000 LCPs incorporating vertical error and the LCPs calculated without incorporating vertical error has identified that the probability of an LCP incorporating vertical error from north-to-south having an accuracy value less or equal to the accuracy of the LCP without incorporating vertical error is 0.507 (i.e. 507 of the 1,001 LCPs calculated via Monte Carlo simulation incorporating vertical error have a value less or equal to the accuracy of the LCP without incorporating vertical error) (Figure 7: A). This is in contrast to the accuracy of the LCP from south-to-north, in which the probability of an LCP incorporating vertical error from south-to-north having an accuracy value less or equal to the accuracy of the LCP without incorporating vertical error is 0.003 (Figure 7: B).

4.3.3Visualising the uncertainty of Least Cost Paths incorporating vertical error
The incorporation of vertical error when calculating LCPs via Monte Carlo simulation has identified regions within the landscape where there is a higher probability of an LCP crossing a cell, with these regions therefore more likely to contain the "true" location of the LCPor in other words, the uncertainty around the "true" location of the least cost path is reduced (Figure 8). Both probabilistic LCPs from north-to-south and south-tonorth show that movement within the southern section of study area is more likely to occur along the ridge, following the High Street Roman road (with particular emphasis on the notable bend). Where the probabilistic LCPs differ is in the northern section of the study area, whereby the uncertainty in the LCPs from north-to-south incorporating vertical error follows a single corridor (Figure 8: A). This is in contrast to the LCPs from southto-north, in which two corridors of similar probability of containing the "true" location of the least cost path can be identified (Figure 8: B).

Discussion
The present case study shows significant differences in the accuracy of the LCP results dependent on whether vertical error is incorporated within the LCP modelling process. Whilst the LCPs calculated without incorporating vertical error identified that the LCP from north-to-south is less accurate (607.68m greater maximum distance) than the LCP from south-to-north, the use of Monte Carlo simulation as a method for incorporating vertical error into the modelling process and propagating its effects on the LCP results has shown that the accuracy of the LCP from south-to-north is statistically unlikely (p-value = 0.003). Therefore, the accuracy of the LCP from south-to-north without incorporating vertical error is not representative of the LCP population. More plainly, this suggests that the LCP from south-to-north without vertical error does not reflect the most probable "true" location of the LCP, but rather the statistically unlikely LCP produced from a single realisation of the DEM. Thus, archaeological interpretations based on this single LCP would be susceptible the change if the LCP was produced using the same DEM with vertical error incorporated. Furthermore, the bimodal density distribution of the accuracy of the LCPs from south-to-north incorporating vertical error suggests that although 1,000 realisations of the LCP incorporating vertical error ensures that the estimated population distribution of the accuracy of the LCP from south-to-north is stable, the results do not converge to a unimodal distribution with an estimated sample mean that is representative of the population mean. The lack of convergence of the LCPs is apparent from the probabilistic LCP from south-to-north, whereby two corridors of similar probability of containing the most probable "true" location of the LCP are identified. Due to this, archaeological interpretations should be based on the likelihood that either corridor is the "true" location of the LCP. More importantly, however, it suggests that the input data and parameters, for example the DEM used or the number of neighbours used in the LCP calculation, does not result in an output where the "true" location of the LCP can be differentiated from another location which shares a similar probability of being the "true" location of the LCP. Therefore, in order to be able to discern the "true" location of the LCP, it is recommended that the different input data and parameters are used in the LCP calculation and tested to see whether the LCP accuracy converges to a unimodal distribution. This is in contrast to the LCP without incorporating vertical error from north-to-south, whereby the accuracy of the LCP without incorporating vertical error is representative of the LCP population (p-value = 0.508). Therefore, archaeological interpretations based on this LCP result will reflect the LCP population and be less susceptible to change due to the effects of vertical error in the DEM.
The high probability of the LCPs from both north-to-south and south-to-north following one single corridor that closely aligns with the location of the High Street Roman road suggests that the Roman road in the southern section of the study area was selected to minimise the difficulty of moving up and down slope, irrespective of the direction of travel. This is in contrast to the probabilistic LCPs from north-to-south and southto-north in the northern section, whereby the LCPs are less likely to follow the High Street Roman road along the ridge, and instead more likely to follow the lower elevation towards the Roman fort of Brougham in the north. Although this is more apparent in the probabilistic LCPs from north-to-south due to less uncertainty on the most probable "true" location of the LCP, the probabilistic LCPs from south-to-north identify two corridors of similar likelihood: one following the lower elevation towards the Roman fort of Brougham; and the other more closely aligned with the High Street Roman road along the ridge. As there is a similar likelihood of the "true" location of the LCP following the High Street Roman road, it may be suggested that the desire to minimise the difficulty of moving up and down slope when moving south-to-north influenced the location of the High Street Roman road in the northern section. However, as the accuracy of the LCP incorporating vertical error does not converge to a unimodal distribution with an estimated mean that represents the LCP population, the corridor of LCPs more closely aligned with the High Street Roman road along the ridge may be caused by the input data and parameters used in the LCP calculation, rather than reflecting the "true" location of the LCP. Therefore, it is not possible to conclude whether the location of the High Street Roman road in the northern section was selected to minimise the difficulty of moving up and down slope.
The most significant contribution of this work lies in the explicit incorporation of the vertical error of the DEM into the LCP results, and how this uncertainty propagates through the LCP modelling process. By incorporating uncertainty into the LCP analysis via Monte Carlo simulation the accuracy of the LCPs can be represented probabilistically, allowing for the identification of LCP accuracy results that are statistically unlikely, and are therefore not representative of the LCP population with random error account for. By ensuring that the LCP accuracy results are stable and represent the LCP population, rather than a single, perhaps, statistically unlikely, LCP realisation, the archaeological interpretations based on the LCP results will be less susceptible to change if the LCP modelling process is run with a slightly modified, but equally likely, realisation of the true elevation surface. Lastly, this research has demonstrated that probabilistic LCPs can be used as a method for visualising and communicating the uncertainty when incorporating vertical error on the LCP results, as well as being used to identify the most probable "true" location of the LCP.

Conclusion
Although LCP analysis has been used extensively in the modelling of past peoples" movement in the landscape, the impacts of vertical error within the digital elevation model on LCP results has not been fully assessed. This research proposes the use of Monte Carlo simulation as a method for incorporating and propagating the effects of vertical error on the LCP results. Through the use of Monte Carlo simulation, the stability of the LCP results can be statistically evaluated as to whether it is representative of the LCP population. Furthermore, Monte Carlo simulation allows for the LCPs to be viewed as probabilistic LCPs, whereby the probability of an LCP crossing a cell can be calculated, and provides a graphical representation of the uncertainty in the LCP calculation due to vertical error, as well as identifying the most probable location of the "true" least cost path.
In the case study, which postdictively assessed the location of the High Street Roman road through the use of probabilistic LCPs, this research suggests that the location of the Roman road in the southern section was selected to minimise the difficulty of moving up and down slope, irrespective of whether movement occurred north-to-south or south-to-north. In contrast, the northern section of the Roman road along the ridge does not follow the most probable "true" location of the LCP from north-to-south, with the probabilistic LCPs following the lower elevation to the east of the Roman road. Therefore, it can be suggested that the location of the Roman road in the northern section was not selected to minimise the difficulty of moving up and down slope from north-to-south. Similarly, the probabilistic LCPs from south-to-north identify two corridors of similar likelihood of containing the "true" location of the LCP: one following the Roman road along the ridge; and the other following the lower elevation to the east of the Roman road. Due to this, it is not possible to discern whether the location of the Roman road in the northern section was selected to minimise the difficulty of moving up and down slope from north-to-south. However, it has been suggested that the two corridors of similar likelihood of containing the "true" location of the LCP may be caused by the input data and parameters used in the LCP calculation. Therefore, this research suggests that different input data and parameters are used in the LCP calculation in order to identify the most likely "true" location of the LCP.

Declaration of Competing Interest
The author declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.