Geospatial applications to landslide riskscape development: a modeling approach to quantify landslide riskscapes in the Colorado Front Range, USA

Riskscapes are interdisciplinary concepts that integrate multiple facets of physical, environmental, and social components in a spatial and temporal context. While the notion of risk is well documented for landslides, riskscapes are a novel approach in the natural hazard and spatial assessment studies. This term, ‘riskscape’, is described in terms of parameters required and quantication methodological approaches. Geographic Information Systems (GIS) or geospatial methods are an appropriate tool to dene the development of these riskscape quantication methods. A weighted sum overlay model for a riskscape is developed with three weighted approaches using GIS to measure the strength of spatial relationships across a regional landscape in Colorado, focused on landslide susceptibility modeling in the riskscape context. Binary riskscapes resulted in a limited understanding of the impact of features related to landslide riskscapes, but both ranked and human-factor weighted riskscape models provided more details to inform policy and plan for response to landslide events. Clustering measures using spatial-autocorrelation tools revealed that riskscape outputs are clustered and can further be used to identify areas of increased risk due to landslides in emerging population-growth areas. In conclusion, ranked and human-factor riskscape models are developed and can support decision-making and prioritization for response deployment based on landslide susceptibility criteria to focus resources on areas of interaction between landslide risk and social factors.


Introduction
Riskscapes for natural hazards are conceptual models that integrate the physical attributes and vulnerability dimensions of risk (Müller-Mahn 2013) within a spatial context. Landslides, commonly de ned as mass movements of earthen material down a slope (Varnes 1958;Varnes 1978) present hazards for many regions and increasingly interact with human-built environmental factors, becoming more threatening, damaging, and costly. Development of riskscape methods and applications to landslides contributes to a regional framework to address gaps in human/hazard interaction studies. A riskscape approach supports mitigation of hazard conditions by identifying areas of increased susceptibility to landslides and includes consideration of urban and human impact factors, allowing for better planning decisions and response prioritization.
Riskscapes as the intersection of natural hazard features in their environment can be assessed more fully by applying spatial autocorrelation measurements to geospatial methods. GIS methods for creating a landslide riskscape are based on GIS methods applied to natural-hazard-assessment and landslide-susceptibility-modeling literature. Many authors have demonstrated the e cacy of various geostatistical and geospatial methods for creating landslide risk zones depending on a variety of environmental factors such as slope, aspect, curvature, lithology, proximity to hydrologic and seismic features, soil characteristics, and land cover/vegetation (Dai and Lee 2002; Lee and Choi 2004;Huang et al. 2017;Zêzere et al. 2017;Carrara et al. 1991;Legorreta Paulin et al. 2014;Magliulo 2008;Huabin et al. 2005;Guzzetti 2006) Reichenbach et al. (2018) provides a review of these landslide models and statistical methodological approaches highlighting key contributors to the eld, including Carrara (1988Carrara ( -2008 and Guzzetti ( -2006. From a GIS historical perspective, Curry (1995, p. 1007) in Pickles (1995) re ected on the placement of GIS and its capabilities in the geographic ethical realm, stating " the necessity of seeing a 'realm'", like GIS, as constituted of sets of interlocking and overlapping patterns of actions". This statement, though regarding the placement of GIS in the literature, applies the notion of overlapping elements in space and how GIS has a role in determining these "patterns of actions". These patterns can lead to the discussion of riskscapes. The addition of riskscapes and focus on affected regions, in other words, the locations where the results matter, augments these GIS models with additional parameters and analytical processes. The goal of this study is to develop a model for GIS methods for riskscaping landslides for the regional study area of Boulder and Larimer Counties, Colorado.
Riskscapes are composite models, and represent the intersection of theoretical perspectives. Places of risk (Hewitt 1971) and risk as similar to landscapes (Cutter 1993) developed risk as a locational and pattern-driven feature. Müller-Mahn (2013) approached riskscapes as a social-theory structure, relating risk to the interaction of humans and environment, focusing on perspectives of riskscape participants, those who occupy the space where riskscapes occur. Applied riskscapes have explored multiple hazards assessments, using risk equations to calculate risk of natural hazard events and damages (Schmidt et al. 2011;Schmidt et al. 2007;Morello-Frosch, Pastor, and Sadd 2001;Jenerette et al. 2011). The RiskScape Program, a GIS application developed in New Zealand, is a module-based platform that creates cost loss estimates for modeled hazard events (RiskScape 2016). This program applies GIS processing to a riskscape platform, but focuses on the economic factors, not on the spatial relationships.
Traditional landslide risk assessment and susceptibility models lack regionally based spatial analyses that incorporate urban factors at a broad regional scale. Tobler's First Law of Geography states, "everything is related to everything else, but near things are more related than distant things" (Tobler 1970, p. 236). Applying this law as an analytical and modeling platform, GIS can be used to develop a riskscape framework for evaluating regional distributions of landslide probability and susceptibility with measurements for spatial autocorrelation and clustering, an applied extension of this First Law of Geography. Geospatial analysis tools operationalize the riskscape through integration of data from diverse sources and across sectors (i.e., urban areas, geology, infrastructure, population) and should factor in proximity as well.
Riskscape parameters extend natural hazard studies to shift the focus from the physical attributes of events to the affected or impacted communities. The region and occupants of the region (community) or stakeholders (externally located but interested parties responsible for supporting response or policy making) are included and the physical attributes are evaluated from a human dimension perspective. For example, a resident whose home is damaged by a landslide is part of the riskscape and needs to be included when evaluating regional support by agencies responsible for safety, transportation, funding response, rst responders, scienti c researchers, and other agencies. Riskscape factors include an additional classi cation step for modeling by using human geography and urban classi cations. These factors include urban classi cation (urban area or non-urban area), road presence, buildings, and potentially, response agencies, municipalities, population statistics, jurisdictions, or regions.

Objective of Study
This study developed a riskscape method for landslides using GIS models to evaluate the applicability of spatiallyweighted models to landslide riskscapes. Using frequency ratios and GIS weighted overlays, a landslide riskscape model was created for two study areas in Colorado. Spatial autocorrelation was calculated using Global Moran's i and Local Anselin Moran's i to test the strength of the spatial relationships. A quantitative riskscape model can be developed to apply analytical models that incorporate human factors and anthropogenic de nitions lacking in traditional landslide susceptibility models. Landslide models in Colorado need to include spatial and human elements and landslide susceptibility factors to better determine areas of potential impact. Areas identi ed as having potential higher impact can be used to develop policy and mitigation/preparedness strategies. The expansion of social construction frameworks to spatially-quanti able riskscapes creates a model and methodological approach that can cohesively describe and quantify landslide riskscapes and their relationship with spatially co-located features.

Study Area
Boulder and Larimer Counties, located in the Front Range region of Colorado, are mixed urban and non-urban counties with a history of landslide activity. Landslides are prevalent due to the mountainous terrain and periodic intense precipitation events, as well as interactions with wildland re burn-scar exposures. Boulder and Larimer Counties ( Fig. 1) have recent histories of landslide events and populations present in urban and rural spaces. These counties have been identi ed by the Colorado Geological Survey as the highest priority for mapping debris ows at the county scale based on their post-2013 Front Range ood-event evaluation (McCoy 2016). Each county region has urban and non-urban areas, elevation changes, differing geologic lithology, hydrologic networks, human and urban factors, and a recent history of landslide events.
Boulder County ranges in elevation from 1490 meters to over 4300 meters above sea level and has an area of approximately 1916 square kilometers (BoulderCounty.org 2019). Larimer county has a similar elevation pro le, from around 1450 meters to over 4100 meters above sea level and extends 6837 square kilometers (Larimer.org 2019).
Geologically, the counties range from shale deposits in the front range plains to the metamorphic units of the Rocky Mountains. Land-use/land-cover classi cations are primarily forested types in the mountains to developed and agricultural classes on the plains.
Boulder County population is estimated at 322,854 and Larimer County population is estimated at 343,853 (Affairs 2017). The population has grown signi cantly in recent years, with approximately 30,000 additional residents in Boulder County between 2010 and 2018 (Review 2020) and approximately 50,000 additional residents in Larimer County in the same timeframe. These increases in population in the riskscape area can lead to increased risk to humans and property. Population growth in the wildland-urban interface (WUI) needs to be considered as a contributing factor to riskscapes. WUI areas, the intersection between urban areas and wildland vegetation areas, are important spatial features indicating the relationship between human and natural features. Populations present in these areas are more at risk to wildland hazard events (Radeloff et al. 2018).

Materials And Methods 3 Materials and Methods
Quanti able approaches to landslide risk modeling incorporate elements related to the hazard, frequency probability, and impacted features (elements, consequences, or vulnerability). Multiple studies outline examples of risk equations for landslide events or natural hazards, including Bell and Glade (2004); Kappes, Keiler, Elverfeldt, and Glade (2012); Schmidt et al. (2011), andvan Westen (2013). These interpretations of risk equations [often given as risk = hazard x element x consequence, from Bell and Glade (2004)] are based on the hazard event, though they approach risk differently through the emphasis on different criteria. Vulnerability to a hazard is a component of several but not all of the approaches. Areal impact, or the region of the impact is only speci ed in one study, Kappes et al. (2012). Loss as a function of vulnerability, expressed as direct or indirect economic costs was also only found in one approach, van Westen (2013). The value of the spatial component is not measured consistently and does not use the geographic distribution as a measurable feature. To operationalize these spatial elements, a geospatial process was developed to integrate spatial aspects into the risk equation, and create a riskscape for landslides based on the in uential factors (e.g., physical features including slope, elevation, aspect, land cover, lithology, soil type, proximity to faults and hydrologic features and social consequence factors including transportation networks, urban area classi cation, population, and buildings).

Data
Data for the riskscape analysis are available from public sources (including US Geological Survey, Colorado Geological Survey, US Department of Agriculture, Colorado Department of Transportation, Colorado Department of Local Affairs, US Census Bureau, Boulder County, and Larimer County). Data are divided into classi cations of physical hazard factors, infrastructure elements, and human factors. Table 1 lists the sources for publicly-available data used in the analysis.

Data preparation
Datasets were preprocessed before analysis to create compatible layers of information. A common framework for projection (NAD83 UTM Zone 13 Colorado North), scale, and area were applied, and a frequency ratio calculation was performed to create a weighting scale (Table 2). Esri's ArcGIS version 10.6.1 and ArcGIS Pro version 2.3 were used to perform the spatial data analysis and to build spatial processing models. A data preparation work ow is shown in Fig. 2. Raster data cell size was 10 meters. Classi cations were based on Natural Breaks (Jenks). Probabilistic likelihood ratios were calculated from the percentages of landslides compared to the percentage of each class of the datasets present within the study area, using the formula from Pourghasemi et al. (2014)

Landslide Riskscape Class Development Methods
To develop the riskscapes for landslides in Boulder and Larimer Counties, 16 factors (lithologic units, distance to faults, soil type, land cover type, elevation in meters, slope in degrees, slope aspect and curvature, areas near owlines (rivers and streams), areas near waterbodies, urban area type classi cation, census population by block, distance to rail lines, distance to airports, distance to roads, and building footprints) were classi ed by intersections with landslide susceptibility areas using Zonal Histogram in ArcGIS Spatial Analyst. Factors were classi ed based on the presence or absence of landslide susceptibility cells within each dataset. The riskscapes were created using three models: the binary weighting (landslide present/landslide absent); ranked weighting based on frequency calculations ranging from 0 to 5 weight; and, a human factor weighing, where buildings, urban classes, and population were all weighted 5 with other factors weighted based on their frequency values. Three riskscapes were developed based on geospatial modeling techniques. Dorshow (2012) described a GIS predictive modeling process for archaeology that was applied to riskscapes and natural hazards. The weighting scheme introduced by Dorshow (2012) included equal weights for all of the models for one iteration of the analysis, and raw data inputs using map algebra, a geospatial raster-cell-value based calculation that creates additive layers based on input criteria scales (e.g. 1 for low risk, 5 for high risk), for the comparative model.
Nominal datasets (geology, soil, and landcover type) were normalized and data were converted to a weighted scale based on the probabilistic frequencies of the features coincident with landslide susceptibility cells. Data classes removed from the weighing models for datasets that had data replicated in other datasets (water bodies from soil and geology types)). Binary riskscape models use the presence/absence of a feature to generate the weighting percentages for an assessment of the features. This model includes only the location of a landslide susceptibility cell in the 10meter grid as a weighted factor. The ranked model used a scaled presence, with weights emphasizing the relative importance of a type class based on the probabilistic frequency ratio, or frequency of occurrence. The human factor model applied the highest ranking to all human (census population) and built environment features (urban class and buildings) regardless of the historic landslide presence or susceptibility ranking. Weights are listed in Table 2. 3.3 Spatial autocorrelation factors and clustering method Two methods of spatial autocorrelation measurements were used, Global Moran's i and Anselin Local Moran's i, to assess the strength of the spatial relationships and clustering within the riskscapes, The Moran index from the Global Moran's i tool ranges from − 1 which indicates a complete dispersion of features, to 0, a complete randomness of features, to + 1, a complete clustering of features. This index can be used to determine the amount of spatial distribution within each feature. Anselin Local Moran's i creates a mapping output of clustered areas based on a local neighborhood mathematical function, showing surrounding cells' relationships. Both clustering and cluster-outlier autocorrelation tools were performed on multipart (classed) rasters from the binary, ranked, and human factor input riskscape models. Global methods apply to all features within the dataset; local methods apply neighboring features to the analysis. The functions of the spatial autocorrelation tools are described in Table 3.
Global Moran's i spatial autocorrelation was calculated using three approaches, inverse distance, inverse distance squared, and zone of indifference. These methods differ in how the neighborhood weighting scheme in calculated but all are based on the fundamental rule of Tobler's law (Tobler, 1970), namely that proximity directly in uences clustering and values within features. Inverse distance without xed band zones (as used here) generates a local neighborhood based on creating a zone large enough that each feature has a minimum of one neighbor. Using the row standardization function can reduce bias introduced by aggregating polygonal data by summing the neighboring feature weights from each row (Esri, 2020). The inverse-distance-squared tool is based on the same concept but by squaring the distance function, it creates a steeper distance decay outside of the neighbor zone. Zone of indifference weights the cells within the distance-boundary setting equally and reduces those outside of the boundary based on the distance decay function (Esri, 2020).
The Local Moran's i spatial autocorrelation (Anselin Local Moran's i cluster outlier) was used to test for the spatial distribution of clusters within a dataset and generated z-scores, p-values, and output maps of clustered and outlier data. The tool created results for high and low value clusters or outliers. The map output indicates zones of clustering of high and low values where the cells are similar in rankings (hot and cold spots), and zones where outliers occur, or high or low values are surrounded by the other type of value (high and low outliers).

Binary riskscape results
The binary riskscape ranged from three to thirteen factors for the combined value of the weighted sum factors. Boulder County's binary weighted-sum riskscape based on ve classes had the medium class (7-9 factors) account for 79% of the areal extent but over 88% of the landslide susceptibility cells. For individual classes of factors, the class with eight total binary riskscape factors accounted for the highest probabilistic likelihood ratio at 1.58, with over 41% of the landslide susceptibility cells and 26% of the surface area. Larimer County had a concentration for the binary riskscape, with almost 90% of the area classi ed as medium, and over 94% of the landslide susceptibility cells occurred in the class shown in Fig. 3. The highest individual class for landslide probabilistic frequency ratio was class 8 (with eight factors), at 1.39, with over 37% of the landslide susceptibility cells and 27% of the surface area. Table 4 shows the binary riskscape occurrences.

Ranked riskscape results
The ranked weighted-sum riskscape generated a larger distribution due to the increased number of factor classes. The ranked weighted sum riskscape was classi ed into 9 classes using manual classi cation based on rounded Natural Breaks (Jenks) shown in Table 5.
Both counties showed that the higher factor classes had more prevelent landslide cells and higher probabilistic likelihood ratios. Larimer County had very high ratios in the highest classes, but these classes had very few landslide susceptibility cells. The results shown in Table 5 demonstrated that while areas of the highest risk represent a smaller areal extent of each study area, their probabilistic likelihood ratio was higher, indicating that the likelihood of the risk exceeds the base percentage of the areal coverage.
Using an equal-interval method for classi cation shows a different approach to viewing the data. This classi cation used integers grouped by ves to display data in an easier to read format. The ranked riskscape using equal intervals showed in Larimer County that the two classes accounting for most of the landslides are class 25-30 factors with 44% of the areal extent and 45% of the landslide susceptibility cells, and class 30-35 factors with 21% of the surface areal extent and over 34% of the landslide susceptibility cells. For Boulder County, classes with 25-30, 30-35, and 35-40 also had the highest probabilistic likelihood ratios, 1.32, 2.22, and 3.23 respectively.

Human riskscape results
The human factor weighted-sum riskscape values ranged from 7 to 51 total weighted factors for Boulder County and from 6.5 to 49 total weighted factors for Larimer County, classi ed into 9 classes (shown in Table 6) using manual classi cation based on rounded Natural Breaks (Jenks). Both counties show the highest probabilistic likelihood ratio in the higher classes, though fewer landslide cells and less areal extent occur in these classes. Human-factor landslide riskscapes are concentrated in classes 5, 7, 8, and 9 for Boulder County and classes 6, 7, 8, and 9 for Larimer County, showing the middle and upper ranges of factors contribute to the highest riskscape occurrence (Fig. 5).

Binary landslide riskscape spatial autocorrelation results
For Boulder County's binary weighted-sum riskscape using all classes Inverse Distance, row standardization and default threshold distances, the p-value was 0.00, and the z-score was 682.01. These show that the data are clustered, and that there is a statistically signi cant low (less than 10%) chance of this being random. The Moran's index value of 0.52 indicated a clustering of values. The Larimer County binary riskscape showed a p-value of 0.00 and a z-score of 835.08.
The Moran's index value of 0.47 indicated a clustering of values. These values indicate that the binary riskscape for Larimer County is also clustered. However, these results are not reliable as the binary riskscape only contained 14 total classes less than a standard sampling minimum. Therefore, it is recommended that the other ranked and human-factor riskscape results be used to determine signi cance. These results demonstrated that this type of classi cation shows clustering, but is not the most appropriate scale to measure autocorrelation.

Ranked landslide riskscape spatial autocorrelation results
For ranked riskscapes using all classes Inverse Distance, row standardization and default threshold distances, both counties indicated clustering for riskscape surface. Boulder County's ranked riskscape z-score was 2.06 and the p-value was 0.04, indicating the distributions were signi cant. The Moran index value of 0.080681 indicates there is measurable clustering in the features. For Larimer County, the ranked riskscape z-score was 11.65 and the p-value was 0.00, showing that the results were signi cant and had a low chance of being randomly clustered. The Moran index value of 0.62 shows that there is clustering in the features.

Human-factor landslide riskscape spatial autocorrelation results
Boulder County's human-factor riskscape z-score 15.22 was and the p-value was 0.00, indicating a low probability of the distributions being random. The Moran index value of 0.61 indicated there is clustering in the features. For Larimer County, the human-factor riskscape z-score was18.47 and the p-value was 0.00, showing that the results were signi cant and had a low chance of being randomly clustered. The Moran index value of 0.57 shows that there is clustering in the features. Figure 6 shows the ranked and human factor riskscape Anselin's local Moran's i clustering.

Discussion
These riskscape and spatial autocorrelation results demonstrate the applicability of a riskscape approach to addressing risk in a regional environment. The different approaches in riskscape modeling show that the classi ed binary approach, while the most basic, does not yield su cient accuracy over the region to be used extensively in mitigation strategies. While a suitability ranking forms an intuitive scale for reading the output, using only three classes did not contribute enough information and distinction amongst the data classes to be applicable to modeling response. Ranked and human-factor weighted-sum riskscapes show more detail in the data and are more appropriate for use in response planning. Urban areas are highly weighted in the human-factor studies, and the mapping output demonstrates the clustering of these urban areas. The ranked distribution in Fig. 8, map A shows the in uence of physical and environmental factors more strongly with lighter colors representing higher classes and numbers of factors. The human-factor map (Fig. 7, map B) shows the urban areas emphasized.
By including the urbanization factors, accounting for populations that may exist in the WUI and more rural areas, more detailed models can be developed to protect exurban populations from landslide damage, or present areas that must be included in planning for response. Extending future planning and mitigation efforts to incorporate riskscape models supports decision-making and prioritization for response to landslide hazard events. Figure 8 shows the human factor riskscapes coincident with WUI medium and high rankings.
Clustering analysis for the landslide riskscapes also emphasizes the need to develop mitigation strategies based on built-environment and human factors. When evaluating spatial features, the strength of distributions can be measured using spatial autocorrelation and clustering measurements. Further, these relationships can be measured spatially for their strength and correlations by using spatial autocorrelation and clustering tools. Spatial autocorrelation measures the amount of clustering within a dataset, and generates p-values, and z-scores to measure if these results are statistically signi cant. Spatial autocorrelation is not frequently expressed in landslide susceptibility studies. In traditional landslide susceptibility studies, these spatial relationships are often ignored (Erener andDüzgün 2010), andCatani, Tofani, andLagomarsino (2016) state "more quantitative information on the spatial variability of the MFDs [magnitude frequency distribution] of landslides in the geographical space is needed, a topic almost totally lacking in the relevant literature." Global and local spatial autocorrelation tools were used to assess an urban hazard environment (termed a hazardscape) for Mexicali, Mexico, using multiple types of hazards (Ley-García, Denegri de Dios, and Ortega Villa 2015). Catani, Tofani, and Lagomarsino (2016) applied spatial autocorrelation to landslide hazards to determine volume and dimensional relationships. Due to the spatial nature of riskscapes, and the value of evaluating spatial distribution of risks, spatial autocorrelation is a tool that can support the assessment of the data for independence or relative clustering. Clustering demonstrates the spatial relationships amongst the riskscape classes, quantifying the strength of the riskscape distribution.
Limiting factors are noted in the data. The clustered nature of urban settlements may in uence the amount of clustering seen in the spatial autocorrelation and clustering measurements. This effect is mitigated in part by the other features used to measure the riskscape such as the physical environmental features. However, human and builtenvironment features exhibit clustering which should be accounted for in the analysis. Classi cation of riskscape factors is also an area for improvement. Several classi cation methods exist, including natural breaks (Jenks), used here, as well as equal-interval and standard deviation. Depending on the content of the data, datasets applied, and number of factors, other classi cations methods may be appropriate.
Similar to the Global Moran's i method, the binary riskscape values were not effectively mapped given the lower feature count of 14 values. The ranked riskscape shows the emergence of patterns, with urban areas showing more high-low clustering in Larimer County in particular. The human factor riskscape provides more detailed clustering, with more urban and high elevation low-high and low-low clustering.

Conclusions
Riskscapes are transactions within space-time events, integrating quantitative and probabilistic approaches expanding across a regional landscape. Riskscapes incorporate elements of physical and human geography with natural hazard events and risks. This interface between disciplines strengthens planning capabilities and measures capacity to respond by identifying priority areas for preparedness and mitigation strategy implementations.
The riskscape incorporates the human factors and infrastructure data into the spatial model to evaluate the spatial relationships with the sensitive factors. Humans and the built environment are part of the landscape and in uence the sur cial changes. Spatial autocorrelation, the measure of the strength of these relationships between data, can be created. For a riskscape, the spatial relationships are measured between all of the factors in the region, including the hazard factors as well as the vulnerable elements and consequence data (human factors and built environment infrastructure).
This riskscape approach is based on using a spatial integration approach that expands on the use of risk equations for natural hazards in the context of landslides. Additional studies with different datasets or changes to scale (state-wide, municipality level) may yield results that place this regional approach in context.

Declarations
Funding: Not applicable.

Con icts
of interest/competing interests: Not applicable. This study is part of a student doctoral dissertation study.

Authors' contribution:
All authors contributed to this study conception and design. Material preparation, data collection and analysis were performed by Heather Hicks. The rst draft of the manuscript was written by Heather Hicks and all authors commented on previous versions of the manuscript. All authors read and approved the nal manuscript.