Landslide hazard and susceptibility maps derived from satellite and remote sensing data using limit equilibrium analysis and machine learning model

Landslide susceptibility mapping and landslide hazard mapping are approaches used to assess the potential for landslides and predict the occurrence of landslides, respectively. We evaluated and tested a limit equilibrium approach to produce a local-scale, multi-temporal geographic information system-based landslide hazard map that utilized satellite soil moisture data, soil strength and hydrologic data, and a high-resolution (1.5 m) LiDAR-derived digital elevation map. The final multi-temporal landslide hazard map was validated temporally and spatially using four study sites at known landslide locations and failure dates. The resulting product correctly indicated low factor of safety values at the study sites on the dates the landslide occurred. Also, we produced a regional-scale landslide susceptibility map using a logistic regression machine learning model using 15 variables derived from the geomorphology, soil properties, and land-cover data. The area under the curve of the receiver operating characteristic curve was used for the accuracy of the model, which yielded a success rate of 0.84. We show that using publicly available data, a multi-temporal landslide hazard map can be created that will produce a close-to-real-time landslide predictive map. The landslide hazard map provides an understanding into the evolution of landslide development temporally and spatially, whereas the landslide susceptibility map indicates the probability of landslides occurring at specific locations. When used in tandem, the two mapping models are complementary to each other. Specifically, the landslide susceptibility mapping identifies the area most susceptible to landslides, while the landslide hazard mapping predicts when landslide may occur within the identified susceptible area.


Introduction
Landslides are common and dangerous natural hazards that occur worldwide, often causing severe direct impacts on human lives, public, and private properties, and critical infrastructure (Klose et al. 2014). Direct costs of landslides in Kentucky are conservatively estimated to be between $10 and $20 million annually (Crawford 2014;Crawford and Bryson 2017). Therefore, it is imperative to identify the landslide susceptible areas to avoid or mitigate the possible damage. Landslides occur when the shear stress along a failure plane within the geologic materials of a slope exceeds the shear strength of the material. Landslides are commonly triggered by intense short period rainfall or prolonged rainfall, earthquakes, or human activities. Landslide susceptibility mapping (LSM), which describes the spatial distribution for the probability of landslide occurrence in each area according to the geographical sett, is considered a common countermeasure for mitigating the effects of landslides (Huang and Zhao 2018;Merghadi et al. 2020). More simply put, susceptibility can be defined as the probability of spatial occurrence of slope failures, given a set of geologic and climatological conditions (Guzzetti et al. 2005).
Several researchers (e.g., Guzzetti et al. 2005;Alvioli et al. 2016; Reinchenbach 2018) have suggested a variety of approaches and methods for landslide susceptibility in different geological and climatic settings. Within the last several years, there is an increasing preference toward machine learning methods (Di Napoli et al. 2020;Pourghasemi et al. 2021;Tehrani et al. 2022). The more common methods for landslide susceptibility modeling include logistic regression, neural network analysis, data-overlay, index-based, and weight of evidence analyses (Reichenbach et al. 2018). While LSM provides the spatial distribution of landslide probability, landslide hazard mapping (LHM) gives the probability that a landslide of a given magnitude will occur in each period and in a given area (Segoni et al. 2018). LHM presents the "when" and "where" landslides are likely to occur, combining spatial and temporal distribution of the landslide susceptibility (Guzzetti et al. 2005).
Rainfall thresholds are the most widely used triggering factor in landslide hazard assessment and early warning tools (Segoni et al. 2018;Monsieurs et al. 2019). However, in many cases, early warnings based solely on antecedent rainfall is not adequate because antecedent soil moisture conditions play a crucial role in the initiation of landslides (Zhuo et al. 2019). Some researchers (Zhuo et al. 2019;Wicki et al. 2020;Guzzetti et al. 2020) have investigated the soil moisture threshold for triggering landslides directly, accounting for the temporal evolution of soil moisture before and after slope movement. These studies mainly used machine learning models to investigate the correlation between the soil moisture threshold and the onset of landslides. However, these types of approaches do not consider the statics and fundamental mechanics associated with landslides (Di Napoli et al. 2020;Wicki et al. 2020;Guzzetti et al. 2020). Approaches that utilize fundamental mechanics, such as limit equilibrium, are required to better understand landslide occurrences. However, the use of a limit equilibrium equation to create and verify LHM at known locations and known dates is very minimal. The aim of this paper is to produce local-scale landslide hazard mapping (LHM) using a limit equilibrium approach (Lu and Godt 2008) that utilized publicly available satellite and remote sensing data, and to compare this mapping to LSM that is produced using a machine learning technique. A LHM model provides insights into the temporal and spatial evolution of landslides, while LSM models the probability and location of the landslides. This paper presents a two-part approach, where satellite and local site geologic data were combined to create a LHM model and a LSM model. A machine learning technique was used to create the LSM using geomorphic variables derived from an aerial LiDARbased 1.5-m digital elevation model (DEM), soil property data from Natural Resources Conservation Service (NRCS) Web Soil Survey (WSS), and land-cover data in northern Kentucky area where the Kentucky Geology Survey (KGS) has an extensive landslide inventory. The produced LSM was verified against the known landslide occurrences yielding an accuracy rate of 84.1%. The LSM is based on relatively static parameters at failure conditions (i.e., the onset of a landslide). Therefore, the LSM represents the probability of the failure at the worst case of scenario.
A limit equilibrium factor of safety equation for a hillslope with transient infiltration conditions (Lu and Godt 2008) was used to produce the LHM with 15.2-m spatial resolution using soil moisture data from the National Aeronautics and Space Administration (NASA) Soil Moisture Active Passive (SMAP) mission, slope angles derived from 1.5-m LiDAR DEM, and soil physical properties from WSS. The LHM was verified temporally and spatially against the landslides reported in the KGS landslide inventory. The LHM on the date with the lowest factor of safety values was compared to the LSM produced by a machine learning approach to assess how the factor of safety values compared to the probability of the landslide susceptibility map.

Area geology
Several landslide sites in northern Kentucky comprised the study area for this investigation. The north-central area of Kentucky is characterized by weathered limestone bedrock of the Ordovician geologic period that has been pushed toward the crest of the Cincinnati Arch and is often exposed at the ground surface (McGrain 1983). This area is defined as Outer Bluegrass physiographic region of Kentucky. The defining geologic elements of the Outer Bluegrass region consist of late Ordovician and Silurian-age limestones, dolomites, and shales (McDowell 1986), deep valleys formed by erodible shale that exists above non-erodible rock types, rich soils for agricultural uses, and gentle rolling hills formed by slopes ranging from 20 to 30%. Topographic relief averages approximately 150 m, ranging from steep slopes (i.e., slope angles between 20 and 35°) along the Ohio River, gently sloping uplands (i.e., slope angles < 6°) and broad dissected valleys. Shaly bedrock in the region weathers easily and produces thin to thick, clayey colluvial soils. Landslides typically occur within the colluvium or along the colluvial-bedrock contact (Crawford and Andrews 2012). Geotechnical reports from highway projects in northern Kentucky show that the average depth to the bedrock in the hillslope soil is about 4.6 m.

Site locations
Landslide data used for this study were obtained from the Kentucky Geological Survey (KGS) Landslide Inventory (https:// kgs. uky. edu/ kgsmap/ helpfi les/ lands lide_ help. shtm). The landslides reported in the inventory are primarily shallow colluvial slides. The inventory contains information about the types, dates, locations, and sizes of the landslides. Northern and eastern Kentucky are regions where landslides are most prevalent in Kentucky. Thus, this study concentrated on these regions of the State. Campbell County and Kenton County were arbitrarily chosen to represent northern Kentucky, and Johnson

3
County and Pike County were arbitrarily chosen to represent eastern Kentucky. A portion of the landslide sites in Campbell County and Kenton County was used as training points for the machine learning model to create the LSM model. Other sites in Campbell County not used as training points were used as validation sites. For the LHM development, two sites in Campbell County (identified in the KGS Landslide Inventory as Site 6501 and Site 6294) and two sites in Johnson County and Pike County (identified in the KGS Landslide Inventory as Site 6396 and Site 6458, respectively) were used to validate the results of LHM. Figure 1 presents the locations of the landslide study sites, relative to the State (shown in the inset).
The terrain information used to develop Fig. 1 was obtained from the 1.5-m LiDAR DEM in the study area that was flown in 2011 and 2012. Sites 6501 and 6294 (Fig. 1b) are in are about 180 m apart. The landslides occurred on 3/30/2018 and 2/04/2016, respectively. The slope angles in the vicinity of the landslide points range between 5 and 38°, with the slope surface partially covered with trees. The landslide runouts at these sites resulted in debris that blocked the road downslope of the sites. Site 6396 (Fig. 1c) located in Johnson County was reported on 1/24/2017. The slope angles in the area ranged from 38° to 58° along the roadway downslope of the site. Site 6458 (Fig. 1d) located in Pike County was reported on 2/10/2018. The slope angle was approximately 63° with sporadic rock outcropping with little to no overburden in some areas.

Landslide susceptibility map: machine learning analysis
The general LSM was created using the logistic regression machine learning model to analyze landslide site data listed in the KGS Landslide Inventory. The logistic regression machine learning model was selected for this study based on the high prediction accuracy reported in previous studies (Kalantar et al. 2018;Kadavi et al. 2019;Nhu et al. 2020;Crawford et al. 2021). The machine learning analysis was implemented using the Scikitlearn library in Python (Pedregosa et al. 2011). Scikit-learn is an open-source machine learning library that supports supervised and unsupervised learning. It also provides various tools for model fitting, data preprocessing, model selection and evaluation, and many other utilities. A Geographic Information System (GIS) model was used to prepare and extract features such as geomorphic variables, soil physical variables, and land-cover variables, for the machine learning analysis. These features are independent variables and are used to predict target variables, known as Landslide Conditioning Factors (LCFs). The LCFs were compiled into database as training points for the machine learning analysis to predict landslide susceptibility. The KGS Landslide Inventory listed 234 known landslide occurrences in Campbell and Kenton counties in northern Kentucky region, including some data entries for the reported failure dates. All 234 landslide sites were used to extract LCF features for the machine learning analysis. The machine learning model was trained using 75% of the sites and 25% was used in testing the data for validation.

Variables and statistics compilation used in the machine learning analysis
Landslide occurrence is influenced by bedrock geology, geological structure, hillslope morphology, soil type and thickness, and hydrogeological conditions (Xu et al. 2018). We used 15 variables in the machine learning analysis from three distinct categories including soil property features, geomorphic features, and land-cover information. Table 1 shows all used features with their definitions. The soil property features play important roles in determining soil strength and mechanical behavior. Figure 2 shows the soil property features and their range of values extracted from the NRCS WSS (https:// webso ilsur vey. sc. egov. usda. gov/ App/ HomeP age. htm) that include percent sand, percent silt, percent clay, saturated hydraulic conductivity, available water capacity, one third bar water content, plasticity index, and liquid limit, all of which were used as features in the logistic regression analysis. In addition, land-cover data from of 2016 National Land Cover Database (NLCD) product suite (www. mrlc. gov) (Fig. 2i) were used in the analysis. Geomorphic features were used as LCF features in landslide susceptibility modeling in several studies (e.g., Reichenbach et al. 2018;Crawford et al. 2021). The geomorphic features used in this analysis were slope, aspect, curvature, elevation, roughness, and plan curvature. All these geomorphic features were derived from 1.5-m LiDAR DEM in ArcGIS and shown in Fig. 3. The 1.5-m LiDAR DEM was taken from the Kentucky Elevation Data and Aerial Photography Program (https:// kyfro mabove. ky. gov) and is publicly accessible through an open-source portal.
The extent of each landslide was assumed to be contained within a circular buffer area with a 45.7-m radius from the centroid point of the landslide. This buffer size was used as a mask layer to extract feature statistics for each landslide. Therefore, the spatial resolution of the LSM was 91 m. Figure 4 shows how the mask buffer areas used to extract each feature were aligned with the landslide extent in a close view. Most of the buffer masks fell within the extent of the landslide boundary. Timilsina et al. (2014) observed that a buffer polygon that represents most of the landslide extent was superior to a single point in accounting for variability in landslide characteristics. Crawford et al. (2021) evaluated the performance of several buffer sizes for predicting landslides in eastern Kentucky using the MATLAB Classification Learner application (https:// www. mathw orks. com/ help/ stats/ class ifica tionl earner-app. html). These researchers found a buffer size of 45.7 m performed best in evaluating landslide susceptibility. Given that the landslides in this current study were in the same region as those investigated by Clay percentage in the soil layer Saturated hydraulic conductivity The ease with which pores in a saturated soil transmit water (micrometers per second) Available water capacity Quantity of water that the soil can store for use by plants (centimeters of water per centimeter of soil layer) One third bar water content Amount of soil water retained at a tension of 1/3 bar (volumetric percentage of the whole soil) Plasticity index Numerical difference between the liquid limit and plastic limit of the soil. It is the range of water content in which a soil exhibits the characteristics of a plastic solid (percent) Liquid limit Water content, on a percent by weight basis, of the soil (passing #40 sieve) at which the soil changes from a plastic to a liquid state (percent)

Land-cover information Definition
Land-cover data Land-cover characteristic (e.g., grassland, shrub, deciduous forest) 1 3 Crawford et al. (2021), the 45.7-m radius buffer size was assumed to be most appropriate for used in this study to extract the features. Basic statistical parameters were used to capture the variation of the distribution of the data in each buffer area for each LCF features. The statistics used in the analysis were the mean, standard deviation, variance, skewness, and coefficient of variance. Descriptive statistics (e.g., mean, range, standard deviation) of elevation and slope in the buffer area are better predictors of the presence (or absence) of landslides than the same indices computed for the single DEM cells (Carrara et al. 1991;Alvioli et al. 2016;Reichenbach et al. 2018). The statistics used for each feature in the buffer areas for the landslide and non-landslide are calculated in ArcGIS and shown in Table 2.
There were 15 independent variables with five statistics for each variable except for NLCD. For the NLCD data, only the median value was used because the spatial resolution the NLCD dataset was 30 m, which was smaller than the buffer size. Therefore, a total of 71 features were used in the machine learning model as training points. There were some dependencies between the features such as elevation variance feature and elevation standard deviation feature. However, the 71 features were narrowed down to an optimal number  of features that produced a higher accuracy of the result as determined by the statistical skill assessment methods used in the analysis.
For the data preparation, all 234 landslide occurrences in northern Kentucky were assigned the class value of 1 for landslide, and all 71 features for each landslide were extracted using the buffer mask area. An equal number of points that represented nonlandslide areas were randomly selected and assigned the class of 0 for non-landslide. All 71 features for the 234 non-landslide areas were extracted the same way, using the same size buffer mask area as the landslide areas. The non-landslide areas were further manually inspected to confirm that there were no landslides in the non-landslide buffer areas. An equal number of landslides (landslide = 1) and non-landslides (non-landslide = 0) was required for an equal class-distribution ratio, which helped to eliminate class bias (Gupta et al. 2019).  Table 2 Statistics for each of the features in a buffer area for logistic regression model x i , grid cell value; Md, median value; N, the total number of samples in the training or testing dataset in each feature for the buffer area

Logistic regression model for the creation of the LSM
Logistic regression is a statistics-based linear model used to model the probability of the existence of a certain class. Logistic regression is also known in the literature as logit regression, maximum-entropy classification (MaxEnt), or the log-linear classifier. In this model, the probabilities describing the possible outcomes of a single trial are modeled using a logistic function. The dependent variables or indicator variables, which were predicted using the logistic regression model, were value of 1 (landslide) or value of 0 (nonlandslide). The probability of absence of landslide (0) or presence of landslide (1) was expressed as Eq. 1: where P is the cumulative estimated output probability of a landslide occurrence (confined between 0 and 1), z was the weighed linear combination of independent variables (ranging from − ∞ to + ∞) that could be expressed by Eq. 2 as the sum of constants (Kadavi et al. 2019): where p∕(1 − p) is the corresponding odds or the likelihood ratio, 0 is the constant intercept; V i ( i = 1, 2, 3, ...., n ) are the independent variables (e.g., curvature standard deviation, slope variance) and i ( i = 1, 2, 3, ...., n ) are the coefficient estimates of the model. The coefficients in Eq. 2 express the effects of the predictor independent variables on the relative risk of being a landslide or not a landslide (0 or 1). The relative risk increases or decreases with each value of the independent variable V i (i.e., the rate of change) in logodds as V changes (Crawford et al. 2021). Consequently, higher z-values suggest a P-value closer to 1, which indicates the presence of a landslide. Conversely, lower z-values suggest a P-value closer to 0, which indicates the absence of a landslide. The i coefficients of the model were estimated using a cost function as shown in the works of Lee and Liu (2003). The threshold value of P = 0.5 in Eq. 1 was used as a decision boundary such that a P-value of 0.5 or greater was classified as Class 1 (landslide), and a P-value less than 0.5 was classified as Class 0 (non-landslide). We assumed that the variables were not normally distributed or did not have linear relationships (Nandi et al. 2016). Therefore, the logistic regression analysis was well-suited for this study because the primary unknowns were the relationships among the variables, which tends to imply nonlinearity (Crawford et al. 2021).
The initial result including all 71 features (independent variables) was developed and analyzed in the logistic regression model using the Scikit-learn library in Python (Pedregosa et al. 2011). The optimal number of features were determined using a univariate feature selection method. The univariate feature selection method is a preprocessing step that works by selecting the features that produced the highest accuracy of the model based on univariate statistical tests. In the univariate selection process, each feature of the indicator variable (1, 0) was compared to determine if there were any statistically significant relationship between the variables. This method is also known as analysis of variance (ANOVA), which analyzes the relationship between one feature and the indicator at a time without the consideration of the other features. The univariate selection method yields test scores for each feature, and the features with top scores were selected. Figure 5 shows the optimal number of features with their corresponding accuracy percentage (given as crossvalidation scores in the figure) using the univariate feature selection method. Based on the univariate statistical test on all the training data, 18 features produced the highest accuracy and were subsequently used in the logistic regression model. The accuracy for the model increased to an optimal of 77% at 18 features and decreased to an approximate accuracy value of 73% after 30 features. Table 3 shows the 18 features that produced the highest model accuracy using univariate selection method with their corresponding ranking selection scores. The selection score was calculated using chi-squared statistics between each feature and indicator value in the univariate selection tool. The selection can be used to select the number of features with the highest values for the test chi-squared statistic from feature sets relative to the indicator value.

LSM model performance and validation
The logistic regression model for the LSM was validated using the Receiver Operating Characteristic (ROC) metric to evaluate the classifier output quality (Landslide = 1, Nonlandslide = 0). The ROC is a practical method to visualize the performance of the model. In the validation process, we used 25% of the landslide inventory points (59 landslide points and 59 non-landslide points) as test data. The area under the curve (AUC) of the ROC provided the accuracy of the classifier for the test points. The ROC was created by plotting "sensitivity" as a function of "specificity". Specificity and Sensitivity were defined as, where TP is true positive i.e., the number of landslides (Class 1) predicted accurately; TN is true negative i.e., the number of non-landslides (Class 0) predicted correctly; FP is false positive-the number of non-landslides predicted as landslides, whereas FN is false negative-the number of landslides predicted as non-landslides. The specificity indicates the  proportion of the negative class that was correctly classified, and the sensitivity indicates the proportion of the positive class that was correctly classified. Figure 6 shows the AUC for the model, which was calculated in the test dataset as 0.841. The aforementioned machine learning model was used to produce landslide susceptibility mapping on a large-scale area such as Campbell County. ArcGIS was used to extract the selected 18 features (see Table 3) from the spatial data, and to enumerate them in the tables that were required for the machine learning model. The landslide susceptibility result from the machine learning model was inserted back into ArcGIS to produce the LSM. Figure 7 shows county-wide LSM for Campbell County using the result of the machine learning model. Crawford et al. (2021) developed a landslide susceptibility classification system based on standard deviation from the mean of landslide occurrence probability. We used a similar approach for this current study. The mean probability value of landslide occurrences for our study was 0.634 with a standard deviation of 0.186. One standard deviation from the mean was 0.448, and two standard deviations from the mean were 0.262. We assumed rounded values for break points. High susceptibility was values greater than the mean (approximately 0.6), high-to-moderate susceptibility was within one standard deviation from the mean (0.41 to 0.6), moderate susceptibility was within two standard deviations from the mean (0.31-0.4), and low susceptibility was values below two standard deviations from the mean (< 0.3). Table 4 shows the susceptibility levels applied to the 131 reported landslides in Campbell County. The table shows that approximately 87% of the reported landslides in the county (approximately 114 landslides) were identified as being in areas of moderate-to-high susceptibility and areas of high susceptibility. Campbell County has a total area of 413 km 2 . In terms of a percentage of the total area shown in Fig. 7, 53.6% of the area was classified as low probability, 20.1% of the area as moderate probability, 19.9% as high-to-moderate probability, and 6.3% of the area as high probability. Figure 8 shows the validation of the LSM using the actual landslide areal extents in a close view in which the extent of seven landslide occurrences in the area 3.5 km by 2 km are all within or crossing areas classified as high moderate to high.
Because the statistical-based LSM is created by static features, the LSM provides estimated probabilities of landslide susceptibility only in spatial terms. LSM describes the relative likelihood of future landslides occurrences based solely on statistical methods using slope variables such as those used in this study (see Table 1) from the past landslide events. This implies that future landslide events will be more probable to occur under the conditions in which the past landslide events occurred. The mapped landslide susceptibility provides useful tools for avoiding or reducing potential damages from the landslides. In contrast, LHM predicts not only the spatial component of the landslide susceptibility, but also temporal aspect of the landslide events, rendering the LHM time-aware dynamic map.  1 3 4 Landslide hazard mapping: factor of safety equation As was previously mentioned, the LSM provides an estimate of "where" a landslide is most likely to occur, whereas the landslide hazard mapping (LHM) provides an estimate of "when" a landslide will most likely occur within a region susceptible to landslides. We estimated the "when" landslides occurred using the limit equilibrium factor of safety equation presented by Lu and Godt (2008). Limit equilibrium factor of safety equations are based on driving forces in a slope being in equilibrium with the resisting forces. Lu and Godt (2008) suggested that the resisting forces in a hillslope system were functions of the hydrologic conditions (i.e., volumetric water content and matric suction) in the slope and those conditions dynamically varied with depth in the vadose zone. The implication of the Lu and Godt (2008) suggestion is the factor of safety of a hillslope will vary as a function of depth with dynamic variations of the hydrologic conditions. Further, the occurrence of landslides in colluvial soils roughly correlate to the seasonal variations in soil water content and matric suction (Zhuo et al. 2019). Therefore, we created a LHM model using the Lu and Godt (2008) limit equilibrium factor of safety equation that varies with varying water content. The efficacy of the LHM model was tested using the four landslide sites used to validate the LSM model.

Factor of safety equation: creation of dynamic landslide hazard map
The Lu and Godt (2008) equation presents the factor of safety for subaerial infinite slopes in the unsaturated saturated soil conditions. The Lu and Godt (2008) equation is given as the sum of three distinct contributions to slope strength; the frictional strength component, where FS is the factor of safety; ′ is the soil friction angle; c ′ is the effective soil cohesion; is the slope angle; is the soil total unit weight; H ss is the depth to bedrock; s is the suction stress, a characteristic function of the soil that describes the inter-particle stresses resulting from the wetting and drying of the soil. The suction stress is given as, where S e is the effective degree of saturation = − r s − r ; is the volumetric water content; r is the residual volumetric water content; s is the saturated volumetric water content; s is the soil matric suction, which is defined in this study using the van Genuchten (1980) soil water characteristic curve (SWCC) model given as, where u a is the pore air pressure; u w is the pore water pressure; is a fitting parameter reflecting the air entry value; n is a fitting parameter related to the inflection point of the SWCC; m is fitting parameter related to the curvature of the SWCC near the residual point, m = (n − 1)∕n.
The shear strength parameters (i.e., c ′ and ′ ) for the four study sites used in the LHM model were acquired from various geotechnical reports obtained near the site locations. The geotechnical reports were obtained from the Kentucky Transportation Cabinet Geotechnical Project Report Database (https:// kgs. uKent ucky. edu/ kgsweb/ kentu cky/ search. asp) and from estimates made using data from the NRCS WSS database. The geotechnical reports indicated that the soils in the northern Kentucky area consisted predominantly of colluvium (i.e., silty clay). When not directly given, the soil friction angles were estimated using the Wood (1990) equation and the plasticity index extracted from the WSS database given as, where PI is the plasticity index. The effective soil cohesion, c ′ , for the colluvium was found to be negligible. Therefore, the second term in Eq. 5 dropped out. Figure 9 shows the estimated map of the soil friction angles, ′ derived from Eq. 8 and the plasticity index data in the vicinity of study sites 6294 and 6501. The soil friction angle, ′ , ranged between 27.6° and 33.6° with the lower end of this range dominating the study area in Campbell County, Kentucky.
The spatial resolution of the LHM was set as 15.2 m (approximately 50 ft). This spatial resolution was chosen to detect landslide occurrences that fell between local-scale (grid size of 5 m or less) landslide occurrences and regional-scale (grid size of 20 m to 30 m) landslide events (Kakavas and Nikolakopoulos 2021). Thus, the 15.2-m resolution represents an approximately average grid size between the local and regional scales. The slope angle map was created using the 1.5-m LiDAR DEM in ArcGIS. Therefore, the resolution of the slope map was 1.5 m. Consequently, the slope map was resampled to a 15.2-m cell grid using a bilinear resampling technique. The bilinear resampling technique performs a bilinear interpolation and determines the new value of a cell based on a weighted distance average of the four nearest input cell centers. The bilinear resampling technique is useful for continuous data and will cause some smoothing of the data (ArcGIS Pro manual, https:// pro. arcgis. com). Figure 10 shows the resampled slope angles used in Eq. 5, which ranged from 0 to 57° at the study sites. Depth to bedrock, H ss , was assumed to be 4.5 m in the study region as a typical depth to the bedrock. This value was taken from the geotechnical reports in the study area and used as the slip surface depth in Eq. 5. Bittelli et al. (2012) pointed out that the infinite slope model often assumes that the failure plane is coincident with the soil-bedrock interface. The unit weight of the soil, , was assumed to be 18.8 kN/m 3 . A summary of the soil shear strength and slope parameters used for the study sites is given in Table 5. All the parameters in Eq. 5 are constant except for s = suction stress, which is a function of hydrologic parameters such as volumetric water content and matric suction. These parameters enabled the LHM to be dynamic.

Estimation of matric suction
As was previously discussed, soil matric suction was defined using the van Genuchten (1980) soil water characteristic curve (SWCC) model. The van Genuchten (1980) model parameters s , r , , n for the study area were estimated using the Rosetta online tool (https:// www. handb ook60. org/ roset ta). These parameters were estimated according to the Schaap et al. (2001) pedotransfer function, which utilizes the percentages of sand, silt, clay, and bulk density, all of which were extracted in ArcGIS from the NRCS WSS soil database. These data are shown in Table 6. Figure 11 represents the sand, silt, and clay percentage used in estimating the hydrologic parameters for Eq. 5. However, if the estimated s from the Rosetta tool was less than the highest value of obtained from volumetric water content obtained through remote sensing, then s was replaced by the highest value of . It is observed from Fig. 11 that the study area mostly consists of silt and clayey silt soil.

Volumetric water content satellite data acquisition
Volumetric water content data from the National Aeronautics and Space Administration (NASA) Earth and Precipitation Satellite Missions were used for a time span of six weeks, three weeks prior to the reported landslide failure dates and three weeks post the landslide failure dates for each study area. Time-series data for satellite-measured soil moisture were acquired from the NASA Soil Moisture Active Passive (SMAP) Earth satellite mission. The SMAP data were accessed and acquired using the Application for Extracting and Exploring Analysis Ready Samples (AppEEARS) tool (https:// lpdaac. usgs. gov/ tools/ appee ars/). SMAP provides four data products, one of which is the Level 4 (SMAP L4_SM) product. This product is a model-derived value-added product obtained by merging SMAP observations with estimates from a land surface model in a data assimilation system. The land surface model component of the assimilation system is driven with observations-based meteorological forcing data, including precipitation, which is the most important driver for soil moisture (Reichle et al. 2011). The model-derived product produces three-hourly estimates of surface and root zone soil moistures (to a depth of 100 cm) at a 9-km gridded spatial resolution with a data availability latency of 7 to 14 days (Chan et al. 2016). The SMAP L4_SM product was used for this study. The multiple three-hourly soil moistures reported over a 24-h period were averaged to provide a daily value of soil moisture. Figure 12 shows the SMAP data in the northern Kentucky area temporally and spatially. Figure 12a shows 9-km spatial resolution grids in Campbell County, over which the soil moisture data are uniform. The volumetric water content ranges from 0.42 to 0.51 in Campbell County (Fig. 12a) from January 14, 2016 to February 25, 2016 during which the failure occurred. Figure 12b shows the variation of soil moisture change with time at Site 6294. The data were reported three weeks before and three weeks after the failure date.

Performance and validation of landslide hazard map
The soil and hydrologic parameters that remain constant ( , ′ , s , r , , n , m , H ss ) and the temporally varying daily volumetric water content at each study site were used in Eq. 5 to generate the daily factor of safety values for the given investigative period (i.e., three weeks prior to the landslide, and three weeks post the landslide). The factor of safety values were validated using the four study sites to assess whether the LHM correctly mapped the landslides spatially and temporally. The LHM covered 300 m by 300-m area with a spatial resolution of 15.2 m at each study site. The validation was made by overlaying the known failure area with the created LHM temporally and spatially at each site. Table 7 shows the factor of safety values a week before the failure date, failure date, and a week after the failure date at each study site. Overall, the landslide hazard map indicated the areas where landslide occurred as having factor of safety values ranging from 1.15 -1.73. We observed that these factor of safety values exceed a value of 1.0 (i.e., the value at incipient failure). This exceedance is most likely due to geomorphic and geotechnical behavior that is too complicated to be captured by a simple subaerial infinite slope model (Griffiths et al. 2011). However, the general trend shows that the failures occurred as the factor of safety values reached their lowest levels. Therefore, the factor of safety values were assumed to be relative to the dates failure were observed rather than being absolute.

Validation of model performance at sites in which the LHM model performed optimally
Site 6396 in eastern Kentucky and Site 6294 in northern Kentucky showed similar results that identified the landslides correctly temporally and spatially. The LHM derived from Eq. 5 for Site 6396 and Site 6294 indicated that only areas around these sites showed a low factor of safety on the date the actual failures occurred. The LHM was symbolized with different colors to illustrate the areas with different factor of safety values by the break points based on factor of safety values; values lower than the factor of safety value on the failure date (approximately 1.15) for each site (red), values between 1.151 and 1.75 (orange), values between 1.751 and 2.5 (yellow), values between 2.51 and 3.5 (light green), and values ≥ 3.51 (dark green). Figure 13 shows the factor of safety evolution for Site 6396 a week before the failure was observed (Fig. 13a), the failure date (Fig. 13b), and a week after the failure (Fig. 13c). As shown in Fig. 13b, only the area around the failure location was correctly identified with the red color grid indicating the low factor safety of 1.15 on the day of failure, and everywhere else, the factor of safety is greater than 1.15. The LHM for a week before the failure (Fig. 13a) shows that the factor of safety is more than 1.15 everywhere in the map,  Fig. 13c a week after the failure capturing the evolution of the landslide area. Figure 13d illustrates the factor of safety values at Site 6396 over a time span of six weeks covering three weeks prior and three weeks after the failure. The factor of safety reaches 1.15 on the failure date which was the lowest value during the time span (Fig. 13d) at the Site 6396. Figure 14 shows the performance of the LHM at Site 6294. Similar to the analysis performed at Site 6396, the factor of safety values were assumed to be relative to the factor of safety on the date of failure rather than being absolute. For Site 6294, the red areas in Fig. 14 correspond to the factor of safety value of 1.29 (rounded to 1.3 in Fig. 14) at the time of failure. Also, like Site 6396, Fig. 14 illustrates the evolution of the factor of safety values before failure, at failure, and after failure. Figures 13d and 14d illustrate that the factor of safety values dropped to 1.15 and 1.29 when the failure occurred, respectively, and both plots correctly showed the time of failure within the six-week observation period. Figure 14d shows that the factor of safety value at Site 6294 began to increase after the initial observed failure. Subsequently, the factor of safety decreased to a value of 0.5. The decrease in the factor of safety corresponds to the soil moisture increasing to the highest value 10 days after the failure date. Equation 5 dictates that the factor of safety will decrease with increasing water content. The factor of safety values below the relative factor  Figure 15 shows that the factor of safety was 1.73 at Site 6501 on the failure date, which was much higher than the factors of safety observed on the failure dates at the other validation sites. Also, a factor of safety of 1.73 is unreasonable to use when evaluating slope failure, given that very stable slopes have factor of safety values greater than or equal to 1.5 (Lu and Godt 2008). Therefore, the relative factor of safety value representing failure was set as 1.3 to compare with Site 6294 (i.e., the other Campbell County site). Several red grid cells (i.e., factor of safety less than 1.3) are observed near Site 6294 (see Fig. 15b). However, these areas failed in 2016. According to the notes reported in the KGS Landslide Inventory, the area around Site 6294 was consequently remediated and stabilized. The slope angles and the soil properties are most likely no longer similar to the data obtained from the NRCS WSS soil database or from the KYTC geotechnical reports. Therefore, the red grids are most likely not representative of the conditions at Site 6294 for the time shown in Fig. 15b.  Figure 15c shows red areas northeast of Site 6501 a week after failure. These areas may be reflective of actual failures, but these areas are in a wooded area up a hill, well away from a highway. There were no reported landslides for these areas in the KGS Landslide Inventory. Therefore, landslides in areas northeast of Site 6501 were not verified. Figure 15d shows the factor of safety at the Site 6501 was 1.73 on the failure date but continued to decrease to 1.26 five days after the reported failure date. While this observed continued decrease in the factor of safety may be indicative of an increase in the post-failure soil moisture, the overall LHM performance at failure and one week after failure is not consistent with reported behavior in the study area during the observed times. Figure 16 shows the LHM for Site 6458. The figure correctly shows the failure location with factor of safety of 1.27 on the failure date, which was assumed to be the relative failure factor of safety at the site. However, Fig. 16b shows more areas with similar or significantly less factor of safety values on the failure date. Figure 16c indicates that a week after failure a significant portion of the study area was estimated to have factor of safety values less than the relative failure factor of safety value of 1.27. The 1.5-m LiDAR DEM indicated that many of the areas around Site 6458 have very steep slope angles (greater than 45°), which are generally too steep to support colluvium of much depth (Glover 1999). Therefore, these steep areas were assumed to be rock outcroppings. We argue that the mixture of rock outcropping and colluvial slopes in the study area resulted in soil material properties that were not accurately reflected in the data obtained from the NRCS WSS soil database. In addition, the steep slopes imply that the actual depth to the bedrock was shallower than 4.5 m. Deviations from the soil strength data and soil depths from the values initially used in Eq. 5 most likely resulted in in situ factor of safety values being higher than are calculated from Eq. 5. The possible difference between the in situ factor of safety values and the factor of safety values estimated from Eq. 5 would explain the abundance of red areas observed in the LHM on the failure date and one week after the failure date (Fig. 16b, c, respectively). Steep slopes and shallow soil cover would also promote an increase in surface runoff and would result in higher estimates of soil moisture due to higher surface reflectance (i.e., the basis for satellite-derived soil moisture estimates). Thus, the argument of differing geomorphic characteristics than were originally assumed may also explain the observation of Fig. 16d that the factor of safety value continued to decrease three weeks after the failure date indicating the soil moisture continued to increase after the failure.

Validation of model performance at sites in which the LHM model performed slightly less than optimal
The implication of the results presented in Figs. 15 and 16 is that Eq. 5 may not be valid for areas comprised mostly of steep rocky slopes with little to no soil cover. Regardless of this implication, the LHM provided insights into the temporal and spatial evolution of factor of safety values. Thus, providing a means to assess changing landslide hazards in near real time.

Comparison of the landslide hazard map to the landslide susceptibility map
Comparison of the LHM and the LSM was performed to assess the spatial accuracy and variances of both models for the area surrounding Site 6501 on the date the lowest factor of safety occurred (4/4/2018). In addition, comparison between the LHM and the LSM was performed to further assess the predictive performance of the models for the area surrounding Site 6501 at a date other than the date of failure. The date was arbitrarily chosen as two weeks before the reported failure date (3/20/2018).

Visual comparison of the LSM and LHM for the area surrounding site 6501
For the quantitative comparison of the two models, the spatial resolution of the mapping models was required to be the same. Therefore, the spatial resolution of LSM was resampled to convert the 91.4-m spatial resolution to a 15.2-m spatial resolution to match the 15.2-m resolution of the LHM. The quantitative comparison between the probability map (i.e., LSM) and factor of safety map (i.e., LHM) was performed assuming a factor of safety values less than 1.75 as likely to be a landslide. The 1.75 value represents the highest factor of safety value calculated for the study sites on the failure dates (see Table 7). For the LSM, a probability value of 0.40 and higher was assumed to be a landslide because 87% of actual landslides in the inventory belonged to this category (see Table 4). Figure 17 shows the comparison of the mapping models. Figure 17a Table 8 compares the LSM and LHM models quantitatively using a confusion matrix. A confusion matrix is a table that describes the performance of a classification model on a set of test data for which the true values are known. In this study, we used the confusion matrix to compare the model classes (landslide versus non-landslide) on a cell-to-cell basis. In the confusion matrix, the classifications of landslide occurrence versus non-landslide occurrence for LSM are read from top to bottom. The classifications of landslide occurrence versus non-landslide occurrence for LHM are read from left to right. As was previously mentioned, the criteria for evaluating a cell as an occurrence of landslide (i.e., landslide = 1) was a factor of safety value less than or equal to 1.75 for the LHM model, and a probability value greater than or equal to 0.4 for the LSM model. Thus, on the day of the lowest factor of safety (4/4/2018), the confusion  Table 8 shows that on the date of the lowest factor of safety, 1633 cells on the maps were classified as landslides by both LHM and LSM, and 4726 cells were classified as non-landslides by both models out of 9508 total cells on each map. Conversely, the confusion matrix shows 977 cells in the LHM were classified as landslides, while the LSM showed these cells as classified as non-landslides. The confusion matrix shows 2172 cells were classified as landslides by the LSM, whereas those cells were classified as non-landslides by the LHM. In the area surrounding Site 6501, two weeks before the date of the lowest factor of the safety (03/20/2018), 57 cells on the maps were classified as landslides by both the LHM and the LSM. Also, 5683 cells were classified as nonlandslides by both models. However, two weeks before the date of the lowest factor of the safety, the LSM classified 3748 cells as landslides, whereas those cells were classified as non-landslides by the LHM.

Comparison of the LSM and LHM for the area surrounding site 6501 using a confusion matrix
The confusion matrix must be evaluated in the context of the landslides reported in the KGS landslide inventory. The inventory indicated that on the date of the lowest factor of safety, no additional landslides were reported in the area surrounding Site 6501. Therefore, cells in the mapping models, other than those immediately adjacent to the location of Site 6501, classified as landslides are considered false positive. The higher number of false positives imply a more spatially conservative and possibly less accurate mapping model. However, true comparisons of accuracy between the LSM and the LHM must also include the visual observations of landslide occurrence at the Site 6501 shown in Fig. 17. From these criteria, the LSM yields a less accurate predictive model than the LHM. This assessment of accuracy can be further supported from the comparisons of the area surrounding Site 6501 two weeks before the lowest factor of safety. The KGS landslide inventory reported no landslides in this area during this time. Therefore, it is expected that the number of cells classified as landslides will be low. For the LHM, this is in fact the case. However, because the LSM is predicated on static factors, it indicates the same number of false positives as before.

Discussion
Although comparative analyses were performed for the LSM and LHM models, true direct comparisons are not strictly valid. The two mapping models serve different functions. The LSM identifies various degrees of probabilities in which areas are susceptible to landslides based on a host of static factors (see Table 3). Thus, the LSM is not a temporally predictive model. In contrast, the LHM is a temporally predictive tool that allows the user to predict the likelihood of a landslide occurring given dynamic changes to the hydrological state. Although the LHM yields spatial predictions of landslides, the model does not provide a probabilistic likelihood of an area being susceptible to landslide. The two models are directly comparative only from the standpoint that landslide occurrence is loosely predicted by the LSM model in specific areas that are highly susceptible to landslides, compared to landslide occurrence predicted by the LHM at a given location, for given conditions. In actuality, the mapping models are complementary to each other. Specifically, the LSM identifies the area most susceptible to landslides, while the LHM predicts when landslides may occur within the identified susceptible area. Uncertainty in model parameter evaluation is an important cause of mismatch between simulated and observed distributions of landslide occurrence (Burton et al. 1998). Therefore, when performing predictive analysis with physics-based models, spatial variability and uncertainties in ground conditions must be considered (Bittelli et al. 2012). The LHM model specifically relied on model parameters that were generalized over specific area. A limitation of this study was the soil strength and hydrologic properties used in the Eq. 5 to create the LHM were estimated using data reported in the NRCS WSS soil database and various KYTC geotechnical reports. These data may not accurately reflect the soil spatial variance at the targeted areas. For example, the internal friction angles used in the analysis for Site 6458 and Site 6396 were taken from the geotechnical reports for the project located approximately 3.2 km away from the study site. Although these friction angles produced a LHM model that matched the location and the time of the failure at the study sites, it is quite reasonable to assume these friction angles are not valid at all locations within the study area. The limitation of the spatial resolution of the input data also pertain to the soil moisture. The soil moisture data from SMAP L4_SM has a spatial resolution of 9 km, and thus, the soil moisture was assumed to be uniform across the entire 9-km grid cell. Variations in soil moisture within the 9-km grid due to abrupt soil changes would not be captured. Another limitation for this study was that the depth of the bedrock was taken from the average depth reported in select geotechnical engineering reports for transportation projects in northern Kentucky. The actual depth to the bedrock at the study sites was not in situ measurement. Future research must investigate ways to obtain these data at higher resolutions. In addition, future research must answer the question, "what is the minimum resolution required to produced desirable results." Regardless, the framework presented in this study can be generalized to develop landslide susceptibility and predictions models with a high degree of confidence, using satellite and remote sensing data that are publicly available.

Conclusion
We created the LSM using a logistic regression machine learning model. Features for the logistic regression model were developed from variables that included six geomorphic variables extracted from 1.5-m LiDAR DEM, eight soil variables extracted from the NRCS WSS soil database, and a land-cover variable from the NLCD product suite. The result of the LSM model produced an AUC of 0.841. A LHM model was created using a subaerial infinite slope factor of safety equation for unsaturated soil conditions (Lu and Godt 2008). The forcing function of the factor of safety equation was satellite soil moisture data from SMAP L4_SM. The geomorphic characteristics of the slopes were obtained from 1.5-m LiDAR DEM, and the soil strength and hydrologic data were derived from the NRCS WSS soil database. We validated the results of the LHM using four study sites with known landslides and the failure dates. The results showed the LHM correctly indicated low factor of safety values at the study sites temporally and spatially. However, some additional areas within the study regions for Site 6548 and Site 6501 estimated the occurrences of landslides at times and locations not identified in the KGS Landslide Inventory. We argue that these areas were comprised of steep slope angles (greater than 45°) and were most likely characterized by rock outcroppings and thin soil cover. The uncertainty of the estimated friction angle, the depth to the bedrock, and the soil moisture data, all of which contributed to the false indication of the landslide areas in the LHM for those sites.
We show that using publicly available data, a multi-temporal landslide hazard map can be created that will produce a close-to-real-time landslide predictive map. The landslide hazard map derived from the limit equilibrium analysis tells us the evolution of the landslide development temporally and spatially, whereas the landslide susceptibility map derived statistically indicates the probability of landslides occurring at specific locations.

Data availability
The datasets and codes generated and used in the analyses of this current study can be found at the following reference: Bryson, L.S., Dashbold, M., Crawford, M.M., (2021). GeoDatabase and Modeling Code used for Landslide Hazard and Susceptibility Mapping in Eastern and Northern Kentucky: UKnowledge Civil Engineering Research Data, https:// doi. org/ 10. 13023/ kbrw-5x87.