Landslide susceptibility in sensitive clay in eastern Canada: some practical considerations and results in development of an improved model

P.E. Quinn

BGC Engineering Inc.

Email: pquinn@bgcengineering.ca

ABSTRACT

This paper describes the development of an improved and expanded landslide susceptibility model for eastern Ontario and southern Quebec, where the clay plains of the former Champlain Sea are affected by large landslides in sensitive clay. The work uses the weights of evidence method, and considers twelve different geospatial themes, of which seven are included in the final model: soil type; overburden thickness; flow accumulation; elevation; stream sinuosity; local relief; and, slope angle. The paper describes some practical considerations in the current modelling work, with emphasis on a number of key aspects, including: identification of target hazard types; selection of analytical and predictive study areas; selection of geospatial thematic data for analysis; image processing and analysis of both discrete and continuous valued raster data; and, combination of thematic weights for development of a meaningful susceptibility model. The new model has significantly better predictive power than the older model, and covers the entire area. It does not fully address earthquake-related landslide hazards, and would benefit from expansion of the digital landslide inventory.

KEY WORDS

Landslide, susceptibility, weights of evidence, sensitive clay, Champlain Sea

1 INTRODUCTION

This paper is an extension of previous work by Quinn et al. (2010), which presented a landslide susceptibility map for the western half of the former Champlain Sea in eastern Canada. The present work upgrades the model with support from additional data and extends the model to cover the entire Canadian footprint of the former Champlain Sea. A number of technical challenges in the bivariate statistical mapping method are discussed, and the upgraded model is presented with a discussion of model validation, and comment on its applicability and limitations. The new model has improved predictive power and is found to be broadly applicable across the complete study area. Subsequent work will attempt to extend this new susceptibility model to a regional scale qualitative risk map for large landslides in sensitive clay.

Large landslides in sensitive clay are widespread in the lowlands of eastern Ontario and southern Quebec, occurring in fine-grained marine sediments within the past limits of the post-glacial Champlain Sea. This paper describes the development of a landslide susceptibility map within the former Champlain Sea, which is shown approximately in Figure 1. That Figure also shows the outline of drainage basins containing watercourses that drain into the study area, as well as the rectangular outline of a digital database of existing large landslides in NTS 31H.

Figure 1. Study area: former Champlain Sea and outline of associated watersheds, with landslide inventory area (NTS 31H) shown.

This paper begins with a review of the analytical methodology, and summary of the previous modelling work. The bulk of the paper describes the current modelling work, with emphasis on a number of key aspects, including: identification of target hazard types; selection of analytical and predictive study areas; selection of geospatial thematic data for analysis; image processing and analysis of both discrete and continuous valued raster data; and, combination of thematic weights for development of a meaningful susceptibility model. The discussion of the current modelling is followed by a brief discussion of limitations and possible improvements to the model.

2 METHODOLOGY – WEIGHTS OF EVIDENCE APPROACH

2.1 General Review of the Method

This paper uses a bivariate statistical method, the weights of evidence method, first described by Bonham-Carter et al. (1989) for use in modelling mineral resources, and since adapted by others for use in landslide susceptibility modelling (e.g. Dahal et al. 2008). The specific application of the method for this hazard and study area was described by Quinn et al. (2010). A brief simplified review of the method follows.

The weights of evidence approach develops a susceptibility model based on the spatial relationships between landslides and terrain features from a variety of different thematic maps, each considered to have some relationship to landslide incidence (i.e. associated with a conditioning factor or triggering factor). For each thematic map (e.g. soil type), weights are calculated for each individual factor (e.g. clay, or sand and gravel for soil type), and a weight map is developed for that theme. All thematic weight maps are then added together to develop an overall susceptibility map.

The weight is an indication of the likelihood of encountering a landslide when the specific factor is present. The weights, Wi, are calculated as the natural logarithm of a ratio of spatial probabilities, as follows:

[1] W_i = ln[P{F_i|L}/P{F_i|L_bar}]

where F_i represents the presence of a specific factor, and L and L_bar represent presence and absence of a landslide, respectively. Note that the weight can have either a positive or negative value, with positive values indicating a positive correlation between factor presence and landslide occurrence.

The spatial probabilities in Equation 1 are obtained by totalling the number of raster pixels where landslides are present or absent and a specific factor (e.g. clay, sand or bedrock) is also present. Equation [1] can be computed in GIS using raster algebra to calculate the following, using “clay” as an example:

[2] W_i = ln[(A1/A2)/(A3/A4)]

where

A1 = area of landslides in clay (calculated as the number of pixels corresponding to both landslides and clay);

A2 = total area of all landslides (number of pixels containing landslides);

A3 = area of clay with no landslides; and

A4 = total area of no landslides.

3 PREVIOUS WORK

The initial development of a landslide susceptibility model for part of the study area was presented by Quinn et al. (2010). The previous modelling work examined eight different geospatial themes for potential inclusion in the model, including: soil type; overburden thickness; land use; elevation above sea level; flow accumulation in the drainage network; bedrock; slope angle; and, slope aspect. The final model incorporated the first five of these themes, and the latter three were excluded for various reasons. Figure 2 illustrates the conceptual combination of five different thematic weight maps to produce the overall susceptibility model. The previous susceptibility model, which divided the map area into three broad susceptibility categories (i.e. low, moderate and high), is shown in Figure 3.

Figure 2. Combination of thematic weight maps to produce original susceptibility model.

Figure 3. Original landslide susceptibility model for western half of study area.

A number of potential modifications were suggested by Quinn et al. (2010) to improve the previous susceptibility model at the regional scale, including:

– Expansion of the landslide inventory to include a more substantial area north of the Saint Lawrence River (e.g. NTS 31I);

– Consideration of additional spatial data, including: variation of landslide frequency from north to south; regional relief; river bank height; and, stream gradient and tortuosity; and

– Expansion of the model to include the eastern half of the former Champlain Sea area.

The current paper addresses improvements to the model related to the latter two points: consideration of additional spatial data, and expansion of the model to cover the whole study area.

4 MODELLING – CURRENT WORK

4.1 Introduction

The weights of evidence method, being a statistical approach, can theoretically be implemented without introducing bias into the analysis. However, there are a large number of important practical decisions that must be made throughout the process, each having some effect on the outcome of the analysis, and thus possibly introducing bias. The following sections describe the major steps in the analysis, including the following: identification and selection of applicable hazard types; determination of extent of the study area; choice of thematic data for consideration in the model; image processing and analysis; selection of themes for final use in the model; model validation and predictive power testing; and, model representation. The author’s experiences with, and approaches to, various key decision points are discussed where necessary in the following sections.

4.2 Selection of hazard types

A study area may be affected by several different geological hazards. Different hazard types are likely to be associated with different terrain features. For example, rockfall, rockslides and rock topples may each be expected to be strongly, but variably, controlled by bedrock type, geological structure and slope angle, whereas earth slides may have very little association with rock characteristics and a different relationship with slope angle. If multiple landslide hazards are present in the study area, they should be modelled separately, or at least be examined separately to test whether the calculated weights apply to all hazards.

The present study focuses on a single hazard type: large landslides in sensitive clay. These are not the only important landslides in the study area, and not all large landslides in the study area have the same cause(s) or mechanism(s). Small, local rotational earth slides are far more common than large landslides, occurring in the same general settings (i.e. along the numerous creeks and rivers dissecting the gentle pains). Large landslides in sensitive clay occur in at least two, possibly three, distinct types, categorized by Quinn et al. (2011) as flows, spreads and flake slides. While spreads are shown by those authors to be far more common than flows in the study area, and flakes are rare or absent in Canada, their spatial distributions appear to be indistinguishable, therefore they are lumped together into the single category of “large landslides in sensitive clay.”

The analysis discussed in this paper assumes natural triggering by river erosion, and this has influenced the selection of thematic data for use in the model. Perhaps 20 % of large landslides in eastern Canada are triggered by human activity, and the rest are generally considered to be triggered by river erosion (Lebuis et al. 1983). However, a number of landslides in the study area, including all of the largest known landslides (see for example Aylsworth and Lawrence, 2003; Desjardins, 1980; Eden, 1967; Hodgson, 1927; and, Leggett and LaSalle, 1978), are known to have been triggered by major earthquakes. Earthquake triggered large landslides in sensitive clay might be properly treated as a separate hazard type, distinct from erosion triggered landslides, for separate consideration in any subsequent risk analysis. However, it is not possible to distinguish triggering mechanism from air photo analysis, therefore the digital inventory of landslides within NTS 31H includes all interpreted large landslides in sensitive, with no distinction of trigger, whether anthropogenic, erosion or seismic. Small rotational landslides, being too small to resolve in the medium scale air photo study, and not included, nor are other landslide types, such as rock slides, which may be locally present within the study area where local geological conditions are anomalous in comparison with the typical gentle clay plains.

4.3 Study area

4.3.1 General

The study area can be divided into two distinct sections: an analytical area, within which the spatial relationships between landslides and other themes are calculated; and, a predictive area, which comprises the rest of the study area outside the analytical area. Both should be selected deliberately, with the objective of obtaining some degree of similarity in geology, geomorphology and physiography. Both sections should contain roughly the same variability in the geospatial thematic data; if a certain attribute widely present in the predictive area is missing in the analytical area – say a major soil type – then the model will not include what may be an important factor in landslide occurrence.

4.3.2 Analytical area and landslide inventory

Development of the susceptibility model relies on statistical comparison of specific geospatial features with the locations of existing landslides. The digital landslide inventory used in this work was reported by Quinn et al. (2007). That paper describes the mapping methodology and limitations associated with the data. A subset of the western half of the former Champlain Sea was selected for detailed air photo review: NTS map unit 31H, which is contained between 45 and 46 degrees north latitude, and 72 and 74 degrees west longitude. This area is shown as a rectangle in Figure 1. The landslide database includes a total of 1264 probable large landslide features. Figure 4 shows all of the interpreted landslides.

Figure 4. Landslide inventory in NTS 31H map area.

Computations of weights from equations [1] and [2] are spatially weighted, and thus selection of the boundary and associated area of the analytical study area can influence the model results. The landslide inventory covers a rectangular area (i.e. NTS 1H) that is largely comprised of gentle clay plains, but also includes rugged hills of the Canadian Shield to the northwest and Appalachians to the southeast. A subset of the inventory area was selected for analysis to simplify the work: the analytical area for the current work is a large proportion of the NTS 31H landslide inventory area, as shown in Figure 5.

Figure 5. Analytical area for landslide susceptibility modelling.

The extent of the analytical study area was selected to isolate the physiographic region affected by large landslides in sensitive clay. The areas outside the limits of marine deposition to the northwest and southeast contain greater relief, steeper slopes, higher elevations, different soil types and different drainage patterns than within the clay plains of the former Champlain Sea. For example, within the clay plains, elevation ranges from sea level to about 200 m above sea level, and outside the clay plains, it ranges from 200 m to roughly 1500 m, within NTS 31H. Large landslides in sensitive clay do not occur at those higher elevations, since marine soils are not present. It is important to include conditions where landslides are rare or absent in the analytical study area so that the model will predict absence where warranted. However, these areas should not be over-represented, which might dilute the statistical value of model results. The analytical study area includes some land with steep slopes and higher elevations at the fringes.

It may be noted that the size of the analytical area influenced the final selection of model parameters; earthquake activity and bedrock type, both expected to influence landslide distribution, did not vary sufficiently within the analytical area to yield a meaningful contribution to the susceptibility model. Expansion of the landslide inventory to yield a larger analytical area would arguably improve the model’s predictive power by allowing inclusion of these other important factors, but would come at the cost of significant manual effort in air photo interpretation.

4.3.3 Predictive area

Susceptibility modelling can be extended outside the analytical study area to include other areas with similar terrain characteristics. The success of the model in predicting landslide incidence outside the analytical study area depends on the degree of similarity between the two sections of the study area. In the present work, the gentle clay plains of the former Champlain Sea are assumed to be relatively homogeneous in terms of physiography and geology. This assumption is believed to be true at an appropriate scale of analysis, but is expected to be invalid at larger scales, as local variability comes into greater focus. For example, the surficial soils within the limits of the former Champlain Sea are all, generally speaking, marine soils, having been deposited in a marine environment, and a significant proportion of these soils are fine-grained, quiet water deposits (i.e. silt and clay). However, there were numerous rivers and streams draining into the Champlain Sea, including the proto-Ottawa river that drained much of eastern North America and may have been larger than any modern river in the world. There are therefore numerous deposits of coarser marine soils encroaching into the silts and clays, being locally present as abandoned deltas or other landforms.

The predictive area for the present work can be taken as either the area of the former Champlain Sea or as the combined drainage basins, which were used in the generation of inferred stream network and associated features, including flow accumulation, gradient and sinuosity of streams in the drainage network. Either choice is possible because the model predicts negligible landslide susceptibility outside the Champlain Sea boundary, due to the rugged terrain and absence of marine sediments, consistent with expectations. The outline of the predictive study area is shown in Figure 6.

Figure 6. Predictive area for landslide susceptibility modelling.

The predictive strength of the susceptibility model within the predictive area outside the analytical area would be improved by expansion of the landslide inventory to a larger area, or by shrinking the predictive area, potentially to coincide with the analytical area. Differences in predictive power between the analytical and predictive areas are discussed later in the paper.

4.4 Geospatial themes considered for the present susceptibility model

The model is generated from the addition of weight maps obtained from different geospatial thematic data. The selection of themes for potential inclusion in the model must be done deliberately, with the intention of selecting thematic data that will yield a useful correlation with landslide incidence or absence. In order to be selected for use in the analysis, thematic data should meet the following criteria:

– there should be some expectation of a correlation between some aspects of the theme and either presence or absence of landslides;

– the thematic data should have spatial coverage of all of the analytical study area, and most, but preferably all, of the predictive study area;

– the data should be available at a consistent scale, and mapped according to a consistent methodology; and

– the data should be generated by a single mapper/analyst, to ensure consistency in meaning of map units.

Additionally, it is preferable for the different themes to be statistically independent, in as much as possible, since the weights of evidence method assumes independence, which allows the linear combination of calculated individual weights. This independence criterion is nearly impossible to meet in practice, and where the data have strong correlations, it may be better to use multivariate statistical methods. However, in the author’s experience, weak correlations between the different geospatial themes have little impact on the outcome of the analysis, and it is generally possible to anticipate where stronger correlations will exist. Use of strongly correlated themes can have the effect of “double counting” a specific physical attribute in the model, and in the author’s experience this either has no net impact to model performance, or results in a slight degradation of predictive power. Some care is necessary to select themes with relatively low correlations, and to de-select themes from the model where warranted.

4.5 Image processing and analysis

4.5.1 General

The analytical method used in this work compares landslide presence/absence with other geospatial raster data within a GIS framework. Landslides in this study are represented as raster data with a singular value, with the rasterized landslide footprints having been converted from vector data (i.e. polygons). While not necessary, it is generally most convenient to convert all geospatial thematic data into raster form, for ease of GIS analysis. Some data (e.g. digital elevation models, or DEMs, and their related interpretations) are originated in raster form. Other data (e.g. bedrock, surficial geology) are more typically available in vector form – most often polygons – and must be converted to raster format for analysis. These latter vector-based themes yield rasters with a limited number of discrete values (e.g. eight distinct soil types). Data originating in raster form more often include a continuous spectrum of values (e.g. elevation ranging continuously from 6.572 to 321.297 m). These two different data types are treated differently in the analysis, as discussed in subsequent sub-sections.

Resolution of the various raster data requires consideration. Use of identical resolution rasters using the same spatial extents and grid for all themes simplifies the analysis and reduces the potential for errors; however, this is rarely possible without substantial manipulation, and a resulting risk of loss of data integrity. Data originating as rasters will come with a pre-determined grid size, and these will vary between data sources. It is possible to re-sample each raster to yield the same grid size, and this can be a useful step where it makes sense (e.g. re-sampling a 100 m by 100 m grid to a 20 m by 20 m grid to be consistent with the other data). Where re-sampling will lead to reduction of data quality (e.g. smoothing a 5 m by 5 m grid to 100 m by 100 m), some thought should be given to the desired resolution of the susceptibility model before deciding to re-sample. The effective resolution of the model will be no greater than that of the lowest resolution weight map. Data originating as polygons can be rasterized at any desired resolution; here, the important consideration is to balance the competing demands of desired resolution and required data storage and computational effort.

Twelve different geospatial themes have been considered for use in the susceptibility model. These are classified as either discrete valued or continuous valued data. The discrete valued data include: land use; bedrock type; and, soil type. The continuous valued data include: LANDSAT imagery; stream gradient; seismic hazard; overburden thickness; stream flow accumulation; elevation; stream sinuosity; local relief; and, slope angle. Each theme is discussed individually in the following sub-sections.

4.5.2 Discrete valued data

Discrete valued raster data, which are generally derived from vector data, are treated differently than continuous data. Both data types require a certain amount of pre-processing in preparation for analysis, and this step may require some trial and error. The need for manipulation of discrete valued rasters is generally less than that required for continuous valued data.

Vector (polygon) data of use in susceptibility modelling, like bedrock or surficial geology mapping, usually have numerous attributes assigned to each polygon, yielding choices in the selection of appropriate attributes for discretization in raster form for analysis. The choice of appropriate attribute should be made on the basis that it will have an expected relationship with landslide occurrence, and that its different values will have varying relationships with landslide occurrence. Further, the number of attribute values should be large enough that the model yields some variation in calculated weights (i.e. only one or two attribute values will add little value to the model), but small enough that the results are statistically meaningful. Using soils mapping as an example, one should choose an attribute that discriminates types on the basis of some property related to mechanical behaviour, like texture or grain size. Soils data used by the author in recent work for a Caribbean Island was subdivided into seven categories on the basis of genesis, and these were further subdivided into 59 soil types, based on local names. Use of seven categories in that work was sufficient to yield meaningful results; attempts to use the 59 sub-categories generated weights that could be characterized as statistical noise, rather than meaningful relationships of use in a susceptibility model.

Discrete valued raster data may also be considered for grouping of similar features, either to reduce the number to a manageable and meaningful level, or to eliminate issues with features that have no association with landslides. If a given soil or rock type is not associated with any landslides, the calculated weight value is negative infinity, which has no practical value in the model. With this situation, the modeller has a small number of choices: the attribute can be assigned some arbitrary low (negative) weight, either lower than all other calculated weights, or equal to the lowest calculated weight; or, the attribute can be lumped with another soil or rock type that has a very low weight, resulting in a smaller number of attributes. This latter approach is valid when there are multiple material types with similar character; for example, in this work the available soils mapping has “alluvial terraces” and “erosional alluvial terraces,” which are expected by the author to consist of similar materials, in similar geomorphological settings.

Land use was included in the previous susceptibility model reported by Quinn et al. (2010), since land use mapping is expected to provide a better local association with soil type than provided in the available surficial geology mapping. Land use was once again considered, using data interpreted from Natural Resources Canada (1999). This theme was found to have negligible impact to the overall model, given the addition of other new data, thus its inclusion would result in unnecessary additional complexity. A statistical relationship was shown to exist between land use and incidence of landslides in NTS 31H; however, when the calculated weight was included with the overall susceptibility model, it neither improved nor degraded the model’s predictive power, and was thus considered redundant, and was de-selected from the final model.

Bedrock had been considered and set aside from the previous model, but was reconsidered for possible inclusion in the expanded and updated model. Bedrock data were interpreted from Wheeler et al. (1997). This theme was again set aside because not all rock types in the overall study area are represented in the analytical study area, and over 90 % of the analytical study area contains unclassified sedimentary rocks. No meaningful statistical relationships can be derived due to the lack of variability in mapped bedrock type.

Table 1. Soil types and calculated weight factors.

Large landslides within the study area are expected to be associated with the presence of sensitive clay, which will exist where fine grained marine deposits have been leached by freshwater. Surficial soil types have been interpreted from Natural Resources Canada (1994), which presents a surficial geology map for Canada at 1:1,000,000 scale. The study area includes eight different soil types, classified on the basis of genesis: alluvial terraces; erosional alluvial terraces; lacustrine deposits; bog, wetland; ice contact sand and gravel; marine deposits; till, moraine; and, till, veneer. These are shown in Figure 7. Preliminary analysis suggested the two types of alluvial terrace could be grouped, as could wetlands and lacustrine deposits, and the two types of till. Table 1 shows the relative proportion of the analytical study area occupied by each soil type, along with the corresponding proportion of landslide area. The soil weight map is shown as an example in Figure 8. Intuitively, where the landslide proportion is higher than the relative proportion of area, landslides are more common than average for the study area, so one would expect a positive weight. Where the relative landslide density is lower than average, one would expect a negative weight. The alluvial terraces and marine deposits have the highest weights, both being about 0.7. The similarity suggests that sensitive clay is likely widely distributed in both of these soil units; therefore the alluvial terraces, being alluvial deposits on terraces adjacent to abandoned major rivers, are likely underlain in places by marine soils. The very low negative weight associated with wetland, lacustrine and till deposits indicates that sensitive clays are largely absent in these map units, as expected. The weight factor for ice contact sand and gravel is very close to zero, indicating little relationship with landslides.

Figure 7. Soil types in the study area.

4.5.3 Continuous valued data

Continuous raster data also require a certain amount of pre-processing for the weights of evidence analysis to be as efficient and effective as possible. Each weight map should have a relatively small number of statistically significant weight categories, each corresponding to distinct ranges of values from the continuous raster. The pre-processing required to yield meaningful results may proceed in an iterative fashion. There are at least three ways the continuous rasters can be pre-processed: divide into a number of intervals of some assumed engineering significance (e.g. divide slope angle into 5 or 10 degree increments); divide the map into approximately equal area segments; or, divide the map into segments with approximately equal proportion of landslides. Each of these three approaches can be valid, but some care is necessary to choose the appropriate method that avoids statistical weakness in the resulting weights. Regardless of which approach is taken, it is also necessary to decide how many segments the map will be divided into. In the author’s experience, there should be at least four or five segments, and generally no more than 15 or 20. Selection of the final number of map segments, which determines the number of distinct weights, may be refined over the course of the analysis. For example, one might start with 20 equal-area elevation slices, and the initial analysis may yield groups of similar weights, suggesting grouping the 20 slices into a smaller number of groups of slices. Some of these considerations are discussed in light of specific data examples in the following paragraphs.

Three continuous valued data sets considered for use in the model have been eliminated from the final model, including LANDSAT imagery, stream gradient, and seismicity. Six continuous value data sets have been examined and retained for use in the final model: overburden thickness; stream flow accumulation; elevation; stream sinuosity; local relief; and, slope angle. The latter four continuous valued data sets were all derived from the same source: the Level 1 DEM available from Natural Resources Canada (2007). The overburden thickness data were obtained from Natural Resources Canada (2004) in vector (polygon) format and converted to raster format. Unlike most polygon data, the overburden thickness data have an effectively continuous value scale, with values ranging from 0 to 168 m, in increments of 1 m.

LANDSAT satellite imagery were obtained with coverage of parts of NTS 31H from Natural Resources Canada (2008). These images contain nine bands: a panchromatic band with 15 m pixel size; six multi-spectral bands with 30 m pixel size; and, two thermal infrared bands with 60 m pixel size. Multi-spectral satellite data are used widely for remote sensing, including interpretation of geology and land use mapping, as examples. LANDSAT data have been used by other landslide investigators to interpret indices that may correlate with landslide incidence, including for example the normalized difference water index (NDWI) and normalized difference vegetation index (NDVI). The individual LANDSAT bands, or interpretive indices, which are a combination of selected bands, may have a relationship with landslide incidence, and this was investigated for the present work. The analysis considered the relationship between landslides, each individual band, and the NDVI, where NDVI = (infrared – red)/(infrared + red), as well as several other combinations of individual bands. The strongest relationships were found to be with the thermal infrared bands, and the strength of correlation was similar to that observed for land use data, discussed in the previous section. When included in the overall susceptibility model, LANDSAT data resulted in a marginal improvement in predictive power; however, the effort required to obtain and manipulate data from consistent time periods for the entire study area was not justified by the slight model improvement, and so this theme was excluded from the final model.

An interpreted stream network was obtained from the DEM using hydrology tools available in GIS. Several different attributes were interpreted for this stream network, including stream flow accumulation, gradient, and sinuosity. Flow accumulation and sinuosity were selected for inclusion in the model, and are discussed later. Bjerrum et al. (1969) showed that for a study area in Norway, the occurrence of large landslides in sensitive clay is linked with river behaviour: landslides tend to occur along actively downcutting younger rivers that have not yet reached a stable profile; more mature rivers, having gentler gradients for the same flow as younger rivers, have a more stable thalweg profile, with less active downcutting and fewer landslides. It was therefore expected that an analysis of stream gradient, perhaps in conjunction with stream flow, would show a relationship with landslide incidence. An analysis of landslide incidence in comparison with stream gradient did not yield meaningful relationships, and this is believed to be due to the fact that the interpreted stream network was comprised primarily of short stream segments. Short stream segments result in significant error in interpreted gradient, due to the resolution of the DEM (i.e. 23 m N-S by 11 to 16 m E-W). A 50 m long stream segment, a fairly common length, would traverse about three pixels in the DEM, each with discrete elevation values. Depending on local relief, the calculated gradient could thus be substantially in error, and since the difference in gradient between young, aggressive and old, stable rivers is relatively small, this uncertainty is believed to mask meaningful differences that may exist. It may also be recalled that the expected relationship would be with some combination of gradient and stream flow, not simply gradient. Such combined effects would be better revealed in a multi-variate analysis; the weights of evidence method is not suited to obtain such combined relationships. In any event, the variability in inferred gradient, being dominated by errors associated with DEM resolution, masked any relationship with landslide incidence, and so stream gradient was not included in the final susceptibility model.

One more continuous valued data set with expected relationship with landslides was considered and eliminated from the model: seismic activity. A number of large landslides, including perhaps all of the largest landslides in the study area, are known or believed to have been triggered by major earthquakes. This includes several landslides in the Ottawa area dated about 4500 years before present (Aylsworth and Lawrence, 2003), two large landslides in the Shawinigan area associated with the 1663 Charlevoix earthquake (Desjardins, 1980), a very large landslide at Saint Jean Vianney associated with the 1663 Charlevoix earthquake (Leggett and LaSalle, 1978), and numerous other landslides triggered by the 1663 Charlevoix earthquake (Filion et al. 1991). Figure 9 shows the distribution of seismic hazard in eastern Canada, in this case represented by the peak ground acceleration (PGA) associated with the 1 in 475 year earthquake, or annual probability of 0.0021, roughly equal to a 10 % probability of exceedence in 50 years. These values have been obtained from Halchuk and Adams (2008). This Figure also shows the spatial distribution of recorded earthquake activity in the area. There is widely distributed seismic activity north of the Ottawa River in southwestern Quebec, and a concentration of seismic activity near the Gulf of Saint Lawrence, east of Quebec City. Note that there may be a lower threshold for earthquake triggering of large landslides in the study area. Lefebvre et al. (1992) studied landslides associated with the M5.9 Saguenay earthquake of 1988 and showed the event triggered only local rotational failures, no large landslides. Aylsworth and Lawrence (2003) suggest a threshold of M5.9-6.0 for triggering large landslides in sensitive clay in eastern Canada.

Figure 9. Seismic hazard in eastern Canada: PGA for 475 year earthquake plus recorded earthquakes.

There are several challenges with attempting to relate seismic activity to landslide incidence. First, the analytical study area (i.e. part of NTS 31H) does not contain the full range of seismic hazard values encountered across the predictive study area, and in particular does not contain values from the top of the range, above 0.4 g. Therefore, the resulting weight factors would not properly address areas of high seismic hazard. Second, landslides triggered by major earthquakes can occur a long distance from the epicentre. For example, the large landslides in Shawinigan triggered by the 1663 Charlevoix earthquake occurred approximately 250 km from the earthquake. Therefore, spatial comparison of landslides with local seismic hazard may be of limited importance. Third, the analysis of earthquake triggering should only consider landslides known to be associated with earthquakes. There is unfortunately no way to distinguish between river erosion and earthquake as trigger for the landslides in the digital inventory in NTS 31H. Therefore, any spatial analysis would necessarily compare seismic hazard with incidence of landslides with any trigger, of which the majority are expected to not be associated with earthquake activity. Hence, the relationships between seismic activity and landslides are not evident in the weights of evidence analysis, and therefore seismic hazard has not been included in the susceptibility model. It should be noted, however, that large earthquakes are known to cause large landslides in the study area, and so this exclusion of seismic hazard should be taken as a limitation of the susceptibility model.

Six continuous valued themes have been included in the model – overburden thickness, plus five interpretations from the DEM: flow accumulation, elevation, stream sinuosity, local relief, and slope angle. Overburden thickness data were obtained from Natural Resources Canada (2004), and were treated in the model in nearly the same way as described by Quinn et al. (2010), with similar outcome to the analysis. The map was initially reclassified into 20 slices of approximately equal area, and these were subsequently grouped, on the basis of similar calculated weights, into five categories. The category limits, map area and landslide area proportions, and calculated weights are shown in Table 2. The expected relationships between overburden thickness and landslide occurrence have been described previously by Quinn (2009) and Quinn et al. (2010). The interpreted weights agree with the expected outcome, with high positive weights for intermediate thickness (i.e. > 17.8 m and < 30.8 m), low positive weights for very high or somewhat lower thickness, and negative weights for low thickness, less than 9.8 m.

Table 2. Factor data and calculated weights for continuous valued raster themes,

The expected relationship between flow accumulation and landslide occurrence has also been described previously by Quinn (2009) and Quinn et al. (2010). Table 2 presents the category data and calculated weights for fifteen flow accumulation slices. The map was originally divided into ten slices of approximately equal area; however, roughly 60 % of the landslide area was contained within one of the slices, therefore that slice was further subdivided into four smaller slices, and the next adjacent decile slice was also further subdivided, into two slices of about 5 % map area. With this theme, the slices could be re-grouped based on similar weight values, for ease of analysis, but they have been left as is to allow some discussion of the choices in analysis. As already noted, two of the original ten equal area slices have been further subdivided to provide better resolution in areas with a high concentration of landslide incidence. It should be evident that it would also be possible to group several sets of slices, and there are many possible permutations. Potential candidates for grouping into larger slices include: slices 1 to 3; slices 4 to 7 or 4 to 9; and, slices 11 to 14. In the author’s experience, all of these potential groupings would have minimal impact to the resulting susceptibility model, and would thus be valid means of simplifying the analysis; however, one would be advised to experiment with the analysis to test whether such groupings might deteriorate the model. Landslide incidence tends to be infrequent along streams with low flow, and become progressively more frequent with increases in flow accumulation, as expected. The weight map for flow accumulation is presented as an illustrative example in Figure 10.

Figure 10. Flow accumulation weight map.

The DEM was divided into 20 slices of approximately equal area, and following preliminary analysis, these were re-grouped into four categories. The relevant data and calculated weights are presented in Table 2. Landslides are most common at elevations less than 34 m above sea level, becoming progressively less common with increases in elevation, and very uncommon at elevations greater than 147 m.

Stream sinuosity, as used in this paper, is defined as the ratio of the total length along a stream segment to the direct (i.e. shortest) length between its endpoints. This gives some indication of the amount of curvature along the stream segment. A perfectly straight stream segment will have a sinuosity equal to 1. The relevant data and calculated weights are presented in Table 2. Landslides tend to be infrequent along streams with low sinuosity, and more common along streams with elevated sinuosity.

Landslide incidence is expected to correlate with riverbank height, with very few large landslides occurring where bank height is below some threshold. Obtaining riverbank height directly from the DEM is not possible, since its measurement is subjective, depending on correct choice of the toe and crest, which cannot be automated. An alternative measure – local relief – was used as a proxy for bank height. In this work, local relief is calculated as the difference between the lowest and highest elevations within a set neighbourhood, which was taken as a 500 m square. Table 2 presents the local relief data and calculated weights. It can be seen from the calculated weights that landslides are less common than average where local relief is less than 8 m, increasing thereafter to a peak frequency for local relief up to about 20 m, and thereafter decreasing, with landslides being uncommon for local relief greater than 40 m. Where calculated relief is low or moderate (i.e. < 30 m), it likely reflects a good approximation of bank height in the gentle lowland plains. Where local relief is greater than about 30 m, this may be more indicative of presence of rugged hills in the Appalachians or Canadian Shield above the gentle plains, and absence of marine soils, hence low landslide frequency.

The final theme examined for the model is slope angle. This theme had been considered for the original model, but was previously set aside. The original decision to exclude slope angle was in part due to the lower resolution of the previous DEM, and in part due to the observed lack of correlation between riverbank slope angle and landslide occurrence. The current work relies on higher resolution topographic data, hence finer interpretation of topographic characteristics, including slope angle. Further, the analysis reveals relationships between slope angle and landslide incidence, although not the same kind of relationships typically observed in other work, with landslides usually being strongly associated with moderately steep slope angles. The thematic data and calculated weights are presented in Table 2. A review of the calculated weights suggests that the nine slices of slope angle could be grouped into three categories: < 1.91, < 8.75, and > 8.75 degrees. Landslides are relatively uncommon for slope angle less than 1.91 degrees, relatively common for intermediate slope angles greater than 1.91 and less than 8.75 degrees, and uncommon for slope angles greater than 8.75 degrees. It may be noted that most of the study area consists of gently sloping plains. Steep slopes are locally present along riverbanks, but are far more common in the rugged uplands to the northwest and southeast. Therefore, while some of the landslides’ areas coincide with relatively steep riverbanks, most of the landslides’ footprints extend across gently sloping land, and most of the steeper land is in areas that are not prone to large landslides in sensitive clay, due to the absence of marine soils.

It may be expected that local relief and slope angle will have significant correlation, thus the inclusion of both themes in the susceptibility model could have the effect of “double counting” a specific terrain feature. The decision was made to include both themes after trial and error confirmed the model would be improved, if only marginally, by their mutual inclusion, at little incremental effort.

4.6 Combined weights: landslide susceptibility model

The landslide susceptibility model was obtained by direct addition in GIS, using raster algebra, of the weight maps from the following seven themes: soil type; overburden thickness; flow accumulation; elevation; stream sinuosity; local relief; and, slope angle. While each of the input weight maps contains a small number of distinct weight values, the resulting combined map is effectively a continuous valued raster, with combined weights ranging in small increments from -8.7 to +8.3. A weight around zero implies landslide density equal to the average within the analytical study area; extreme high and low values represent the highest and lowest landslide densities encountered within the study area.

It is necessary to manipulate the combined weight map further to obtain a meaningful model, of use in practical applications. The combined weights raster was first reclassified into twenty equal area slices, and the total landslide area within each slice was then obtained. The combined weights raster was also reclassified into twenty slices with roughly equal proportion of landslides. These two different approaches provide complementary insights into the distribution of landslides at the lower, and higher density ends of the spectrum, respectively. In the case of equal map area slices, the first ten slices (i.e. half of the map) contain roughly 3 % of the landslides, and the final slice (i.e. 5 % of the map area) contains roughly 54 % of the landslides. In the case of the equal landslide area slices, the first slice, containing 5 % of landslide area, covers about 69 % of the map, and the last ten slices (i.e. 50 % of landslide area) cover only 3 % of the map. A combination of results from these two comparisons can be used to show the predictive power of the map, which is illustrated in Figure 11. In this Figure, a “perfect” susceptibility map would isolate only existing landslides, or about 0.1 % of the study area, as being susceptible to landslide activity. This is shown conceptually as the upper dashed line, which rises to 100 % of landslides area within a very small proportion of map area. By contrast, a completely useless susceptibility model, with zero correlation between computed susceptibility and landslide occurrence, can be modelled by a straight line with slope equal to unity. The current and previous susceptibility models’ performance are also shown in the Figure, and it can be seen that the updated model is markedly superior to the original model, with 90 % of landslides being contained within less than 20 % of the map area, as compared with the same proportion of landslides being contained within about 45 % of the map area in the older model. The new model is thus significantly better in isolating areas of landslide activity.

Figure 11. Model predictive power within NTS 31H.

Table 3. Susceptibility values for the landslide susceptibility map.

The combined weights susceptibility model has been reclassified to four broad qualitative categories, as detailed in Table 3: low; low to moderate; moderate to high; and, high. The category limits were chosen on the basis capturing specific percentages of landslide area, as detailed in the Table (i.e. 20, 30, 30 and 20 %). These qualitative landslide susceptibility categories provide a relative indication of the spatial probability of landslides. Table 3 shows the average landslide density within each of the susceptibility categories, as well as the relative density, which here means the ratio of landslide density in the category area with the average density across the study area. The low and high categories, for example, have landslide densities equal to 0.227 and 36.568 times the overall average, respectively. Therefore landslides are more than 100 times more common in high versus low susceptibility areas of the map. The final susceptibility model is shown for the analytical study area in Figure 12, and for the complete study area in Figure 13.

Figure 12. Combined susceptibility map for NTS 31H area.

Figure 13. Landslide susceptibility map for complete study area.

A numeric meaning can be assigned to the four qualitative susceptibility categories through Monte Carlo analysis. Figure 14 shows the distribution of some of 100,000 random points distributed within the NTS 31H landslide inventory area. Each of these points was associated with the underlying susceptibility category, and distance to the nearest existing large landslide was calculated in GIS. From this analysis, it is possible to determine the probability of encountering existing large landslides within any given distance, given knowledge of the susceptibility category. The resulting statistics are illustrated in Figure 15. From the illustrated probability distributions, one can interpret, for example, that at a random point on the map located in a high susceptibility area, there is approximately 60 % probability of encountering an existing large landslide within 1000 m, or 80 % probability of encountering a large landslide within 2000 m.

The statistics illustrated in Figure 15 apply within NTS 31H, but it is important to recall that other geographic, geologic or geomorphic factors may exist elsewhere in the Saint Lawrence Lowlands outside NTS 31H. It is necessary to test the model validity outside NTS 31H prior to using it outside the analytical study area. This has been done both qualitatively and qualitatively. The former approach involved a comparison of documented landslide locations with the susceptibility model, and the quantitative method involved air photo analysis.

Figure 14. Random points for Monte Carlo analysis of susceptibility categories.

Figure 15. Statistical meaning of susceptibility categories.

Figure 16 shows the locations of a large number of landslides in sensitive clay that have been documented in the literature by various authors, as compiled by Quinn (2009). Virtually all of the documented large landslides overlap areas of either high or moderate to high susceptibility. A small number of exceptions occur, primarily along the abandoned banks of the proto-Ottawa River east of Ottawa, where a number of large landslides, believed to have been triggered by major earthquakes, have occurred. It may also be noted that some areas with indicated elevated susceptibility contain no existing landslides, suggesting the model overpredicts landslide incidence in some areas. This issue is less significant than for the previous susceptibility model, as evidenced by the improved predictive performance illustrated in Figure 11. The inclusion of additional thematic data in the model, notably local relief (i.e. a proxy for riverbank height) has eliminated some areas of overprediction; for example, the previous model was not able to distinguish between deeply incised channels (where landslide potential would be higher) from shallower channels with similar flow accumulation.

Figure 16. Locations of documented large landslides as recorded in the literature.

Figure 17 shows 174 pseudo-random points outside NTS 31H used for quantitative model validation. Air photos were reviewed for each of these points to determine the proximity of the nearest recognizable existing large landslide. The compiled results of these observations are shown in Figure 18. The results for categories 2 and 3 (i.e. low to moderate, and moderate to high) showed no evident difference, and have been grouped together. It can be seen that landslides appear to be more common outside NTS 31H. This may be partially explained by the higher seismicity outside NTS 31H, although other factors could also be significant contributors. The general patterns, and qualitative meaning of susceptibility categories, are preserved. This is consistent with prior observations by Quinn (2009). The map can therefore be used outside NTS 31H with the proviso that users recognize landslide frequency may be higher, on average, than observed within NTS 31H.

Figure 17. Pseudo-random points used to validate the model quantitatively outside NTS 31H.

Figure 18. Validation results for locations outside NTS 31H.

5 DISCUSSION: LIMITATIONS AND POTENTIAL IMPROVEMENTS

The updated susceptibility model improves substantially on the previous version reported by Quinn et al. (2010), offering greater predictive power, and expanding to cover the whole area of concern for large landslides in sensitive clay. However, there remain some important limitations that need to be taken into consideration in any future use of the model. First, the model does not address the important linkage between seismic hazard and landslide incidence. Second, the model’s predictive power is significantly better within the analytical study area than elsewhere, since the relationships are calculated within, and extrapolated outside, the analytical study area.

Both weaknesses in the model can be addressed by expanding the landslide inventory to cover the whole study area and re-doing the weights of evidence analysis. While this would be a relatively significant effort, requiring perhaps one person-year for an experienced terrain analyst to expand the landslide database, it would lead to substantially more confidence in the model output. Further, expansion of the landslide inventory to cover the whole study area would allow direct integration of seismic hazard as an input in the model, effectively addressing that important component of landslide hazard and risk.

6 CONCLUSIONS

This paper has described the process of updating and expanding a regional landslide susceptibility model for eastern Ontario and southern Quebec, with description of some important decisions necessary in the analysis, as well as discussion of model limitations and potential improvements. The new susceptibility model involves seven different factors, including: soil type; overburden thickness; flow accumulation; elevation; stream sinuosity; local relief; and, slope angle. The work also included five other geospatial themes, which were excluded from the final model, including: land use; bedrock; LANDSAT remote sensing imagery; stream gradient; and, seismicity. This updated combination of weighted thematic maps led to a new susceptibility model with significantly improved predictive powers, and the model has been expanded spatially to cover the full footprint of the former Champlain Sea where large landslides in sensitive clay may occur.

The analytical method used to develop the susceptibility model – the weights of evidence method – is a bivariate statistical method that can theoretically avoid bias from pre-conceived judgement. However, a large number of practical decisions are necessary in execution of the analysis, and each such decision can introduce bias and affect the outcome. This paper has described different approaches taken at various stages of the analysis, with some commentary on important considerations affecting decisions in the analysis.

While the model appears quite robust, on the basis of both qualitative and quantitative validation, it is known to be inadequate in relation to earthquake-triggered large landslides, that can be expected to occur in the identified elevated susceptibility area, but that may also occur along abandoned banks of the proto-Ottawa and Saint Lawrence Rivers, where susceptibility levels are not always shown as elevated. The model could be improved again in a substantial way by extending the landslide inventory to cover most or all of the Saint Lawrence Lowlands and re-do the analysis.

7 ACKNOWLEDGEMENTS

The analytical work discussed in this paper was funded by Transport Canada, CN Rail and CP Rail. Geospatial analysis was conducted by Neil Ripley, a geospatial specialist with BGC Engineering Inc.

8 REFERENCES

Aylsworth, J.M., and Lawrence, D.E. 2003. Earthquake-induced landsliding east of Ottawa: A contribution to the Ottawa Valley Landslide Project, In: Proceedings of Geohazards 2003, 57-64.

Bjerrum, L., Loken, T., Heiberg, S., and Foster, R. 1969. A field study of factors responsible for quick clay slides. Proceedings of the 7th International Conference on Soil Mechanics and Foundation Engineering, 531-540.

Bonham-Carter, G.F., Agterberg, F.P. and Wright, D.F. 1989. Weights of evidence modelling: a new approach to mapping mineral potential. Geological Survey of Canada Paper 89–9: 171–183.

Dahal, R.K., Hasegawa, S., Nonomura, A., Yamanaka, M., Masuda, T. and Nishino, K. 2008. GIS-based weights-of-evidence modelling of rainfall-induced landslides in small catchments for landslide susceptibility mapping. Environmental Geology, 54: 311-324.

Desjardins, R. 1980. Tremblements de terre et glissements de terrain: Corrélation entre des datations au 14C et des données historiques à Shawinigan, Québec. Géographie physique Quaternaire, 34(3): 359-362.

Eden, W.J. 1967. Buried soil profile under apron of an earth flow. Geological Society of America Bulletin, 78: 1184-1184.

Filion, L., F. Quinty, and Bégin, . 1991. A chronology of landslide activity in the valley of Rivière du Gouffre, Charlevoix, Quebec. Canadian Journal of Earth Sciences. 28: 250-256.

Fransham, P.B., and Gadd, N.R. 1977. Geological and geomorphological controls of landslides in Ottawa Valley, Ontario. Canadian Geotechnical Journal, 14: 531-539.

Halchuk, S., and Adams, J. 2008. Fourth generation seismic hazard maps of Canada: maps and grid values to be used with the 2005 National Building Code of Canada. Geol. Survey of Can. Open File 5813, 32 pgs.

Hodgson, E.A. 1927. The marine clays of eastern Canada and their relation to earthquake hazards. The Journal of the Royal Astronomical Society of Canada, XXI (7): 257-264.

Lebuis, J., Robert, J.-M., and Rissmann, P. 1983. Regional mapping of landslide hazard in Quebec. In: Symposium on Slopes on Soft Clays, Swedish Geotechnical Institute Report No. 17, Linkoping, 205-262.

Lefebvre, G., D. Leboeuf, P. Hornych, L. Tanguay, 1992, “Slope failures associated with the 1988 Saguenay earthquake, Quebec, Canada,” CGJ, Vol. 29, No. 1, pp. 117-130

Leggett, R.F., and LaSalle, P. 1978. Soil studies at Shipshaw, Quebec: 1941 and 1969. Canadian Geotechnical Journal, 15: 556-564.

Morin, L.-G. 1947. La coulee d’argile de Saint-Louis (Comté de Richelieu). Le naturaliste Canadien, 74(5-6): 125-143.

Natural Resources Canada. 1994. Surficial Materials of Canada, Map 1880a, Geological Survey of Canada, Terrain Sciences Division, Scale 1:5,000,000. Available http://gsc.nrcan.gc.ca/urbgeo/ [Accessed September 2006].

Natural Resources Canada. 1999. Canada Land Inventory Level II UTM Digital Data. The Canada Centre for Remote Sensing, Geological Survey of Canada, Scale 1:250,000. Available http://geogratis.cgdi.gc.ca/CLI/frames.html [Accessed September 2008].

Natural Resources Canada. 2004. Urban Geology of the National Capital Area – Drift Thickness, Geological Survey of Canada, Terrain Sciences Division. Available http://gsc.nrcan.gc.ca/urbgeo/ [Accessed September 2006].

Natural Resources Canada. 2007. Canada Digital Elevation Data, Level 1. The Centre for Topographic Information, Scale 1:50,000. Available http://www.geobase.ca/geobase/en/data/cded/description.html [Accessed October 2010].

Natural Resources Canada. 2008. Landsat 7 orthorectified imagery over Canada, Level 1. The Centre for Topographic Information. Available http://www.geobase.ca/geobase/en/data/imagery/index.html [Accessed October 2010].

Quinn, P.E. 2009. Large landslides in sensitive clay in eastern Canada and the associated hazard and risk to linear infrastructure. PhD Thesis, Department of Geological Sciences and Geological Engineering, Queen’s University.

Quinn, P.E., Hutchinson, D.J., Diederichs, M.S., Rowe, R.K., Harrap, R., and Alvarez, J. 2007. A digital inventory of landslides in Champlain clay, Proceedings of the 60th Canadian Geotechnical Conference, Ottawa, 713-720.

Quinn, P.E., Hutchinson, D.J. Diederichs, M.S., and Rowe, R.K. 2010. Regional scale landslide susceptibility mapping using the weights of evidence approach: an example applied to linear infrastructure. Canadian Geotechnical Journal. 47: 905-927.

Quinn, P.E., Hutchinson, D.J. Diederichs, M.S., and Rowe, R.K. 2011. Characteristics of large landslides in sensitive clay in relation to susceptibility, hazard and risk. Canadian Geotechnical Journal. 48: 1212-1232.

Wheeler, J.O., Hoffman, P.F., Card, K.D., Davidson, A., Sanford, B.V., Okulitch, A.V., and Roest, W.R. 1997. Geological Map of Canada, Map D1860A, Geological Survey of Canada.