Spatial modeling of multi-hazard threat to cultural heritage sites

Cultural heritage is the foundation upon which global and historical values are based on. It connects us to the legacy left by our ancestors and identifies who we are as part of the modern society. Globally and specifically in the northeastern Romania, the landscape where cultural heritage sites were built on is constantly evolving due to mass wasting processes. Among these processes, landslide and gullies can disrupt the gravitational equilibrium directly or around these sites, threatening their very existence and our capacity to pass them on to future generations. Because landsliding and gullying are stochastic processes, the use of spatial statistics has often been employed to map locations at risk. In this work, we make use of advanced spatial Bayesian statistics to model landslide and gully erosion susceptibilities, separately. And, we ultimately combine these two outputs into a unified multi-hazard susceptibility model which we cross with the known cultural heritage sites in a study area close to the city of Iaşi, in Romania. Specifically, we implement a Bayesian version of a Generalized Additive Model (GAM) which assumes that the two separate landslide and gully presence/absence distributions to behave according to a Bernoulli probability distribution. Contrary to common practices in the literature, the two susceptibility models both feature fixed and random effects, including covariates acting at a latent level. We do this to also capture the unexplained but spatially coherent distribution of properties not directly included in the model. As for the properties directly expressed as covariates, our GAM features terrain attributes obtained from a LIDAR survey, in addition to land use and soil layers. The two single models outstandingly perform (AUC > 0.9) both during the calibration and validation phases. This modeling procedure ensures that the probability of occurrence of the two mass wasting processes under consideration is well estimated and therefore can be used to reliably plan conservation practices for local cultural heritage sites deemed at risk.


Introduction
Natural hazards do not only pose a threat to infrastructure and human lives (Bell and Glade, 2004) but also to the heritage left by our ancestors (Nicu, 2017a). Heritage sites are a fundamental part of our ancestry and they are a testimony of those that preceded us and the cultural values they built; therefore being a symbol for our identities to cling on. Global cultural heritage sites are not only being threatened by anthropic (e.g., Aykan, 2018) but also by climate change effects (e.g., Fatorić and Seekamp, 2017). At an international level, constant efforts are made to detect (Tapete and Cigna, 2019), monitor (Agapiou et al., 2020), and assess (Nicu, 2016) mass wasting processes for cultural heritage management. Other efforts are aimed towards the improvement of adaptation measures (Guzman et al., 2020), sustainable development (Guo et al., erosion is also a threat (Kincey et al., 2017) because the intense action of gully formation and development through time can mobilize large volumes of soil and damage those sites (Nicu, 2018b).
Of course other natural agents can pose an equal or even greater threat. For instance, weathering (Kuchitsu et al., 2000), tectonic stress (Topal and Doyuran, 1997), earthquakes (Okamura et al., 2013), rock falls (Mineo and Pappalardo, 2020), sea level rise (Marzeion and Levermann, 2014), tsunami (Daly and Rahmayati, 2012), freeze-thaw processes (Grossi et al., 2007) and floods (Ortiz et al., 2016) are known to have caused and will cause damage to cultural heritage (Mäntyniemi et al., 2004). However, the study of these phenomena and their implications are usually implemented in a phisically-based framework, involving geotechnical and engineering solutions (e.g., Pappalardo et al., 2016). As for geomorphological processes such as landslides and gully erosion, the scientific community commonly relies on statistical models to estimate the probability of occurrence (susceptibility, Rahmati et al., 2019;Reichenbach et al., 2018) of each of these processes over space (e.g., Arabameri et al., 2019;Lombardo et al., 2016b). And, on the basis of this information, territorial management institutions can make decisions to protect the heritage sites at risk (Tarragüel et al., 2012). The literature on susceptibility modeling of geomorphological hazards is rich and it vastly grew together with technological advancements since its early years (Brabb et al., 1972), from direct geomorphological mapping (e.g., Verstappen, 1983), to inventory analysis (e.g., DeGraff and Canuti, 1988), heuristic (e.g., Leoni et al., 2009), deterministic (e.g., Bout et al., 2018), statistic (e.g., Frattini et al., 2010) and data mining (e.g., Lombardo et al., 2015) approaches. The last two types nowadays correspond to the vast majority of approaches in susceptibility studies. However, on the one hand the data mining applications are growing in numbers and variety -several articles and relative comparisons are published every time a new, even slightly, different algorithm is proposed from the machine learning community (e.g., Felicsimo et al., 2013;Arabameri et al., 2020) -, while on the other hand, the statistically-based literature is much more stable in terms of adopted methods. Studies are commonly carried out via Generalized Linear Models (Budimir et al., 2015) by assuming that the distribution of a given geomorphological process is explained via the Bernoulli probability distribution (also referred to as Logistic Regression) (Chang et al., 2017;Lombardo and Mai, 2018). Nevertheless, several arguments are still under scientific debate. Among these, we particularly address three of them: 1. The choice of mapping unit is currently a relevant topic, with the landslide community increasingly leaning towards a slope unit partition (Carrara et al., 1995), especially after automatic tools for their calculations have been made freely accessible (e.g., Alvioli et al., 2016). Conversely, gully erosion studies are almost unanimously based on grid-cells (e.g., Arabameri et al., 2018), with very few cases adopting unique condition units (e.g., Conoscenti et al., 2013). Therefore, testing and proving that other mapping units can be used is a topic of particular geoscientific relevance. 2. Even among the Logistic Regression studies, more and more articles have been published by using an extension to the simple linear case, the Generalized Additive Model framework (GAM; e.g., Goetz et al., 2015;Gómez Gutiérrez et al., 2009). The GAM approach is much more flexible than its linear counterpart and allows for the inclusion of any sort of nonlinear effects (Brenning, 2008). Among these effects, Lombardo et al. (2018a) has recently shown that it is possible to capture residual effects which are not expressly represented in the data other than the latent level. 3. Another scientific topic of particular interest consists of the integration of more than one hazard phenomenon into the spatial prediction (e.g., Chen et al., 2016;Pourghasemi et al., 2020) This work investigates the three topics mentioned above. Specifically, we aim at testing a slope unit partition to separately explain both landslide and gully erosion occurrences in a study area located in the north-eastern part of Romania where cultural sites of Neolithic age are widely spread  and vulnerable to natural hazards (Asăndulesei et al., 2020). While doing so, we test a Bayesian version of a binomial GAM to predict presence-absence cases for the two geomorphological processes, by featuring linear properties as well as nonlinear ones (e.g., slope steepness), and also expressed at the latent level (e.g., effects not directly expressed or defined as covariates, see Bakka et al., 2018, and Section 3.3). Ultimately, we combine the two separate susceptibilities into a unified susceptibility map and intersect it together with the heritage sites. As a result, we propose a modeling tool that not only can assess locations where the two separate processes are likely to occur but that can also indicate likely heritage sites at risk from both geomorphological processes.
The article is structured as follows: §2 introduces the study area and the three inventories namely, heritage sites, landslides and gully heads; §3 describes mapping units, covariates and modeling strategy; in §4 we report the results which are discussed in §5; and in §6 we summarize our concluding remarks.

Study area
Bahluieţ river basin is located in the north-eastern part of Romania ( Fig. 1a) and has a surface of 550km 2 . Administratively speaking, belongs to the Iaşi county. The catchment is part of the Moldavian Plateau. Due to its main geological features, Bessarabian deposits of Sarmatian age, the area is highly susceptible to landslides and gully erosion (Nicu and Asăndulesei, 2018).
This underlying geology represents the base for initiation and development of geomorphological processes (landslides and gullies). Bessarabian deposits mainly correspond to poorly consolidated clay marls with sand intrusions and oolithic limestones. These deposits are friable with a simple touch and because of the fine granulation can be weathered and mobilized by hydrological agents with relative ease. Field experience has shown that landslides primarily trigger in those poorly consolidated deposits when slope steepness is high. Conversely, our field experience suggest that gullying is favored by lithology, steep slopes, lack of vegetation and deforestation. Within the area, gullies develop with considerable depths ranging between 20 and 25 m.
Climatically speaking, the area is exposed to continental temperate conditions with alternate heavy rainfall and periods of draught. Precipitations range between 500 and 700 mm/year and the average annual temperature is between 8.3 9.6 °C. The dominant land use in the area consists of arable lands and pastures. The area is intensively used for agriculture and animal husbandry; these represent the main economic activities in a dominant rural area (Romanescu and Nicu, 2014;Nicu, 2018aNicu, , 2018b. Detailed description of the study area can be found in Romanescu and Nicu (2014); Nicu (2018b,a), Nicu and Asăndulesei (2018). The area is located in one of the most favourable landscape in Eastern Europe for the development of Neolithic civilisation (Cucuteni culture). As shown by Nicu et al. (2019), the area has a tremendous potential for hosting Neolithic sites still undiscovered. Within the basin, 107 Neolithic sites are currently known (see Fig. 1b). This inventory was compiled from the national databases (National Archaeological Registry -RAN, Institute of Cultural Memory -CIMEC and the NAtional Heritage Institute -INP), Archaeological Registry of Iaşi County (Chirica and Tanasachi, 1985), and field trips with archaeologists.
Both landslides and gully inventories were made initially using LiDAR data and then confirming them via field investigations. Overall, 764 landslides (predominantly classified as debris slides, sensu Hungr et al., 2014) and 1042 gullies were manually digitised, whose spatial distributions are shown in Fig. 1, panels c and d. As for the corresponding size characteristics, these are shown in Fig. 2 where landslides are shown to be much larger and gullies are shown to be much more elongated, as expected. Previous studies have argued that gullies and landslides partially overlap in the area (Nicu, 2018b); this may imply that they could have a faster evolution in the future.

Slope unit and presence/absence assignment
Mapping units subdivide study areas into geographic polygonal structures for which any predictive model makes estimates. As a result, the choice of a mapping unit regulates the spatial scale of a given susceptibility model. The geoscientific literature reports several mapping units, although in terms of number of applications, they essentially boil down to two types, grid-cells and Slope Units (Reichenbach et al., 2018). On the one hand, grid-cells allow for a fine and regular spatial partition but they are inevitably sensitive to all the uncertainties in the geomorphological mapping and the process itself. Conversely, Slope Units (SUs) subdivide the geographic space into a coarser and irregular spatial partition that is closely linked to the geomorphological process under consideration; but, they force the user to make some assumptions on how to summarize the covariate distribution from the grid-cell to the SU level itself (e.g., Castro Camilo et al., 2017).
The dilemma related to the use of the most appropriate mapping unit between grid-cells and SUs starts in the late 90's since the first SU appearance (Carrara et al., 1995). And, it has been investigated in several studies since then (e.g., Van Den Eeckhaut et al., 2009). In the landslide community, both grid-cells and SU have become more and more common. However, in gully erosion studies, grid-cell are the standard with no SU ever appeared to the best of our knowledge. For this reason, here we test whether SU can be a meaningful spatial partition even in the context of gully erosion susceptibility studies. More specifically, to combine landslide and gully erosion susceptibility in a single map, we had to choose for a common spatial partition, which we opted to be made into SUs. By using the software r.slopeunits (Alvioli et al., 2016), we have tested three different parameterizations (keeping the flow accumulation and circular variance fixed, and changing three combinations of the minimum SU area) whose result are not reported here. And, we selected a SU partition with a total number of 4978 polygons, with a mean surface of 0.1 km 2 and a 95% confidence interval of 0.43km 2 . The presence-absence status for landslides and gully heads respectively, was assigned with the following criterion: if the landslide or gully heads locations intersected a given slope unit, then we assigned a Fig. 1. Panel (a) shows the broad geographic context and location of the study area. Panels (b), (c) and (d) summarize the spatial distribution of Neolithic sites, landslides and gully heads, respectively. The three inventories are reported in terms of their density over space, this being computed via a 2D density kernel with a radius of 5 km. Black points are the locations of each element in the three inventories, associated with each panel. Panel (e) shows an example of a buried (along the black solid line) Neolithic site threatened by geomorphic processes acting along the slope.
presence condition (or 1) and vice-versa for the absence case (or 0).

Covariates
We start from a large set of morphometric and thematic covariates because by planning to implement a latent field in the model (see §3.3), we aimed to avoid capturing effects for which we could have a direct interpretation. All morphometric covariates originate from a LIDAR survey of the study area, from which a high resolution DEM was built and resampled at 5 m. The only morphometric exception corresponds to the SU areal extent. Thematic properties correspond to the CORINE Land Use/Cover (e.g., Feranec et al., 2007) and a soil map from local pedological studies at 1:10,000 scale. Both the Land Use and Soil classification are detailed in Appendix 1. In addition to those, we computed the Euclidean distance from each streamline to the centroid of the squared lattice coinciding with the 5 m DEM resolution.
Also, during fieldwork activities we noted that gully erosion almost systematically occurred in the proximity of the paths taken by local sheep herding freely through the landscape. Our interpretation was that the constant disturbance of the sheep may lead to incise the soil and compact it, hence locally modifying the hydraulic properties. We thought of mapping the actual paths but the resolution of satellite scenes nor available orthophotos was not sufficient to see clearly the paths. Overgrazing is known to contribute to the development of gully erosion in Iran (Shahbazi et al., 2017) and South Africa (Boardman, 2014). Therefore, as an alternative, we thought of mapping the locations where sheep are held and assumed that, if our observation is valid, then a regression model would estimate a negative relation between gully occurrences and the distance from sheepfolds. As a result, we also computed the Euclidean distance from each sheepfold to the centroid of the 5 m lattice.
Below we report the list of covariates (acronyms in italic) we used during the modeling phase: 1. Distance to sheephold (Dist2Sheep) 2. Distance to stream (Dist2Stream) 3. Eastness and Northness (e.g., Lombardo et al., 2018b) 4. Elevation 5. Planar (PLC) and profile (PRC) curvatures (Heerdegen and Beran, 1982) 6. Relative Slope Position (RSP) (Böhner and Selige, 2006) 7. Slope steepness (Slope) (Zevenbergen and Thorne, 1987) 8. Slope Unit area (SU area) (e.g., Amato et al., 2019) 9. Topographic Position Index (TPI) (Guisan et al., 1999) 10. Topographic Wetness Index (TWI) (Beven and Kirkby, 1979) 11. Soil classes in % per SU 12. Land Use classes in % per SU Among these, we chose Dist2Stream because water flowing along channels can undercut slopes leading to landslides and provide the essential water for soil erosion, hence for gullying. Eastness and Northness are known to be a proxy for strata attitude, thus they can represent a vital information for debris slides. And, they can also control exposition to the sunlight, hence to wetting and drying cycles that could contribute to shrink and swell for clays. In turn, these can open up cracks along the soil profile where gullies can develop. The Elevation can be a proxy for precipitation. Where topographic barriers exist, rainfall can be confined either on one or the other side of the divide. As for the two Curvatures, these have been shown to focus or diverge overland flows (see, Ohlmacher, 2007). The SU area can control the type of landslides. For instance, if the SU is small and narrow, then it is geomorphologically unlikely that shallow L. Lombardo, et al. Engineering Geology 277 (2020) 105776 debris slides can develop. The TWI is a property that should indicate the terrain tendency of retaining water, hence could provide relevant information for both processes, whereas the TPI can explain the location across the landscape where the two geomorphic processes can take place. Because of the larger polygonal structure of the SU partition, we adopted two criteria to summarize the distribution of the covariate set described above: for continuous properties we computed the mean and standard deviation values of the grid-cell distribution per SU. For categorical properties we computed the ratio between the extent of the given class within a SU and the extent of the SU itself (then expressed in percentage).
It is important to note that some of the covariates do not share the same resolution. This is an issue which is omnipresent in any susceptibility model and will remain present until technology would allow us to spatially co-register all environmental properties. Out of all the mapping units proposed in the literature, non-co-registered properties may be a problem specifically for grid-cells, and especially when positional errors affect the landslide inventories (Steger et al., 2016). However, for a SU partition, properties with different resolutions have a much more limited effect, if any at all. In fact, the SU act as a smoother over these problems because the user needs to scale (generally) up the resolution of the covariates to the SU resolution.
Finally, before using the covariates in the modeling phase, we also adjust their respective distribution in the same range. We do so by rescaling each covariate with mean zero unit variance (Paulson, 1942), a procedure which ensures that all the regression coefficients will be expressed in the same scale; enabling the comparison of their influence with respect to the whole multivariate model.

Binomial generalized additive modeling with INLA
A binomial GAM offers the possibility to model linear (or fixed) and nonlinear (or random) effects representing the functional relation between a set of explanatory (see §3.2) and dependent (presence/absence of landslides and gullies per SU, separately) variables (e.g., Brenning, 2005). Among the random effects one can include categorical covariates, or properties represented by discrete classes each one independent from the others, and model them as iid (or independent and identically distributed) effects. Additionally, one can model ordinal covariates, or properties represented by discrete classes which retain an ordinal structure (hence from small to large corresponding values per class). Ordinal covariates can be modeled accounting for adjacent-class-dependence by using a random walk function of the first order (RW1) between class estimations . Furthermore, one can model covariates expressed at a latent level, either in space, in time or both (Lombardo et al., 2019b). This concept in the geomorphological literature has been theoretically introduced in (Lombardo et al., 2018a). The idea behind it, and in the context of the present work, assumes that the covariates we commonly adopt in geomorphology may not express the whole variability in the landslide or gully distribution over space. In other words, there is virtually always an unexplained but still spatially structured component which is hardly accounted for in landslide susceptibility studies. This component or missed covariate or Latent Spatial Effect (LSE), can reflect the spatial signal of the trigger for event-specific inventories (Guzzetti et al., 2012) or any other unaccountable effects for which we do not have explicit data, in case of geomorphological inventories.
We note here that we opted to model both landslide and gully erosion susceptibility (separately) by using a LSE, three ordinal effects (SU Area, Slopeμ and TWIμ), whereas we considered all remaining covariates as fixed effects.
The INLA package (Rue and Held, 2005;Rue et al., 2009) of the programming language R (Team, 2000) offers a rapid and precise estimation for Bayesian statistics featuring latent Gaussian models (Rue et al., 2017). In case of landslide susceptibility models we can summarize the present models (more details are provided in Lombardo et al., 2019aLombardo et al., , 2019c, as follows: where, η is the logit link, P is the susceptibility either considered for landslides or for gullies, β 0 is the global intercept, β j are the regression coefficients estimated for each of the z j covariates, f SU Area , f Slope μ , f TWIμ are the three ordinal properties we selected and f LSE is the latent field we compute by using a Besag model (Blangiardo and Cameletti, 2015).
There are entire books dedicated to the calculations of latent covariates, both over space and time and a detailed explanation can be found for instance in Krainski et al. (2018). Here we try to provide a quick numerical and graphical explanation directed to a geo-scientific readership.
Any model, be it physically-or statistically-based produces a distribution of residuals between the prediction estimates and the observation. This concept is valid also for susceptibility models. In advanced spatial statistical applications, these residuals are evaluated over space. And, whether they are non-randomly geographically distributed their signal can be captured and re-integrated in the modeling procedure as a covariate acting at a latent level (e.g., Bakka et al., 2018). Simply computing the residuals and re-introducing them would inevitably lead to overfitting issues and misinterpretation of the results; for this reason, solutions have been proposed inspired by concepts such as smoothed residuals (Baddeley et al., 2005). In our case, being our space partitioned into SUs, we opted to implement a latent spatial effect which is driven as a function of neighboring mapping units. In practice, we compute a moving window through space which returns the average residual per SU from all the adjacent SUs (which share a border in the corresponding shapefile) and adds for each SU a random error component (ε). In this way, one can retrieve the effect of unexplained covariates and drive modeling performance strength from its inclusion in the model itself. This LSE approach is commonly referred as Besag (Gómez-Rubio, 2020) and it is assumed to be normally -N f~(0, ) LSE LSE -and multivariate normally -f LSE ~ mvn LSEdistributed.
To drive the residual spatial averaging window, we computed the adjacency matrix (Condon et al., 2002), which in our case is essentially a sparse matrix 4978 × 4978 (the squared number of SUs) with ones (or adjacency links) reported for SU that share a boundary. This is graphically summarized in Fig. 3.
Ultimately, we implemented a 10-fold cross-validation (CV) scheme for landslides and gullies separately. Each randomly selected subset is constrained to be sampled only once for a total of 10 subsets containing 10% of the dataset for validation (the complementary 90% is used to fit). We note here that because of the one-time sampling constraint, the integration of the 10 predicted subsets corresponds to the whole study area, whose susceptibility can be plotted in a fully predicted map. To estimate the model perfomance we used Receiver Operating Characteristic curves and their Area Under the Curve. More specifically, we used the AUC classification proposed by (Hosmer and Lemeshow, 2000).

Modeling performance
Based on modeling performances, both landslide and gully susceptibility models show outstanding (Hosmer and Lemeshow, 2000) predictive skills. The performance of the 10-fold cross-validation scheme is summarized in Fig. 4, where we report the Receiver Operating Characteristic (ROC) curves and their Area Under the Curve (AUC). As mentioned above, their predictive skill is quite high, with CV landslide susceptibility having a median AUC of around 0.91 and the same for gullies with a median value of 0.97. We recall here that these are predictive performances for which even the LSE is being re-estimated for each CV replicate. As for the robustness of the CV scheme, the inter-quartile distance in the AUC distributions is less than 0.2 AUC in both cases, indicating strong stability even for randomly extracted subsets. Such high predictive skills are surprising and one may intuitively think that the inclusion of the LSE is boosting the performance. However, during additional tests we have run, the median AUC of a CV landslide susceptibility without LSE (unreported results) is approximately 0.89 and the same for the gully case is approximately 0.86 (unreported results). This implies that for the landslide case, analogous estimates even in a simpler context can be achieved. As for the gully model, the latent component of the model retrieved more spatial information than it did for the landslide case. This could be due to the absence of a unknown but spatially-structured covariate in the gully susceptibility model, for which we could not account for in our initial covariate set. Fig. 5 shows the posterior distribution of covariates estimated to be significant at least in one of the two susceptibility cases under consideration. What stands out the most is that only four covariates are significant in the landslide susceptibility model whereas the gully one has a much larger number of covariates for which the model is at least 95% certain of their role. And yet, the difference in performance is negligible between the two, which implies that a gully erosion susceptibility model expressed at the SU scale, requires a much larger number of properties to explain the spatial distribution of gully occurrences. Conversely, the landslide case is well determined and the presence/absence distribution is well explained even with less predisposing factors.  . 4. Panel (a) shows 10 ROC curves for both landslide and gully susceptibility models. Each ROC curve correspond to one of the 10-fold cross validations we performed. Panel (b) summarizes the corresponding AUC distribution across the 10-fold CV scheme, for both landslides and gullies.

Random effects
Here we separately show the random effects of ordinal and latent covariates. Fig. 6 reports the ordinal nonlinear effects of SU area, Slopeμ and TWIμ. Two out of three covariates appear to be significant through their entire distribution for the landslide case. More specifically, the SU area mean effect increases almost as a logarithmic function where a negative regression coefficient is estimated for small SUs and a rapid increase in the regression coefficients can be seen up to SUs with an extent of about 0.8km 2 after which, the effect remains essentially the same. This is even more exacerbated for the Slopeμ case, where the nonlinear function appears more sigmoidal in shape, from small to large steepness values. As for the TWIμ, the effect is mostly nonsignificant although the mean contribution to the model is shown to be positive for small TWI values after which, it rapidly decays to negligible effects.
The situation for the gully susceptibility model is quite different. The SU area, on average, does not contribute to the model with the exception of very small SU for which the contribution is clearly negative. The Slopeμ effect is negative for the left and right tail of mean steepness distribution per SU, whereas intermediate mean slope values per SU between 5° and 10°. As for the TWIμ the whole distribution is not only nonsignificant but the mean value of the regression coefficient per class aligns with zero, which implies a negligible effect to the multivariate scheme.
As for the LSE computed for the two processes, Fig. 7 shows their posterior distribution in map form, highlighting some relevant differences. The LSE for the landslide susceptibility model is much smoother than its gully counterpart. And, its posterior distribution is confined between two relatively small minimum and maximum values, compared to the LSE for gully. This leads to two considerations. The landslide susceptibility is well explained even just by using fixed and ordinal random effects. And, that the gully counterpart is missing instead a much larger spatially-structured information. This does not come as a surprise because SU have a well established literature in the landslide context whereas this is the first attempt for the gully case.
Notably, the two LSE have significantly different distributions both in amplitude and geographic patterns. This is shown in Fig. 8 where the minima and maxima of the two dimensional space between the two LSE are very different as well as their respective relationship. This plot implies that the two geomorphological processes have very different residual distributions and also that the two processes themselves may be independent from each other, at least in the present study area.

Multi-hazard susceptibility mapping
Our modeling framework led to estimate two separate susceptibility models, one for landslide and one for gully erosion occurrence. This is geographically reported in Fig. 9, panels a and b, where the mean posterior probability spectrum of the two geomorphological processes has been reclassified according to 3 quartiles (hence mimicking the central box of a boxplot).
To complete the information on the posterior distribution of the susceptibility estimates, we also report the uncertainty around the mean susceptibility for the two cases, this being shown in Fig. 9, panels c and d. For a susceptibility model to be useful for any spatial planners or stakeholders it should reliably estimate the left and right tails of the probability distribution. This is a requirement because it is upon these two extremes that decisions are usually made, whereas the central portion of the probability spectrum can exhibit a larger variation (Rossi et al., 2010;Lombardo et al., 2016a). In other words the space defined between the mean susceptibility and its uncertainty should be  . 6. Random effects for SU area, Slopeμ and TWIμ. The distribution for each nonlinear function is summarized with its mean (solid line) and 95% credible interval (transparent background). Blue and red are colour-codes for landslide and gully susceptibility models, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) distributed in a bell shape. This is shown in Fig. 10 where the two susceptibility models well exhibit small uncertainties for the extremes of the probability range and much larger variations in between. Also, the gully erosion model shows a larger variability in the susceptibility estimates compared to its landslide counterpart.
Ultimately, we combined the two reclassified mean susceptibility maps in a unified multi-hazard one. This is done by taking the four susceptibility classes for landslides (L1 to L4) and gullies (G1 to G4) and by computing each of the 16 possible combinations -or using the combine function in any GIS environment-. As a result, we obtained 16 multi-hazard classes, from very low (L1-G1) to very high (L4-G4) susceptibility estimates in both cases. This is shown in Fig. 9e in map form, where the highest combined susceptibility mostly affects the southern sector of the study area; and it transitions to the lowest combined susceptibility northward with intermediate combined susceptibility conditions confined to the north-west corner.

Heritage sites at risk
We recall here that our goal was to assess the exposure of cultural heritage sites to the combined susceptibility of landslides and gully occurrences. For this reason, we have added an additional step where we intersect the combined susceptibility shown in Fig. 9e together with known locations of Neolithic sites in the study area. The result of this procedure is shown in Fig. 11a, where the spatial distribution of Neolithic sites is plotted with a palette corresponding to the combined susceptibility classes. Here they appear quite scattered with very few highly-susceptible clusters shown in the south-and north-western sectors. A graphical summary of the number of sites under threat is also reported in Fig. 11b, where, for instance, 12 Neolithic sites fall in the highly susceptible class for both processes (S4-G4).
As shown in Nicu (2018a) and Nicu and Asăndulesei (2018), gully erosion and landslides are a real threat for cultural heritage in the northeastern part of Romania. Currently, there are no mitigation measures or future management plans for what concerns the landslides in the area, except for the stone gabions erected in order to stabilise the north-western slopes of Oilor Hill (see, Nicu, 2018b). As for gully erosion mitigation measures, the most notable example is the one of Baiceni-Cucuteni gully. In the upper part of this gully, trees were planted and 14 concrete barriers were built to stop the erosion process. As a result, the erosion has considerably decreased (see, . Notably, there are no fast moving landslides within the study area. And, the way different types of geomorphic processes affect cultural heritage, both buried and built, has been discussed in (Nicu, 2017a). There, the Fig. 7. Posterior distribution of the Latent Spatial Effects for the landslide (panels a and c) and gully (panels b and d) susceptibility models. Panel (a) and (b) report the posterior mean of the two latent fields, whereas panels (c) and (d) report the 95% credible interval estimated as the difference between the posterior 97.5 and 2.5 LSE percentiles. Fig. 8. 2D density scatterplot of the posterior mean for the LSE belonging to the landslide (x-axis) and gully erosion (y-axis) susceptibility models.
author argued that, for this specific area, gully erosion is more destructive than landslides. (Nicu, 2017a) explains that in case of gullies, the archaeological remains are washed away and carried even for long distances through the small river valleys, while in the case of landslides, they are moved together with the landslide body towards the base of the hill. In such a complex environment, multi-hazard approaches represent a significant step forward in the prioritisation of risk reduction, especially for geo-archaeological applications (Sevieri et al., 2020).

Discussions
The modeling workflow we propose here is a unique case in the geoscientific literature for several reasons. Firstly, the vast majority of the community working with statistically-based susceptibility models implements a frequentist version of GLM. Conversely, here we do not only extend the study to the Bayesian case, but we also adopt a GAM, with random effect of ordinal and latent spatial nature. Secondly, we model more than one susceptibility in the same study, and we combine the two probability spectra to assess vulnerable Neolithic sites. We do this by adopting a SU partition which was tested so far only for landslides. We stress here that SUs may still not be the best spatial partition for gully erosion occurrence because gullies can occur in relatively flat conditions where the SU calculation itself struggles to find homogeneous aspect conditions. In this work, we tried to minimize this effect by generating high-resolution SUs; we recall here they have a mean extent of 0.1km 2 and a 95% confidence interval of 0.43km 2 .
Few interesting points should also be stressed in terms of covariate contributions. The landslide susceptibility model required a much fewer number of fixed effects (only four estimated to be significant) to produce outstanding results. Conversely, the gully erosion susceptibility estimated eleven fixed effects as significant. Conversely, ordinal covariates show an inverted situation where the landslide susceptibility Fig. 9. Panel (a) shows the mean landslide susceptibility together with its 95% credible interval in panel (c). Similarly, the mean gully erosion susceptibility is shown in panel (b) and its 95% level of uncertainty in panel (d). For the two mean susceptibilities (a) and (b), we have reclassified the probability spectra into four classes namely, low (S1-G1), moderately low (S2-G2), moderately high (S3-G3) and high susceptibility (S4-G4). This has been obtained by slicing the two probability distributions according to the three quartiles (τ = 0.25, 0.5 and 0.75). Panel (e) shows the combination of the four probability classes for landslides and gullies into a single multi-hazard susceptibility map. model draws explanatory strength from the three random effects, and the gully susceptibility only shows significant contributions from the Slopeμ. This is inverted again for the two latent spatial fields, where the landslide susceptibility model is quite stable and highly performing with or without the LSE. As for the gully erosion case, the LSE has a much larger contribution over the final probability estimates. Moreover, we tested our field observation that gully head started in close proximity to sheep tracks. We show here an example photo taken during mapping campaigns (see Fig. 12).
Moreover, even in the same statistical modeling setting, an ideal situation would have directly featured the distance to the sheep tracks. However, because we lacked high resolution satellite images, we could only compute the distance from sheepfolds, which we recognize to be an approximation of the phenomenon we observed. Nevertheless, the regression coefficient of the distance to sheepfolds appears to be significant and among the largest contributors (in absolute value) for the gully susceptibility model. And, it appears to be non-significant and with a much smaller coefficient in the landslide case. Notably, despite the significant contribution of this covariate to the gully erosion model, we do not claim this to be a cause for gullying. To make such a claim, further geotechnical tests are required to demonstrate significant hydrological and mechanical changes in soils compacted due to the continuous passage of sheep herds.
We recall here that the two geomorphological processes have been modeled separately. However, the literature suggests that landslide and gullies may influence each other. For instance, Lucà et al. (2011) explains that debris flows can occur in gullies thanks to the volume of entrainable sediments (Bovis and Jakob, 1999). They also mention that landslides can affect gully headcuts and banks steeper that the adjacent slopes (VanDine and Bovis, 2002). Also, Kukemilks and Saks (2013) notes that landslides themselves can form on the banks of newly developed gullies showing high correlation one another. For this reason, we are also testing more complex joint probability models where the two processes can be estimated together rather than separately. For instance, one can initially use one of the two processes as a reference, let us assume here to be the landslides. By adopting a suitable likelihood the spatial probability of landslide occurrences can be estimated. Then, a second likelihood can be used to measure the necessary deviation from the landslides to spatially estimate gullies (see, Hessellund et al., 2019). We see this as the direction to take in the future because joint probability models would ensure that reciprocal spatial influences would be properly accounted for. In the present state of the geoscientific literature, the implementation of such models is yet to be tested, likely because of their added complexity, and recent examples still consider multiple processes separately (Cama et al., 2020).
From the heritage sites perspective, from a total sample of 107 Neolithic sites, our combined susceptibility highlighted: i) 12 sites falling in the class with high probability for both processes (S4-G4), ii) 6 sites falling in the class with high probability for landslides and moderately high probability for gullies (S4-G3), iii) 7 sites falling in the class with high probability for gullies and moderately high probability for landslides (S3-G4), and iv) 16 sites falling in the class with moderately high probability for both processes (S3-G3). The full list of endangered cultural sites is included in Appendix 2.
Out of the 12 Neolithic sites that fall in the class with high probability for both processes (S4-G4), one (La Cruce in Fundoaia) is featured both in the List of Historical Monuments (LMI) and National Archaeological Registry (RAN). Out of the 16 Neolithic sites that fall in the class with moderately high probability for both processes (S3-G3), one site is listed in the LMI and two in RAN. Moreover, out of the 6 sites falling in the class with high probability for landslides and moderately high probability for gullies (S4-G3), one site is listed in RAN; and out of the 7 sites falling in the class with high probability for gullies and moderately high probability for landslides (S3-G4), one is listed both in RAN and LMI. This means that five sites of national importance may be exposed to both landslides and gully erosion. However, even if the remaining sites (102) are not listed as monuments or geo-archaeological sites of great importance at the national level, they are still culturally relevant as they give to the landscape a higher historical value. In fact, cherishing the cultural heritage gives to the local communities a sense of identity and belonging.

Conclusion
Our combined susceptibility to landslide and gullies offers an integrated view of geomorphological threats to cultural heritage sites in the study area close to the city of Iaşi (Romania). Usually research outputs rarely find a direct translation into real-world applications. The added value of this work resides in the multi-disciplinary components namely, geomorphology, statistics and anthropology, which make the direct translation of the results a natural consequence of our research. In fact, we are in the process of contacting local administrations to inform them of our findings and update the heritage site status in the Romanian registry.
More generally, this type of integrated assessment is hardly available both in the geomorphological and anthropological communities. We obtained it by using advanced spatial Bayesian statistical models which offer a higher flexibility, interpretability and performance compared to common approaches reported in the literature. However, two remarks must be made: 1. Because gullies are so active in the region, temporal spatio-temporal variations in the landscape itself as well as the gully and landslide occurrences can be expected. For this reason, one of the limitations of the present contribution is related to the pure spatial connotation of the model we have built. Therefore, we envision to temporally map the evolution of two geomorphic processes to produce two multi-temporal inventories. This database would support more appropriate space-time models necessary in such a dynamic landscape. 2. We combined the two susceptibility estimates only at the end of the modeling procedure. And, if from the geomorphological perspective this is already at the research forefront, from the statistical perspective further improvements may still be required. This improvement corresponds to joint probability models, where landslides and gully locations would be modeled together and their combined susceptibility would be the natural modeling output. Joint probability models are much more demanding in terms of computational resources when compared to what is required to compute the two susceptibilities separately. For this reason, we envision the investments of future efforts to build a tool based on joint probability for a fully-integrated and unified multi-hazard model.

Declaration of Competing Interest
None.

Acknowledgement
We are grateful to the Prut-Barlad Water Administration who kindly provided the LIDAR coverage of the study area making it possible to create a detailed inventory of both landslides and gully occurrences. L. Lombardo, et al. Engineering Geology 277 (2020) 105776 Appendix B. Endangered heritage sites Table B1 List of cultural heritage sites for which the combined landslide and gully erosion susceptibility indicates highly prone conditions for at least one of the two geomorphological processes. The projected coordinate system of Easting and Northing corresponds to the "Dealul Piscului 1970 Stereo 70"; Class is the combined landslide-gully susceptibility (see Fig. 9); LMI/RAN codes are the label for which some of the sites are reported in the Romanian archive.  Table B2 List of cultural heritage sites for which the combined landslide and gully erosion susceptibility indicates moderately high prone conditions for both geomorphological processes. The projected coordinate system of Easting and Northing corresponds to the "Dealul Piscului 1970 Stereo 70"; Class is the combined landslide-gully susceptibility (see Fig. 9); LMI/RAN codes are the label for which some of the sites are reported in the Romanian archive.