Landslide susceptibility near highways is increased by one order of magnitude in the Andes of southern Ecuador , Loja province

Introduction Conclusions References


Introduction
Anthropogenic denudation has increased exponentially over the last decades and now exceeds natural denudation by several orders of magnitude (Cendrero et al., 2006).Accordingly, human activity has also become a major factor in increasing landslide susceptibility.It affects the intensity and frequency of landsliding globally through global warming and enhanced precipitation (Crozier, 2010), and locally especially through land cover changes (Glade, 2003) and constructions (Dikau et al., 1996).Deforestation and roads in particular seem to enhance the likelihood of mass movement occurrence (Glade, 2003;Goetz et al., 2011).Increased landslide activity in the vicinity of populated areas and transportation infrastructure may lead to human catastrophes (Fassin and Vasquez, 2005) and vast economic damage (Vranken et al., 2013).Tropical mountain ecosystems are especially sensitive and respond quickly to any change (Vanacker et al., 2007).In the Andes of southern Ecuador, accelerated human pressure has increased sediment yield manifold (Molina et al., 2008), and areas adjacent to highways have proved to be most susceptible to landsliding (Muenchow et al., 2012).
Scientific research may, therefore, contribute to preventing such catastrophes by providing policy makers with maps of high landslide potential and information regarding possible implications of construction along hillsides in highmountain areas.Inadequate design of drainage systems and the mechanical destabilization by undercutting and overloading are typical problems related to mountain roads in de-veloping countries but which also affect unimproved roads in developed countries (Sidle and Ochiai, 2006).Empiricalstatistical models of landslide susceptibility are an efficient and effective means of identifying areas of high susceptibility to landsliding, as well as quantitatively measuring the increase in the odds of landslide occurrence that is associated with known or hypothesized risk factors such as highways while accounting for environmental confounding factors (Brenning, 2012a).Among the challenges encountered in the application of these models is the possible presence of nonlinear relationships and the need to deal with spatial autocorrelation, especially when trying to provide confidence intervals for the estimated odds ratios (Brenning, 2012a(Brenning, , 2005;;Goetz et al., 2011;Vorpahl et al., 2012).
The objective of this study is to estimate the difference in landslide susceptibility in close proximity versus greater distance to paved interurban highways in the tropical Andes of southern Ecuador along a climatic gradient.For this purpose, data from a landslide inventory compiled within a 300 m corridor along these highways were analyzed using generalized additive models (GAM) and generalized linear models (GLM) in order to empirically estimate the effect of highways on landslide initiation while accounting for possible topographic, climatic, and geological confounders.The dependence of the relationship between highway distance and landslide initiation on other environmental factors is further examined.

Study area
Our study area is the 300 m buffer on both sides of paved interurban two-lane highways of the Loja and the Zamora-Chinchipe provinces in southern Ecuador (highways Troncal de la Sierra E35 and Transversal Sur E50 (Fig. 1).The buffer zone comprises 88 km 2 within a bounding box of 51 km × 47 km.
The Cordillera Real runs midway through the study area, and as a climate divide it strongly shapes mean annual precipitation patterns within the study area (Beck et al., 2008).Tropical easterlies are forced to ascend the eastern escarpment which results in values > 6000 mm (Rollenbeck and Bendix, 2011).After passing the main ridge, massive föhn walls develop, which leads to dry conditions west of the main divide.The area around Catamayo in particular is characterized by a dry climate (< 400 mm rainfall per year; Rollenbeck and Bendix, 2011).Vegetation patterns reflect this strong climatic gradient.Tropical dry forest formations dominate the west, and páramo formations and tropical mountain and cloud forests the east (Muenchow et al., 2012;Peters et al., 2010).
The area's surficial geology is furthermore divided into two main units, namely metamorphic rocks and sedimentary The numbers along the street refer to the corresponding geological unit (1: unconsolidated rocks; 2: sedimentary rocks; 3: volcanic rocks; 4: metamorphic rocks; 5: plutonic rocks).The area of the detailed map (lower right panel) will be used as a sample area for the visualization of a predictive map in Fig. 5. Precipitation data are taken from the study of Rollenbeck and Bendix (2001).
rocks.The strike of the metamorphic rock is north-south with a slight dip to the east.Moreover, the metamorphic rocks exhibit an orthogonal joint set, and are frequently interspersed with layers of highly weathered phyllite and clay schist (Beck et al., 2008;Muenchow et al., 2012).Metamorphic rocks prevail in large parts of the study area along the highways.The inter-Andean Sierra, in contrast, served as a sediment trap, resulting in conglomeratic and sandstone formations (Litherland et al., 1994;Beck et al., 2008).These sedimentary rocks almost always exhibit a horizontal layering in conjunction with an orthogonal joint set (Muenchow et al., 2012).
Land use is visible throughout the study area, with some differentiation according to local climate (Beck et al., 2008).Many hillslopes have been deforested and converted into pas-ture, while the fertile valleys of the inter-Andean Sierra are used to grow, for example, bananas, coffee, and sugar cane (Pohle et al., 2013;Rodríguez et al., 2013).Podocarpus National Park is a nearby major protected area; however, it only overlaps with our study area locally in the surroundings of the El Tiro pass (Fig. 1).The smaller Reserva Biológica San Francisco (located in the surroundings of the Estación Científica San Francisco in Fig. 1) is near the road but presents land use and deforestation in proximity to the highway.

Data
All landslides occurring within a 300 m corridor along the highways were mapped by M. Schwinn during several months of field work in 2010 as well as from an orthorectified aerial photograph of the year 2000 (scale 1 : 5000, resolution of orthoimage 1 m × 1 m; data source: E. Jordan and L. Ungerechts, Düsseldorf).The following attributes were recorded for all field-mapped landslides: type of movement (classification is in accordance with Dikau et al., 1996), material type (soil, debris, rock), and state of activity (Cruden and Varnes, 1996).These classifications are available for all landslides observed along the highways in the field (843 movements) but not for the other slides within the 300 m buffer.All movements were subsequently digitized as polygons from the aerial photograph.
Landslide initiation points were manually digitized by selecting a point in the central part of the uppermost portion of each landslide polygon in order to represent its detachment zone (Muenchow et al., 2012;Goetz et al., 2011).Overall, 2185 landslide initiation points are available within the 300 m buffer for this study.Since only 47 field-mapped landslides were not visible in the aerial photographs, and given the high image resolution, we consider the inventory to be complete for landslides > 100 m 2 and to present no bias toward areas that were examined more intensively in the field.
We used a photogrammetrically derived digital elevation model (DEM; 10 m × 10 m resolution; data source: E. Jordan and L. Ungerechts, Düsseldorf) to derive several terrain attributes that serve as proxies for landslide-controlling processes.SAGA GIS 2.1.0and the RSAGA package for the statistical software R were used for all GIS operations (Conrad, 2006;Brenning, 2008;R Development Core Team, 2014).Local slope angle (in • ), plan and profile curvature (rad m −1 ; positive values represent a convex shape), and slope aspect were calculated based on local polynomial approximations according to Zevenbergen and Thorne (1987).Sine and cosine transformations were applied to slope aspect in order to express this circular variable using two independent variables that represent north-south and east-west exposure components (Brenning and Trombotto, 2006).Upslope contributing area (in m 2 ) and its average slope angle (in • ; referred to as catchment slope angle) were derived using the multiple-flowdirection algorithm (Quinn et al., 1991), the former being transformed logarithmically (to the base 10) to reduce skew-ness.Upslope contributing area serves as a possible proxy for soil moisture and soil depth, while catchment slope angle may be interpreted as a proxy for destabilizing forces upslope from a location.While the highways in the study area locally modify hillslope geometry through undercutting, these local slope modifications are not visible in the DEM or any of the derived terrain attributes used here.Rollenbeck and Bendix (2011) compiled a mean annual precipitation raster for the study area.By blending weather radar data and meteorological field observations, they reconstructed both the altitudinal and the longitudinal precipitation gradient apparent in the study area (Fig. 1).This is of great value as rainfall is a common trigger of landslides in the area (Muenchow et al., 2012), and previously published precipitation maps (National Weather Service INAHMI; Hijmans et al., 2005) were strongly biased toward lower precipitation levels and failed to represent the study area's complex precipitation patterns.
The geological maps of Loja and Gonzanamá (scale: 1 : 100 000; data source: Mapa Geológico del Ecuador map sheets 56 and 57 published by the Instituto Geográfico Militar of Ecuador, 1975) were digitized to provide information on underlying bedrock types at a general level.Five geological units are distinguished in this study: metamorphic rocks -phyllite, diorite, muscovite, gneiss, quartzite, graphite (Zamora series, Variscan orogeny, Paleozoic); plutonic rocks -granite (Variscan orogeny, Paleozoic); volcanic rocks -andesite, basalt (Laramide orogeny, Tahuin series); sedimentary rocks -conglomerate, sandstone (geological formations Quillollaco, San Cayetano, Trigal, and Loma Blanca; Oligocene to Pliocene); and unconsolidated sediments -fluvial gravel (Holocene).Metamorphic rock is the most abundant rock class in the study area and was therefore used as the reference class in the statistical analysis.With regard to a possible tectonic influence, earthquakes can be ruled out as a major triggering factor in our study area (Muenchow et al., 2012).
Spatial data representing different land use and vegetation types were not available.Areas that fall within urban settlements (according to field observations in combination with the available aerial photograph) were excluded from further analyses.Grazing is common along the roads throughout the study area.Urban and agricultural areas are mostly confined to less inclined areas, and pristine areas to Podocarpus National Park, which only overlaps with our road buffer in the surroundings of the El Tiro pass (Fig. 1; Peters et al., 2010).
Euclidean distances to highways were calculated using standard GIS operations.Given the positional accuracy of highway geometry data relative to other thematic data, we consider areas within up to 50 m of distance from highways to be potentially directly influenced by changes in hillslope geometry and hydrology due to highway construction and maintenance, while areas between 200 and 300 m away from the highways will serve as control areas that are not directly or at least substantially less influenced by the highway.

Statistical analyses
In this study we follow the recommendations of Brenning (2012a) and Goetz et al. (2011), who emphasize the suitability of the GAM and the GLM or logistic regression model for landslide susceptibility modeling compared to alternative approaches such as weights of evidence or machine learning techniques.A GAM and a GLM were therefore used to characterize the empirical relationships between topographic, climatic, and geological predictor variables as well as highway distance and landslide occurrence in the study area.The logistic GLM, or logistic regression, linearly models the logit, i.e., the logarithm of the odds of landslide occurrence, as a function of l linear predictors x 1 , . . ., x L (Hosmer et al., 2013): While the GLM is a well-established tool for landslide susceptibility modeling (e.g., Ohlmacher and Davis, 2003;Atkinson and Massari, 2011), linearity is unrealistic in many environmental modeling situations, and it may limit predictive performance.The GAM as a semi-parametric, nonlinear extension of the GLM has therefore been proposed as a more flexible alternative to the GLM (Goetz et al., 2011).
The GAM replaces (all or some of) the linear terms of the GLM with nonlinear transformation functions s i (Hastie and Tibshirani, 1990): The transformation functions are typically based on spline smoothers whose flexibility can be adjusted using, for example, the Akaike information criterion (AIC) or generalized cross-validation (GCV) procedures (Wood and Augustin, 2002;Hastie and Tibshirani, 1990).
In this study, the GAM and GLM are both used to provide alternative assessments of patterns of landslide occurrence and, in particular, the ratio of the odds of landslide occurrence in proximity of versus distance to the highways in the study area.The GLM was fitted by iteratively weighted least squares and confidence intervals obtained by profiling.To test the null hypothesis that a coefficient equals 0 against the two-sided alternative, a χ 2 likelihood ratio test was used for quantitative predictors, and z tests for the indicator variables representing geological units.The GAM was fitted with spline-type variable transformations of two equivalent degrees of freedom.The GAM implementation in the "gam" package and GLM implementation in the "stats" and "MASS" packages of R version 3.1.1were used (Hastie, 2013;R Development Core Team, 2014;Venables and Ripley, 2002).
In addition, GAM models that incorporate interaction terms of distance to highway and other quantitative predictors were fitted in order to explore a possible additional differentiation of highway-related effects.The interaction term was represented by a bivariate loess smoother based on firstdegree polynomials and using a span of 0.5 (Hastie, 2013).These smoothers estimate, for any combination of values of the two predictors involved, its contribution to the logit, while accounting for the other predictors in the model.Since the bivariate loess smoother implemented in the gam package uses an isotropic kernel, the variable interacting with distance to highway was linearly rescaled so that 95 % of its values fell within a 0-300 value range that is comparable to the range of values of the distance variable.Possible interactions with geological units were examined by fitting separate models for each geological class.
Since none of these GAM and GLM implementations account for possible spatial autocorrelation, confidence intervals for GLM coefficients may be biased, and the GAM may overfit to the training data (Atkinson and Massari, 2011;Brenning, 2005).In order to obtain alternative estimates of the effect size of distance to highway and its confidence interval, a non-overlapping spatial block bootstrap was applied, which accounts for spatial autocorrelation.The bootstrap is a resampling-based estimation method that mimics the process of drawing a new random sample from a population by drawing (with replacement) from the available data set.This general procedure needs to be modified in the case of dependent data (Davison et al., 2003).To account for possible spatial dependence, resampling in this study was performed at the level of sub-regions or "blocks", which were obtained by mmeans clustering (m = 100) of point coordinates.Bootstrap resampling of m out of m blocks was repeated 1000 times.A GAM and GLM were fitted for each bootstrap data set in order to obtain odds ratios of landslide occurrence near the highway (25 m distance) versus distant from the highway (200 m distance) while controlling for all other environmental variables included in the model.Bootstrap mean and percentile-based 95 % confidence intervals were obtained for GAM and GLM.
While predictive modeling was not the primary goal of this study, the predictive performance of the GAM and GLM was assessed by spatial cross-validation estimation (Brenning, 2012a) using the "sperrorest" package (Brenning, 2012b).In r-repeated k-fold spatial cross-validation, the study region is partitioned into k (here: k = 10) disjoint sub-regions using k-means clustering of coordinates in this study (Ruß and Brenning, 2010).One sub-region at a time serves as a test set, the remaining (k − 1) sub-regions being used as the training set.This is repeated for each partition and repeated r times (here: r = 100) to obtain results that are independent of a particular partitioning.The area under the receiver operating characteristic (ROC) curve, or AUROC, was used to assess a model's ability to discriminate landslide initiation points ver- The final learning sample in this study consisted of 2106 out of the 2185 mapped landslide initiation points, and 4177 additional randomly selected points within the study area's 300 m buffer that were located outside of landslide polygons.To avoid bias due to possibly more detailed landslide information near the highway, 13 landslides < 100 m 2 were excluded from the analysis, and landslides located in urban areas or with missing data in one of the predictors were furthermore omitted.
Prior to statistical modeling, descriptive and exploratory analyses of univariate relationships between predictors and response were performed using correlation coefficients and AUROC.

Descriptive and exploratory data analysis
Correlations among predictors were mostly weak to moderate (absolute value < 0.5) according to both Pearson's and Spearman's coefficient, with the exception of local and catchment slope angle with a Pearson's (Spearman's) correlation coefficient of 0.73 (0.72).Precipitation correlated with elevation, slope angle, and catchment slope in the 0.30 to 0.40 s, and plan and profile curvature as well as logarithmic contributing area partly showed correlations in that range as well (absolute values for Pearson's and/or Spearman's coefficient).All other correlations were weaker.
Of the quantitative predictors, distance to highway, local slope angle, catchment slope, and elevation had the strongest univariate discriminatory power according to the ROC analysis (AUROC ≥ 0.58; distance to highway: AUROC = 0.81; Table 1).

Statistical modeling
Based on the models without interactions, the estimated odds ratios describing the landslide susceptibility differences between areas near the highway (at 25 m distance) and distant from it (at 200 m distance) while accounting for all other predictors were of the order of 18-21 according to all estimation procedures, with lower 95 % confidence bounds > 13 (Table 2).The GLM's OR estimate was sensitive to the choice of a distance value for the "distant from highway" category because its assumed linear relationship tends to extend the strong decrease in odds between 0 and 150-200 m distance toward greater distances.Parametrically derived confidence intervals from the GLM were substantially narrower than spatial bootstrap confidence intervals (interval width 6.5 versus 12.9).
Estimated GLM coefficients and effect sizes for "meaningful" increments indicate a predominance of distance to highway over other predictors (Table 3).Elevation, topographic attributes (except slope aspect), and geology were additional important predictors, each with an odds ratio of the order of 2 for meaningful increments in the predictor's values.In contrast, mean annual precipitation was unrelated to landslide initiation.The fitted GAM displayed nonlinear relationships, especially for distance to highway and slope angle (Fig. 3).Susceptibility to landslide initiation decreased steadily over the first ∼ 150 m of distance to highway, where it started to level off.
The explored nonlinear interaction terms of distance to highway with each of the other quantitative predictors suggest that highway-related effects on landslide initiation vary by less than a factor of 2 depending on the values of interacting variables, according to both bootstrap and parametric estimates (Table 4).However, with the exception of the interaction with plan curvature and upslope contributing area, variations of this magnitude can be explained by random sampling variability alone.According to these results, the highway-related odds ratio is 54 % greater on divergent slopes (positive plan curvature) than on convergent slopes (Fig. 4), and 56 % greater where the upslope contributing area is 500 m 2 compared to 5000 m 2 .Differences in highway-related effects between the geological units seem to be more pronounced.However, these are also subject to greater uncertainty due to smaller subsamples and substantial differences between bootstrap and parametric estimates.Highway-related effects appear to be enhanced -and possibly strongly so -in units 13 (Laramide andesite and basalt) and 2 (Holocene fluvial gravel; Table 4).
The GAM-derived landslide susceptibility map in Fig. 5 highlights highway-related effects as well as local topographic modifications.The overall ability of the GAM and GLM models to discriminate landslide initiation points versus stable locations is very good, with slightly improved performances for the GAM (spatial cross-validation AU-ROC 0.853; training set estimation: 0.866) compared to the GLM (0.838/0.850).Training set estimates of AUROC are slightly higher and therefore overoptimistic.

Highway-related landslide hazards: empirical findings
Our results indicate that landslide hazard was strongly increased in close proximity to mountain highways in the Andes of southern Ecuador compared to control areas at ∼ 150-300 m distance.The estimated odds ratio of landslide initiation at 25 m distance versus 200 m distance was 18-21, with lower 95 % confidence bounds > 13 in all analyses while accounting for several topographic, climatic, and geological confounders, but without interaction terms.Spatial bootstrap estimation using the GAM supports the higher odds ratio estimate of 21.2 (95 % confidence interval: 15.5-25.3).This odds ratio furthermore appears to vary to some extent depending on plan curvature, upslope contributing area, and geological unit, according to the analysis of interactions.The estimated increase in landslide hazard near roads -interurban highways in a developing country -is comparable to increases encountered in other studies in proximity to unimproved roads, largely harvest roads.Sidle and Ochiai (2006) indicate an increase by 1 order of magnitude compared to clearcuts and 2 orders of magnitude compared to undisturbed forest land in humid temperate climates, while landslide occurrence near forest roads in the Oregon Coast Range was only doubled to tripled (Miller and Burnett, 2007).In the humid tropical mountains of Puerto Rico, landslide erosion along highways was 5 times higher in proximity (< 85 m) to highways compared to adjacent forests (Sidle and Ochiai, 2006;Larsen and Parks, 1997).In the humid Andes of south- ern Ecuador, in a smaller area that overlaps with a highway segment of the present study, landslide-related material mobilization rates in a human-influenced area were, on average, more than twice as high as in the surrounding natural tropical mountain rainforests (Muenchow et al., 2012).In addition, the mobilization rate in close proximity to the highway was increased by a factor of 2-4 compared to the humaninfluenced background, which was partly related to increased landslide frequency but also larger (and deeper) landslides near the highway (Muenchow et al., 2012).The present study shows that comparable increases in landslide hazard in proximity to an interurban highway in a developing country by at least 1 order of magnitude can, overall, be generalized across a variety of environmental conditions found in this study region.This contrasts with results from the Indian Himalayas, where areas near a highway only showed an increase in landslide hazard of the order of 50 % (Das et al., 2012).
Our results furthermore suggest that highway effects extend up to a distance of ∼ 150 m from the highway.This distance seems to be greater than the influence distance of up to 100 m observed by Larsen and Parks (1997) for Puerto Rican highways, or of up to 50 m for logging roads on Vancouver Island (Goetz et al., 2014).Our observations suggest that this may be due to the growth of landslides as a consequence of successive reactivation.
Nevertheless, the present findings also have limitations related to study design and data quality.In this study, we can be confident that landslide inventory quality was constant throughout the 300 m corridor along the highway regardless of field-assisted mapping or mapping based on 1 m × 1 m aerial photographs.A bias in estimates of highway effects would have been introduced if landslide mapping based on aerial photographs had been less complete than field-based mapping, which was focused on areas near the road.This was not the case in our study.Direct causal mechanisms related to highway construction and design may furthermore be confounded, to a limited extent, with a possibly higher intensity of land use in close proximity to the highways.However, we do not have direct evidence that any of these possible effects are of substantial concern in this particular study region.

Highway-related landslide hazards: causal mechanisms
While a variety of causal mechanism can contribute to an increase in landslide incidence near highways, it has been pointed out that highways in developing countries are particularly vulnerable due to often poor engineering design (Sidle and Ochiai, 2006;Sidle et al., 2006).
In addition to the effects that highways have by undercutting or mechanically overloading hillslopes, the lack of a drainage system in this study area has to be pointed out as an additional feature that may reduce slope stability.
Slides associated with the highway in our study area are also more susceptible to reactivation than slides located at greater distance from the road (Muenchow et al., 2012).Since the highways of the study area have existed in their present design since at least the beginning of the 1960s (Beck et al., 2008), we suggest that the landslide population in proximity to the highway has been growing for half a century, Table 3. Results of generalized linear modeling: coefficient estimates with standard errors (in parentheses), and odds ratios with 95 % confidence intervals for meaningful increments of the predictors.Boldface indicates significant tests at the 5 % significance level of the null hypothesis of the true coefficient being equal to zero against the two-sided alternative.Spatial autocorrelation and multiplicity of tests are unaccounted for.

Variable
Coefficient with only a low rate of landslide "disappearance" due to regrowth of vegetation or, to a small extent, highway repairs and engineering design.While land use in the area is in general believed to increase persistence and reactivation of landslides (Muenchow et al., 2012;Richter, 2009), we believe that this effect is enhanced in particular by the effects of highways, thus limiting the landscape's ability to recover from landslide disturbances.However, the Ecuadorian government modernized the highway between Loja and Zamora during the study period with basic slope-stabilizing and drainage-facilitating measures.Hence, our study provides a baseline to assess the effectiveness of these efforts in future studies.

Topographic and climatic controls on landslides
Overall, areas near the highway, on or at the foot of steep hillslopes, at higher elevations and with small upslope contributing areas in metamorphic and granite bedrock tended to be more likely to initiate landslides.Generally, differences among geological units exerted only a relatively minor influence on landslide occurrence.The slightly increased susceptibility in the metamorphic areas might be due to interspersing layers of highly weathered phyllites and clay schists and a slight dip to the east (Beck et al., 2008;Litherland et al., 1994).Overall, however, geomechanical differences between the general geological units used here appear to result in very minor variations in landslide susceptibility relative to the overwhelming highway effect.Nevertheless, highway-related effects are likely enhanced in the geological units that are comprised of volcanic rocks (Laramide andesite and basalt) and unconsolidated material (Holocene fluvial gravel).The geological units furthermore consist of a variety of subunits, which may obscure more prominent differences among geomechanically more distinct subunits.
As expected for this area (Muenchow et al., 2012;Vorpahl et al., 2012Vorpahl et al., , 2013;;Richter, 2009), slope steepness was an important predictor of slope failure in terms of both local slope angle and mean slope angle of the upslope contributing area.Steep upslope contributing areas may serve as a proxy for mechanical destabilization due to overloading.
Landslides furthermore tended to occur in locations with smaller upslope contributing areas, i.e., near the local ridges.This may seem counterintuitive as hillslope hydrology would suggest that soil saturation is more likely to occur at locations with larger contributing areas (Montgomery and Dietrich, 1994).However, while landslide initiation points are mapped in the uppermost portion of the observed landslide area, additional locations further downslope from the observed initiation points may also be unstable but are not mapped as initiation points.This may result in the observed tendency of mapped landslide initiation points to exhibit smaller upslope contributing areas.In addition, larger upslope contributing areas > 1 ha often correspond to drainage channels, which are less steep and are often associated with surface runoff.
Numerous studies have shown that precipitation (Guzzetti et al., 2008), especially if exceeding certain thresholds (Caine, 1980;Giannecchini, 2006;Jibson, 1989), is one of the most common triggers of (shallow) landslides (Aleotti, 2004) in mountainous regions, including our study region (Muenchow et al., 2012;Vorpahl et al., 2012).Accordingly, the negligible influence of (mean annual) precipitation in our landslide susceptibility models may come as a surprise, especially when considering that the majority of the mapped landslides are shallow debris slides.However, despite the enormous differences in mean annual precipitation across the study area, this variable may not be a good proxy for the frequency of potentially landslide-triggering extreme rainfall events in this area since critical rainfall intensities occur even in the drier parts up to once per year on average (Muenchow et al., 2012).Mean annual precipitation may also be a proxy for land use differences.In addition, the persistence of landslide scars in the landscape may blur the expected rainfall effect since landslides in the drier area are thought to persist longer than landslides in the more humid part, where faster regrowth of vegetation occurs (Peters et al., 2010).Thus, multi-temporal landslide inventories would be needed to account for these differences in persistence in the landscape.
Elevation was another important predictor, which may be a proxy for altitudinal differences in vegetation and land use as well as local trends in landslide-triggering rainfall events that are not captured by the available precipitation data.However, given the substantial shift in climatic regimes in the westeast direction in our study area and the lack of detailed land use and land cover data, separating these effects and their possible interactions should be the object of future studies.

Statistical methodology
The approach of modeling landslide-environment relationships using the GAM and estimating odds ratios using a spatial bootstrap was preferred in this study over the use of the simpler GLM, as well as over attempting to model spatial autocorrelation in a parametric framework.
The GAM is able to account for nonlinear relationships (e.g., Goetz et al., 2011), which may or may not be strong, but which cannot be ruled out in advance in a heterogeneous environment.While the performance difference between GLM and GAM appears to be small, the inability of the GLM to represent nonlinear relationships may produce an important local bias in predicted landslide susceptibility, for example in close proximity to the road, where the distance decay is greatest.The estimation of odds ratios with the GLM for the nonlinear distance variable is furthermore more sensitive to the choice of the control distance that is thought of as being unaffected by the highway, e.g., 200 or 250 m distance.A logarithmic transformation of distance to highway as performed by Muenchow et al. (2012), by contrast, may obscure the transition between highway-influenced and unaffected distances.The comfortable size of the present data set also allowed us to fit this more flexible model type without overfitting to the training data.
Spatial autocorrelation may violate the independence assumption underlying GLMs and GAMs, which in turn would invalidate statistical hypothesis tests and render confidence intervals for model coefficients and derived odds ratios invalid (Dormann et al., 2007).Extensions of the GLM (Venables and Ripley, 2002) and GAM (Wood, 2006) that incorporate parametric representations of residual spatial autocorrelation are, especially for large data sets, often computationally very intensive, and results can depend on the particular implementation as well as the choice of a specific approximation of the model's likelihood function (Venables and Ripley, 2002).The spatial bootstrap, although itself subject to the choice of a suitable spatial block size, was therefore selected in this study as a nonparametric alternative that is transparent and can be combined with complex models such as the GAM.The comparison of bootstrap and ML-based confidence limits for the odds ratio of distance to highway in the GLM suggests that a non-spatial parametric estimation would substantially underestimate the margin of error (interval width 6.5 versus 12.9; Table 2).

Conclusions
Landslide susceptibility was found to be increased by more than 1 order of magnitude in close vicinity to paved interurban highways in the Andes of southern Ecuador.This overwhelming influence of highways, which fades at about 150 m distance, persists along a strong climatic gradient as well as throughout areas with metamorphic as well as sedimentary rock types.Topographic factors are of secondary importance (odds ratios < 2) in modifying highway effects, while road influence appears to be enhanced in geological units with Holocene gravel and Laramide andesite/basalt.Further research is needed to determine the role that land use may play in contributing to or modifying road-related effects.
Model predictions identifying locations that are most susceptible to landsliding can be instrumental in planning mitigation measures in a cost-effective way.The present study may furthermore serve as a baseline for assessing the effectiveness of resulting improvements to engineering design.

Figure 1 .
Figure1.Overview of the study area with mean annual precipitation patterns (top panel), and its location in southern Ecuador (lower left panel).Highways Troncal de la Sierra E35 and Transversal Sur E50 extend in the north-south and east-west direction, respectively.The numbers along the street refer to the corresponding geological unit (1: unconsolidated rocks; 2: sedimentary rocks; 3: volcanic rocks; 4: metamorphic rocks; 5: plutonic rocks).The area of the detailed map (lower right panel) will be used as a sample area for the visualization of a predictive map in Fig.5.Precipitation data are taken from the study ofRollenbeck and Bendix (2001).

Figure 2 .
Figure 2. Landslides occurring along the investigated highways.(a) Typical landslides of the wet metamorphic part of the study area in the east.(b) Typical landslides of the semi-arid, conglomeratic part of the study area in the west.(c) Highway destroyed by landsliding.(d) A highway is cleared from a recent landslide occurrence.

Figure 3 .
Figure 3. Transformation functions of the generalized additive model (GAM) without interaction term.

Figure 4 .
Figure 4. Odds ratio of landslide initiation versus distance to road based on (a) GAMs with interaction between distance and plan curvature, and (b) separate GAMs for the different geological units.Odds ratios are relative to the odds at a distance of 200 m (black dot) at a straight plan curvature and in geological unit 11, respectively.Grey dashed line: GAM without interaction term for comparison.

Figure 5 .
Figure 5. Landslide susceptibility index maps for a portion of the study area (Fig. 1) based on the GAM.

Table 1 .
Descriptive statistics of the data set used for statistical modeling.

Table 2 .
Estimated odds ratios and their 95 % confidence intervals for landslide occurrence distant from roads (200 m distance) versus near the road (25 m distance) using the GAM and GLM.Estimation based on the spatial block bootstrap and the standard parametric approach.
* Parametric confidence intervals not available for the GAM.

Table 4 .
Dependence of the GAM-derived odds ratio of landslide occurrence near versus distant from highway (200 m versus 25 m) on other predictors, expressed as the ratio of odds ratios (with 95 % bootstrap confidence intervals).As an example, the odds ratio expressing the highway effect is 52 % greater at a plan curvature of 0.01 compared to a plan curvature of −0.01 rad m −1 , based on parametric estimation.Estimates whose confidence intervals do not contain the ratio of 1 are printed in boldface.