The Forest Practices Division of the Department of Natural Resources is mandated to regulate timber harvest operations in Washington State. It was recognized that the existing screen for shallow landslides was not sufficient and that several viable GIS-based models for slope stability were currently in use by timberland managers and scientists. However, no information existed that compared the effectiveness of the

various models. Prior to creating a new screen of modeled slope stability, a comparative analysis of several quantitative GIS models was undertaken to determine which model was most effective in predicting shallow landslide potential.

Introduction

Over the past several years, a variety of entities have developed GIS-based models for shallow-rapid slope stability. These models, however, have not been rigorously compared or adapted for statewide application to management and regulation of forest lands. This report briefly describes the methods, results, and conclusions of our comparative analysis. This test was conducted under contract to the Washington Timber/Fish/Wildlife (T/F/W) Program (i.e., a cooperative group of regulatory, tribal, environmental, and industrial sponsors who collectively makes recommendations to the Washington Forest Practices Board (WFPB) on matters related to forest management; T/F/W, 1992) and Washington Forest Protection Association (WFPA), as a precursor to developing the statewide slope-stability screen required by the WFPB.

During the course of this study, our focus expanded from evaluating models for use in regulatory watershed analyses and routine forest management, to include an assessment of their potential as statewide landslide-screening tools. This shift was driven primarily by the T/F/W negotiations and the resulting commitments of the legislature to promote the development of a statewide screen. We are developing a similar test for watersheds in each of the distinct geomorphic provinces in eastern Washington, as groundwork for creating a statewide screen of shallow landsliding. This test should help determine whether any of these GIS-based models can accommodate the geology and climatic regimes east of the Cascades Range.

People interested in land management in the Pacific Northwest historically have possessed limited means for evaluating landslide potential where activities are proposed. Existing information on site characteristics and failure potential typically has been confined to small geographic areas (e.g., 20 km² or less) in which landslide inventories, geomorphic research, or semi-empirical stability analyses have been conducted. More recently, private landowners and natural-resource agencies in Washington State have initiated a regulatory form of watershed analysis (WFPB, 1995) for specific landscape units (i.e., Watershed Administrative Units (WAUs), usually less than 200 km² or 78 mi² in size), in which landslide inventories are developed largely with the aid of aerial photographs and limited field reconnaissance. Landslide assessments in only about 60 of the 764 Watershed Administrative Units, however, have been finalized and approved by the state during the last seven years (Washington Department of Natural Resources (WDNR), 1999). Furthermore, incomplete and often imprecisely mapped state soil surveys and their slope-failure ratings still constitute the main source of information used by state regulatory foresters to evaluate management proposals in areas outside of those where reliable landslide assessments have been performed.

GIS-based slope stability models can be useful to managers for screening potential landslide areas and determining where land-use or habitat-restoration activities should be concentrated, to regulators for determining whether environmental checklists or impact statements are required, and to analysts for developing preliminary hazard-zonation maps. Isolated tests of GIS-based models in the Pacific Northwest have suggested that preliminary landslide-failure or hazard-zonations maps can provide more accurate slope-stability information than customarily can be interpreted from topographic, geologic, or soil maps alone (e.g., Shaw and Johnson, 1995; Montgomery et al., 1998).

Description of Test Models

Three GIS-driven models have been selected for this evaluation, based on their current availability, potential for adaptation to management decision-making, and/or use by T/F/W cooperators in field applications or previous tests of model performance. They are the current statewide soil-stability screen, maintained by the WDNR and herein labeled SOILS; the shallow landslide model of Montgomery and Dietrich (1994), nicknamed SHALSTAB by its authors; and the shallow landslide model of Shaw and Johnson (1995), herein referred to as SMORPH.

The three selected models have a number of elements in common. They use geographic information systems (GIS) to couple DEM data with assumptions regarding topographic attributes that influence slope destabilization and with algorithms for calculating slope stability. Whereas the SHALSTAB and SMORPH models assume that topographic relief (i.e., hillslope gradient) and form (i.e., slope curvature) are the principal driving factors in promoting shallow landslides, the SOILS screen assumes that only gradient is a critical variable. These assumptions derive from previous studies suggesting that shallow landslides occur most often above a threshold gradient and in topographic convergences where shallow subsurface flow concentrates, such as hollows and channelized depressions, with consequent effects on soil moisture and strength (e.g., Dietrich and Dunne, 1978; Swanson et al., 1981; Swanson and Fredriksen, 1982; Sidle et al., 1985; Montgomery and Dietrich, 1994). This simplifying assumption permits a number of key slope-stability factors to be treated implicitly, including substrate type, bedrock structure, rainfall duration and intensity, soil depth, soil conductivity and strength, plant transpiration, root strength, and subsurface drainage properties.

In addition, each model is limited similarly by the accuracy of the DEM data; that is, these models are only as good as the DEMs on which they are based. Much of western Washington is mapped with DEMs at a 10-meter resolution. For regions in which DEMs are available only on a 30-meter grid, however, all models suffer correspondingly in their precision and accuracy.

The three model differ primarily in the sophistication with which independent physical parameters affecting slope stability are addressed. The SOILS screen relies on hillslope gradient and soil type to rate slope-stability potential (Table 1)( WDNR, 1988). The SMORPH model explicitly treats gradient and slope curvature, while the SHALSTAB model treats these topographic attributes as well as several key soil physical and hydrological properties. From the standpoint of practical application, there are advantages and disadvantages to each approach. Simpler models in which key influencing factors are treated implicitly can be employed readily (i.e., with little to no data collection) and for larger geographic areas. The level of site-specific accuracy, however, might be reduced by assuming static or invariant hydrologic and geomorphic conditions, and by extrapolating local data on soil and hydrologic properties to the basin or regional scale. The advantage of explicitly treating parameters such as rainfall, subsurface hydrology, and soil properties is that the model might identify patterns of potentially unstable ground at a higher resolution. Consequently, such models are useful for predicting site conditions in the local area for which the input data apply. Conversely, employing local data might limit the ability of the model to predict accurately the spatial distribution of unstable slopes at a landscape scale. This approach also requires considerably more data collection in the field. Some factors, for example subsurface hydrologic and soil strength properties, might be very difficult to analyze and measure due to their spatial and temporal variations and their complex physical interactions.

A number of other models were considered but not chosen for this comparative test because of availability and software-development issues. (Wu and Sidle, 1995, Wu and Abdel-Latif (1995, 1997), Pack et al. (1998), and Earth Systems Institute, (pers. comm.). Other methods were too site-specific to be applied over large geographic areas, as required of a watershed analysis or statewide landslide screen (e.g., LISA and DLISA; Hammond et al., 1992). For a general review of analytical methods other than GIS-based modeling, see literature reviews in papers by Montgomery and Dietrich (1994) and Wu and Sidle (1995).

Methods

Study areas and landslide data

We chose eight areas in western Washington for this comparative test. The test basins range in size from 81 km² to 331 km² (Table 4). Existing Watershed Administrative Units (WAUs) were used as the test-basin boundaries, wherever possible. WAUs, defined for the purposes of regulatory watershed analysis, typically follow major drainage divides; the larger-order river systems, however, may be divided into several WAUs to limit the watershed analyses to a maximum acreage that reasonably could be assessed in the limited time period permitted by law (WFPB, 1995). Hence, some of our test basins comprise only the upper or mid- sections of a major river system (e.g., Chehalis Headwaters WAU, Middle Hoh WAU). Preference was given to those WAUs with recently completed watershed analyses, to utilize existing databases and to take advantage of the standardized format of data collecting used in this regulatory process.

We attempted to include at least one test watershed in each of the major geologic provinces in western Washington (Table 4; Thorsen, 1978). Parent materials range from glacial till/outwash and lightly metamorphosed sediments to volcanics and igneous intrusives. Test basins also vary in topographic relief (i.e., lowest to highest elevation points) from 818m., in the Chehalis Headwaters basin, to 1941m. in the Jordan-Boulder basin. Our intent was to examine model performance in areas with different combinations of relief and parent materials, as a means for exploring model versatility and the feasibility of using each model as a management tool in diverse topographic and geologic settings. The eight test basins contain a total of 2524 known landslides (Table 4), including shallow and deep-seated landslides (i.e., earthflows). We retained data on deep-seated landslides (e.g., earthflows) in the test database to evaluate the ability of each model to predict shallow landslide features that often are superimposed on more areally extensive earthflows.

Existing digital landslide inventories were acquired from the appropriate landowners in the test basins where watershed analyses had been performed (Table 5). Where inventories were not current or were spatially incomplete (i.e., original inventories covered only portions of the test area), we conducted aerial-photograph and field surveys to fill in data gaps. Aerial-photo series extended from the mid-1940's through 1996, in most instances. All inventories were updated chronologically to include, at a minimum, the most recent storm event known to have triggered widespread landsliding throughout Washington State (i.e., the high-intensity, long-duration storm of February, 1996; Gerstel, 1996). In addition, most inventories were checked in the field to verify database accuracy (e.g., landslide type, location, size). Road-related failures were retained in the test database, to evaluate the theory (e.g., Montgomery et al., 1998) that their locations are governed largely by hillslope gradient and topographic convergence. Standardized field data-forms were designed similar to the those used in the mass wasting assessment of the regulatory watershed analysis (WFPB, 1997, Appendix A). Newly identified landslides were mapped on to 1:24,000 scale topographic maps and then digitized into the GIS (Arc/Info^tm, version 7.2, for UNIX on a Solaris platform), coded, and edit-checked for positional and tabular accuracy.

In some cases, we updated the landslide inventories to include small landslides (i.e., less than 100m²) that might have been omitted due to time and mapping-resolution limitations that customarily constrain the regulatory watershed-analysis process. We increased the number of recorded landslides on these inventories by an average 12%, during our field and aerial-photo verifications of the databases. In the Upper East Fork Lewis River watershed, for example, our reanalysis of the GIS landslide-inventory cover maintained by the USFS resulted in a 70% increase in the number of recorded landslides. Hence, the watershed-analysis-derived landslide inventories really only provide a lower limit on the number of landslides present during the time period evaluated by the analyst (i.e., typically coinciding with the aerial-photo record). Consequently, landslide inventories were used here only as a common basis for comparing model abilities to predict known contemporary landslides, recognizing that other shallow landslides have been overlooked or perhaps no longer can be discerned in the field and photo records due to such obscuring factors as vegetation regrowth.

All inventory data were projected into Washington State Plane, south zone, North American Datum 1927. Having all data in the same projection allowed us to easily incorporate other existing data (e.g. hydrography, transportation), as well as provide a uniform projection from which to work.

We encountered a number of problems with existing landslide data while updating and verifying mass-wasting inventories from the completed, regulatory watershed analyses. These included incorrect basemaps on which landslides were recorded, as well as incorrectly mapped landslides. Discrepancies between USDI Geological Survey (USGS) topographic maps and basemaps created from GIS for use in watershed analysis typically included differences in topographic-contour delineations and stream-channel positions. Keying landslide locations to these features on USGS topographic maps, for example, apparently cause a positional offset when data are transferred to GIS DEM-based topography. A number of mapping errors also appeared to be related to inaccurate transfer of field data onto basemaps or incorrect digitizing from basemaps. In the Sol Duc watershed, for example, we determined from a reassessment of aerial photographs that several landslides were mapped in tributaries adjacent to the ones in which they actually exist. Hence, we remapped and redigitized landslides wherever we encountered such discrepancies during field or aerial-photo verification.

Another common mapping problem is related to landslide size. Mapping techniques used by analysts ranged from representing landslides as a point or symbol (e.g., circle) to delineating slides as polygons of finite area. The latter technique also included a range of mapping styles, from mapping the failure scarp separately to delineating the entire portion of slope involved in landsliding (e.g., some combination of the contributing area, initiation point, transport zone, debris-flow runout track, and depositional area), generally accompanied by little or no explanation of mapping style. In addition, landslide mapping is prone to some amount of inaccuracy, given that data are transferred between a number of different media (e.g., photos, maps, digital databases) with varying levels of resolution and precision, and often between different workers (e.g., field technicians, analysts, cartographers).

To address problems of mapped landslide location and size, we created a buffer around landslides mapped as points or symbols, or polygons smaller than 100m². The buffer, mapped as a polygon of radius 15m. (50 ft.) around the presumed center of the landslide feature, assured that landslides registered in a 100m² DEM grid cell when inventory data were compared with GIS model output. In many cases, landslide scarps and bodies were remapped, during aerial-photo and field verification of the existing databases, to exclude associated features (e.g., contributing areas and debris-flow runout tracks). The landslide polygons then were joined with the buffered landslide points to create a single coverage of mapped landslides. The polygon and buffer method also served to extend the mapped landslide area by an amount slightly larger than a DEM 10-m. grid cell, to account in part for imperfectly aligned digital landslide-inventory data and DEM topography.

Landslide hazard-zonation maps, created as a product of regulatory watershed analyses, were employed in this study to evaluate the ability of GIS models to predict areas considered by field analysts to have a potential for instability. Hazard-zonation maps produced via the regulatory watershed-analysis process (i.e., Mass-Wasting Map Unit maps; WFPB, 1997, Appendix A) typically delineate areas of presumed low, moderate, and high potential for landsliding and delivery of debris to downstream (or downslope) areas with sensitive public resources. Digital hazard-zonation maps were available in only four of the eight test basins (Table 5).

The principal dilemma faced with hazard-zonation maps is mapping resolution. Watershed analysts appear to use two styles of mapping: fine-scale and broad-brush. Fine-scale mappers delineate map units in detail, attempting to include in a high-hazard polygon only those slopes a high probability of shallow landsliding and to exclude any stable ground (e.g., the ridge lines between hollows in steep, dissected terrain). Given that such resolution can be intractable on 1:24,000 scale maps, another mapping option is to include the entire area in a generic mapping unit and explain in the report text how to differentiate high and low hazard zones on the ground. These broad-brush techniques promote Type II mapping errors, in which more area is included in a high-hazard unit than likely would fail.

GIS model calibration and database development

The SOILS screen required no adjustments to be employed in this study, and in fact cannot be adjusted to accommodate any new information, including altered soil classifications or gradient classes, without significant revamping of the GIS cover. The digital soils database for federal lands, maintained by the USDA Forest Service on the Internet, was merged with that maintained for state and private lands by the WDNR (1988). Nonetheless, six of eight test basins had incomplete digital soil covers (Table 5), due largely to gaps in soils-layer coverage on federal property. For statistical analysis of comparisons between the digital landslide inventories and soils slope-stability cover in these test basins, an existing landslide was given a "no data" value where the soils cover was lacking.

The SMORPH model was calibrated in each test basin with its respective landslide-inventory data to adjust the critical slope classes and their hazard-rating designations in the gradient-curvature matrix (Table 2). A slope map derived from the DEMs was intersected with the landslide inventory to determine the maximum gradient found in each landslide polygon. A curve of maximum gradient versus cumulative frequency percent was created with the lowest gradient at which a landslide occurred being used to determine the lower class limit of the moderate hazard rating. The lower class limit of the high hazard rating was established at a value for which 15% of the landslides occurred (Table 6), to guarantee a model-prediction rate of at least 85% of observed landslides.

For consistency with other published tests of the SHALSTAB model (e.g., Montgomery et al., 1998), we used the following soil-property values: soil depth (z) = 1.0m; soil bulk density (_s) = 2000 kg/m³; internal friction angle () = 33; effective cohesion (C') = 2 kN/m²; and transmissivity (T) = 65 m²/day. These values were selected by Montgomery et al. (1998) based on extensive field measurements in a small catchment in coastal Oregon (Montgomery et al., 1997), and the authors felt that they gave reasonable results for their test watersheds in western Washington, including the Chehalis Headwaters WAU that we also use as a test basin. We then compared predictions of unstable-slope potential for the range of angles and effective cohesions set internally in the model to yield a standard range of outputs (i.e., default parameters; = 33and 45, and c'= 0, 2, 5, 8, 15 kN/m²), to evaluate the effect of modifying these parameters. In section 4.2 of this paper, we discuss the sensitivity of model output to variations in input values.

Comparing SHALSTAB with the other GIS models required that we reduce all model outputs to a common denominator. SMORPH and the SOILS screen yield output in terms of management hazard ratings (e.g., low, moderate, high), in which the more subjective determination of what constitutes "hazard" and "risk" previously has been made in the policy arena. For example, the SMORPH slope matrix is calibrated with landslide-inventory and hazard-zonation databases created during regulatory watershed analyses for which definitions of hazard and risk have been set by T/F/W policy and WFPB regulations (WFPB, 1995, Chapter 222-22 WAC). Likewise, the SOILS screen hazard designations are derived from unstable-slope ratings in the state soil surveys. In the absence of another mechanism for converting all model outputs to the same units of measure, we therefore elected to assign hazard ratings to the SHALSTAB model output values of predicted critical rainfall, by using rainfall intensity and duration as the diagnostic criteria.

Given that SHALSTAB model output is expressed as rainfall in mm/day, we created "precipitation rules" for each test basin by clipping the two-year, 24-hour storm isohyte data (WDNR-GIS; Miller et al., 1973) and computing the minimum, maximum, and mean precipitation values for each basin. A high hazard rating was given to each DEM grid cell in which the predicted critical-rainfall value fell in the model-defined Q_c -stability class occupied by the mean precipitation value calculated for that basin (Table 7). A high rating was also given to any predicted Q_c less than the minimum two-year, 24-hour calculated precipitation. A moderate hazard rating was assigned to a DEM cell in which the critical rainfall value occupied the Q_c -stability class corresponding to the maximum calculated precipitation. A low hazard rating was assigned to all other Q_c stability classes. See Table 7 for the precipitation rules and slope-stability hazards created for each test basin.

The two-year, 24-hour recurrence interval was chosen as the precipitation regime for which data were readily available and which yielded the most conservative estimate of failure potential. The SHALSTAB model is configured such that the less frequent rainfall event yields a greater percentage of the basin area predicted to fail (Montgomery and Dietrich, 1994). Theoretically, then, a higher-intensity storm event characteristic of a longer recurrence interval, and/or a longer-duration rainfall, would result in greater spatial distribution of potential shallow landslides.

This method of assigning management criteria to SHALSTAB output was chosen in the absence of established techniques or direction provided by the authors (e.g., see discussion of management applications in Montgomery et al., 1998). A preferred approach might be to adjust the model in each test basin by using measured values of input parameters (e.g., soil transmissivity, bulk density, cohesion, internal friction angle), and calibrating predicted distributions of slope stability with observed landslide inventories and/or associated hazard-zonation maps in which management criteria have been assigned (i.e., similar to the approach used by SMORPH). Adjusting input parameters in the current version of the SHALSTAB model is problematic, given the relative paucity of soil-property data and the current lack of published algorithms for modelling stochastic elements or calibrating them from landslide inventories. Obtaining sufficient soil-parameter samples to adequately describe their spatial variability also could be intractable or prohibitively expensive for creating a landscape or regional GIS cover of predicted slope stability.

Calibrating model output with landslide-potential ratings from hazard-zonation maps is problematic. We found, for example, that hazard map units with different management designations (e.g., low and high) might contain DEM grid cells with the same range of Q_c - slope stability class values (e.g., 2 through 7; see Table 3), making it difficult to segregate the eight model-output classes into discrete management categories of low, moderate, and high. Calibrating model outputs solely on the basis of landslide inventories also can be misleading because, as discussed previously, they typically represent only contemporary rates of shallow landsliding, thus conceivably underestimating the density of potential landslide sites. Landslide density commonly has been a key factor in assigning management criteria to hazard-potential map polygons created from inventories (e.g., WFPB, 1997).

The precipitation rules imposed by this study make a number of assumptions, not the least of which is steady-state throughflow of subsurface water. The SHALSTAB model, however, is founded on the assumption of steady-state rainfall, constant transmissivity, and spatially uniform soil saturation (Montgomery and Dietrich, 1994). Hence, the steady-state precipitation rules are consistent with these assumptions.

Results and Discussion

We evaluated the performance of each model by using the GIS to intersect the updated, digital landslide inventories and hazard-zonation maps with model predictions of slope stability. For each model, output was expressed in terms of management criteria (i.e., low, moderate, high "hazard"), as described in the report section 3.0, so that model performances could be compared directly. We statistically analyzed the following, as a measure of the performance of each model:

(1) intersection of the digital landslide inventory with model predictions of hazard potential, expressed as the number of incorrectly identified landslides per total number of landslides in each test basin (i.e., Type I model errors);

(2) intersection of the hazard-zonation maps with model predictions of hazard potential, given as the percent probability that the model predicts a low landslide potential where it is likely that landslides have occurred or will occur (i.e., Type I model errors); and,

(3) intersection as in (2) but expressed as the percent probability that the model predicts the potential for landslides where they are not likely to occur (i.e., Type II model errors).

Inventories of known existing landslides and maps of hazard potential often are used in different management contexts. For that reason, we calculated Type I errors first by intersecting model outputs with the landslide inventories, to evaluate the ability of each model to predict the spatial distribution of existing landslides. We then computed Type I errors associated with comparing model outputs and hazard-zonation maps, to assess model abilities to predict the spatial distribution of existing and potential slope instability. Given that landslide inventories typically provide only a minimum estimate of contemporary landslide rates, the hazard-zonation maps theoretically yield a more complete view of the spatial distribution of past, present, and potential future landslide occurrences.

Table 8 lists, for each model, the number of incorrectly identified landslides per total number of landslides in each test basin (i.e., Type I errors). We assumed that an existing landslide was identified incorrectly if all DEM grid cells overlapping the landslide polygon or its 15m. (50 ft.) buffer (e.g., see report section 3.1) were coded by the model as having a low potential (hazard) for shallow landsliding. Conversely, an existing landslide was assumed to be identified correctly if any overlapping DEM grid cell was predicted to have a moderate or high potential (hazard) for landsliding. DEM cells with no data entry in the SOILS screen (i.e., missing soil-survey data) were coded as an incorrect identification, to account statistically for the incomplete nature of the data coverage. For this test, the SHALSTAB model was run using default parameters = 33 and C = 2 kN/m² and assuming that the two-year 24-hour storm recurrence interval is a reasonable criterion for assigning hazard-potential ratings to the model output (i.e., see report section 3.3).

A principal assumption of the model comparative tests is that predictions of landslide probability densities can be compared even though the GIS covers contain known mapping artifacts (e.g., elevation banding), as described in section 3.2. Given that model predictions of slope stability are evaluated using the same DEMs and landslide databases, the model outputs could be evaluated relative to one another. However, computed statistics (e.g., average number of landslides incorrectly identified by each model) should be viewed as estimates rather than absolute values, because the errors in model predictions associated with database noise (e.g., DEM elevation banding, field-mapping accuracy and resolution).

Table 8 indicates that the SOILS screen did not identify 32% of the total known landslides in all eight test basins, whereas the SMORPH and SHALSTAB models misidentifed 3% and 8%, respectively. The significantly higher percentage of landslides missed by the SOILS screen can be attributed to the lack or near lack of soil-survey data for two of the test basins (i.e., the North Fork Stilliguamish and Upper East Fork Lewis watersheds; see Table 5), given that missing data were coded as undetected landslides for the purposes of comparing model performances (see report section 3.1). Where the SOILS screen was complete (e.g., Morton and Chehalis Headwaters watersheds), however, it misidentified a significantly higher percentage of landslides than the other two models (e.g., for the Chehalis Headwaters watershed, 32% versus 2% each for the SMORPH and SHALSTAB models).

In the Olympic Peninsula test basins, the SOILS screen misidentified more landslides than SMORPH but fewer than SHALSTAB (e.g., in the Hoh watershed, 67 versus 53 and 84, respectively). The fact that these were the only basins for which 30-m. DEMs were used was ruled out as a likely cause. In other test basins for which model results were compared using both 10-m. and 30-m. DEMs, there was no change in the ordering of models based on their predictive accuracy, although the relative magnitudes of predicted landslide occurrence (i.e., number of correctly identified existing landslides) differed between 10-m. and 30-m. DEM test results for each model. Hence, the seemingly better performance of the SOILS screen might be explained by at least two compounding factors. One is that, for the portions of the test basins in which soils data exist, the SOILS screen classes 68% of the Sol Duc and 84% of the Hoh basin terrain as potentially unstable or very unstable, so that the majority of the landscape and its associated landslides fall within the high-hazard-potential category. Although this result lends the appearance that the SOILS screen more closely reflects the spatial distribution of known landslides than does SHALSTAB, it also tends to over-predict significantly the percent of watershed area predicted by field-derived, hazard-zonation maps to be potentially unstable (see further discussion of the SOILS screen in this paper section).

Another compounding factor is that the SOILS screen and SMORPH model consider hillslopes as being potentially unstable at gradients somewhat lower than the threshold gradient defined in the SHALSTAB model. In the latter model, slopes are considered unconditionally stable when tan tan [1 - (_w/_s)] which, for = 33 and _s = 2000 kg/m³, means any slopes less than 18 (32.5%). Field evidence suggests that non-road-related shallow landslides have occurred in this region on slopes closer to 25% (e.g., Shaw and Johnson, 1995; D. Parks, WDNR, pers. comm.), particularly in gently sloped, groundwater-seepage areas whose downslope margins coincide with the top of steep, inner-gorge slopes, which are quite common in this terrain. Hence, the SHALSTAB model has the potential for under-predicting the spatial distribution of unstable ground on hillslopes with gradients less than the threshold value set internally by the model.

The SMORPH model predicted an average of 22 times fewer Type I errors than the SOILS screen and five times fewer than the SHALSTAB model. The greatest discrepancy in SMORPH and SHALSTAB model predictions occurred in the Hazel watershed (1% versus 32% incorrectly identified; Table 8). Given that the Hazel watershed is dominated by deep-seated landslides in thick glacial deposits (Table 4), we expected the predictive capability of both models to diminish correspondingly, with respect to locating earthflow-influenced topography. It appeared, however, that SMORPH was better able to distinguish the local slope and curvature of numerous shallow-landslide headscarps superimposed on the larger earthflows. Hence, the polygons representing deep-seated failures effectively were identified by SMORPH predictions of high hazard potential on the basis of these smaller secondary features.

This variation in results might be explained by the manner in which the two models identify "hazard" potential in adjoining DEM grid cells. The SMORPH model analyzes variations in topographic relief between adjacent cells based on their relative steepness and curvature, then assigns a value according to the slope matrix (Table 2); hence, the model can discern topographic changes between a flatter upslope cell and a steeper downslope cell (i.e., a landslide headwall). On the other hand, the SHALSTAB model can smooth (i.e., not detect) subtle variations in topographic relief at the DEM-cell scale, by assigning a given flow tube a Q_c value depending on the flow across its upper boundary (i.e., variable "a" in Equation 2) from upslope contributing areas, which, in turn, is governed by the way in which flow is dispersed from that contributing area to any one of a number of downslope grid cells. Hence, if the upslope contributing area has a lower gradient and requires a relatively higher water flux to create "wet" soils, then a relatively steeper cell downslope (e.g., a landslide headwall) might not be predicted to fail until the same "wetness" is achieved. Hence, the grid cell downslope of the contributing area is given a lower slope-stability rating, whereas SMORPH assigns a higher value based solely on topographic factors.

Although Table 8 indicates that SMORPH yielded 43% fewer Type I errors in predicting known landslide occurrences than SHALSTAB (Table 8), we wanted to evaluate whether these differences in model performance, based on a comparison in eight watersheds, were significant statistically. We used a non-parametric test for non-normally distributed, small, independent samples to evaluate the hypothesis that there is no difference in the average performance of the SMORPH (SM) and SHALSTAB (SH) models, in terms of their ability to predict the spatial distribution of known landslides. The null hypothesis is that the means (µ) of the population of Type I errors for each model are equal when only eight independent samples (i.e., test basins) exist; H₀: µ_SM = µ_SH. Equality of means was tested with the Wilcoxon rank-sum statistic for two populations (Walpole, 1974; MathSoft 1998), in which the null hypothesis was true if:

Pr [W w = (a - n(n+1)/2)] > ,

where Pr is the probability distribution, W is the test statistic, a is the smaller of the summed ranks for each population, n is the number of observations corresponding to a, and = 0.01, 0.05 is the level of significance. Table 9 indicates that the test statistic is significant at a confidence level of 95%, permitting rejection of the null hypothesis, which suggests that the models differ somewhat in their ability to predict known landslide distributions; that is, µ_SM < µ_SH. However, the test statistic proved insignificant at the 99% confidence level (Table 9), allowing acceptance of the null hypothesis and implying that the difference in model predictive capability is relatively small. A similar statistical comparison of SHALSTAB and the SOILS screen indicated that the test statistic was significant at the 99% confidence level, implying that the screen and model are considerably different in their ability to predict existing landslide distributions.

Tables 10 and 11, respectively, give the estimated Type I and Type II model errors for the SMORPH and SHALSTAB model based on a comparison of model output with hazard-zonation maps. Error distributions were not computed for the SOILS screen, given that soils-survey data were complete in only two of the test basins, neither of which had usable hazard-zonation maps. Type I errors were calculated, for each model in each test basin, by intersecting the low-hazard DEM cells predicted by the model with the moderate- and/or high- hazard map units produced via watershed analysis (i.e., incorporating all map units intersecting with known landslides in the GIS inventory layer). This database intersection was expressed numerically as a percentage of model-predicted, low-hazard areas (in km²) overlapping field-mapped hazard areas. Type II errors similarly were analyzed by intersecting the high-hazard cells predicted by the model with the low-hazard map units and computing respective areas. These estimates were made for the four basins in which we had access to complete, digitized, hazard-zonation maps. To facilitate comparison (see Table 10 and 11), the percent error for each model (A/M) in each basin was normalized by the basin area in a given hazard class (A) divided by the total A for all four basins (T), that is: E = (A/M)(A/T).

Analysis of Type I error estimates with respect to hazard-zonation maps indicates that the SMORPH and SHALSTAB models similarly under-predict the percent area of hazard map units determined to be of moderate and/or high failure potential, by an average 6% and 5%, respectively. Using the Wilcoxon rank-sum statistic for two populations, as described previously, the computed test statistic proved insignificant at the 95% confidence level (Table 9), implying that the models perform similarly in predicting areas of relatively low hazard potential inside mapped landslide-hazard areas.

Whether the observed discrepancies between model predictions and hazard-zonation map units represent true " Type 1 errors" in the statistical sense is debatable, given that three of the four hazard-zonation maps (i.e., Jordan-Boulder, Hazel, and Sol Duc River) were drawn using broad map polygons that incorporated both unstable slopes and intervening stable ground. In the Jordan-Boulder basin, for example, hazard-zonation units intentionally were drawn to include potential landslide sites (e.g., hollows, groundwater seeps, inner gorges) and intervening divergent topography (e.g., ridge lines) because it was not possible to delineate them on 1:24,000 scale maps (Coho, 1997). Hence, the GIS-based models might discriminate, more accurately than the hazard-zonation maps, the topographic features potentially influencing shallow landslide initiation in finely dissected terrain.

As a test of the influence of mapping resolution on hazard zonation maps, we intentionally created the hazard-zonation map units in the East Fork Lewis test basin with as fine a resolution as possible on 1:24,000 scale maps. This allowed us to compare model predictions with two different scales of hazard-map resolution. Type I "errors" generated by SMORPH and SHALSTAB decreased substantially, from 14% and 9% for the Jordan-Boulder basin, respectively, to 1% and 2% for the East Fork Lewis basin (Table 10; values normalized as described previously). One implication of this result is that GIS-based model predictions of slope-stability potential could be used advantageously by analysts in drawing hazard-zonation maps with higher resolution than demonstrated.

Table 11 shows the distribution of Type II errors generated by the SMORPH and SHALSTAB models, based on comparisons with hazard-zonation maps. As in Table 10, error values are given as normalized relative percent areas. Calculated error estimates for each of the test basins suggest that SMORPH over-predicts the percent area of hazard-zonation map units designated as high landslide potential, by an average amount slightly less than predicted by SHALSTAB (i.e., 3% versus 7%, respectively). In all four test basins, SHALSTAB tended to over-predict, by a factor of two greater than SMORPH, the spatial distribution of high-hazard areas observed on hazard-zonation maps. With respect to the East Fork Lewis basin, which we believe was mapped fairly carefully for the purposes of this study, some amount of model over-prediction (i.e., 16% for SMORPH and 43% for SHALSTAB) might be true Type II errors. That is, the models likely do over-predict observed spatial patterns of slope-stability potential, as can be discerned from observed spatial patterns of existing and potential landslides. Particularly in the case of SHALSTAB, however, some portion of this over-prediction might be an artifact of the manner in which hazard-potential criteria were derived (i.e., the Q_c - slope stability classes assigned by precipitation rules to be included in the high-hazard management designation), as discussed previously.

To evaluate the potential for model use in a management context, we developed a ranking scheme to quantify model performance and a number of other comparative criteria. We employed a statistical method for ranking models in terms of their ability to correctly and incorrectly identify known, existing shallow landslides. A numeric value was assigned to each of the possible database-intersect outcomes:

For example, an existing landslide was considered to be identified by a particular model if any superimposed DEM grid cell was coded "high hazard" (p = 0) or "moderate hazard" (p = 1). The assigned values for all correctly and incorrectly identified landslides in each of the eight test basins were added to yield a cumulative score for each model, which then was normalized by the total number of landslides in each basin. Where landslides occurred in areas for which the soils survey data were missing, the SOILS screen grid cells were given a score of p = 2. These normalized scores then were added to a score sheet including results of other tested criteria, as will be described in report section 5.0.

Table 12 shows the results of this ranked test. SHALSTAB gained approximately twice as many points as SMORPH, reflected in the normalized cumulative scores (i.e., 1.9 versus 0.8, respectively). The SOILS screen received a significantly higher score (i.e., 6.7) than the other two models, due in part to the partial or total absence of soils-survey data in most test basins. SHALSTAB received a greater cumulative score than SMORPH, largely due to more frequent intersections of identified landslide polygons with model-predicted low and moderate hazards. Some of the discrepancy theoretically could be attributed to our assignment of management criteria via the precipitation rules.

At the outset of this study, we posed the following questions with regard to model performance: (1) How do model predictions of shallow landsliding compare with existing landslide inventories and hazard-zonation maps?; and, (2) How do model predictions compare with respect to each other? In summary, test statistics imply that the SMORPH and SHALSTAB models predict fairly well the spatial distribution of known existing landslides in the eight test basins (i.e., error frequency of 3% and 8%, respectively). These models, in general, also compare favorably with maps of shallow-landslide potential produced via watershed analyses (i.e., 6% and 5% Type I errors, respectively; and 3% and 7% Type II errors, respectively). The SOILS screen performed least well, missing 32% of the known existing landslides (i.e., Type I errors) and providing an incomplete cover of a substantial percentage of western Washington terrain (e.g., full data coverage existed in only two of the eight test basins). Test statistics also indicated that the mean differences in predictive model capability between the SOILS screen and either model were statistically significant, whereas the mean differences between SMORPH and SHALSTAB were marginally significant statistically. Hence, we conclude that the SOILS screen is comparatively less accurate and certainly less complete than the two tested models. While the average differences in predictive capability of SMORPH and SHALSTAB were not great, the former model tended to produce slightly fewer Type I and II errors. Contingent on the appropriateness of the precipitation-rule algorithm used to calibrate the SHALSTAB model, we conclude that SMORPH is slightly more accurate than SHALSTAB in predicting existing and potential landslides as represented in our updated landslide-inventory and hazard-zonation-map databases.

Acknowledgments

This study was funded by the Washington Department of Natural Resources, the Washington Forest Protection Association, and the Cooperative Research and Effectiveness Monitoring Committee (CMER) of the Washington Timber/Fish/Wildlife program. We thank Weyerhaeuser Company, the Tulalip Tribe, Murray-Pacific Corporation, and WDNR for providing data and local knowledge of test watersheds. Landslide and soils data for the Lewis watershed were obtained, in part, from the Web site maintained by the U.S.D.A. Forest Service, Gifford Pinchot National Forest. We thank David Montgomery, Tien Wu, and David Johnson for valuable input and assistance in compiling and running their models. We also thank Matthew Brunengo, Kate Sullivan, Daniel Miller, Tien Wu, George Pess, and Paul Kennard for valuable input in the early phases of this project. In addition, we thank Harvey Greenberg for assistance in making the SHALSTAB model operable on the WDNR GIS system; Daniel Miller for sharing his model and many critical discussions; Andria Villines and Elizabeth Freeman for assistance with statistical analyses; and Venice Goetz for her field assistance in verifying landslide data.

Table 1. Criteria for determining slope stability from the SOILS data.