Setting of parameters

In this section we discuss issues of parameterization of the model, the effects of the drainage area algorithm on the results, and the effects of grid size. In the next section (Testing) we examine how to evaluate the performance of the model in explaining observed patterns of shallow landsliding. This testing is important in assigning relative risk when the model is used for landuse prescriptions, a subject we discuss in a separate section.

As described above there are four parameters that could be evaluated in the use of equation (7): tanf, rs, T and q. The first three are soil properties and the last is effective steady state precipitation. Each of these parameters also varies spatially, but so far in most applications we have elected to assign a single value to an entire landscape (see Tang and Montgomery, 1995, for an exception). The three soil properties on average also vary between different landscapes. In the temperate rainforests of Coastal Oregon, for example, our data suggest that the wet soil bulk density is about 1600 t/m3 and the friction angle is in the mid 30's (Schroeder and Alto, 1983), whereas in parts of the California coast, the wet bulk density is about 2000 and the friction angle is in the 40's (Reneau et al., 1984). Root strength is an important contributor to overall strength, but as mentioned earlier it was eliminated in order to simplify parameterization of the model.

If cohesion is not considered, we have found it useful to set the friction angle equal to 45 degrees, and not let it vary between landscapes. This accomplishes two things: 1) it is a high value that reduces the overall area of potential instability relative to that which would be predicted with the more common lower values (say mid-30 degrees) and to some extent this makes up for the lack of cohesion in the problem (by making slopes as steep as 20 to 27 degrees stable) and 2) by holding it constant it not longer needs to be parameterized and different landscapes can be compared - a particularly useful exercise when using this model in a coarse screen to identify watersheds in need of watershed analysis.

The low range of bulk density that is likely to be encountered in the field has a small effect on the predicted pattern of slope instability when cohesion is neglected. In Oregon and Washington, we normally use about 1600 kg/m3.

Hence, by plotting the model outcome in terms of the log(q/T), and fixing friction angle at 45 degrees and the bulk density between 1600 and 2000 kg/m3, the model is parameter free.

In field application of SHALSTAB, the important things to measure are the drainage area to the landslide scar, the width of the scar (we normally use the crown width) and the local ground slope at the failure site. Surveys that fail to collect these data at individual sites can not be used to test SHALSTAB with field observations alone and would have limited value in slope stability analysis and comparison of failure conditions between sites. A great value of field measurementsis that the data can be used directly to test the model (without having to rely on often inaccurate digital elevation data). Of course, field mapping of landslide scar location alone is useful for testing SHALSTAB but without field observations it is often difficult to ascertain whether landslide sites that are not predicted as unstable by SHALSTAB are incorrectly predicted because of poor digital topography.

Use as a digital terrain model

Although equation (7) is a simple algebraic expression of the coupled subsurface flow and slope stability model that defines SHALSTAB, there is not a unique procedure for applying this model in a digital terrain framework. Surprisingly, the value of sinq, tanq, a, and b, may in fact vary considerably with choice of procedure. The best procedure we have found is that used in so-called contour based digital terrain models (i.e. Moore et al, 1988). Unfortunately, such methods are very difficult to use over large areas, and as mentioned earlier, we found it necessary to shift to grid-based methods. Here we describe how the current version of SHALSTAB performs these calculations, and mention some alternatives that could be used. The approach we selected was based on extensive testing that showed it was essentially free of grid artifact, i.e., results did not depend on the orientation of the topography relative to the grid system. This was considered essential for a model that was directed at site specific analysis of relative slope stability. Recent proposals for topographic analysis that are more mechanistic than our approach but that claim to be relatively free of grid artifacts, i.e. Tarboton (1997) warrant consideration.

slope determination

The use of SHALSTAB as a digital terrain model requires calculation of local slope for each grid cell. We have explored various algorithms to perform this calculation, including those available through ARC/INFO. In the end, we elected to use a program developed for our own use out of desire to minimize grid artifacts. Comparisons with the ARC/INFO SLOPE function in the GRID module give relatively minor differences for the cases we have tried and, in practice, this function has been regularly used by the Washington group.

The local slope is estimated as the geometric mean of the four directions that slope can be calculated across a given cell using the surrounding eight cells. The two directions normal to the cell have a spacing of 2 times the grid size, whereas the slopes diagonally across the cell have a distance that is 2.83 times the grid size. This procedure smoothes the landscape relative to taking the local maximum slope. Both our method and SLOPE in ARC/INFO, and most other such slope calculations, use the eight cell approach which means that the slope is being determined over a distance of at least 2 grid lengths (for 30 m grids this is 60 m, for example). This is important to appreciate when attempting to translate the digital terrain based calculation to the field.

We make no distinction between the slope gradient used in the sinq term and in the tanq term. A case can be made that these could, and perhaps should be different because sinq represents the total head gradient driving shallow subsurface flow, whereas tanq is the slope of failure plane (presumed parallel to the ground surface). We will discuss this further under the "alternatives" subsection below.

drainage area

Drainage area to each grid cell is equal to the component of every cell upslope that contributes some of its area to the cell of interest. One approach is to use the maximum fall direction out of each grid cell. In this case the fall path out of every cell is determined and it is then possible to count the number of grids that drain into each cell. The number of contributing cells times the area of each grid cell determines the drainage area to a cell. This way of calculating drainage area introduces huge artifacts, in which the drainage area to a point on the landscape may depend strongly on where it is relative to the orientation of the grid system.

To avoid this artifact, we employ a multiple-direction algorithm rather than maximum fall method of distributing area (as do others, i.e. Quinn, et al., 1991; Costa-Cabral and Burgess, 1994; and Tarboton, 1997). One can think of each cell as having its area and all the area draining into it to hand-off downslope. The total area the grid cell has to hand-off is divided in proportion to the local gradient out of the cell to each adjacent and diagonal cell. To illustrate, suppose there are three cells downslope of cell A. The B and C cells are adjacent to A and D is on a diagonal. Let the gradient to cell B be 0.27 and to C be 0.36 and to the D it is 0.32. Then the total area that cell A has (all upslope drainage area plus its own drainage area) is distributed to B as 0.27/(0.27+036+0.32) times the total area, to C as 0.36/(0.27+036+0.32) times the total drainage area and to D as 0.32/(0.27+036+0.32) times the total drainage area. Extensive testing has shown that this approach gives results that are independent of orientation of the topography relative to the grid. In SHALSTAB the proportion of slope in each direction is first calculated. Then starting at a low point in the topography, the contributing line is followed to the divide and then the area to the point at the bottom is calculated. This process is repeated for all cells. The specific catchment area, a/b, in our model is the total drainage area for each cell divided by the cell width.

Alternative approaches for topographic quantification.

We have made several choices about calculating a, b, sinq and tanq. First, we calculate drainage area in a manner that testing demonstrated has no grid artifact (i.e. the results do not change if the topography is rotated relative to the grid system and the model is rerun). As Tarboton (1997) correctly points out, the method we use (which is like the Quinn et al. method) is "dispersive" in that cells well outside the natural flow boundaries end up contributing to the drainage area of the cell of interest. This happens because any cell downslope of a given cell receives some proportion of the upslope drainage area, hence even on a planar slope oriented parallel to the grid system, each cell will hand-off to three lower cells, which in turn hand off to three more lower cells and so on. But, importantly, the amount area hand-off is progressively smaller, and each cell that hands-off some area receives some in return. In our tests, the cells end up receiving the correct amount of total drainage area and that area is smaller than the plan area of all the contributing cells upslope (many of which would lie outside the actual drainage area boundary because of this dispersive effect). Hence, while the algorithm we use would give incorrect results if one were to use it, say, for tracing contaminants down a hill, it does give reasonable estimates of the total drainage area to each cell.

Having then defined the drainage area of each cell, we then calculate the geometric mean of the local slope and assume that this slope is appropriate both for flow discharge and slope instability calculations - importantly we recognize that it is sinq that must be used as the flow routing slope. We assume that all the received water- based on the drainage area-then flows out according to this slope and across a unit cell width. In effect, then we find out how much water (based on drainage area) a cell must discharge and calculate its outflow according to a mean slope out one side of the cell.

This method is a compromise. It effectively avoids grid artifact because it permits outflow out only one cell boundary, but by doing so we do not define the actual flow path down the grid system, as the mean slope we calculate will rarely be pointed directly out only one cell boundary. The discharge we effectively calculate, and the resulting ratio of saturation, h/z, however, should still be a reasonable estimate for this cell. This compromise is possible because we use a steady state model in which the details of flow paths are not important. Furthermore, we avoid the large artifact which results from calculating drainage area according to just the maximum fall path.

It would perhaps be more consistent with how we route the drainage area to use an approach like that proposed by Quinn et al., 1991, to calculate b and the hydrologic slope, sinq. They set b equal to the sum of length of outflow planes it crosses (weighted by whether the path is diagonal or cardinal to the grid system), and calculate the slope as mean of the flow-width weighted local slope in each downslope direction. Simple geometric reasoning and our own testing shows that any algorithm that calculates b based on the sum of flow directions will be influenced by the orientation of the topography relative to the grid system, and that is why we have avoided using it. This can be easily understood from a simple example: a cell in a narrow valley bottom parallel to the grid orientation will slope only in the cardinal direction along a grid system, whereas a cell in a valley that cuts 30 degrees to the grid system may have three or more cells downslope of it. This is unavoidable.

The Quinn et al. method, however, may more effectively capture the effects of flow divergence on ridges than does just the drainage area algorithm we use. Preliminary tests show that the number of outflow directions (and, therefore, the total length of b) varies with topographic setting with the number typically being 5 or more on ridges, and 3 or less in valley axes. In both the 2 m grid Mettman Ridge case and the 30 m Greenleaf case shown in this report, the mean number of flow directions was 4. The length of b must vary with the number of flow directions, and hence a b calculated in a manner similar to that proposed by Quinn et al. would reduce a/b on ridges and planar slopes relative to that found in valley axes. This will have the effect of requiring a larger q/T for instability on ridges and planar slopes, and therefore reduce the overall predicted potential instability of the landscape. The local slope method for calculating sinq while consistent with how the drainage area is allocated, may be influenced by important local grid artifacts that averaging over nine cells avoids.

Pack and Tarboton (1997) use the Tarboton (1997) procedure in their solution of equation (7). The Tarboton model allows only two outflow directions and uses a local maximum slope calculation. Although Tarboton reports that his model is free of grid artifacts, it is not entirely obvious from the examples he shows. An interesting test would be to compare his results for different orientations of the same topography relative to a constant grid orientation.

As we stated at the introduction to this section, there is no unique and most correct way to covert equation (7) into a digital terrain model. Our approach is a compromise focused on avoiding grid artifacts. It is worth exploring whether sacrificing some to grid artifacts gains significantly in model performance.

Importantly, two things should always be stated in applying SHALSTAB: what methods were used to calculate the topographic attributes and what efforts were taken to map landslides to calibrate the hazard rating interpretation.

Topographic source and grid size

SHALSTAB relative hazard rating depends on local slope and topographic convergence (a/b). It is well-known that the coarser the grid size representing the landscape in a digital terrain model, the gentler the local slopes (e.g. Zhang and Montgomery, 1994). It is also not difficult to imagine, that similarly, the larger the grid, the less convergent the topography, hence local concentrated areas of high a/b like that shown in Figure 8 disappear. Hence, coarse grids produce relative smooth (lacking in fine-scale ridge and valley topography) landscapes.

For the United States, the most readily available digital elevation data are the 30 m quadrangle data. In many cases, these data were not derived from digitizing elevations every 30 m, but instead were generated from digitized profiles spaced at intervals farther than 30 m in which the 30 m grids were generated from extrapolation from widely spaced profiles. Topography generated from 30 m data fail to capture the fine scale ridge and valley topography that often dominates shallow landslide location. Figure 15 shows a comparison between SHALSTAB predictions for a portion of the Greenleaf quadrangle in the Oregon Coast Range in which 30 m USGS data (Figure 15a) and 6 m data derived from photogrammetric-based topographic maps created for the Bureau of Land Management (BLM). The percentage of the landscape in the BLM-based map in the < -3.1 value of log(q/T) is about twice (7.7 versus 3.3%) that in the coarse USGS-based map. Although the percentage of the landscape in the < -2.8 category is similar between the two data sets (12% for the 6 m BLM data and 15% for the USGS 30 m data), spatial pattern differs in important ways. In the finer scale BLM data, the low log(q/T) values (high slope instability potential) are more concentrated in steep valleys, rather than spread out across the landscape. These high slope instability areas follow distinct, but subtle hollows which are absent on the USGS topographic maps. Patches of "chronic" instability (where the local slope is ³ 45°) appear in the finer resolution topographic map (and our field experience in Oregon tells us these steep areas really do occur). Hence, the finer resolution topographic map increases the number of cells of highest instability, but more importantly the map enables delineating these high instability sites into distinct areas rather than as broad zones across the hillside. Furthermore, some landslide mapping suggests that with the higher resolution BLM topographic base, the majority of the landslides occur in the log(q/T) of <-2.8 which is only 15% of the landscape (as compared to <-2.5 and 37% of the land using the 30 m data). Hence the higher resolution topography, while delineating steeper more convergent areas, should also more accurately pin point that areas most prone to shallow slope instability.

Another comparison provides further evidence of the value of fine scale topography. Figure 16 shows the SHALSTAB predictions for 10 m grided data derived from 7.5' digitized contour lines and for 2 m grided data for the same area derived from airborne laser altimetry. As in the previous comparison, the major difference is that the pattern of relative slope stability is much more strongly defined by local ridge and valley topography in the finer resolution data. Comparison with landslide locations (Figure 2b) reveals that the minimum log (q/T) for all the scars in the 2 m case is -2.8 or smaller, whereas for the 10 m case here it is -2.5 , hence a much smaller area would be in the highest hazard zone in the higher resolution data (should the highest hazard condition be one of including all mapped landslides). While the 10 m data are much better than the 30 m, they still miss the finest scale topography in highly dissected landscapes.