Surge-type glaciers experience cyclic flow instabilities characterized by alternating periods of slow and fast flow. The geographical distribution of surge-type glaciers has been shown to be distinctly non-random in that they are clustered in some regions yet are completely absent from others. In order to identify factors that influence glacier surging, a number of environmental and glacial attributes were examined for an area in the Karakoram Himalaya, Central Asia. A new GIS-based glacier inventory was produced using a combination of ASTER and Landsat remotely sensed images and paper maps. A total of 150 glaciers were digitized, with 19 of these (12.6%) being classified as of surge-type. Attribute data for 10 glacial and environmental attributes were recorded either during the digitization phase or extracted automatically from the GIS. Simple data visualization techniques revealed a positive correlation between glacier surging and glacier length, area, perimeter size, average width, debris cover, and orientation. The use of univariate logit regression analysis showed that length, area, perimeter, average width, and the heaviest debris cover class showed significant correlation with surging. Multivariate logit regression techniques were employed to show that length, area, average width, and debris cover were all multicollinear, with the strongest statistically modeled relationship using the variable perimeter size. The significance of glacier perimeter on surging may be explained by an increased availability of avalanche-fed snow and debris material which may act as a mass balance proxy. The findings that glacier size (in particular length and perimeter) is most strongly related to surging are consistent with the findings of studies in a number of different regions.

## Introduction

Glacier surging is an internally regulated flow instability characterized by cycles of slow and fast flow, which are accompanied by cycles of alternating negative and positive mass balance. Surge-type glaciers tend to experience long periods of quiescence or almost stagnant ice flow, typically 15–100 years, interrupted by shorter surge periods of 1–10 years where flow increases by as much as 100 times, often resulting in terminus advances of up to several kilometers in a single year (Meier and Post, 1969). Advances by surging glaciers and ice sheets can have catastrophic consequences including glacier lake outburst floods, or jökulhaups (Bruce et al., 1987), and potentially cause significant flood events at nearby communities or hydroelectric installations. Surges by larger ice bodies may also cause the speed up of iceberg release and potentially the blocking of navigable shipping routes or fuel exploration projects (Dowdeswell, 1989). Alley (1991) proposed that surges of the Laurentide ice sheet may have caused the repetitive glaciomarine sediment layers found in Pleistocene deep ocean cores that are known as Heinrich layers. Less than 1% of the world' s total glacier population have been classified as being of surge-type. Despite their rarity, surge-type glaciers are of crucial importance for the study of ice flow dynamics and the nature of fast glacier flow. The identification of surge-type glaciers is also important in order to attempt to predict hazardous events and to remove obscured or absent proxy signals from the long term climate record (Jóhannesson and Sigurðsson, 1998).

The Karakoram region of the Himalayas lies to the northwest of the greater Himalayan Mountains bordering the area of southwest China and the northern part of Kashmir (northern India/Pakistan) ( Fig.1). The glaciers and ice fields of the Karakoram make up the dominant freshwater store in what is generally an arid and drought-prone region (Hewitt, 1997). Meltwater from glaciers and ice fields of the Karakoram constitute the most important contribution to the flow of the Yarkand and Indus rivers, which affect the livelihood of an estimated 130 million people (Hewitt, 1998). The Karakoram Himalaya is an area with a recognized history of glacial surge events (Paterson, 1994). The most recent inventory of this region outlines in detail a number of surge events (Hewitt, 1997) usually described by eyewitness accounts or examination of historical aerial photography archives. A number of recent papers have summarized the nature of surge events at individual glaciers, most prominently being the Bualtar (Gardner and Hewitt, 1990) and Chiring (Hewitt, 1998) glaciers.

There is at present no consensus view as to why some glaciers surge and others do not. Observational evidence during the active phase of surging has led to the view that fast flow results from motion at the glacier bed rather than internal ice deformation (Kamb et al., 1985). This has led to the formation of two opposing views of surge inception based upon the facilitation of fast flow by water trapped at pressure beneath the glacier (Paterson, 1994). These two views have emerged based upon whether the glacier is underlain by soft permeable sediments or hard impermeable (possibly frozen) material or bedrock. The soft-bed view is that surges are promoted by an increase in the deformation of sediment underlying the glacier (Boulton and Jones, 1979; Clarke et al., 1984; Fowler et al., 2001; Truffer et al., 2000, 2001). Build-up in the reservoir zone and steepening during the quiescent phase causes an increase in basal shear stress that is known to be a factor in till deformation (Clarke, 1987). If deformation is uneven, the build-up of pockets of pressurized water may ensue, potentially destabilizing the sediment and triggering a surge (Murray and Dowdeswell, 1992). Clarke et al. (1984) have also proposed the likelihood of this destabilization occurring during the collapse of subglacial sediment drainage paths. The hard-bed; theory is based upon an assumption that fast flow is promoted by pressurized water at the bed (for example, in linked cavities Kamb, 1987). Kamb proposed that the initiation of surge events takes place during a switch between an efficient channelized water drainage system and the less efficient linked cavity system. This theory may be tested more easily than the soft-bed approach as it does not require specific access to the ice-bed interface, rather it predicts that surge-type glaciers are likely to have preferentially lower surface slopes (Kamb, 1987). Direct observations at the bed of Variegated and Black Rapids glaciers, Alaska, have shown instabilities in their drainage systems overlying soft beds and conclude that the key issue in surging is the relationship between subglacial till and internal hydraulics (Harrison and Post, 2003).

Surge behavior may be characterized by longer and slower “Svalbard-type” or shorter and faster “Alaska-type” events (Murray et al., 2003), which have been identified in different geographical areas, in glaciers of a range of sizes, shapes, and forms and are thought to be caused by different processes. Surging glaciers tend not to be evenly distributed around the world's glaciated regions rather, their spatial distribution has been observed to be noticeably non-uniform on both local and global scales (Post, 1969; Clarke et al., 1986). This results in concentrations or clusters of surge-type glaciers occurring in some regions yet being completely absent from others. Such clusters have been observed in the glacier populations of western North America, Yukon Territory, Svalbard, Canada, Iceland, Greenland, and parts of Central Asia. Research has taken place at regional scales using a variety of observational and analytical methods in order to identify the controls on glacier surging. A number of studies have utilized qualitative observations (e.g. Post, 1969) and univariate statistical methods (Clarke et al., 1986; Hamilton, 1992; Hamilton and Dowdeswell, 1996). However, studies have shown that univariate (independent) statistical models are unable to adequately explain the distribution of surging and that multivariate (combined) models are more appropriate. These types of analyses have so far been restricted to statistical samples of glaciers in the Yukon Territory (Clarke et al., 1986, 1991), Svalbard (Jiskoot et al., 1998, 2000), Iceland (Hayes, 2001), and East Greenland (Jiskoot et al., 2003). In this paper we apply multivariate models to a population of surge-type glaciers in the Karakoram Himalaya.

## Methodology

The data used in this analysis are derived from a new GIS-based glacier inventory constructed using a combination of satellite image interpretation, existing inventory information (e.g. Hewitt, 1997), recently published reported surges (e.g. Diolaiuti et al., 2003), and paper maps. Delineation of glacier outlines was performed using ice flow and topographic data derived from four remotely sensed images. Regional trekking maps and the inventory of Hewitt (1997) were used to create a base map of the geographic area, which was then used to order the appropriate georeferenced satellite images to cover the extent of the glacierized area. The images chosen consisted of three Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) scenes each covering an area of 60 × 60km and one Landsat scene of total size 185 × 185km. This corresponds to pixel sizes in the visible near infrared band of 15 and 30m, respectively. ASTER images of this area are available under the Global Land Ice Measurements from Space (GLIMS) project wherein large amounts of glacier-covered Earth have been tracked and recorded by the instrument aboard the NASA satellite Terra's Earth Observing System (EOS). ASTER images have been used previously to identify and analyze glacial surface features in similar areas (e.g. Wessels et al., 2002). The Landsat 7 Enhanced Thematic Mapper (ETM+)) instrument has also become a popular primary data source for glaciological research partly due to the availability of imagery and large spatial coverage (e.g. Dowdeswell and Williams, 1997; Winther et al., 1999; Bindschadler et al., 2001). Images from late summer were used where available, as surface features diagnostic of surging are most easily recognized on bare glacier ice (Dowdeswell and Williams, 1997). Delineation of glacier boundaries and flow lines took place using ASTER images viewed in the visible near infra-red bands 1, 2, and 3N (0.45–0.69µm), and the Landsat image was viewed in the higher resolution panchromatic band 8 (0.52–0.90µm). Gaussian and linear contrast stretching techniques were employed in the delineation of boundary positions. Elevation and surface slope information were derived from Shuttle Radar Topography Mission (SRTM) global land cover data obtained from the Global Land Cover Facility at the University of Maryland, interpolated to a 50-m-resolution digital elevation model.

Glaciers in the study were assigned attribute information throughout the digitization process, most important being a surge-index or a likelihood that a particular glacier is of surge-type. This attribute value assigns a numerical value between 0 and 3 to the probability of surge behavior at each of the glaciers based on the strength of evidence ranging from no evidence of surging (0), ambiguous morphological evidence (1), a likely surge-type glacier based on strong morphological evidence (2) (see Fig.2), and glaciers with a previously observed or reported surge (3) (after Clarke et al., 1986). Our definition of surge-type means that either one or a number of separate surges have been recorded or observed at a glacier, or that morphological characteristics show enough strong evidence to assume that a surge has happened in the past and may happen again. Unfortunately, it is not possible to examine individual surge dynamics using our approach. Surge-type glaciers (index values 2 and 3) were identified from published data sources (Hewitt, 1969, 1997; Gardner and Hewitt, 1990; Diolaiuti et al., 2003) and satellite image interpretation. The following check list was applied to each of the digitized glaciers in order to apply a systematic and reliable methodology for designating surge-index:

Does direct evidence for surge(s) exist within the literature?

Are surface morphological features present such as looped or contorted moraines that are diagnostic of surge behavior (Meier and Post, 1969?) ( Fig.2b).

Are other surface features present that are potentially diagnostic of surging, e.g. potholes or surface lakes (Sturm, 1987?)

Are significant areas of transverse crevassing present? These may be caused by longitudinal extension stresses placed on the glacier, potentially during an active surge phase.

Are pro-glacial features present that are possibly diagnostic of surging, e.g. Usherbreen-type concentric moraine circles (Lefauconnier and Hagen, 1991), mud, or landslides (Gardner and Hewitt, 1990?)

Although some of the features listed above are thought to be diagnostic for surging, they may not be unequivocal, and must therefore be examined within the context of surge evolution (Jiskoot et al., 2000). Uncertainties are inherent in the identification and classification of surge-type glaciers, and for this reason only glaciers with previously reported surges (point 1) and those with very distinct morphological evidence (point 2) were classified as surge-type (surge-index 2 or 3). This is because potholes, surface lakes, and crevassing may be indicative of a number of ice dynamic processes, not just surging, and as such should not be used in isolation to identify surge-type glaciers. A population of 150 glaciers were digitized of which in total 19 (12.6%) were positively identified as being of surge-type ( Table 1). Table 2 provides some descriptive statistics of the resulting inventory. Although only 12.6% of the total population of glaciers in the sample are identified as being surge-type, these glaciers comprise 22.2% of the total glaciated area.

## Table 1

Surge-type glaciers in the Karakoram Himalaya.

## Table 2

Descriptive statistics of glaciers in the Karakoram Himalaya database.

For all glaciers a set of 10 different glacial and environmental variables were generated to test against the likelihood of being classified as surge-type. These variables included glacier length, area, perimeter, complexity, average width, aspect (orientation), degree of debris cover, maximum and minimum surface elevation, and mean surface slope. Following the construction of the GIS inventory, information for the continuous variables glacier length, area, average width, perimeter, and complexity was generated automatically. Glacier length data (center flowline from backwall to terminus) were compiled through measurement of center flowlines added during digitization. Values for glacier complexity and average width were derived by dividing values for perimeter by area and area by length, respectively. Complexity refers to the extent of perimeter in relation to length a complex glacier therefore may have a large number of tributary trunks contributing to a large perimeter with a short main ice flow. Where the statistical relationship between surging and the continuous variables was non-linear, log transformations were employed. The categorical debris cover class was assessed on a scale of 15 according to the amount of debris material observed on the glacier surface during the image analysis stage. The scale used was: (1) very little or no debris (05 of ice surface covered), (2) small amounts of debris material present (515), (3) moderate debris cover (1525), (4) heavy debris cover (2550), and (5) extremely heavy debris cover (more than 50 of ice surface covered). Maximum and minimum glacier elevation data were extracted from DEM cell values at both ends of the digitized flowlines. Mean surface slope (or slope index) data for each glacier were derived by subtracting the minimum surface elevation from the maximum elevation, then dividing by glacier length.

### Multivariate Modelling: Logit Regression

Attempts to understand many processes in the natural world are often aided by modeling the relationships between different components of their systems. The aim of multivariate analyses are to explore the strength of explanatory relationships between a single response variable and a number of predictor variables. However, traditional linear models such as multivariate regression rely on the assumption that variables are independent and normally distributed. Linear models, which are based on ordinary least squares (OLS) regression, are also only applicable to continuous response variables. For our data analysis, logit regression models, fitted as generalized linear models (GLMs), were used. Rather than linking observed response variables and model-predicted values linearly, generalized linear modeling (Aitkin et al., 1989) defines a link function which allows predicted values to be transformed onto a variety of distributions including a binary response. The response and explanatory variable may be either continuous, categorical, or a combination of the two. This is important for certain types of geographical data, which may both be strongly multicollinear and include both continuous and categorical data types. Both univariate linear regression and multivariate logit models have been used previously as analytical and predictive tools to characterize the effects of glacial and environmental characteristics on the incidence of surging (Atkinson et al., 1998; Jiskoot et al., 1998, 2000, 2003).

Multivariate logit models were fitted as GLMs in the generalized linear interactive modeling (GLIM) system (Aitkin et al., 1989; Francis et al., 1993). The logit regression model was fitted in order to examine the strength of statistical relationships between environmental and glacial variables and surging. Presence or absence of surging was modeled as a discrete binary process whereby 0 = “absence” (non surge-type glaciers) and 1 = “presence” (surge-type glaciers). The binary response was then linked to our set of continuous and categorical explanatory variables to be assessed while controlling for the influence of the remaining variables. In order to do this we calculated the probability of occurrence of an event, modeled in terms of the log odds of improvement, called the logit, and expressed with the general form of

where*P*is the probability of a surge event occurring associated with a given observation

_{i}*i*. The logit model is used as a link function in the GLM to produce the logit regression model, given as where

*α*is the base estimate,

*β*are the

_{k}*n*coefficients to be calculated in the model, and

*X*are the

_{ik}*n*values for each of the explanatory variables (Atkinson et al., 1998).

The logit transformation (Equation 1) fitted in the GLM has the effect of calculating a multiple linear regression of the log-odds ratio of each of the explanatory variables. The log-odds ratio can be described in terms of the logarithm of the odds that an individual is a “presence” (is surge-type). The GLM is fitted in a stepwise fashion; that is, bivariate models are fitted between the response variable and each explanatory variable, as the next most significant variable is added to the model at each further step. This means that the effect of each new variable on the model response is generated while controlling for the collinearity effects of previously included variables.

Categorical variables are fitted in GLIM as model “factors”. A baseline is generated from the parameter estimate of the first factor from which successive parameter estimates from each added factor are given as differences from the baseline (Jiskoot et al., 1998). Categorical variables may be aggregated into fewer classes in order to model the significance of groups of categories as the combined effects of closely related groups of variables may be more significant than those of individually modeled categories. For example, the aspect/orientation variable was aggregated into two clusters (N, NW, and S, and NE, E, SE, SW, and W) based upon descriptive results. The results for debris cover were aggregated by comparing modeled significance levels resulting from a range of different combinations of two aggregate classes based on the five original categories. The aggregate combination with the greatest significant reduction in model deviance resulted from two groups: (1) classes 1–4 as a baseline (approximately 0–50% cover) and (2) debris class 5 (extremely heavy debris cover, more than 50% of ice surface covered).

Model significance is determined by examining the reduction (if any) in model deviance as each new variable is added. Model “deviance” is a measure that cannot be used as a direct goodness-of-fit indicator as one cannot compare fitted probabilities and the binary observations (Collett, 1991. However, the difference (reduction) in deviance between models is useful as an indicator of the significance of each of the explanatory variables. The most important explanatory variables are tested by comparing the reduction in deviance from the null model (the simplest possible model where no explanatory variables are included) with the associated loss of degrees of freedom at a given confidence level (in this case 95%). Each variable is added to the model in a stepwise manner allowing assessment of the effects of individual variables and combinations of and interactions between variables. A term is tested by comparison with standard errors so that, at 95% significance level, a variable is significant if its coefficient is greater than around twice the standard error (Atkinson et al., 1998).

## Results

Simple data visualization techniques such as histograms and scatterplots were utilized in the first instance to identify some of the features that distinguish surge-type from non surge-type glaciers. Identification of these features at an early stage is useful in order to aid the reduction of variables during logit regression modeling. Non-normally distributed continuous variables such as length, area, perimeter, and complexity were dealt with using a logarithmic transformation. As mean and standard deviation may be used as measures of central tendency only when referring to normally distributed data sets, median values were calculated from the data prior to the log transformations.

### Glacier Length

Glaciers in the full data set range in length from approximately 1.57km to 56.8km, with a median length of 7.5km. Figure 3 shows that glacier length is log normally distributed and that surge-type glaciers appear to be longer than normal glaciers. The range of length of normal glaciers is 1.57 to 56.8km with a median of 7.3km while surge-type glaciers range from 4.2 to 49.9km with a median length of 13.6km.

### Glacier Area

Figure 4 shows the area distribution of normal and surge-type glaciers. The range of glacier areas in the full data set extends from 6.1 to 1331km^{2} . The range of normal glaciers is also 6.1 to 1331km^{2} , while surge-type glaciers vary in area from 391 to 1248km^{2} . The median area for normal glaciers is 68km^{2} compared with 164km^{2} for surge-type glaciers. Figure 4 shows that area is log normally distributed, and that surge-type glaciers are more likely to have larger total areas than normal glaciers.

### Glacier Perimeter

Glacier perimeter has a direct relationship with area; i.e. glaciers with large areas will have correspondingly large perimeters. The frequency distribution of perimeter values for surge- and non surge-type glaciers is therefore similar to those of the variable area. This suggests that surge-type glaciers are more likely to have large perimeters. The range of perimeter values for normal glaciers is 4.4 to 141.9km compared with a range of 16.0 to 155.9km for surge-type glaciers. The median values are, respectively, 24.6 and 41.7km for normal and surge-type.

### Complexity

The complexity data are normally distributed and as such central tendency may be measured with the mean. Mean complexity for the total population of glaciers is 3.38, with values of 3.35 for normal glaciers and 3.55 for surge-type glaciers. It may be concluded from this that surge-type glaciers are more likely to be “complex.”

### Average Width

The average width variable is calculated by dividing glacier area by length. The data for average width are normally distributed with a mean of 1.1km. The mean for normal glaciers is 1.1km, while surge-type glaciers have an average width of 1.4km. This would suggest that glaciers with greater average widths are more likely to be of surge-type. As width is related to length and area this is a result that we might expect. Later stages within the logit analysis will indicate how much these variables are related and which is the most significant.

### Aspect/glacier Orientation

Figure 5 shows the orientation of normal and surge-type glaciers expressed as percentages of the total glacier population. A higher percentage of normal glaciers are oriented in the N/NE or SW plane compared to surge-type glaciers. A total of 28% of the normal glaciers are oriented either NE or SW in comparison to just 10% of surge-type glaciers. The predominant direction of orientation for surge-type glaciers appears to be N, NW, and S with a total of 63% of glaciers oriented in one of these directions. This compares with 42% of normal glaciers oriented in the same three directions. These results appear to suggest that surge-type glaciers are unlikely to be oriented in the NE–SW plane but are more likely to be oriented along the N/NW–S plane.

### Debris Cover

The frequency distribution of glaciers belonging to each debris cover class are shown in Figure 6. The distribution appears to show that surge-type glaciers are more likely to experience heavier surface debris cover than normal glaciers. Of the 131 normal glaciers in the inventory, 81% are included in the first three debris cover classes (very little or no debris, small amounts, or moderate cover). The percentage of surge-type glaciers in the first three classes is 47%; 16 of the 19 identified surge-type glaciers (84%) lie within the last three classes: moderate, heavy, or extremely heavy debris cover.

### Elevation and Slope

Both maximum and minimum glacier elevation variables are normally distributed and thus do not require a log transformation. Mean maximum elevation of all glaciers is 5136.5m a.s.l. (above sea level) compared with 5123.7 for normal glaciers and 5220.1 for surge-type glaciers. The mean minimum elevations are 3788.9m a.s.l. for all glaciers in the sample, 3761.1 and 3638.7 for normal and surge-type glaciers, respectively. Surface slope index (mean slope per glacier) values range from 0.02 (gently sloping glaciers) to 0.67 (the steepest glaciers) with the mean of the total sample being 0.17. The averages for normal and surge-type glaciers are 0.17 and 0.14. Although this may suggest that surge-type glaciers have less steep slopes, this calculation is strongly influenced by the presence of an outlier among the surge-type glaciers, the steeply sloping Baltbar glacier with a slope index value of 0.57.

## Logit Analysis Results

Logit analysis was carried out in two distinct phases. Firstly the relationship between each individual variable and the incidence of glacier surging was explored using univariate regression models. These models alone are not sufficient to explain the various causes of glacier surging as they may suffer from collinearity effects. This condition occurs where variables are so highly related (for example, glacier length may be predicted by area) that it is difficult to calculate reliable estimates of their individual regression coefficients. Collinearity will not affect the ability of the regression equation to predict the response rather, problems occur if the user wishes to estimate the individual contribution of each explanatory variable. The second phase of analysis therefore involves fitting the most significant univariate models to the multivariate logit model in a stepwise manner. If two variables that are diagnostic for surging are unrelated, then they will not be rejected by the logit model. If they are significant and not related, then they can be added to the model as further explanatory variables.

### Univariate Model Results

Individual logit regression models are provided in Table 3. Four of the eight subclasses of continuous variables provided significant results at a confidence level of 95%. Significant results were achieved for the explanatory variables log length, log area, log perimeter, and log average width. Of these variables, glacier perimeter results in the greatest reduction in model deviance (10.02) from the null model. The parameter estimate of 1.7 indicates that the log-odds of a glacier being surge-type increase with increasing glacier perimeter. Log-transformed glacier length reduces the null model by 8.69 and with a positive parameter estimate of 1.01 is also shown to be strongly related to surging. Log area reduces the null model by 8.49 but with a parameter estimate of 0.661 is not as strongly related to surging as the variables perimeter or length. The average width variable is also a significant term as its univariate model results in a reduction in deviance from the null of 4.578. The strength of this relation (estimate 1.039) shows that the odds of a glacier being surge-type also increase with average width, a result that is expected given the relationship between average width and length and area.

## Table 3

Univariate logit regression models (each explanatory variable is fitted separately).

The aspect/orientation variable, aggregated into two classes based on descriptive results ( Fig.5), resulted in a reduction to the null model of just 0.14 and with neither class reporting parameter estimates greater than two standard errors is shown not to be significantly related to the incidence of surging. The second aggregated debris cover class (containing just those glaciers experiencing heavy debris cover of more than 50% of the ice surface covered) shows a reduction in deviance from the null model of 8.87 with a parameter estimate of 2.639. This indicates that heavy debris cover is significantly and positively related to surging. The elevation variables and mean surface slope show parameter estimates less than half of the modeled standard error and are therefore not significantly related to the incidence of surging.

### Multivariate Model Results

The results of the multivariate logit modeling are shown in Table 4. Univariate analysis showed that the explanatory variable that results in the greatest reduction in model deviance was glacier perimeter. It is this variable that is carried forward to the multivariate analysis. Each of the remaining significant explanatory variables (log length, log area, log average width, and debris cover) were added in a stepwise manner to the log perimeter values in the logit model. The results show that none of the variables added to log perimeter in the multivariate model result in a significant model improvement. The logit model identifies interconnectivity between the variables and has rejected them as additional explanatory factors for an optimum model. The correlation between the “size” variables glacier length, perimeter, area, and average width is perhaps unsurprising given the relationships among them. Longer glaciers are likely to have larger area, longer perimeters, and greater average widths, for example. Univariate analyses showed that these variables may provide a good explanation of surge potential on their own but (as the multivariate model has shown) due to their interdependencies cannot be used as individual parameters in a multivariate logit model.

## Table 4

Results of the multivariate logit regression analysis (^{a}Log debris cover variables aggregated into two classes as outlined in methodology).

The categorical variable class debris cover is also rejected by the multivariate logit model. Debris class 5 was significant in a univariate model but is not following addition to the multivariate model. This means that the multivariate model is detecting collinearity and therefore there is a link between perimeter size and heavy debris cover. Altering model factor parameters to control for different aggregated combinations of the original debris classes failed to bring any changes in the significance of these results.

## Discussion

The most significant problem when conducting large-scale glacier research such as this is the availability and quality of database information. Inventory data quality and quantity is of critical importance when deriving a model response using univariate and multivariate statistical techniques. Glacier inventory data may suffer from inaccuracy or inadequacy (missing data). It may also be the case that data are available in an incompatible format that does not allow for direct comparison with other data types. The ability of multivariate logit modeling to utilize a variety of statistical data types solves this problem to a certain extent. Data quality, however, plays a vital role. Misclassification of glacier outlines and more importantly, of surge-index will significantly reduce the validity of the results. For this reason, a rigorous methodological approach to digitization and identification of surging was applied. Glaciers with ambiguous surge evidence such as transverse crevassing and small surface features were not classified as being surge-type while only those with previously recorded surge evidence and very strong morphological evidence were included. In order to obtain a statistically representative sample of glaciers in the Karakoram, satellite images were used which corresponded with the spatial extent of the four most detailed maps of the region. The three ASTER and one Landsat images covered the majority of this area. Problems occurring when an individual glacier was found to be on or overlapping the edge of an image were avoided by excluding these incomplete glaciers from the inventory. Some glaciers to the eastern extent of the region around the Siachen glacier complex are missing from the inventory due to unavailable or poor quality satellite imagery. Inclusion of further surge-related parameters that are currently unavailable such as accurate geological and lithological data may have improved the analysis.

Our results showed that glacier length, area, perimeter, average width, and debris cover are statistically related to glacier surging. Simple data visualization shows that surge-type glaciers in the Karakoram tend to be longer, larger in area, wider, have longer perimeters, and are more likely to be heavily covered in debris. They also appear to be more likely to be oriented in N, NW, or S directions compared to normal glaciers, which most commonly face N, NE, S, or SW. Univariate analyses allow the strength of these relationships to be tested with the results showing that log-transformed perimeter offers the most significant explanatory variable with a reduction in modeled deviance of 10.02 at a confidence level of 95%. Length, area, and debris cover are also shown to have a significant effect on the incidence of glacier surging, with reduction in modeled deviance from the null of 8.69, 8.49, and 8.87, respectively. Our findings show that mean glacier surface slope appears to have little effect on the incidence of surging.

Our results show the advantage of using a multivariate logit regression analysis approach. Logit analysis can deal with the condition of collinearity which was shown to be evident in several of the size-derived explanatory variables. The use of generalized linear modeling (GLIM) also allowed the combined analysis of continuous and categorical variables leading to the rejection from the logit model of the debris cover variable, implying a link between perimeter and heavy debris cover. Hewitt (1997) identified the fact that many of the surge-type glaciers in the Karakoram are fed by avalanche and debris-fall material. This could explain the positive correlation between surging and heavy debris shown in our results. However, the rejection of debris cover by the logit model suggests that debris cover and perimeter are statistically related (collinear) and therefore cannot be used in combination to provide an optimal logit model of glacier surging.

Gardner and Hewitt (1990) have suggested that the incidence of surging in the Karakoram may be influenced by an unusual build-up of avalanche debris material leading to a complex thermal layering effect as a result of debris-rich ice horizons. It is theorized that this may result in a build-up of deformable bed sediments and/or unstable transitions between thawed and frozen bed conditions. Our results support the link between surging and avalanche-feed material however, we believe that the effect of avalanche material is more likely to act as a mass balance proxy providing increased snow, ice, and debris material to the mass input of the glacier rather than directly affecting processes at the ice-bed interface.

After perimeter, glacier length is the variable most strongly correlated with surging, a conclusion which is consistent with those reached by Clarke et al. (1986), Hamilton and Dowdeswell (1996), Atkinson et al. (1998), and Jiskoot et al. (1998). Hewitt (1997) acknowledged the importance of length as a control on surge-type glaciers but conceded a lack of data from his study to compare to others. This study explores a glacier population approximately the same as his, and shows results that support the previously published research for a number of different glaciated regions.

The relationship between length and surging has been explored by Jiskoot et al. (2000) who suggested a link between subglacially derived sediment and surge behavior. Longer glaciers have more time to erode particles, therefore bed grain size would be expected to be finer over the same geology. As fine-grained sediments may be more permeable and flow instability is directly related to permeability and till thickness, it is possible that long glaciers overlaying fine-grained beds will be more likely to surge (Boulton, 1996). Longer glaciers may also overlie one or more geological or lithological boundary resulting in a change in substrate that may have implications for the stability of the system. Longer glaciers are more sensitive to longitudinal variations between the upper and lower zones and will therefore experience larger stress gradients which could lead to damming effects and the development of trigger zones for surges (Robin and Weertman, 1973). The separation of cold and warm ice in the thermal regime of polythermal-type glaciers, a factor which is related to their size, has also been proposed as a causal factor in surging (Murray et al., 2000).

The geology of the Karakoram is dominated by thrust faulting with the central band of the glaciated region being composed of a mix of both igneous and sedimentary units. Hewitt (1997) concludes that no distinctive relationship of surging to lithology, as indicated in other regions, has yet been identified. Most of the glaciers cross two or more geological formations (Hewitt, 1997) and are distributed within highly active tectonic zones that display globally extreme rates of uplift and denudation (Searle, 1991). It is possible that these factors play some role in the dominant surge mechanism. However, the geological information required to test this hypothesis was not available.

## Summary and Conclusions

We have derived a new GIS-based inventory of glaciers of the Karakoram Himalaya. Environmental and glacial attributes derived from this inventory have been tested for the incidence of surging. We have displayed the advantage of multivariate logit regression over conventional regression techniques by testing a number of different data types against surging while controlling for interdependencies and the confounding effects of collinear variables. Analysis of the results allows us to derive the following conclusions:

Nineteen of a total inventory of 150 glaciers in the Karakoram were identified as surge-type from published material and satellite image interpretation. This represents 12.6% of the glaciers in the region (compared to the concentration of surging glaciers worldwide, thought to be around 1%) and 22.2% of the total glaciated area.

Simple data visualization techniques showed that surge-type glaciers are more likely to be longer, larger in area, wider, have longer perimeters, and be heavily covered in debris. They are also more likely to be oriented in the N/NW–S plane and unlikely to be oriented in the NE–SW plane.

Univariate regression shows that log-transformed values for length, area, perimeter, average width, and debris cover are significant in relation to surging. The models most strongly related to surging are perimeter, heavy debris cover, and length (with a reduction in deviance from the null of 10.02, 8.87, and 8.67, respectively). Glacier maximum and minimum elevations and mean surface slope show no statistically significant correlation with surging.

Multivariate logit modeling shows that a number of continuous and categorical variables added in a stepwise manner to the most significant variable perimeter are collinear (area, length, average width, and heavy debris cover). This interdependency means that they cannot be used in combination with perimeter to provide a single explanatory multivariate model of surging. Therefore, the strongest relationship to surging between any modeled variable or combination of variables is explained by glacier perimeter.

The significance of the “size” variables, in particular glacier length, is consistent with the findings of previous research on the controls on glacier surging in areas including Svalbard and the Yukon Territory.

The relationships shown between surging and: (i) perimeter, and (ii) heavy debris cover support the observations of Gardner and Hewitt (1990). However, we suggest that the surge potential resulting from increased avalanche material may be the result of proxy contributions to glacier mass balance rather than having a direct effect on processes at the bed.

## Acknowledgments

This study was undertaken while the first author was studying the M.Sc. program in Geographical Information Systems at the School of Geography, University of Leeds. Thanks are offered to Stuart Barr and Hester Jiskoot for their advice. ASTER and Landsat ETM data are distributed by the Land Processes Distributed Active Archive Center (LP DAAC), located at the U.S. Geological Survey (USGS) Center for Earth Resources Observation and Science (EROS) ( http://LPDAAC.usgs.gov). Thanks are given to Adrian Luckman for processing SRTM data. This data was made available courtesy of the USGS via the Global Land Cover Facility (GLCF) ( http://www.landcover.org). Suzanne Anderson, Will Harrison, and an anonymous reviewer suggested valuable improvements to the manuscript. The Karakoram Himalaya glacier inventory is available free of charge to interested parties for non-commercial usage. Please contact the corresponding author for further details.

## References Cited

- M. Aitkin, D. Anderson, B. Francis, and J. Hinde . 1989. Statistical Modelling in GLIM Oxford Oxford University Press. Google Scholar
- R. B. Alley 1991. ?Deforming-bed origin for the southern Laurentide till sheets. Journal of Glaciology 37:12567–76. Google Scholar
- P. Atkinson, H. Jiskoot, R. Massari, and T. Murray . 1998. Generalized linear modelling in geomorphology. Earth Surface Processes and Landforms 23:1185–1195. Google Scholar
- R. Bindschadler, J. A. Dowdeswell, D. Hall, and J. G. Winther . 2001. Glaciological applications with Landsat 7 imagery: early assessments. Remote Sensing of Environment 78:150–162. Google Scholar
- G. S. Boulton 1996. Theory of glacial erosion, transport and deposition as a consequence of subglacial sediment deformation. Journal of Glaciology 42:14043–62. Google Scholar
- G. S. Boulton and A. S. Jones . 1979. Stability of temperate ice caps and ice sheets resting on beds of deformable sediment. Journal of Glaciology 24:9029–43. Google Scholar
- R. H. Bruce, G. A. Cabrera, J. C. Leiva, and L. E. Lenzano . 1987. The 1985 surge and ice dam of Glacier Grande del Nevado sel Plomo, Argentina. Journal of Glaciology 33:113131–132. Google Scholar
- G. K. C. Clarke 1987. ;Fast glacier flow ice streams, surging and tidewater glaciers. Journal of Geophysical Research 92:B98835–8841. Google Scholar
- G. K. C. Clarke 1991. Length, width and slope influences on glacier surging. Journal of Glaciology 37:126236–246. Google Scholar
- G. K. C. Clarke, S. G. Collins, and D. E. Thompson . 1984. Flow, thermal structure and subglacial conditions of a surge-type glacier. Canadian Journal of Earth Sciences 21:232–240. Google Scholar
- G. K. C. Clarke, J. P. Schmok, C. S. L. Ommanney, and S. G. Collins . 1986. Characteristics of surge-type glaciers. Journal of Geophysical Research 91:B77165–7180. Google Scholar
- D. Collett 1991. Modelling Binary Data London Chapman and Hall. Google Scholar
- G. Diolaiuti, M. Pecci, and C. Smiraglia . 2003. Liligo Glacier (Karakoram): a reconstruction of the recent history of a surge-type glacier. Annals of Glaciology 36:1168–172. Google Scholar
- J. A. Dowdeswell 1989. On the nature of Svalbard glaciers. Journal of Glaciology 35:120224–234. Google Scholar
- J. A. Dowdeswell and M. Williams . 1997. Surge-type glaciers in the Russian High Arctic identified from digital satellite imagery. Journal of Glaciology 43:145489–494. Google Scholar
- A. C. Fowler, T. Murray, and F. S. L. Ng . 2001. Thermally controlled glacier surging. Journal of Glaciology 47:159527–538. Google Scholar
- B. Francis, M. Green, and C. Payne . 1993. GLIM4: The Statistical System for Generalised Linear Interactive Modelling Oxford Clarendon Press. 82. Google Scholar
- J. S. Gardner and K. Hewitt . 1990. A surge of the Bualtar Glacier, Karakoram Range, Pakistan: a possible landslide trigger. Journal of Glaciology 36:123159–162. Google Scholar
- G. S. Hamilton 1992. Investigations of surge-type glaciers in Svalbard. Ph.D. dissertation. Scott Polar Research Institute, Cambridge. Google Scholar
- G. S. Hamilton and J. A. Dowdeswell . 1996. Controls on glacier surging in Svalbard. Journal of Glaciology 42:140157–168. Google Scholar
- W. D. Harrison and A. S. Post . 2003. ?How much do we really know about glacier surging. Annals of Glaciology 36:1–6. Google Scholar
- K. Hayes 2001. Controls on the incidence of glacier surging. Ph.D. dissertation. University of Leeds. Google Scholar
- K. Hewitt 1969. Glacier surges in the Karakoram Himalaya, Central Asia. Canadian Journal of Earth Sciences 6:1009–1018. Google Scholar
- K. Hewitt 1997. Recent glacier surges in the Karakoram Himalaya, South Central Asia. American Geophysical Union ( _http://www.agu.org/eoselec/97106e.html). Google Scholar
- K. Hewitt 1998. Catastrophic landslides and their effects on the Upper Indus streams, Karakoram Himalaya, Northern Pakistan. Geomorphology 26:1–347–80. Google Scholar
- H. Jiskoot, P. Boyle, and T. Murray . 1998. The incidence of glacier surging in Svalbard: evidence from multivariate statistics. Computers and Geosciences 24:4387–399. Google Scholar
- H. Jiskoot, T. Murray, and P. Boyle . 2000. Controls on the distribution of surge-type glaciers in Svalbard. Journal of Glaciology 46:154412–422. Google Scholar
- H. Jiskoot, T. Murray, and A. Luckman . 2003. Surge potential and drainage-basin characteristics in East Greenland. Annals of Glaciology 36:142–148. Google Scholar
- T. Jóhannesson and O. Sigurðsson . 1998. Interpretation of glacier variations in Iceland. Jókull 45:27–34. Google Scholar
- B. Kamb 1987. Glacier surge mechanisms based on linked cavity configuration of the basal water conduit system. Journal of Geophysical Research 92:B99083–9100. Google Scholar
- B. Kamb, C. F. Raymond, W. D. Harrison, H. Engelhardt, K. A. Echelmeyer, N. Humphrey, M. M. Brugman, and T. Pfeffer . 1985. Glacier surge mechanism: 19821983 surge of Variegated Glacier, Alaska. Science 227:4686469–479. Google Scholar
- B. Lefauconnier and J. O. Hagen . 1991. Surging and calving glaciers in east Svalbard. Norsk Polarinst. Meddelelser 116:130. Google Scholar
- M. F. Meier and A. Post . 1969. ?What are glacier surges. Canadian Journal of Earth Sciences 6:807–817. Google Scholar
- T. Murray and J. A. Dowdeswell . 1992. Water throughflow and the physical effects of deformation on sedimentary glacier beds. Journal of Geophysical Research 97:B68993–9002. Google Scholar
- T. Murray, G. W. Stuart, P. J. Miller, J. Woodward, A. M. Smith, P. R. Porter, and H. Jiskoot . 2000. Glacier surge propagation by thermal evolution at the bed. Journal of Geophysical Research 105:B613491–13507. Google Scholar
- T. Murray, T. Strozzi, A. Luckman, H. Jiskoot, and P. Christakos . 2003. ?Is there a single surge mechanism Contrasts in dynamics between glacier surges in Svalbard and other regions. Journal of Geophysical Research 108:B5 No.2237. Google Scholar
- W. S. B. Paterson 1994. The Physics of Glaciers,
*3rd edition*Oxford Pergamon Press. 480. Google Scholar - A. Post 1969. Distribution of surging glaciers in Western North America. Journal of Glaciology 8:53229–240. Google Scholar
- Gde Q. Robin and J. Weertman . 1973. Cyclic surging of glaciers. Journal of Glaciology 12:643–17. Google Scholar
- M. P. Searle 1991. Geology and Tectonics of the Karakoram Mountains New York Wiley. Google Scholar
- M. Sturm 1987. Observations on the distribution and characteristics of potholes on surging glaciers. Journal of Geophysical Research 92:B99015–9022. Google Scholar
- M. Truffer, W. D. Harrison, and K. A. Echelmeyer . 2000. Glacier motion dominated by processes deep in underlying till. Journal of Glaciology 46:153213–221. Google Scholar
- M. Truffer, K. A. Echelmeyer, and W. D. Harrison . 2001. Implications of till deformation on glacier dynamics. Journal of Glaciology 47:156123–134. Google Scholar
- R. L. Wessels, J. S. Kargel, and H. H. Kieffer . 2002. ASTER measurement of supraglacial lakes in the Mount Everest region of the Himalaya. Annals of Glaciology 34:399–408. Google Scholar
- J. G. Winther, S. Gerland, J. B. Ørbaek, B. Ivanov, A. Blanco, and J. Bioke . 1999. Spectral reflectance of melting snow in a high Arctic watershed on Svalbard: some implications for optical remote sensing studies. Hydrological Processes 13:–11122033–2049. Google Scholar