Forest harvesting is associated with increased landslide activity during an extreme rainstorm on Vancouver Island , Canada

Introduction Conclusions References


Introduction
Landslides triggered by heavy rainfall are a common occurrence on the Pacific coast of British Columbia (BC), Canada.The shallow soils and heavy precipitation associated with the mountainous terrain of the Pacific Northwest coast create favourable conditions for frequent landslides (Schwab, 1983;Guthrie et al., 2010a).Additionally, commercial forest management practices, including timber harvesting and forest road construction, contribute to greater landslide frequency than observed in natural forest areas (Jakob, 2000;Slaymaker, 2000;Guthrie, 2002).There is also general concern about a possible increase in landslide activity in response to timber harvesting and deforestation in different areas worldwide (Rickli and Graf, 2009;Turner et al., 2010;Imaizumi and Sidle, 2012;Muenchow et al., 2012).
As regulated by the provincial government of British Columbia, those involved in timber harvesting must ensure that "the primary forest activity does not cause a landslide that has a material adverse effect on forest resource values" (Wise et al., 2004).Furthermore, timber harvesters are required to provide a "terrain stability hazard map" outlining potentially unstable terrain where landslides may occur.This type of map is often the result of a susceptibility model that predicts the likelihood of landslide occurrence for any given location (Dai et al., 2002).
There are a variety of methods for landslide susceptibility modelling available (Guzzetti et al., 1999).Traditionally in British Columbia terrain mapping approaches have been applied to classify landslide susceptibility using quantitative methods based on the presence and density of landslides in a terrain unit, i.e. in categorized units of terrain with similar slope, surficial material and slope morphology (Howes and Kenk, 1997;Rollerson et al., 2002;Chatwin, 2005;Guthrie, 2005;Schwab and Geertsma, 2010).Typically, the terrain units are represented as polygons.These terrain polygons are interpreted for stability and may over-or underestimate actual probabilities, particularly where the interpretation is qualitative.For a large Pacific coast island, such as Vancouver Island, more detailed approaches utilizing statistical methods to susceptibility modelling for individual cells of a raster data set have also been explored (Chung et al., 2001;Goetz et al., 2011).
The application of statistically based models for landslide susceptibility modelling is common practice (Chung and Fabbri, 1999;Brenning, 2005;Guzzetti et al., 2006;Frattini et al., 2010;Blahut et al., 2010;Sterlacchini et al., 2011;Goetz et al., 2011).Statistical models are typically applied to spatially predict areas susceptible to future landslides.These predictions are based on the assumption that future landslide events will occur under similar preparatory conditions as past events.However, the performance of a susceptibility model relies on the quality of the landslide inventory, the modelling method and on the selection of variables representing preparatory conditions that are used to predict landslides (Demoulin and Chung, 2007).
An advantage to using a statistical model is one's ability to learn more about the landslide controlling processes by creating a model that describes the association of landslides with the preparatory conditions (Dai and Lee, 2002).The knowledge gained can be used to improve predictions wherein the set of preparatory conditions has a well-defined relationship to landslides (Hastie and Tibshirani, 1990;Demoulin and Chung, 2007).In general, models are used to form more effective forest management practices aimed at reducing damages caused by landslides.They analyze the impacts of human activity, such as forest cover conditions (related to timber harvesting) and forest roads, in relationship to other physical (e.g.precipitation) and geomorphological (e.g.topography) influences on landslide occurrence (Turner et al., 2010).Ideally, the relationships between landslides and highly variable temporal and spatial conditions (e.g.precipitation and forest cover) should be investigated based on event-specific models.
In this study, we focus on an extreme weather event on the pacific coast of British Columbia that hit Vancouver Island around 15 November 2006.This storm brought high-intensity rainfall and extreme winds that resulted in over 600 landslides.These landslides encompassed > 5 km 2 and generated > 1.5 × 10 6 m 3 of sediment (Guthrie et al., 2010b).In addition to this landslide inventory, the availability of a high-resolution digital elevation model (DEM) for terrain analysis, cloud-free satellite imagery to capture forest conditions and weather station data associated with the time period leading up to and including this extreme weather event has created an opportunity to explore the controls of rainfall-induced landslides on Vancouver Island using statistical models.
The purpose of this study was to explore the effects of forest harvesting, precipitation, topography and rock type on landslides that were triggered during the November 2006 storm on Vancouver Island.Our methodological approach is based on nonparametric statistical techniques for modelling landslide susceptibility to provide insights into landslide initiation patterns on a regional scale.

Study area
Vancouver Island (31 788 km 2 ) is located off the west coast of mainland British Columbia, Canada (Fig. 1).The mountainous terrain, referred to as the Vancouver Island Ranges, is a sub-range of the Insular Mountains that runs along the Pacific coast and includes the Queen Charlotte Mountains to the north.The island landscape has 2200 m of relief and has been heavily modified during the Pleistocene glaciation, with over-steepened U-shaped valleys, deep coastal fjords and ample post-glacial sediment supply.
The pattern of precipitation is typical of west coast mountain ranges in a temperate climate.The mean annual precipitation ranges from 800 to 1200 mm along the east coast (mountain shadow) sides and increases towards the west coast (windward side) to more than 3000 mm (McKenney et al., 2006).Precipitation is most abundant from early autumn to midwinter, which is generally from October to March.
The mountain ranges of Vancouver Island contribute to various climate conditions that affect the distribution of tree species (Meidinger and Pojar, 1991).The most common tree species is the western hemlock (Tsuga heterophylla (Raf.)Sarg.), which is predominantly located from sea level to 1050 m a.s.l.Red alder (Alnus rubra Bong.), which is known as a pioneer species, is found in areas that have been recently been disturbed (e.g.timber harvesting or landslides; Meidinger and Pojar, 1991).Mountain hemlock (Tsuga mertensiana Bong.), amabilis fir (Pacific silver fir; Abies amabilis (Dougl.Ex Loud.)Dugh) and yellow cedar (Chamaecyparis nootkatensis (D.Don) Spach) are common at higher elevations from 900 to 1800 m a.s.l.The treeline associated with the high mountain peaks of the Insular Mountains is located at 1650 m a.s.l.Western and mountain hemlock in the coastal region are subject to uprooting caused by severe winds, i.e. windthrow.Uprooting is common where there is an impermeable layer in the soil or a high water table that causes trees to be shallow rooted (Ruth, 1964).Windthrow has been observed to contribute to landslide initiation in the Pacific Northwest (Johnson et al., 2000).
In addition to being located throughout most of Vancouver Island, Douglas fir (Pseudotsuga menziesii Mirb.) is also predominantly located in relatively small segments of the rainshadow on the southeast coast (< 150 m a.s.l.).This low elevation area has experienced heavy timber harvesting during the early 20th century, with old growth forests only remaining in parks (Meidinger and Pojar, 1991).Highly productive forestland has been logged since the early 20th century using various techniques.Landslide initiation has increased as a result (Rollerson, 1992;Rollerson et al., 1998;Jakob, 2000;Guthrie, 2002Guthrie, , 2005;;Guthrie and Evans, 2004b;Chatwin, 2005;Pike et al., 2010), but precise estimates of the amount of increase on a landscape scale are still lacking.
Shallow debris flows and debris slides are the most common type of landslides occurring on Vancouver Island.These shallow landslides usually have a maximum depth of 1 m (Guthrie, 2005).The most common triggering mechanism for these landslides is heavy precipitation; however, rapid snowmelt (Guthrie et al., 2010b) and seismic activity are also known triggers of landslides in the area (Rogers, 1980;Van-Dine and Evans, 1992).The debris slides typically occur in till, colluvium and fluvial deposited sediments (Pike et al., 2010).Debris flows are predominantly triggered by an initial failure of a debris slide, which can occur on a channel sidewall or headwall (Brayshaw and Hassan, 2009).Also, flow initiation is much more likely to occur in steep channels than low gradients channels (Brayshaw and Hassan, 2009).

Landslide inventory
The British Columbia Ministry of Environment mapped the inventory of 638 shallow landslides for the winter of 2006-2007 (Fig. 1).These landslides were mapped from 5 m × 5 m resolution SPOT satellite imagery using an automated change detection method comparing scenes from the summer (May to September) of 2006 to the summer of 2007, which is the season when landslides are least frequent (Guthrie et al., 2010b).While the mapped landslides could have occurred anywhere between these dates, it is believed that vast majority occurred during the November 2006 extreme rainstorm (Guthrie et al., 2010b).More detail regarding the mapping and description of these landslides can be found in Guthrie et al. (2010b).
Initiation points digitized from the mapped landslide polygons were used for the subsequent analysis (Fig. 2).The location of the initiation point was manually digitized where the main scarp may be expected.These initiation points were used to approximate the environmental conditions that led to slope failure.Only a single initiation point per landslide was mapped to provide equal treatment of landslides regardless of their size and to reduce the potential effect of spatial autocorrelation in our statistical modelling, which can occur when observations are separated by small distances (Brenning, 2005).The polygons were used to estimate the landslide area.Non-landslide points, which are required for the statistical analysis, were determined by spatial random sampling from the entire study area, excluding the landslide polygons.Overall, 638 observed landslide points and 638 nonlandslide points, a 1 : 1 sampling ratio, was used to conduct this study.

Geostatistically interpolating precipitation
The saturation of forest soils from antecedent precipitation has been identified as an important variable for the prediction of landslide initiation (Crozier, 1999;Jakob and Weatherly, 2003).The relationship of precipitation to the landslide inventory was therefore analyzed using the antecedent precipitation leading up to and including the storm on 15 November 2006(Environment Canada, 2013), which is assumed to have triggered most of the landslides in the inventory (Guthrie et al., 2010b).By observing daily weather station data across the island, it was determined that a 2-week period (2-15 November 2006) would be used to account for these antecedent conditions.These 2 weeks correspond to a period of continuous daily precipitation that preceded the storm.
The actual distribution of rainfall received on the island during the 15 November 2006 storm varied with aspect, elevation, location of storm cells and the wind-driven component of rain.Analysis of meteorological data by Guthrie et al. (2010b) revealed that precipitation varied to more than 200 mm in 24 h.Where rainfall was the only triggering factor, Guthrie et al. (2010b) determined that 88 % of landslides were associated with > 80 mm of rain and 70 % with > 100 mm of rain.In addition, ambient surface temperatures increased by 1-10 • over the same period and, coupled with strong winds, resulted in considerable snowmelt (a rain-onsnow event) for areas were precipitation alone was not responsible for the landslide response.
The weather station data were compiled from a variety of Canadian government sources: Environment Canada (36 stations), BC Ministry of Transportation (3 stations) and BC Ministry of Forests, Lands and Natural Resources Operations (12 stations).Altogether, 51 stations ranging in elevation from 0 m to 580 m a.s.l. were used for the interpolation of precipitation.Three of these stations were located on small islands that are located within 15 km off the east coast of Vancouver Island.Given these sources of precipitation data, hourly rainfall intensity was not included in our analysis because the majority of our weather station data, which came from Environment Canada, did not have hourly measurements.
On Vancouver Island, the precipitation patterns are strongly dependent on the mountainous topography.Therefore, we applied universal kriging to geostatistically interpolate precipitation based on a spatial linear model that ac-counts for a relationship to elevation (Goovaerts, 2000).The association between elevation and precipitation is generally dependent on the spatial scale (Daly et al., 2008).In general, an upscaled elevation model can better capture the effects of air movements through topographic obstacles (Daly et al., 1994;Funk and Michaelsen, 2004;Sharples et al., 2005).Therefore, the spatial scale for precipitation interpolation in this study followed the methods used by Daly et al. (2008), which utilized an upscaled 800 m × 800 m resolution DEM.
The performance of the precipitation interpolation was assessed with leave-one-out cross-validation that estimated the root-mean-square error (RMSE) and mean bias (Goovaerts, 2000).

Representing forest harvesting activities
Forest cover was mapped by applying a supervised maximum likelihood classifier (MLC) of Landsat TM scenes from the summer of 2006 reflecting pre-landslide conditions.The MLC was applied to image data obtained by a tasseled cap transformation (Crist et al., 1986).The training data comprised of 100 (visually interpreted) observations per each desired class.The Landsat scenes were first classified and then mosaicked to cover Vancouver Island.The Landsat data were  TM (earth.google.com).The Landsat imagery used for the classification of forest cover was utilized as ancillary data to ensure the reference conditions were relevant to the time of the classification (i.e. the summer of 2006).The locations for the reference data collection were based on the random spatial sample used for the non-landslide point locations.At each reference location a suitable forest cover class was assigned, e.g.sparse, open, semi-open or closed forest.Additionally, an attribute indicating whether the reference sample was in a location that has been previously logged was included.This attribute was used to determine the proportion of each forest cover class associated with logging.
Since forest harvesting has happened and continues to happen on Vancouver Island, the general forest covers were meant to implicitly relate to forest harvesting activities (see class descriptions in Table 1).The classification we used in this study may not be specific enough to only represent "logging forests".However, the general characteristics of these forest types can be used to relate to forestry practices and investigate the relationships to landslide initiation (Table 1).We based our classification of forest cover on the work of Cohen et al. (1995), who classified different forest stand structures in similar Pacific coast temperate-forest conditions from Landsat TM data.
The classification of forest cover was characterized by comparing the resulting map to BC government forest inventory data.The relative stand density was estimated by crown (canopy) closure values obtained from the vegetation resources inventory (VRI): a geographic data set that represents forest-sampling units with polygons (British Columbia government, 2013a).The VRI represents numerous projected forest conditions for 2013, e.g.stand age, crown closure, dominant forest species.We attempted to represent conditions associated with the year of our study ( 2006) by including only crown closure estimates in the analysis that had a corresponding projected stand age greater than 6 years.
A BC forestry cutblock inventory (British Columbia government, 2013b) was also compared to our forest classes to gain a sense of how well the classification represents forest-harvesting activities.This inventory consists of digitized cutblock polygons across Vancouver Island for harvesting years up to 2013.The cutblock polygons have a corresponding "end of disturbance date" that was used to approximate the recovery period, in years, associated with our forest cover classes.This value was calculated by subtracting the year corresponding to the "end of disturbance date" to the year our landslides occurred, 2006.Consequently, cutblocks harvested after 2006 were not included in the recovery period calculation.The proportion of samples that were located within the cutblock polygons was calculated to highlight the relative association of each forest cover class to forestry activities.Higher proportions indicate a closer association to the forest cutblocks.
Both comparisons to crown closure and cutblock data were based on the sample of landslides and non-landslide points.Since these data sets do not comprehensively cover Vancouver Island, the comparisons to crown closure and recovery period could only be conducted where the data sets overlapped with the landslide sample.Also, due to the lack of completeness and unknown quality of these government inventories, the comparisons to the forest classification only serve to assess the potential of the forest cover classes as a proxy for recovery stage and were not considered for use as predictors for statistical modelling.
The Euclidean distance from forest service roads was calculated to explore the relationship to landslide initiation in the statistical analysis.Since landslides related to service roads in steep terrain can occur in either the fill material or the cut slopes (Jakob, 2000;VanBuskirk et al., 2005), distances upslope and downslope from roads were considered.Also, assuming road-related landslides were initiated in road fill and road cuts, only distances from roads up to 100 m were considered to have a potential influence on landslides, and consequently the distance variable was trimmed by assigning a value of 100 m to all distance observations > 100 m.Forest road data were obtained from the British Columbia Digital Road Atlas (BCDRA), which was last updated in 2007 (British Columbia Government, 2007).Through visual inspection of Landsat scenes overlaid by the BCDRA, it was determined that the forest service roads were best identified based on road surface attributes for "loose" and "rough" surfaces in the BCDRA data set.

Geology and terrain analysis
The geologic conditions were represented with generalized formations of rock type, which included igneous (intrusive and volcanic), metamorphic and sedimentary rock.The steep slopes of the Insular Mountains are supported by mainly coarse-grained intrusive rock made of durable minerals (quartz, feldspars) that are relatively resistant to weathering but are subject to mechanical breakdown (Pike et al., 2010).The finer-grained volcanic rocks (e.g.basalts) are highly susceptible to mechanical and chemical weathering.Landslides can occur in volcanic sequences made up of severely weathered rock material or clay residues (Pike et al., 2010).Landslides on Vancouver Island have been observed to occur more likely in metamorphic rocks types such as schist and amphibolite than in granitic intrusive rocks (Guthrie and Evans, 2004a).Also, increased numbers of landslides on the island have been associated with silt and clay sedimentary complexes (Rollerson et al., 2002).
Terrain analysis of topographic characteristics commonly forms the basis of quantitative landslide modelling (van Westen et al., 2008).Slope angle, elevation, plan curvature, profile curvature and upslope contributing area were used to represent topographic conditions related to landslide occurrence.These conditions were derived from a 0.75 × 0.75 (∼ 20 m × 20 m) resolution Canadian Digital Elevation Data (Canadian Government, 2000) DEM using the standard terrain analysis module implemented in the open-source SAGA GIS (Conrad, 2006).
All of the analysis was performed at a 20 m × 20 m resolution to be consistent with resolution of the DEM that the topographic predictors are derived from.Generally in landslide susceptibility modelling the mapping unit, e.g.pixel size, is based on the DEM resolution (Ayalew and Yamagishi, 2005;Schicker and Moon, 2012).The classified forest cover was resampled from 30 to 20 m spatial resolution using the nearest neighbour method.Proximity to forest service roads was calculated at the 20 m pixel size, and the remaining predictors representing geology and precipitation were also resampled to this resolution in spite of their stronger spatial generalization.

Statistical modelling of landslide susceptibility
The statistical relationship of the potential preparatory conditions to landslide initiation was modelled with the generalized additive model (GAM) with multiple predictor variables.The GAM is a semiparametric extension of generalized linear models, such as logistic regression, that has the flexibility to represent the response's dependence on a predictor variable as either linear or nonlinear (Hastie and Tibshirani, 1990).The advantage of this method over logistic regression is the ability to represent nonlinear effects that are known to exist in many geomorphological analyses (Phillips, 2003;Brenning et al., 2007;Goetz et al., 2011;Muenchow et al., 2012).These supposed nonlinear effects can be modelled with smoothing functions that fit the data without rigid assumptions regarding the dependence on the response (Wood, 2006).A detailed explanation of the numerous smoothing functions can be found in Hastie and Tibshirani (1990) and Wood (2006).In this study, a smoothing function for multiple predictor variables was estimated with penalized thin plate regression splines with the smoothing parameters selected by generalized cross-validation (Wood, 2006).The flexibility of the smoothers was controlled with up to 3 effective degrees of freedom to allow for general modelling of nonlinear terms while also avoiding the potential to overfit.
A binomial GAM can provide a straightforward interpretation of the odds of an event occurring -in our case, a landslide -conditional on different preparatory conditions.The odds of an event that occurs with probability p is p/(1 − p).Odds ratios can be used as a measure of effect size, i.e. the strength of a predictor.The odds ratio is defined as the odds of an event occurring in one type of condition (e.g.open forest cover) divided by the odds of an event occurring in a different type of condition (e.g.closed forest cover).In this study, we utilize the GAM to assess the effect size of the predictor variables while accounting for possible nonlinear relationships.
The association of possible preparatory conditions with landslide initiation was modelled with predictor variables for topography (five quantitative predictors), rock type (one categorical predictor with four classes), forest conditions (one categorical predictor with four classes) and antecedent precipitation accumulation (one quantitative predictor: Table 4).In addition, the interaction between forest cover and precipitation accumulation was explored with a separate model that also included all the other predictor variables and a linear representation of the corresponding interaction term.We in-cluded topography, rock type and precipitation accumulation in our models to account for the differences in slope mechanics and hydrology, while focussing on the association with forest harvesting activities using forest cover and proximity to forest service roads.Odds ratios were estimated from the GAM to assess the empirical effect of these predictor variables on landslide initiation.In general, the effect of a variable was considered strong when the odds ratio was greater than 2, which is equivalent to a 2-fold increase of the odds of an event occurring.
Since there is inherently a dependent data structure in spatial data (Brenning, 2005;Atkinson and Massari, 2011), alternative estimates of the predictors' effects size and confidence limits were obtained by applying a non-overlapping spatial block bootstrap to account for possible spatial autocorrelation within the models (Brenning, 2012).The bootstrap is a resampling-based estimation technique that draws independent resamples (with replacement) from the available data set.Following Brenning et al. (2012), spatial dependencies were accounted for by resampling disjoint subregions or "blocks" that were obtained by m means clustering (m = 100) of the sample (landslide and non-landslide) coordinates.The bootstrap resampling of m out of m blocks was repeated 1000 times, and odds ratios were obtained from a GAM fitted to each bootstrap data set.The bootstrap results were summarized by the mean and percentile-based 95 % confidence intervals for the odds ratios.
Although predictive modelling was not the main objective of this study, we assess performance by spatial crossvalidation estimation to illustrate the general quality of the susceptibility models (Brenning, 2012).This method randomly partitions the samples into k disjoint spatial subsets using k means clustering, trains a model on k − 1 of these subsets, estimates its performance on the remaining subset and repeats this k times until each subset has been used once as a test set.Since the result is dependent on the particular partitioning used, this step is repeated with different random partitions.In this study, a 5-fold spatial cross-validation was repeated 20 times.In each repetition the area under the receiver operating characteristic curve (AUROC; Zweig and Campbell, 1993) and the sensitivity of each model at a specificity of 90 % were estimated; a specificity of 90 % means that (only) 10 % of the non-landslide area is misclassified as susceptible to landslides, creating false-positive predictions.The estimation of sensitivity at this high specificity assesses the ability of a model to spatially predict areas that are highly susceptible to landslide initiation (Brenning, 2012).The differences in model performance estimates were tested using unpaired t tests.
In addition to the statistical modelling, an analysis was conducted to characterize the landslides associated with different forest conditions and the remaining predictors.The density (number of landslides per square kilometre) was calculated for landslides initiated within each forest cover type and within close proximity to forest service roads (< 20 m).
The size distribution corresponding to landslide areas for each forest cover type was presented using a density plot based on Gaussian kernel density estimation.Also, the median and interquartile range (IQR) associated with each predictor was calculated for landslide and non-landslide observations.This univariate analysis was simply exploratory and does not account for other preparatory factors as confounders.

Precipitation distribution
Two-week antecedent precipitation was generally higher in the southwest portions of the island and in the Insular Mountains (Fig. 3a).Precipitation strongly correlated with elevation values from the smoothed DEM with Pearson's correlation coefficient of 0.67 which implies that almost one-half of the variation in precipitation was accounted for by elevation.The residual semivariogram of precipitation after accounting for elevation showed a distance of spatial autocorrelation of ∼ 90 km and almost no nugget effect or sub-scale variation (Fig. 4).Universal kriging predictions consequently combine the overall elevation-related precipitation pattern with local interpolation adjustments in closer proximity to weather stations.The uncertainty of the predicted precipitation values was highest at high elevations and for areas of the island that had a low density of weather stations, which was generally in the north (Fig. 3b).Leave-one-out estimation gave an RMSE of 211 mm and a positive bias of 13 mm, which indicates a slight overestimation of mean precipitation and an acceptable precision considering the high precipitation amounts and large spatial variation (median 481 mm, IQR 613 mm).

Forest cover characteristics
The forest cover classification had an overall accuracy of 93.5 % (±2.2 % with a 95 % confidence interval) and 0.896 kappa coefficient.The open and closed forest classes had the most confusion between the reference data and the classification prediction (Table 2).
Overall, the forest conditions of 13 998 km 2 on Vancouver Island were classified (Fig. 5).The majority (59 %) of the study area was classified as   the closed forest samples were in a location that was logged.The sparse forest samples were least related to logging with a 46 % proportion of samples in logged areas.
The proportion of samples that were located within the forest cutblocks was highest for sparse (33 %) and open forest (31 %) followed by semi-open (16 %) and closed forest (9 %), which indicates that sparse and open forest cover classes had a stronger association to forest cutblocks than the other classes (Table 3; Figs.5c and 6).Additionally, the recovery period estimates for samples that were within the cutblock polygons were lowest for the sparse (2 years) and open forest (11 years) classes (Table 3; Fig. 6).This result reinforces that the sparse and open forest cover types were more associated to relatively recent forest harvesting disturbances than the other classes.The sparse and open forest cover classes were also associated with the lowest crown closure estimates: 20 and 45 % respectively (Table 3; Figs.5d  and 6).The median estimate of crown closure was the same (55 %) for semi-open and closed forest.However, this may be expected because these classes had the greatest classification confusion (Table 2).

Landslide characteristics
Closed   in volcanic rock (62 %), followed by intrusive (28 %), sedimentary (7 %) and metamorphic rock (3 %; Table 4).When looking at the number of landslides initiated in each forest cover class per rock type, volcanic, intrusive and sedimentary rock had a similar distribution to the one for the entire study area.Although closed forest had the largest cumulative area affected by landslides (calculated from the polygons representing the entire landslide area; 2.3 km 2 ) and sparse forest the least (0.43 km 2 ), landslides initiated in sparse forest were associated with the largest landslides in terms of area (median = 8490 m 2 ; Table 5;  According to the univariate analysis, the landslide density in open forest conditions (0.079 landsl./km 2 ) was 2.3 times as great as the density in closed forest conditions (0.035 landsl./km 2 ; Table 5).This was the largest difference between of landslide density amongst the forest conditions.Semi-open forest had a landslide density (0.052 landsl./km 2 ) that was almost equal to the average density found for open and closed forest.
The landslide density in areas adjacent to forest roads (< 20 m) was 30.5 times as high as the general density of landslides for the entire study area.However, regarding this and the other univariate-based findings on the landslide characteristics, topographic and other environmental confounders need to be accounted for in order to confirm the results of this exploratory analysis.

Effect size: odds ratios
Landslide initiation was more likely in close proximity to forest service roads than further away when accounting for confounding variables in the GAM.The estimated odds ratio of landslide initiation at 10 m from the road was 2.1 times as high as at distances of 50 m, all other predictors being accounted for (Fig. 8; Table 6).Thus, while considering influences of topography, forest cover, rock type and precipitation, it appears the classified "forest service roads" were associated with increased landslide initiation up to 50 m from the road (Fig. 8).
While investigating the effects of forest cover on landslide initiation, we found only a strong association with open forest cover and generally weak effects for the remaining forest cover classes, relative to closed forest conditions.The odds of a landslide initiation compared to closed forest cover was estimated to be as high as 2.2 times in open forest, 1.1 times in sparse forest and 1.9 times in semi-open forest, when accounting for the other predictors in the model.
The linear effect of antecedent precipitation accumulation was weak in the model without an interaction term for forest type and precipitation.The estimated odds of landslide initiation at 1500 mm were only 1.3 times as high as the odds at 500 mm, all other predictors being accounted for.However, the model containing a forest-precipitation interaction term revealed substantial differences in the empirical effect of forest type for areas with different antecedent precipitation amounts; conversely, the effect of precipitation depended strongly on forest cover (Fig. 9).While the estimated odds of landslide initiation increased with antecedent precipitation in semi-open and open forest, they were nearly The association with rock type was generally weak.Landslide initiation was more likely in metamorphic rock compared to volcanic rock when accounting for the other predictors (estimated odds ratios of 1.9 for metamorphic rock, 1.5 for sedimentary and 1.1 for intrusive, each relative to volcanic rock).
In terms of topography, the majority of modelled effects on landslide initiation were represented by a nonlinear smoother (Fig. 8).Slope angle had the strongest effect on landslide initiation (Fig. 8).Its nonlinear effect levelled out at slope angles greater than 40 • .The estimated odds of initiation were 35.5 times as high as at a 50 • slope angle than at 10 • .Large upslope contributing areas also had a strong effect on landslide susceptibility (Fig. 8).The estimated odds of initiation were 7.7 times as high as slopes with an upslope contributing area size of 100 000 m 2 than at an area of 1000 m 2 .
Figure 9. Plots illustrating the predicted odds of landslide occurrence as a function of precipitation for each forest cover type.These values were predicted under otherwise average conditions for a GAM that included a linear interaction term for forest cover type and precipitation, as well as partly nonlinear representations of the other confounders (topography, rock type and proximity to forest service road).
The likelihood of initiation was substantially higher where plan curvature, a proxy for subsurface flow conditions, was convergent (Fig. 8).The estimated odds were 4.4 times as high as for convergent plan curvatures (−0.015 rad.per 20 m) than linear (no curvature; 0) surfaces.The odds of initiation at divergent surfaces (0.015 rad.per 20 m) were almost unchanged with an estimated odds ratio of 1.3 compared to linear surfaces.Landslide initiation was also more likely to occur on concave (0.015 rad.20 m −1 ) than convex (−0.015) profile curvatures, which is used to characterize the subsurface acceleration and deceleration of flow down a slope.The odds of initiation on concave slopes were 2.5 times as high as on convex slopes.
A nonlinear effect of elevation was also noticeable above 800 m a.s.l.when accounting for the other predictors (Fig. 8).For example, while the estimated odds of landslide initiation at 300 m was only 1.4 as high as at 800 m, at 800 m a.s.l.they were 1.94 times as high as 1100 m a.s.l. and 8.53 times as high as 1600 m a.s.l.
Overall, the GAMs provided a robust estimate of the effect sizes, which was indicated by the similarity of the odds ratios for the study area sample compared to the bootstrap resamples (Table 6) and the relatively small bootstrap estimated confidence limits.

Model performance and mapping
The predictive performances of the models in this study were very good, achieving similar cross-validation AUROC values of 0.86 (Table 7).Sensitivity at 90 % specificity was 54 % in cross-validation, indicating that the models were able to concentrate around one-half of the landslide incidences in the most susceptible ∼ 10 % of the area.Including an interaction term did not significantly change model performance (unpaired t test, p value = 0.262).
A map of the model prediction outputs has been provided to visualize the modelling effects (Fig. 10).Consistently with the model interpretation, steep slopes, gullies and close proximity to forest roads are slope-scale features that stand out as having relatively high susceptibility to landslide initiation.

The effects of forest cover and service roads on landslide initiation
In this study we explored the regional-scale empirical effects of forest harvesting related activities on landslide initiation.The combined effect of open or semi-open forest cover, interpreted as regrowth areas, and proximity to forest roads was substantial with a 3-4 times increase in the odds of initiation compared to closed forest conditions and distances greater than 50 m from forest service roads, without considering the interaction between precipitation accumulation and forest cover.When this interaction is accounted for, this combined effect was estimated to increase the odds of initiation by 6-9 times at high precipitation accumulation (1500 mm).
Considering that open and semi-open forest types in this landscape can in general be interpreted as stages of forest regrowth after timber harvesting, this result clearly illustrates that conditions related to forest harvesting were associated with an elevated likelihood of landslide initiation during the November 2006 storm.Modelled landslides in the classified open forest type were found to be twice as a high as the odds of initiation in closed forest.These findings are similar to May (2003), Guthrie (2002) and Turner et al. (2010) where mature forest had the lowest landslide densities.Although this study did not explore detailed causal mechanisms of landslide initiation, differences in forest cover are associated with different soil hydrological and mechanical conditions (Sidle, 1992) that may be responsible for the lower likelihood of initiation in closed forest cover.
In particular, forest canopies may increase stability by reducing the intensity of precipitation reaching the surface (Keim and Skaugset, 2003).Interception has been found to capture 30-50 % of annual precipitation in dense forest canopies (Dingman, 1994).Also, evapotranspiration in temperate climates is typically higher for closed forest with ob- served rates that are 5-10 times higher than in exposed soil (Jones, 1997).Deforestation can lead to reduced evapotranspiration, an increase of the amount and intensity precipitation reaching the soil and, as a result, a reduction in the time it takes for the groundwater table to rise, which can lead to slope failure (Johnson et al., 2000).While the protective effect of vegetation may be less pronounced during long-lasting high-intensity precipitation events as observed in this study, we suggest that these effects contribute to both the overall empirical effect of forest cover type and the interaction of forest cover with precipitation in this study.
It was also found that the open forest conditions, representative of young forest, were relatively more susceptible to landslide initiation than sparse forest, representative of recently disturbed forest.This result may empirically support previous findings that there is a lag period between logging activities and increased proneness to landslides (Swanston and Swanson, 1976;Wu et al., 1979;Sidle et al., 2006).Forest harvesting can increase the frequency of landslides during heavy rainfall, with younger forests having more landslides triggered by rainfall events with relatively lower return periods than more developed forests (Imaizumi and Sidle, 2012).A possible explanation for the mechanism involved in this lag period is the deterioration of tree roots, which can increase susceptibility of landslide initiation during high intensity rainfall events (Sidle, 1992).
For the majority of reference data sampled, logging was the most important contributor to the classified forest cover conditions.It was generally observed during the collection of reference data that other site conditions that lead to sparse or open forest were exposed bedrock, thin patches of forest located in high mountain regions and landslide run-out or snow avalanche paths.Thus, although some of the landslides initiated in sparse forest may also be related to timber harvesting activities, other forest disturbances may also play a role in these areas (Guthrie and Evans, 2004b).However, relative to logging disturbances of forest stands in the BC coastal temperate forests, natural disturbances observed at a regional scale are only of a small extent (Pearson, 2010).
Areas adjacent to forest service roads were also found to increase landslide susceptibility.Prior to the November 2006 rainstorm, it was already observed that there were higher densities of landside initiation near (< 20 m) forest service roads for several catchments across Vancouver Is-  land (Guthrie, 2002).Our results suggest that the effect of forest roads on landslides, generalized across Vancouver Island, is up to 50 m from the road.This distance is relatively short compared to the findings of Larsen and Parks (1997) in humid tropical mountains of Puerto Rico and to those of Brenning et al. (2015) in humid Andes of southern Ecuador, who observed the effects of highways on landslides up to 100 and 150 m respectively.An earlier study conducted in the Oregon Coast Ranges, also located in the Pacific Northwest, found most of the road-related landslides initiated during severe winter storms occurred in fill material (Mills, 1997).We believe that this may also be true for road-related landslides on Vancouver Island and would support why we observed the relatively short-distance effect on landslides.In another Coastal Western Hemlock area in the Pacific Northwest of BC, it was observed that most road fill landslides were likely due to loss of strength in foundation soils of older roads built before the 1990s (VanBuskirk et al., 2005).Failure on newer roads built after 1995 was usually associated with inadequate road drainage (VanBuskirk et al., 2005).Broadly speaking, out of all of the factors related to landslide initiation explored in this study, the influence of roads is the most controllable from a forestry practice perspective.Adequate management aimed at improving construction and regular maintenance of road drainage structures is critical to mitigate the influence of forest service roads on landslide occurrence (Slaymaker, 2000;Wemple et al., 2001;VanBuskirk et al., 2005).

Accounting for the effects of possible confounders
We based this analysis on a statistical approach, using the GAM to account for possible confounders (topography, rock type and antecedent precipitation) and potential nonlinear relationships.In general, the topographic predictors (slope, elevation, curvatures and upslope contributing area) were strong at accounting for the spatial variation in this analysis of re-gional landslide initiation.As usual, slope was an outstanding predictor.Slope angle is a well-established spatial predictor of landslides that is essential for all susceptibility modelling, including physically based models (Lee and Min, 2001;Dai and Lee, 2002;Goetz et al., 2011); in addition, previous studies on Vancouver Island have found it an important variable for statistical susceptibility modelling (Chung et al., 2001;Goetz et al., 2011).The modelling also illustrates the regional importance of topographic curvature for landslide initiation.Chances of landslide occurrence were higher for convergent and concave curvatures.Their combined effect resulted in a 10-fold increase in the odds of landslide initiation.These hillslope depressions can be generally described as geomorphic hollows, which have been associated with repeated landslide occurrence in steep terrain (Sidle and Ochiai, 2006).Johnson et al. (2000) also observed in similar study area conditions, just north of Vancouver Island on Prince of Wales Island, that these hillslope depressions have a higher tendency to slope failure.
Several mechanisms may be responsible for linking elevation and landslide initiation in this study.Elevation was used as a surrogate for climatic differences and the associated biogeographic characteristics, such as tree species distribution.We observed that the likelihood of initiation gradually decreased in the subalpine and alpine regions (> 800 m a.s.l.).This climate region has soils that are generally saturated throughout the year with most of the precipitation (20-70 %) falling as snow (Meidinger and Pojar, 1991).We conjecture that the reduced likelihood of landslide initiation may be attributed to the presence of a snowpack developed before the 15 November storm.Snowfall is usually the dominant form of precipitation over 800 m a.s.l. on Vancouver Island (Gutrhie et al., 2010b).However, the role of snowpack and snowmelt on landslide occurrence is not very well known; a better understanding of this relationship may help to improve prediction of landslides (Guthrie et al., 2010b).
A previous study on Vancouver Island by Guthrie and Evans (2004a) found that landslides were 1-4 times more likely to occur in metamorphic complexes containing schists and amphibolites.However, they acknowledged that it was difficult to know the role logging activities had on landslide distribution.In our study, these affects and others (e.g.topography) were accounted for in our GAM models.We similarly observed that the effect of rock type on landslide susceptibility was greatest for metamorphic rock types.Our regional estimate observed that landslides were twice as likely to be initiated in metamorphic rocks than in volcanics, the dominant rock type on the island.Guthrie and Evans (2004a) suggest that the degree of metamorphism associated with complexes containing schists and amphibolites has increased the vulnerability to erosion.

Scope and limitations
It is important to note that we explored the regional characteristics of landslide initiation on Vancouver Island.Thus, the results obtained represent average conditions and characteristics of landslides that occurred during the high-intensity rainfall event.Specific site characteristics of landslide initiation are generally different for catchments across Vancouver Island (Guthrie, 2002).Therefore, we expect that modelled effects on landslide initiation would also vary depending on specific catchment characteristics.The regional modelling approach provides an overview of the underlying effects on landslide initiation to learn more about the corresponding controlling processes, which in turn can be applied to improve model predictions for regional landslide susceptibility.
The landside sampling design may pose some challenges for accurately modelling initiation conditions.The spatial prediction of landslide initiation was modelled based on a single pixel per landslide, ∼ 400 m 2 area, that was mapped in the main scarp; however, initiation of debris flows may also occur in non-head scarp regions.Ideally, the exact point of initiation would be modelled.Additionally, only the pixel associated with initiation was labelled with a particular forest cover type.The surrounding or coarser-scale forest cover conditions were not considered.We observed that many of the landslides in the closed forest cover type occurred at or near the boundary between the other forest types.Therefore, we believe that the effect of sparse forest, open and semiopen forest types relative to closed forest may have been underestimated.
The DEM spatial resolution and quality can influence the prediction of landslide susceptibility (Claessens et al., 2005).The quality of the topographic predictors is generally associated with the quality and accuracy of the DEM and the algorithms used to derive them (Raaflaub and Collins, 2006).Previous works on the effects of DEM resolution on prediction of landslide susceptibility found that better model performance was generally associated with DEM resolutions around 10 m (Tarolli and Tarboton, 2006;Penna et al., 2014).For this study we used a 20 m spatial resolution DEM because it was the only available data set with comprehensive coverage of Vancouver Island.However, in terms of predictive ability, we found that the susceptibility models performed well (AUROCs = 86 %) at this resolution for our regional-scale analysis.
The additive models may have underestimated the effect of antecedent precipitation for a 2-week period on landslide initiation as a result of bias introduced by random measurement error (i.e.regression dilution bias; Carroll et al., 2012).Landslide initiation on Vancouver Island is largely related to high-intensity local precipitation (Guthrie and Evans, 2004a).The antecedent precipitation was interpolated based on sparse and poorly located data; the maximum elevation of the weather stations was only 580 m a.s.l.The RMSE and kriging variance provide an indication of how much local variation we are missing in our estimate of antecedent precipitation distribution.The "true" effect in our additive models would be expected to be larger than the estimated one.The thematic coarseness associated with our representation of geology by rock type may have also contributed to regression dilution.

Conclusions
In this study, the regional effect of forest harvesting activities, topography, antecedent precipitation and rock type on over 600 landslides initiated during a major rainstorm on Vancouver Island, coastal British Columbia, was explored.Using a statistical approach, several factors with a strong effect on landslide initiation were identified.
Forest cover types, as a surrogate for forest harvesting activities, and proximity to forest service roads were found to increase the likelihood of landslide initiation during the major rainfall event.These forestry activities were collectively associated with a 6-to 9-fold increase in the odds of landslide initiation.This finding supports earlier qualitative studies and univariate analysis of landslide factors in the Pacific Northwest that found forestry activities contributed to increased landslide activity.However, the statistical approach in this study improves on previous analysis of landslide factors by estimating the effect of forestry activities, while accounting for possible confounders such as topography, and by considering possible nonlinear relationships in the statistical models (i.e. the GAM).
While it was found that the forest harvesting related variables (forest cover and distance from forest road) significantly contributed to the spatial prediction of the landslide initiation, the effects of the topographic variables (upslope contributing area and concavities in the slope) on this prediction were much stronger.It is an impractical task from a landslide risk management perspective to attempt to mitigate the effects of landslides by controlling topographic conditions all across Vancouver Island.However, the knowledge of areas susceptible to landslide initiation can be used as a guide to direct where mitigation measures should be applied for logging areas.We acknowledge that forest programs and legislation in British Columbia encourage forestry in locations with relatively low landslide susceptibility.Yet, it seems more can be done.In particular, road design, construction and remediation are perhaps the forestry activities that we have the most control over.Our expectation is that this analysis continues to highlight and emphasize the regional effects of forest cover type and service roads on landslide initiation to assist landslide risk management associated with forestry practices.

Figure 1 .
Figure 1.Vancouver Island and the locations of the landslides triggered from heavy rainfall in November 2006.

Figure 2 .
Figure 2. Examples of landslides in the inventory in a variety of forest cover conditions.The white-dotted polygon lines outline the landslide shape, and the solid-yellow circle points show the mapped location of landslide initiation.(a) Debris flows initiated in closed forest and extend into semi-open forest.(b) Debris flow occurring entirely in closed forest.(c) Debris slides initiated adjacent to roads in open and semi-open forest.(d) Debris flow initiated in adjacent to a forest road in open forest and travels through closed forest and semi-open forest conditions.(e) Debris slides initiated in open forest.Bing Maps Aerial © 2013 Microsoft Corporation; data are © 2010 DigitalGlobe Image courtesy of USGS, © 2010 GeoEye and © 2010 the province of British Columbia.

Figure 3 .
Figure 3.A map of interpolated 2-week antecedent precipitation (a) and a map of interpolation uncertainties expressed by the kriging prediction standard deviation (b).

Figure 4 .
Figure 4. Empirical residual semivariogram of 2-week antecedent precipitation after accounting for elevation as well as a fitted spherical semivariogram model.

Figure 5 .
Figure 5.A false-colour composite highlighting forests around Tofino, Vancouver Island, from Landsat TM imagery (a) and the corresponding result from the forest cover classification (b).A BC forest cutblock inventory showing the recovery period relative to 2006 (c).Crown closure (%) from the VRI data set (d).

Figure 6 .
Figure 6.Box-and-whisker plots for characterizing the forest cover classification.(a) Forest recovery years since the last forest disturbance were extracted to sample points that were within the BC forest cutblock inventory (209 samples).(b) Crown closure (%) was extracted to each forest cover type where the sample data and VRI data overlapped (437 samples).The boxes' widths are proportional to the number of observations for each forest cover type.

Figure 7 .
Figure 7. Empirical density function of landslide area for each forest cover type using Gaussian kernel density estimation.

Figure 8 .
Figure 8. Transformation of predictor variables in the generalized additive model without an interaction term.Terms of the form s(predictor) indicate a nonlinear smoothing spline transformation.The effective degrees of freedom (EDF) refer to the flexibility of the smoothers, which were controlled with a maximum 3 EDF.The light grey bounding polygons (and dashed lines for categorical variables) represent confidence bands at a 95 % level.

Figure 10 .
Figure 10.Model-derived landslide susceptibility map for Vancouver Island: overview and (a-e) selected example areas.

Table 1 .
Description for forest cover classes.
obtained from the United States Geological Survey (USGS) Global Visualization Viewer(USGS, 2011).The performance of this classification was assessed from reference data collected by visual interpretation of highresolution DigitalGlobe TM and Spot images from 2003 to 2012 available in Google Earth

Table 2 .
Confusion matrix for forest cover classification.

Table 3 .
Forest cover classification characteristics related to forest cutblock and crown closure inventories.

Table 4 .
Statistical summary of predictor variables for landslide and non-landslide observations within the study area.Only the trimmed distance from forest road variable was included in the modelling. *

Table 5 .
Characteristics of landslide initiation sites with respect to forest cover and forest service roads.

Table 6 .
The effect size of various contrasting values estimated from the two GAM models (with and without an interaction term) for each preparatory condition associated with landslide initiation.The odds ratios were estimated for the complete observations in the study area and using a spatial block bootstrap technique.The mean and percentile-based 95 % confidence limits (C.L.) of the estimated bootstrap odds ratios are presented.

Table 7 .
Median (and IQR)model performance values measured by AUROC and sensitivity at 90 % specificity estimated from spatial cross-validation.