Landslide Susceptibility Mapping in the Vrancea-Buz ă u Seismic Region, Southeast Romania

: This study presents the results of a landslide susceptibility analysis applied to the Vrancea-Buz ă u seismogenic region in the Carpathian Mountains, Romania. The target area is affected by a large diversity of landslide processes. Slopes are made-up of various types of rocks, climatic conditions can be classiﬁed as wet, and the area is a seismically active one. All this contributes to the observed high landslide hazard. The paper analyses the spatial component of the landslide hazard affecting the target area, the regional landslide susceptibility. First, an existing landslide inventory was completed to cover a wider area for the landslide susceptibility analysis. Second, two types of methods are applied, a purely statistical technique, based on correlations between landslide occurrence and local conditions, as well as the simpliﬁed spatial process-based Newmark Displacement analysis. Landslide susceptibility maps have been produced by applying both methods, the second one also allowing us to simulate different scenarios, based on various soil saturation rates and seismic inputs. Furthermore, landslide susceptibility was computed both for the landslide source and runout zones—the ﬁrst providing information about areas where landslides are preferentially triggered and the second indicating where landslides preferentially move along the slope and accumulate. The analysis showed that any of the different methods applied produces reliable maps of landslide susceptibility. However, uncertainties were also outlined as validation is insufﬁcient, especially in the northern area, where only a few landslides could be mapped due to the intense vegetation cover. H.H.; H.-B.H.; H.H.; H.H.; H.H., H.-B.H.; H.-B.H.;


Introduction
Landslides cause catastrophic events every year in Europe [1,2] when affecting a high number of exposed elements [3]. Landslide susceptibility (addressing only the spatial component) and hazard (analysing both spatial and temporal information), as such not targeting the risk/danger component, have been studied at various scales: at a global scale, e.g., by Nadim et al. [4], NASA [5,6] and the World Bank [7]; at the European scale, e.g., by [8][9][10]; at the national scale, e.g., in France, Turkey, Pakistan, Cuba, Costa Rica and Romania by [11][12][13][14][15][16][17]. In general, susceptibility models are considered to be more suited for larger areas evaluated at regional scales than for local scales, where other methods, such as numerical modelling, are better adapted to assess related hazards [9].
The occurrence of a landslide is determined by one or more triggering factors. These triggers (considering those potentially active in the target region of Romania) encompass the effect of water (precipitation and snowmelt), seismic activity or human activity in the landslide-prone area. The most common triggering factor is exceptional heavy rainfall that exceeds the level normally encountered in a given area [18]. Seismicity is another important landslide triggering factor, especially in the Romanian target region [19]. In fact, the area of interest has been exposed to both large historical earthquakes and changing climatic conditions that contributed to extensive slope instabilities [19]. Their combined effect has been considered for the present landslide susceptibility mapping study.
Landslide susceptibility mapping has become very popular since the development of efficient spatial analysis tools implemented on GIS platforms [20]. A milestone in the history of the application of landslide statistical analysis is the work published by Carrara et al. [21]. It provides a critical review of almost all relevant existing statistical methods. Generally, landslide studies are performed according to several hazard assessment models: expert or heuristic models, statistical models [22], statistical and geotechnical models [23,24], spatial multi-criteria evaluation (SMCE) [16], spatially explicit deep learning neural network models [25], machine learning models [26] and some of them also consider seismic microzonation models [27]. Regarding the landslide data input, it has been shown that the performance of a susceptibility model varies slightly after progressive updating and addition of landslide data over time [28]. However, estimates of landslide susceptibility remain largely insensitive to data changes over larger periods of time [28]. Spatial or temporal stratified sampling provides only minor variations in model performance. The results of Ozturk et al. [28] call for further testing of the concept of dynamic susceptibility and its interpretation in data-driven models, especially in the broader context of landslide risk assessment in the context of environmental changes.
Another aspect of scale that is very important for assessing landslide susceptibility is the choice of the appropriate map resolution. According to Schlögel et al. [29], finer resolutions of the chosen digital terrain model do not necessarily lead to higher predictive accuracy in landslide susceptibility mapping (at least below a reasonable resolution of e.g., 30-90 m), while the frequency ratio model seems to be optimal for the coarsest resolutions (i.e., 70, 80 and 90 m, when larger regions covering several watersheds are considered).
Landslide processes of different types represent a widespread geohazard in Romania, where they affect mainly hilly and plateau regions as well as mountain areas developed on flysch formations. Due to the highly active morphodynamics, the region of the Carpathian Curvature can be classified as one of the top landslide occurrence areas in Europe [9,30]. While performing accurate landslide susceptibility assessments over large areas (i.e., including continental territories [31]) remains a challenge, mainly due to the lack of relevant landslide inventories and monitoring, in Romania, recent studies have shown some progress in understanding landslide dynamics at different scales [16,17,32].
At the continental and national scales, given that the influences exerted by landslide conditioning factors are complex and diverse as a function of climate-physiographical parameters, to account for this variety, several authors (e.g., [8,11]) propose a differentiated susceptibility zonation at the regional level based on a previous subdivision of the territory according to the existing physiographical landscape (plains-hills-mountains) and the Köppen's climatic classifications [33]. In Romania, a first attempt to produce a landslide susceptibility map, based on a previous regionalization of the territory, showed that the criteria at the continental level were not suitable for a national scale because of their high generality and a new division of the territory was needed [16]. In this latter study, the territory of Romania was separated into distinct units of analysis based on specific criteria adapted to the national environment. In addition to physiographic and climatic elements (the latter being mainly associated with the different vertical zones), morphostructural divisions were considered to be of significance in imprinting regional peculiarities to the interaction of landslide conditioning factors.
The study area of the present work represents a part of the Buzău river catchment area, situated in the seismically very active Carpathian Curvature region of Romania. The Buzău catchment extends, in its upper and middle parts, over the Buzău Carpathians and Subcarpathians and, in its lower sector, over the Romanian Plain ( Figure 1; Figure 2). The current study is only focused on the upper and middle sectors of the basin, corresponding to the mountain and hilly areas and developing on a NW-SE direction. The aim of this work is to establish landslide density and susceptibility maps according to two different large-scale approaches. The first approach consists of a statistical calculation based on correlations between landslide occurrence and various environmental factors, while the second approach is based on the calculation of Newmark displacements (ND) calculated for earthquakes with magnitudes of 7.2 and 7.8 in combination with different thresholds of water saturation for the landslide formation. The results of the two approaches are discussed in the environmental context of the region, validated and compared. Figure 1. Seismicity map of Southeastern Romania and position of the study area: earthquake events marked in terms of moment magnitude M w (provided by NIEP, downloaded from the ArcGIS online server) and representation of the focal mechanisms of the events between 1943 and 2020 (USGS, 2021). Note, this map actually shows thousands of circles, each one corresponding to an event with a specific magnitude; their overlay locally produces the bluish colors, only the M w > 7 earthquakes appear as individual violet circles.

Geology and Morpho-Seismo-Tectonics
Due to its geological evolution, Romania is divided into two major types of tectonostructural units: the platform regions (Walachian, Dobrogea, Scythian and East European platforms) and the orogene units (the Carpathian orogenic chain and the North Dobrogea orogen) [34] (Figure 2). The intra-continental Carpathian Mountains formed as a result of continental collision by the gradual accretion and exhumation of continental material towards the foreland, driven by a foreland-coupling process. In the East and Southeast Carpathians, the topographic expression, inherited from the Late Cretaceous-Paleogene, was overprinted by the Miocene thrusting and exhumation due to an accelerating subduction. In this part, the exhumation gradually migrated towards the foreland during Pliocene-Quaternary and is still active [34,35]. Considering the above, from a geostructural point of view, the study catchment area overlaps a typical intra-continental plate collisional area with an intense tectonic activity. The litho-stratigraphic conditions are represented by major structures, oriented NE-SW and corresponding in the upper part to the outer flysch unit (Carpathian Mountains), composed of intensely folded and faulted Cretaceous and Paleogene flysch sediments and in the lower part to the Pericarpathian unit (Curvature Subcarpathian hills), made up of the less cohesive Mio-Pliocene molasse deposits. Only the extreme southeast of the study area belongs to the Romanian Plain ( Figure 2). The flysch formations consist of alternations of thick sandstones, shales, limestones with intercalations of marls, clays or black shales, while the molasse formations are built up of a heterogeneous mixture of clays, marls, sands, gravels, sandstones and shales, with intercalations of salt breccia, tuff and gypsum. In terms of morphometry and morphology, the low-to-mid-altitude Carpathians have elevations of 900-1700 m, a relative relief of 500-800 m and steep slopes with inclinations of 15 to 45 • . They are characterized by narrow valleys and continuous ridges. The Subcarpathian region has elevations of up to 900 m, a relative relief of 300-500 m and slope inclinations of 10-30 • , featuring large depressions and valleys, rounded hilltops and slopes almost entirely covered by colluvial deposits. The surface formations along the Buzău floodplain and in the lowland sector are of Quaternary age (alluvial and loessoid deposits).
Romania is a seismogenic province that has experienced violent and disastrous earthquakes mostly concentrated in the Vrancea, Făgăras , , Banat and Maramures , regions. The Vrancea region, to which our study area belongs, is the most active seismic region in Europe with magnitudes M w reaching 7.9 (according to ROMPLUS Catalogue, NIEP) [38]. The relatively recent catastrophic earthquake of 4 March 1977 (magnitude M w = 7.2) caused enormous damages in the Subcarpathian Curvature and the Romanian Plain. The earthquake considerably enhanced the morphodynamics in this area by reactivating old landslides and triggering new ones, some rivers being dammed by the displaced material (Zăbala River) [39,40].

Climate
Romania has a temperate-continental climate. The altitude of the Carpathian range and its central location blocks the circulation of air masses on both sides of the range, while also imposing a climatic vertical zonation [41]. In the study area, three climate zones from those listed in Table 1 are encountered, with mean precipitation amounts of over 500 mm/year. Important spring showers, often overlapping the snowmelt of late spring, and torrential summer rainfalls are responsible for causing landslide occurrences in the region [42]. Considering that the influence of predisposing factors such as lithology, slope and land use on slope stability varies according to climatic conditions [17] made a climatephysiographical zonation map of Romania ( Figure A1). The combination of Köppen's climatic classification (according to the National Meteorological Administration (NMA), 2008) with that of the mountainous zones resulted in a climate-physiographical classification ( Figure A1D). According to the latter, our study area is largely located in zone Z6-mountains under cold climate and, secondarily, towards the SE, in zone Z3-lowlands and hills under cold climate without a dry season and with a warm summer. Two small areas in the NW fall into the Z4 class of lowlands and hills under cold climate with a cold summer.
Climatic variables have a permanent and lasting effect on the behaviour of rock materials. Diurnal, monthly and seasonal variations in the air and soil temperatures, duration and intensity of precipitation, and quantitative variations in precipitation coupled with thermal variations are involved in the initiation of geomorphological processes [37].

Landslide Processes
The Curvature Subcarpathians are characterized by a large number and variety of landslides (shallow and medium-seated earth and debris slides and flows). They often associate with erosional processes, forming large complexes. In the Curvature Carpathians, common processes are high-magnitude deep-seated rock slides, compound and complex landslides, generally being old processes, considered periglacial or immediately postglacial, and showing periodic reactivations [43].
Several studies on landslide susceptibility have targeted the region at the Curvature of the Carpathians [44 -48]. Based on semi-quantitative or quantitative assessments, they have shown that larger areas classified as highly and very highly susceptible to shallow landslides characterize the Subcarpathian hilly unit, whereas in the mountainous sector, these classes are much more restricted and appear mainly along the rivers. As for deepseated landslides, the high and very high susceptibility classes are much more developed in the Carpathians. This situation conclusively reflects the typology of landslides in terms of magnitude [46]. Several landslide risk scenarios have been elaborated according to different factors: landslide-inducing precipitation under present and future climate conditions and earthquakes, such as the one of 4 March 1977 or one having a 100-year return period [49].

Data
One of the most important steps in establishing the susceptibility map is to carry out an inventory of landslides in the study area. This not only allows estimating the importance of the phenomenon in the study area but also to elaborate an updated landslide database. Thus, the inventory allows the preparation and evaluation of landslide density maps for each factor and subsequently to calculate the susceptibility index.
To differentiate between the effects of different parameters on the initiation of mass movements and their effects on the spreading of the material after failure, the landslide area is subdivided into two zones, namely, the initiation (source) zone and the accumulation (indicated as runout) zone [23,50], so as to be able to obtain differentiated results.
We used and completed the databases of Zumpano et al. [45] and Damen et al. [46]. The updating was based on the photo-interpretation of morphological landforms visible on high resolution satellite images (dating from 2009 to 2018) available on the Google Earth platform by digitising polygons of recent landslides still observable at a scale of 1:14,000. Subsequently, the landslide inventory was subdivided it into two parts ( Figure 3), each containing the two landslide areas, source and runout. One data subset, composed of 1982 landslides, with 1029 source zones and 953 runout zones, was introduced in the calculation of the density map, while the other database, composed of 409 landslides with 206 source zones and 203 runout zones, was used for validation.  [45] and Damen et al. [46] and were used for modelling and testing (about 10% of the landslides are reused for validation). Outside the rectangle, the landslides inventoried in the present study were used for validation. The landslide data have been used and are reproduced here with permission from the authors of both papers.
Sixteen parameters, predisposing, aggravating and triggering, were collected, processed and adapted to prepare them for the desired calculation process. While these data come from various sources at different scales (1/200,000, 1/5000, 1/14,000), we are aware that small-scale data may influence the model performance, but larger scale data were not always made available for the studied area. The following parameters were used: slope, geology, roughness, distance to the hydraulic network, topographic energy index (TPI), wetness, curvature, profile curvature, plane curvature, flow accumulation, internal friction angle (ϕ), cohesion (C), rock density (γ), distance to faults, distance to earthquakes and normalized difference vegetation index (NDVI).
To evaluate the dependency between these different explanatory variables, we calculated the correlation matrices for the source and runout (accumulation) areas, using the 'Band Collection Statistics' tool of the ArcMap platform (10.8.1). The tables in the appendix show the correlation matrices of the 16 parameters used (Tables A2 and A3). Due to collinearity reasons, in the case of source areas, we eliminated the factors slope, curvature and the friction angle (ϕ). In the case of runout zones, we removed slope, curvature, cohesion and friction angle. This leaves 13 parameters for the source and 12 parameters for the runout susceptibility analyses.

Statistical Susceptibility Map
To establish the susceptibility maps, we applied the approach of the landslide factor analysis according to Vanacker et al. [23,24] but by including more diverse landslide factors. On the ArcMap platform, we transformed the shapefiles (polygon) of the landslide inventory into raster files (15 m grid). Thus, the landslide pixels inform about the presence of the two landslide components that were mapped: One raster map resulted for the source parts and a second raster map for the runout parts. These two specific landslide raster maps were then used to extract information from the other raster layers featuring morphological, geological, climatic and seismotectonic factors. Finally, we obtained a database of pixels subject to landslides. This allowed us to perform statistical analyses of the factors [24,51].
To calculate the density of landslide source or runout zones in a given class of a parameter (factor) within the catchment, we used the following formula [24,51]: where d i is the normalized landslide density in each class i, N LPC is the number of pixels affected by landslides in each class of the environmental factor, N LPT is the total number of pixels affected by landslides within the environmental factor map, N PT is the total number of pixels within the environmental factor map and N PC is the number of pixels in each class of the environmental factor map. The results of this equation are presented in Table 1 for each environmental factor used. The susceptibility maps for the prediction of source and runout (accumulation) zones were obtained by calculating the geometric mean of the factors used according to the following formula [24]: where Sc is the susceptibility, Sc i the density of factors i and n the number of parameters.

Newmark's Method
The Newmark method [52] is based on a simple model of the sliding of a block on an inclined plane. It aims to calculate the sliding distance, the ND, and is based on a double integration scheme of a time history of acceleration. For this purpose, only values above a certain value, i.e., the critical acceleration and the threshold acceleration necessary to initiate sliding, are taken into account [24]. Several simplified formulae have been proposed that relate the ND to the Arias intensity and the critical acceleration. The formula of Miles and Ho [53] is as follows: where D is the displacement (cm), Ia the Arias intensity (m/s) and a c the critical acceleration (m/s 2 ). The critical acceleration is given by the following formula [53,54]: where g is the gravitation (m/s 2 ), α is the slope ( • ) and FS is the safety factor calculated according to the formula using the simplified Janbu's method (original: in [55]): where C is the cohesion (MPa), ϕ is the angle of internal friction ( • ), t is the total thickness of the layer (m), m is the thickness of the saturated layer, γ w is the unit weight of the material (kg/m 3 ) and γ is the unit weight of water (kg/m 3 ).

Arias' Intensity
The effect of earthquakes on the occurrence of mass movements is not negligible. This effect is introduced by the energy that propagates in the ground during earthquakes inducing complex periodic displacements that can unbalance the slopes. This displacement is measured by Arias [55], who gives the degree of ground shaking. This concept has been applied to the assessment of earthquake-induced landslide risk [53,54].
The seismic correlation of earthquake-induced slip with this intensity measure was first performed by Keefer and Wilson [56][57][58] developed a new generic formula that gives the Ia and the Cumulative Absolute Velocity (CAV): (6) where F InSlab is a dummy variable taking the value 1 for events occurring on the slab and 0 otherwise, M w is the magnitude, R is the distance, H is the depth of the hypocentre and V S30 is the shear velocity (of S-waves) over a depth of 30 m. The coefficients are given by the author in Table 2. Table 2. Repression coefficients for the calculation of Ia (data from Foulser-Piggott and Goda [58]).

Coefficient
Ia (m/s) The earthquakes with a magnitude greater than 7.0 that have struck the region around the studied catchment area since the year 984 (a period of 1036 years) are presented in Figure 4 (from [38]). Inside the catchment area, only earthquakes of magnitude M w = 7.1 have occurred, while in the NE of the catchment one earthquake of magnitude M w = 7.9 (maximum magnitude recorded over a period of 1036 years, in 1802) has also been recorded.  (Table 3).

Evaluation of Model Results
To evaluate both the fitting and the prediction performances of the obtained susceptibility models, we used the ROC (receiver operating characteristic) curve, which measures the performance of a binary classifier that classifies items into two distinct groups based on several factors. In our case, the two classes are represented by the values 1-equivalent to true, i.e., landslide presence pixels, and the value 0-equivalent to false, i.e., landslide absence pixels.
The ROC curve provides both a graphical view and a relevant measure of the performance of a classifier [59,60]. Although several evaluation methods are available in the literature, the ROC curve is widely used to evaluate the outcome of landslide susceptibility maps [60][61][62]. The skill of the model and its discrimination ability are measured by the height of the curve above the diagonal line as well as by its steepness.
As is recommended [63], in addition to testing the statistical performance of the models, in a second phase, the geomorphological meaning and reliability were analysed in relation to the regional characteristics.

Susceptibility Zonation Maps
Two landslide susceptibility zonation maps were obtained: one for the source zones  The evaluation of the fitting performance of susceptibility zones was conducted by calculating the ROC curves related to the calibration landslide dataset, while the prediction performance was tested by calculating the ROC curves in relation to the validation landslide dataset. The evaluation of the density map for the initiation zones (source) shows a ROC curve ( Figure 5D) with an associated AUC of 0.84. The evaluation of the density map for the runout areas shows an upwardly convex ROC curve ( Figure 5D) deviating from the straight line (45 • ), with an area under the curve (AUC) of 0.8. Both AUC values are greater than 0.7, which is, in practice, the value that justifies the validity of a classifier [64]. The ROC curves plotted in relation to the validation inventory ( Figure A2) show AUC values of 0.79 and 0.75 for the source and the runout maps, respectively (values above the threshold for model validation [64]). These results show that the model has a valid and acceptable discrimination power. Indeed, the accuracy of the spatial distribution of susceptibility related to the source zones is higher than that related to the runout zones. This difference can be explained by the uncertainty in digitizing the boundaries of the runout zones, as they are less evident on the satellite images than in the case of source zones, which are clearly visible through the escarpments. This uncertainty could introduce non-runout pixels into the runout zones database and subsequently significantly influence the result.
The map of landslide source susceptibility (based on the landslide source zones; Figure 5A) shows a strong influence of the geological formations. High and very high susceptibility values are rendered for large areas in the Carpathian sector corresponding to the folded and faulted external (Paleogene) flysch units, in particular to the shale-sandstone flysch deposits and secondarily to the sandstone flysch. The southern border of this highsusceptibility area follows the limit between the mountains and the Subcarpathian hills. Such a distribution correlates with a predominance of large dormant landslides (deepseated rock and debris slides and compound or complex landslides) featuring sectors of recent reactivations, partially located in the scarp areas. They are considered to be strongly connected to structural and tectonic characteristics, e.g., conditioned by faults which act as structural weaknesses during earthquakes [30]. The susceptibility values tend to decrease towards the north and northwest, where the intramontainous depression of Intorsura Buzăului lies. The inner Subcarpathian hills, underlain by intensely folded and faulted structures of Miocene molasse deposits, feature moderate landslide initiation susceptibility, although this unit is known to have an extremely large concentration of landslides ranging from slides to earth and debris flows. This situation could rather reflect the distribution of one certain landslide typology in terms of magnitude [46], namely, deep-seated movements, which register a much lower prevalence here than in the mountainous sector, pointing to a larger share held by deep-seated movements as opposed to shallow and medium-seated ones in the calibration landslide source dataset. Low and very low susceptibility areas are yielded for the south and southeast, in correspondence to the Pliocene-Quaternary strip of the outer Subcarpathian, as well as to the large floodplain and lowland areas. The Pliocene-Quaternary homocline relief, built on thick deposits with a high content of sands and gravels, features an interaction of landslides with erosional processes, the latter of which become slightly more predominant towards the exterior.
The susceptibility to landslide spreading (based on the landslide runout zones; Figure 5B) is highest in the median part of the basin and is elongated on a NE-SW direction parallel to the geological structures, located on the shale-sandstone flysch. High values also correspond to other mountain areas, where the lower slopes sectors are especially covered by thick landslide deposits, periodically affected by reactivation due to river deepening or forest cuttings. They are also encountered within the Miocene molasse deposits carved on less cohesive clays, marls or loose sandstones with clayey-marly intercalations, commonly subject to earth flows. Low and very low susceptibilities to landslide runout are registered along the outer Subcarpathian strip, made of Pliocene-Quaternary deposits mainly represented by sands and gravels, as well as in the lowlands.
The difference between the two susceptibility maps ( Figure 5C) shows that the susceptibility to landslide accumulation is higher, i.e., indicative for a prevalence of flow-type of processes, in the inner Subcarpathian hills, where earth and debris flows are abundant. In addition, mountainous areas built of softer rocks (shale flysch, shale-sandstone flysch) are more susceptible to flows (where, indeed, typically mudflows, debris flows, compound and complex landsides are found). The central higher parts of the mountains, made-up of sandstone flysch, are clearly marked by higher source susceptibility values-which is consistent with the development of more massive rock slides featuring large scarp areas. The higher source susceptibility in the upper slope sectors, with massive landslides originating preferentially from the mountain crests, hints at the influence of seismic triggering of slope failures near the mountain tops that are exposed to stronger earthquake shaking.

Newmark Displacement Maps
According to the ROMPLUS catalogue [38], the hypocenter depths of 7.1 magnitude (M w ) earthquakes vary between 100 and 150 km, although most are 150 km deep, as is the case of the M w = 7.9 earthquake, the single biggest magnitude earthquake in the region. To investigate different possible scenarios of these extreme cases, we calculated several models considering the following assumptions: For the whole basin, the thickness (t) of the layer subject to sliding was set to t = 10 m; the saturation rate m = 0 stands for totally dry conditions, m = 0.5 for half saturated conditions (half of the sliding layer, i.e., 5 m, is saturated) and m = 1 for 100% saturation (fully saturated over the 10 m sliding layer thickness). For the scenario calculations we also considered two magnitudes, namely, M w = 7.1 and M w = 8.
Eighteen maps of ND as a function of the saturation rate (m), magnitude and depth of the hypocentre were produced (Table 3). These 18 maps (Figures 6 and A3) show that the most important displacements are related to the 100% saturation (m = 1, Figure 6). The map with the largest ND is the one considering a magnitude M w of 8 and a depth of 100 km with a maximum displacement of 50 cm ( Figure 6B), while the one with the smallest ND is the map assuming a magnitude M w of 7.1 and a depth of 200 km ( Figure A3C). We consider as interesting maps the extreme case of magnitude M w = 8 and depth 100 km (even if, in the catalogue, earthquakes of this order have depths of 150 km) and that of magnitude M w = 7.1 at a 100 km depth, often produced in the basin. To validate these two maps, we performed a statistical analysis by calculating the normalized probability for the runout and source landslide zones separately and then for the two grouped together. This analysis was conducted on the landslide inventory used for validation and also on the landslide inventory used for the calibration of the susceptibility maps (Figures 7 and 8).  The comparison of the 2 maps ( Figure 6A,D) shows that the ND in the M w = 8 map is 10 times greater than that in the M w = 7.1 map. The maximum values are located in the centre of the basin within the Carpathian Curvature, generally located in the concave surfaces of the mountainsides. These can focus the seismic energy and amplify the amplitude (site effect), thus enhancing the potential for mass movement triggering as compared to other locations with convex topographic surfaces.
For the magnitude M w = 7.1 map model ( Figure 6A), the normalised probability of landslides (of any type) is directly related to the ND. It increases proportionally with the ND (Figure 7A,B). The runout zones behave differently from the source ones, as the runout zones show a logarithmic curve, with a steep slope at the beginning followed by a very low slope plateau. The source zones show a stable plateau at the start and then a logarithmic behaviour. For the calibration landslide dataset, the pattern is the same as for the validation landslides, except that, this time, the probability of the runout zones is higher than that of the source zones. The normalised probability of the two grouped inventories ( Figure 7C) illustrates a similar trend between sources and runouts from the 0.1 cm displacement. The general trend for the combined runout and source zones of both databases shows a rapid increase in the normalised probability between displacements 0 cm to 0.2 cm and then another stage with a small increase according to the following logarithmic law ( Figure 7D): The M w = 8 map ( Figure 6D) also has curves similar to the M w = 7.1 map, logarithmically increasing the normalised probability with ND ( Figure 8). For all landslides in both databases, the normalised probability varies according to the following logarithmic law ( Figure 8D): The difference between the normalised probability of the models for the simulated M w = 7.1 and M w = 8 (represented in Figures 7 and 8, respectively) is equivalent to the maximum displacement values of 1.2 cm for the M w = 7.1 scenario and of 30 cm for the M w = 8 scenario. Additionally, according to Equations (7) and (8) for the same displacement (x), we have the relationship between Y 7.1 (normalised probability of the M w = 7.1 model) and the Y 8 (normalised probability of M w = 8 model): This shows that for the same displacement value, the normalised probability of the magnitude M w = 8 model is greater than the normalised probability of the magnitude M w = 7.1. So, the difference of 0.9 magnitude units causes this large influence on the normalised probability.

Discussion and Conclusions
Landslide susceptibility zonation maps are generally calculated based on several environmental factors that authors consider to be the most important according to the climate-morphostructural and geodynamic contexts and the data available at the time of the study. In the literature, there are no standard environmental factors (norms) that must imperatively be used to establish a susceptibility zonation, especially as the calculation approaches are diverse. It is a common practice to search for correlations among these factors in order to reduce their number and eliminate redundant information. For a given model and apart from the mathematical means and techniques of validation which can 'show' satisfactory results (acceptable errors), 'expert' validation in the field is indispensable for the adjustment and correct interpretation of an obtained model.
The comparison of susceptibility zonation maps and ND maps can only be achieved by understanding the meaning and/or physical concept of the magnitude given by each map. The ND maps indicate the maximum displacement of the ground during an earthquake: the larger the displacement, the higher the susceptibility for landslides. However, it is difficult to establish a threshold (minimum value) of ND that can trigger a landslide if it is associated with other factors. Even if we have data on landslides caused by some earthquakes, we do not really know if the displacement generated by these earthquakes is the minimum that was needed to cause the trigging or if there was an additional contribution which was responsible for the observed failures. Even if the ND is a physical quantity, it is based on a statistical calculation, so the interpretation of displacements has a relative aspect. Moreover, even if we have the maximum displacement in several places, those places do not contain any landslide occurrences, which indicates that there are other factors that come into play in addition to the ND. In the case of susceptibility zonation maps, these represent the total influence of environmental factors according to their spatial distributions. This concept, which is different from that of ND maps, allows a high number of factors that are deemed useful to be considered and introduced into the model (in this study we used 16 factors). The comparison of the two types of maps helps combine the advantages of a statistical method with those of a process-based one.
The comparison between the susceptibility zonation maps and the ND maps is carried out firstly by spatial correlation between the four maps (Table 4; Table 5). Secondly, the comparison is made visually by looking at the patterns of the maps according to the geographical areas (Figure 9). The two ND maps (magnitudes M w = 7.1 and M w = 8) are very well correlated with each other, which is logical as the only difference in the calculation of the two maps is the energy of the earthquake. Yet, there are also differences: First, it can be noticed that the landslide source zonation can be better correlated with the ND maps (M w = 7.1 and M w = 8) than the accumulation zonation. The five locations that show a marked difference between the maps of runout susceptibility and ND generated by magnitude M w = 7.1 ( Figure 9B,C) are Dălghiu, NW of Barcani, SE of Pătârlagele, Murgeşti and S of Pârscov. The source zonation shows the lowest susceptibility in these areas while the ND is large, especially at Dălghiu (1.8 to 2.5 cm) ( Figure 9A,C). The highest susceptibilities that correspond to the largest ND are distributed along a SW-NE direction, passing through Nehoiu, Gura Teghii and Lopătari towards Secuiu. Further analysis using unmanned aerial vehicles should focus on these locations in other to better understand the local environmental conditions with up-to-date very-high-resolution imagery and elevation models [65]. The Buzău floodplain in the plain unit shows low values of both ND and susceptibility. To summarize, the zonation of landslide initiation susceptibility ( Figure 9A) shows a distribution more similar to the ND ( Figure 9C), with the difference being especially marked in the mountains (Dălghiu, Intorsura Buzăului and west of Barcani). The direction of this agreement between the two maps correlates well to the main structures. This is consistent with the fact that the landslide source prediction more closely reflects the spatial distribution patterns of deep-seated landslides, many of which show a high probability of being linked to a seismic cause (e.g., large scarps close to the ridge-tops, vicinity to fault lines, accumulation sectors with no apparent interaction with the river network) [43].  Table 5. Correlation matrix for ND maps (for magnitudes M w = 7.1 and M w = 8 with 100-km depth) as well as source and runout susceptibility maps (SS and AS). One uncertainty affecting the output landslide source and accumulation susceptibility patterns stems from the spatial incompleteness and heterogeneity of the landslide inventory used in the analysis. This dataset is mostly constrained to the area mapped by Zumpano et al. [45] and Damen et al. [46], hence it does not cover the entire study area. Completion of the landslide data outside this area was rather difficult. Especially in the northern sector, landslide detection was rendered almost impossible due to the increased surface covered by forests. Since the ratio of areas corresponding to factor classes within the landslide sample area cannot be the same as for the whole study area, calculating landslide densities in relation to entire class areas is suspected to produce biased outputs. This situation applies, for instance, to the lithological class 'sands with marls and clays, locally gravels' present especially in the SE of the study area, for which almost no information on landslides could be extracted from the mapped sample. Nevertheless, statistically, this translated into a 'zero density', thus cancelling the predisposition for landslide runout susceptibility throughout the study area where the respective lithology appears. Future improvements of the predictions are expected if either using a landslide database with a better coverage across the study area or enhancing the representativeness of the landslide sample, restricting the susceptibility computations to the sample area, and finally extrapolating the results to the entire study area.
Additionally, considering the differences in the distribution of shallow versus deepseated landslides, typologically differentiated landslide initiation and runout susceptibility assessments are envisaged for a future work. This could have an impact when comparing the deep-seated landslide source susceptibility map with ND maps, since it was shown that numerous such processes have a morphology showing a potential involvement of the seismic factor [43].
In order to improve the ND map, since the parameter 'm' (water saturation) has a significant influence on the result, different values of 'm' should be assigned to the different lithologies, considering that permeability and water retention vary according to the rock type, which controls the saturation rate. Adding a general 'm' value to all pixels of the study area does not allow to 'introduce' the differentiated effect of water specific for the different lithological formations. This parameter also varies according to the season, so future work will aim to establish separate maps for the wet period (m = 1, typically for the months of March-June and October-December) and for the dry period (m = 0-0.5, corresponding to June-September, January-February), while also taking into account the variation of 'm' with geology.
In addition, the parameter 't' (thickness of the layer, fixed at 10 m in this study), used in the calculation of the safety coefficient (above), is a parameter that can highly vary over the mapped area. However, defining the relative thickness of landslides for different slope conditions, geological conditions, proximities to rivers and vegetation covers is difficult. Possibly, some remote-or close-sensing-based analyses could help establish some statistical links between landslide thicknesses and local factors influencing slope instability.
Moreover, one could also consider that instead of purely comparing the single ND map with the statistical landslide susceptibility zonation maps, it could be interesting to introduce the ND map as one of the environmental factors used for the calculation of the susceptibility zonation map. This could be achieved in the frame of a more detailed future study.
Finally, it might be helpful to introduce the effect of differential weathering caused by climatic conditions as an environmental factor in the landslide susceptibility mapping procedure or via the 't' factor in the ND map calculation. Indeed, the ageing of the rock implies its alteration on the surface, which controls its mechanical weakening and consequently its susceptibility to ground movements. The same lithology does not necessarily have the same physical properties-which evolve over time and change with the slope's exposure (orientation with respect to wind and sun radiation effects) due to variable fracturing, the leaching rate by meteoric water, the presence of groundwater, etc. However, collecting related information over such wide areas represents a major challenge.