New Insights into Factors Controlling Drainage Density in Subarctic Landscapes

ABSTRACT

Fluvial processes play a crucial role in landscape development in periglacial regions. Geomorphological studies have yielded contrasting inferences about the importance of various environmental factors to drainage density. In this paper, 12 explanatory variables, based on relief, vegetation, soil, and rock type, were used to explain the controlling factors of drainage density in an area of 2880 km² in subarctic Finland. Analyses were performed at multiple spatial scales (1, 4, 16, and 64 km²) based on two novel statistical methods, namely generalized linear modeling (GLM) and hierarchical partitioning (HP). The results of the GLM and HP analyses were broadly concordant. Most of the variation in drainage density can be explained by soil and vegetation variables. Drainage density increased with proportion of rock and gravel soils and alpine vegetation, and decreased with peat cover. Variables based on topography and rock type explained only a small amount of the variation in drainage density. The results also confirm the profound influence of the spatial scale on the geomorphologic systems: factors influencing drainage density yield different results depending on the spatial scale at which drainage density is analyzed. The large data sets and the use of rigorous statistical techniques at four different scales add the confidence and generality for the present results.

Introduction

Drainage density (D) has long been recognized as one of the key characteristics of natural terrain. Horton (1932) defined drainage density as the total stream length per area. Mathematically it is expressed as: where L_T is the total length of streams within area A. Drainage density is a measure of how well a watershed is drained by stream channels. According to Goudie (1990), D affords one of the best available quantitative expression of drainage texture—an expression of the relative spacing of channels in fluvially dissected terrain (Gregory and Walling, 1973; Tucker and Bras, 1998). Moreover, D is also related to several other measures of landscape dissection, including valley density (Tucker and Bras, 2000) and channel-head source area (Montgomery and Dietrich, 1989).

The study of drainage density has been a central theme in geomorphology and hydrology and is of considerable conceptual and practical interest (Tarboton et al., 1992; Oguchi, 1997). Previous studies have shown that drainage density is related to climate (Chorley, 1957; Gregory and Gardiner, 1975), vegetation (Melton, 1958), soil and rock properties (Wilson, 1971; Kelson and Wells, 1989), permafrost (McNamara et al., 1999), time (Ruhe, 1952), hypsometric integral (Strahler, 1952), and relative relief (Schumm, 1956; Montgomery and Dietrich, 1989; Oguchi, 1997). The theoretical basis for these observations has been studied within the context of models of dynamic processes (Kirkby, 1987; Montgomery and Dietrich, 1989; Tucker et al., 2001).

Although the research on controlling factors of drainage density has produced several important analyses of the main determinants of drainage density, two troublesome features have not been fully addressed: (multi)collinearity between the explanatory variables and spatial scale. Ignoring these potentially confounding phenomena may result in biased models, which can provide misleading inferences about the determinants in the drainage density.

First, multicollinearity among explanatory variables may result in excluding more plausible variables from multivariate models if other intercorrelated variable(s) happen to explain better, in statistical terms, the variation in response variable (Buckland and Elston, 1993; Mac Nally, 2000). To overcome the problems caused by collinearity, researchers have traditionally used two means: (1) investigating the correlations among explanatory variables and excluding some of the most clearly intercorrelated variables using a stepwise model building procedure (Philippi, 1993), or (2) summarizing the variation in several environmental variables into synthetic composite variables using techniques such as principal components analysis (Quinn and Keough, 2002; Venables and Ripley, 2002). One recent approach to tackling collinearity problems is the method of hierarchical partitioning (HP) (Mac Nally, 1996; Mac Nally, 2004). HP provides estimates of the independent effects of each predictor variable separately, by considering all possible models in a multivariate regression setting. In the current study, I employed stepwise generalized linear modeling (GLM) and hierarchical partitioning (HP) to evaluate the relative importance of several potential factors for the drainage density.

Second, a perennial problem in geomorphological studies is identifying the appropriate scale at which to sample, because spatial patterns detected in any landscape are influenced by factors operating at different scales or at several organizational levels (Turner et al., 1989; Walsh et al., 1998). This has fundamental significance for the study of geomorphological systems, because the distribution patterns and processes that are unique to any range of scales will have unique causes and consequences. Recently, GIS-based approaches have been used on different scales to analyze and model geographical distribution patterns (Walsh et al., 1998). Understanding the factors that influence different scales, however, is limited. Knowledge is lacking on how the different spatial scales affect the understanding of geomorphologic systems. As a consequence, an increasing need exists to evaluate how the analysis and modeling results respond at different spatial scales. In this study, the effect of scale on the determinants of drainage density was evaluated using four different scales (size of the study cell): 1, 4, 16, and 64 km².

Although rivers in a periglacial environment flow only during a few summer months, fluvial processes play a crucial role in landscape development (Williams and Smith, 1989; French, 1996). To my knowledge, however, the effects of environmental variables on drainage density of the periglacial landscape have not been analyzed simultaneously using rigorous statistical methods. Thus, the effects and relative importance of different factors on the drainage density are unclear. Improved knowledge of the factors controlling drainage density in periglacial environments is essential for the development of a thorough understanding of the landscape processes and landforms in cold regions.

This paper illustrates a novel multivariate evaluation of the effects of relief, vegetation, soil, and rock type on the drainage density in a subarctic landscape in northern Finland. A total of 12 environmental variables at 4 different scales were used in a stepwise manner to construct generalized linear models to explain the factors that control drainage density. Using hierarchical partitioning (HP), I also segregated the independent effect of each predictor variable separately to provide further assessments of the potentially causal variables. Specifically, I asked: (1) Which environmental variables correlate most markedly with the drainage density? (2) How great are the independent contributions of the 12 predictors in explaining the drainage density? (3) Are the results of the GLM and HP in concordance? (4) What is the effect of scale on the results?

Periglacial Hydrology

Periglacial and temperate geomorphology and hydrology differ significantly in a number of ways. The various elements of the hydrologic cycle in periglacial regions are modulated in intensity, magnitude, and significance when compared to temperate regions (Williams and Smith, 1989). Winter snow storage is the single most important feature of the hydrological cycle in cold regions. As a consequence, stream flow is intensely seasonal or even ephemeral. Streams are frozen over in winter, and are often frozen solid. Only the largest rivers maintain a reduced flow beneath the ice. Flow regimes of rivers in periglacial environments are extremely highly peaked because of the springtime melting of snow. On smaller rivers, the spring flood may account for as much as 90% of the discharge (McCann et al., 1972). In summer, the runoff response is modulated by the thawing of the active layer.

In most temperate environments, saturation excess is the dominant model of storm runoff generation in the vicinity of channels and channel heads (Dunne, 1978). The main factors that are responsible for this are high-frequency, low-magnitude rainfall distribution, in combination with large values of infiltration capacity and soil hydraulic conductivity, because of abundant vegetation and soil fauna. Under periglacial conditions, vegetation and soil fauna are scarce, soil hydraulic conductivity is, therefore, more likely to be small, and infiltration capacity may become a limiting factor. Hydraulic conductivity and infiltration capacity is further decreased by the presence of permafrost or seasonally frozen soils. Consequently, overland flow can be the dominant runoff-generating process in periglacial environments (Bogaart et al., 2003).

Soil erodibility is, in general, a function of soil texture, vegetation cover, and soil organic matter. Under periglacial conditions, the vegetation cover is sparse, and soil organic matter is often low. Soils, therefore, are easily eroded. Additionally, other periglacial processes reinforce erodibility; e.g. freeze-thaw cycles weaken the soil and riverbank strength (Bogaart et al., 2003). As a consequence, though periglacial streams flow only during a few summer months, flowing water is still capable of great erosional and transporting activity (Williams and Smith, 1989).

Drainage network is not a static variable. During the course of individual storms or spring floods, networks expand and contract. Drainage networks also vary in extent on much longer timescales. The morphometry of valley systems may, therefore, be used as a means of paleohydrological retrodiction (Gardiner, 1983). For example, Bogaart et al. (2003) concluded that channel networks expanded under permafrost conditions and contracted under nonpermafrost conditions in Netherlands during the Pleistocene.

Study Region

The study region is located in the zone of discontinuous permafrost in the northern part of Finnish Lapland (Fig. 1). The cover of the study region is 2880 km². The study region is mostly open highland, where low hills covered with sparse vegetation and gently sloping fells alternate with mires with a shallow layer of peat (Fig. 2). Elevations in the area range from 65 to 641 m a.s.l. with a mean of ca. 330 m (Atlas of Finland, 1986; Seppälä, 1997b; Luoto and Hjort, 2004). The climate of the region is subarctic: the mean annual air temperature was −2.0°C (mean annual min −39.0°C and max +28.4°C) during the period 1962–1990 measured at the Kevo Meteorological Station (69°45′N, 27°01′E; 107 m a.s.l.), which is located 10 km east from the study region. Mean annual precipitation at Kevo was 395 mm in 1962–1990 (Climatological Statistics in Finland 1961–1990, 1991). Like in most subarctic environments, snowmelt in the study region is usually a very rapid process. At the Kevo Meteorological Station, almost 50% of the snowpack may disappear in only 10 days, and in 30 days, on average, 85% of the snowpack is melted away (Dankers, 2002). The river valleys become free of snow in May or occasionally in early June, but this can be several weeks later in the surrounding uplands (Seppälä, 1976).

Figure 1 The location of the study region in northern part of Finland (left), and the study design (right) based on four different spatial scales: 1 km², 4 km², 16 km², and 64 km².

Figure 2 A three-dimensional terrain view (left) and the drainage pattern and catchments (right) of the study region. Drainage networks of the study region were constructed from the digital water database of the National Land Survey of Finland measuring streams shown in blue in topographic maps. The terrain view was produced using the shaded relief surface model derived from a digital elevation model with 25-m resolution (sun angle  =  45°, azimuth  =  315°).

Botanically, the region lies to the north of the northern limit of the continuous Scots pine (Pinus sylvestris) forest in the orohemiarctic zone, with birch (Betula pubescens ssp. czerepanovii) as the dominating tree species (Hustich, 1960; Ahti et al., 1968). Geologically, the study region belongs to the Precambrian granulite complex (Atlas of Finland, 1986; Seppälä, 1997b). Bedrock is characterized by gneissic granite and garnet-bearing paragneiss (granulite) with some felsic metavolcanic rocks (Kaitanen, 1969). Superficial deposits contain podsolized glacial tills covered by younger biogenic deposits on fells (Luoto and Seppälä, 2000). Valleys and other depressions with silty sediments are the sites of extensive peat deposits. Mires, especially palsa mires, cover ca. 5% of the study region, mostly occurring between 360 m and 390 m a.s.l. (Luoto and Seppälä, 2002). In the valleys, glaciofluvial gravel and sand deposits are relatively abundant. Frost heaving and physical weathering of bedrock have created extensive boulder fields (Atlas of Finland, 1986; Seppälä, 1997b; Luoto and Hjort, 2004).

The study region is on the watershed of Tenojoki and Paatsjoki catchments, which drain into the Barents Sea. Hydrologically, the study region is a typical subarctic region dominated by a high discharge peak in spring. The short and extreme spring floods are the result of rapid snowmelt and the lack of significant lake basins that would attenuate flood peaks (Dankers, 2002). The total cover of the lakes and ponds is ca. 51 km² (1.8%). A number of smaller streams characterize the area, and almost half of them are temporary, i.e. they dry during the summer months.

Methods

Drainage Density

Drainage networks of the study region were constructed from the digital water database of the National Land Survey of Finland. The database was created using aerial photographs and field checks on the scale of 1∶10,000. Two main methods have been proposed to determine stream lengths on topographic maps: measuring streams shown in blue, and measuring drainage lines suggested by V-shaped contours. In mountainous landscape the “blue line” method has been criticized because it is inadequate to describe actual stream networks (Morisawa, 1957; Oguchi, 1997). On the other hand, the second method measures valley density rather than channel network density (Montgomery and Dietrich, 1988; Dietrich and Dunne, 1993). An important difference is that valley density is a morphological rather than a hydrological property per se (Tucker et al., 2001). Further, in landscapes with relatively low relief, the determination of drainage networks using the “V-shaped contour” method or automated channel identification using GIS techniques are problematic, because they are sensitive to small errors and noise of digital elevation model (see Tucker et al., 2001). This is particularly true in the study region, which is characterized by extensive peat plateaus and gentle sloping fells with rather small changes in relief. Consequently, I used the blue line method in this study.

Explanatory Variables

The terrain parameters were calculated by Arc/Info‵s GRID from a digital elevation model (DEM). The DEM is a regularly spaced matrix of altitude values with planimetric coordinates. The DEM with 25 m grid size was created by linear interpolation from contour lines (elevation isolines) from 1∶10,000-scale paper copies using TOPOGRID command from Arc/Info (Esri, 1991). Mean altitude, topographical roughness (standard deviation of the altitude), and relative altitude (difference of the highest and lowest point) were calculated directly from the DEM by ZONAL functions.

The proportions of four soil types, as a percentage cover for each study cell, were calculated from a digital soil map using an Arc/Info Grid. The following list describes the variables related to soil types used in the analysis and modeling: (1) peat (peat thickness more than ca. 20 cm), (2) glacigenic deposit (till), (3) sand and gravel (glaciofluvial deposits), and (4) rock terrain (surficial deposits less than ca. 50 cm thick). The proportions of two major bedrocks were calculated from a digital bedrock map: (1) igneous rock (gneissic granite), and (2) metamorphic rock (garnet-bearing paragneiss and felsic metavolcanic rocks). The proportions of three vegetation types were calculated from the CORINE2000 land cover data set. CORINE is a European land cover classification originally introduced by the European Union. Data for land cover are based on supervised interpretation of Landsat ETM images at a resolution of 30 m completed with digital databases (Härmä et al., 2004). The following list describes the vegetation type variables used in the analysis and modeling: (1) alpine (sparsely vegetated alpine heaths and bare ground), (2) bush (dominated by dwarf birch [Betula nana] and willows [Salix sp.]), and (3) forest (dominated by Scots pine [Pinus sylvestris] and fell birch [Betula pubescens ssp. czerepanovii]).

Generalized Linear Modeling

Generalized linear models are mathematical extensions of linear models (McCullagh and Nelder, 1989) that have become increasingly popular recently. Popularity results from the ability of GLM to handle non-linear relationships and different types of statistical distributions, and because it is technically closely related to widely used linear modeling and analysis of variance (ANOVA) (Guisan et al., 2002). Generalized linear models can handle distributions such as the Gaussian, Poisson, Binomial, or Gamma, with respective link functions set e.g. to identity, logarithm, logit, and inverse (Guisan and Zimmermann, 2000). They are, thus, more flexible and better suited for analyzing many different types of distributions that are poorly represented by classical Gaussian distributions (McCullagh and Nelder, 1989; Crawley, 1993). In this paper, I build GLM models using a backward stepwise approach, in which explanatory variables are excluded simply according to the statistical significance (p < 0.05). Because overdispersion was observed in the analyses of data for drainage density at all spatial scales, a quasi-Poisson error distribution and an F ratio test were used in the GLM models (McCullagh and Nelder, 1989; Venables and Ripley, 2002). For each predictive model, a final model was fitted with the significant terms only. Further, I examined the possibility of curvilinear relationships between explanatory and dependent variables by including the quadratic terms of the variables in the models.

Hierarchical Partitioning

In hierarchical partitioning, all possible models for the drainage density were considered in a hierarchy. This process involved computation of the increase in the fit (measured in this case as R ²) of all models with a particular variable compared to the equivalent model without that variable (Mac Nally, 1996, 2002; Quinn and Keough, 2002). For example, for the variable “mean altitude,” improvements in model fit when going from the null model f(Ø) to f(mean altitude), from f(relative altitude) to f(relative altitude + mean altitude)…etc., using all combinations of predictor variables to which “mean altitude” can be added, are averaged to yield the independent explanatory power of this variable. As a result, hierarchical partitioning provides, for each explanatory variable separately, an estimate of the independent contribution with all other variables. Hierarchical partitioning was conducted using the “hier.part package” version 0.5-1 (Mac Nally, 2004), which was run as a part of R statistical package.

Results

Generalized Linear Modeling

The drainage density varied from 0 to 7618 m km⁻² on the scale of 1 km², and from 612 to 2090 m km⁻² on the scale of 64 km², with a mean of 1308 m km⁻² at all studied scales. Separate testing of the relationship between the drainage density and all explanatory variables based on GLM models showed that only 4 variables of the 12 explanatory variables were significantly (p < 0.05) related to the drainage density at all spatial scales, namely mean altitude, peat cover, forest cover, and alpine cover, showing often a quadratic relationship between D (Table 1). Four variables were non-significantly related to D at all spatial scales: relative altitude, till cover, igneous rock cover, and metamorphic rock cover.

Table 1 Summary of the univariate generalized linear modeling (GLM) analyses for the drainage density on the scales of 1, 4, 16, and 64 km². Only those linear and quadratic functions of the variables that account for a statistically significant (p < 0.05) change in the deviance are shown, except linear non-significant terms are shown which have a significant quadratic function. Moreover, only those quadratic functions where the change from linear to quadratic model is significant are listed. Direction of the effect is presented by + and − symbols. The first symbol indicates linear term and second the quadratic term (***  =  p < 0.001, **  =  p < 0.01, *  =  p < 0.05, ns  =  p ≥ 0.05).

The set of variables which together best explained the drainage density in the multivariate GLM models indicated that D generally increased with covers of rock and gravel, and decreased with cover of peat (Table 2, Fig. 3). Some of the variables showed a different direction of the effect at different spatial scales, namely mean altitude and cover of forest. They showed a positive quadratic response at the two finest scales, whereas on the scale of 16 km² mean altitude had a linear negative response and on the scale of 64 km² a positive response, and forest cover a negative quadratic and linear responses on the two largest scales. The GLM models explained 3.7–46.8% (total deviance change) of the variation in the drainage density. A clear positive correlation exists between the deviance change of the GLM models and the spatial scale.

Table 2 Variables in the final GLM models on the scales of 1, 4, 16, and 64 km² for the drainage density, showing estimated parameters and approximate standard error values. GLM models were built using a backward stepwise approach starting with a full model (all variables included), in which the explanatory variables are excluded simply according to the statistical significance (p < 0.05).

Figure 3 The relationship between drainage density (m km⁻²) and (A) cover of peat soils (R ²  =  0.27), (B) cover of alpine vegetation (R ²  =  0.36), (C) mean altitude (R ²  =  0.15), and (D) igneous rock cover (R ²  =  0.00). The solid line is a second order polynomial fit. The data is based on 45 grid squares of 64 km² in size in subarctic Finland.

Explanatory power (proportion of explained deviance) of different variable groups, considered here, showed clear differences (Table 3). In general, soil and vegetation were the single most powerful explanatory variable groups for drainage density, e.g. they explained 23.0 and 33.8% of the variation in the drainage density on the scale of 64 km². In general, topography and rock type were the poorest variable groups explaining the variation in drainage density.

Table 3 Explanatory power (proportion of explained deviance) of different variable groups individually and all variables included in the models for the drainage density on the scales of 1, 4, 16, and 64 km². Number of statistically significant (p < 0.05) variables in the models is shown in the brackets.

Hierarchical Partitioning

The results of the hierarchical partitioning were largely in accordance with those of GLM models (Fig. 4). Five explanatory variables had a substantially greater independent explanatory power than the other predictors included in the analysis for the drainage density. On the two finest scales, covers of forest and shrub, and cover of gravel soils had the highest independent contribution among the explanatory variables. On the two largest scales, variables with the highest independent contributions were covers of forest and alpine vegetation and cover of peat soils. Cover of alpine vegetation had the highest independent contribution among the explanatory variables on the two largest scales, whereas forest cover has the highest contribution on the scale of 4 km² and gravel cover at the finest scale. Variables of topography and rock type had obviously lower independent contributions than other variable groups in HP analyses. The two variable groups with highest explanatory power in GLM models also appeared important in the HP results: variables for soil and vegetation. As a contrast to the GLM results, soil variables, particularly gravel cover, had in HP somewhat higher independent contribution in relation to vegetation variables at the two finest scales.

Figure 4 The independent contributions (given as the percentage of the total explained variance) of the environmental variables for the drainage density at four different spatial scales, as estimated from hierarchical partitioning.

Discussion

Drainage Density

Geomorphological studies have yielded contrasting inferences about the importance of various environmental factors to drainage density. In this study, the results of both the generalized linear modeling and hierarchical partitioning analysis were broadly concordant, with the variables highlighted as significant by the two approaches being the same. Most of the variation in the density of the stream network can be explained by soil and vegetation variables in subarctic Finland. In other words, major soil and vegetation types and stream network density have a rather intercorrelated spatial structure. Grid squares with the highest drainage density had a relatively dispersed distribution pattern in the study region. They were characterized by high covers of gravel and rock soils, and alpine vegetation.

In the study region, numerous parallel channels on the slopes of the fells indicate that other hydrogeomorphological conditions existed during the past. Glaciofluvial action was very intensive during the deglaciation. Stagnant ice produced a huge amount of meltwater that grooved a number of subglacial, marginal, and extramarginal channels into the moraine cover and even into the surface of bedrock (Seppälä, 1971). Marginal and submarginal glaciofluvial channels are predominant on the fell slopes. At the higher altitudes, between the fell summits, occur many overflow channels that may continue as subglacial channels into the valleys. Lateral meltwater channels, incised by glacial meltwaters flowing between the receding ice and the slope, have created a unique drainage network on the slopes of the fells characterized by a stream network running obliquely to the contour lines. Those areas have high drainage density, often more than 5000 m of streams per square kilometer. At the lower altitudes, glacial meltwater channels form fan-like systems at the slopes of the valleys. Most of the glaciofluvial erosion channels, developed in the last deglaciation, control the characteristics of the present-day drainage system.

The drainage density was particularly low in the southern part of the study region with relatively low relief characterized by extensive peat lands and forests, typically less than 1000 m km⁻². The low density of the stream network in mire areas is most likely controlled by the high infiltration capacity and hydraulic conductivity of peat soils (see Williams and Smith, 1989). In contrast, at the high altitudes the soils are often totally or partially frozen during the spring flooding, which increases the overflow and erosion of the uppermost already-thawed layers of soil (see Fig. 5) (Bogaart et al., 2003).

Figure 5 Áhkojohka river valley in the northern part of the study region (69°36′N, 26°18′E, 400 m a.s.l.). The water level has already lowered close to the mid-summer level. Photograph taken by Jan Hjort on 2 July 2002.

Catchment morphology and drainage density are strongly influenced by hillslope processes. For example, Oguchi (1997) reported that drainage density correlates negatively with relative relief in Japanese mountains, whereas Schumm (1956) found a positive correlation between drainage density and the relief ratio. In general, the relationship between drainage density and relief is positive in semiarid, low-relief landscapes, and negative in humid landscapes (Tucker and Bras, 1998). In northern Finland, mean altitude showed a positive relationship with the D, whereas neither relative altitude nor topographical roughness had a clear statistical relationships with the D at any studied scales. Moreover, topographical variables group was the second poorest variable group in the explanatory power.

The geologic control of drainage density is widely understood: drainage density tends to be large in areas underlain by rocks with low infiltration capacity or transmissibility (Gregory and Walling, 1973). In this study, the importance of rock type was negligible and statistically non-significant in all analyses. Igneous and metamorphic rocks of the study region are relatively similar in rock structure, weathering resistance, and infiltration capacity (Kaitanen, 1969; Atlas of Finland, 1990).

The overall message from the results of both statistical methods was similar: the unique effects of topographical and bedrock variables turn out to be rather small in subarctic landscapes of northern Finland. It is plausible, therefore, to assume that drainage density is not related to relief or bedrock in the study region, but depends more on the distribution patterns of soil and vegetation. Further, it is also reasonable to avoid modeling drainage density in terms of single causes. In particular, exploring only relationships between D and relief or bedrock, while ignoring the direct links between other potential explanatory variables and D, may give biased results (see Borcard et al., 1992).

Methodological Aspects

Statistical analysis of drainage density has focused typically on studies of drainage basins or river systems with different sizes and shapes (Oguchi, 1997). A number of studies have related variation in drainage density to different geographical factors, e.g. relief, climate, bedrock, vegetation, and soil. Collinearity between explanatory variables and the problems with spatial scale, however, can hamper detection of key variables that control drainage density. This study utilized two alternative multivariate methods at four different scales to address these difficulties in explaining variation in drainage density, namely generalized linear modeling and hierarchical partitioning.

Generalized linear modeling (McCullagh and Nelder, 1989) and hierarchical partitioning (Mac Nally, 2000) offer appealing alternatives to traditional least-squares-regression techniques when assessing drainage density–environment relationships, particularly if used in a complementary manner. The results of the two modeling approaches applied here, GLM and HP, were broadly concordant. The results of the two approaches, however, differed slightly for some variables; e.g. the relative importance of the cover of gravel was higher in the HP results than in the GLM results. This discrepancy may result from curvilinear relationships that were effectively captured in GLM analyses by including polynomial functions of variables into models, but are more problematic to take into account in hierarchical partitioning (Heikkinen et al., 2004). If the relationships between response and predictor variables are clearly monotonic, hierarchical partitioning provides a powerful method for identifying the variables with the largest independent effects, whereas the importance of polynomial variables cannot be assessed by this method (Chevan and Sutherland, 1991). This is a potential shortcoming of the method, because curvilinear responses, including humped ones, are common in geomorphological systems (Luoto and Hjort, 2005). Under such circumstances hierarchical partitioning is a less powerful method than GLM. On the other hand, a possible weak point in GLM and other regression settings is the potential of collinearity problems during the stepwise selection procedure of predictor variables (Mac Nally, 2000). In such cases, hierarchical partitioning can be used in confirming whether the predictors selected in regression models are among the most likely causal variables (Gibson et al., 2004, Heikkinen et al., 2004).

One important step toward a more general understanding of drainage density might be to emphasize the spatial scale (size of the analysis grid cell) at which D is calculated. Recently, progress in the understanding of geomorphological systems has been achieved by the awareness that patterns and processes are strongly scale dependent (Turner et al., 1989; Walsh et al., 1998). In the present study, the same data sets were analyzed at four spatial scales, varying the area over which D was calculated by only less than two order of magnitude (1–64 km²) and keeping the geographical extension of the study constant. Nevertheless, the relative importance of the independent variables based on HP and the modeling performance of GLM models changed considerably. These results confirm the profound influence of the spatial scale on the results: factors influencing D yield partially different results depending on the spatial scale at which D is calculated and analyzed. The results showed that the density of the stream network is usually influenced by the combined effect of several variables, each operating in a different way on different spatial scales. Hence, it is useful to perform geomorphological and hydrological analyses at different spatial resolutions and extents (Walsh et al., 1998). In this way one can determine the scale at which factors influence geomorphological systems and which factors have scale-invariant and which factors have scale-variant effects.

A clear positive correlation occurred between the size of the study grid and explained variation of drainage density. Model of the finest scale explained less than 5% of the variation of the data, whereas model based on the data calculated at the largest scale explained 47% of the variation in the drainage density, which is approximately double the explanatory power of the second largest scale. This means that the 64-km² grid squares used in this study reveal the variation in cover of the vegetation types rather accurately, and also capture many aspects of the soil and vegetation variation. For example, some of the 64-km² grid squares were totally covered by alpine vegetation and provided opportunities to examine the relationship between variables of vegetation cover and drainage density.

In this study, climatic variables were not included in the analyses, because climatic influences are difficult to assess in mountainous or fell areas. Complex topography produces steep gradients in the local climate, and climate-monitoring stations are very sparsely distributed. In addition, the use of interpolated temperature data, which is normally calculated using the temperature gradient of free atmosphere, is problematic because of the strong winter inversion over the region (Virtanen et al., 1998; Luoto and Seppälä, 2002). In lieu of accurate climatic data at appropriately fine spatial scales, topographical variables are commonly used as surrogate predictors (Franklin, 1995).

The occurrence of permafrost has a profound effect on the generation of runoff and morphometry of the channel network because the presence of permafrost makes the subsoil impermeable. Therefore, overland flow can be the dominant runoff-generating process during periglacial conditions (Bogaart et al., 2003). In general, streamflow response to snowmelt or rainfall is rapid, and the pattern and density of drainage is controlled by spatial and temporal variations in the active layer (Williams and Smith, 1989). Unfortunately, accurate, fine-scale information on permafrost distribution was not available in the study region, thus permafrost occurrence was not used as an explanatory variable in the analyses.

In general, the study region is characterized by two main types of permafrost: palsa hummocks in the mires and permafrost in the bedrock on the barren fell summits (Seppälä, 1997a). Overall, palsa mires are characterized by rather low drainage density, whereas on the fell summits streams are often abundant. At the local scale, palsa hummocks act as a barrier to water movement in the subarctic mires. Thus the local configuration of permafrost in mires has a significant influence on the pattern of the drainage. On the fells, the extent of frozen ground in a river basin is one of the most important factors determining the base flow of streams. Consequently, groundwater flow becomes progressively more important as one moves downhill and permafrost become less continuous (Williams and Smith, 1989). By and large, the patchy occurrence of permafrost in subarctic Finland has a profound effect on the local configuration of the drainage network. Moreover, a large number of other hydrological processes are active in the periglacial regions, such as interhummock flow, icings, and slush flows (Seppälä, 1997a).

Although the regression models uncovered statistically significant relationships between the explanatory variables and the density of stream network, the goodness was relatively modest. More than half of the variation in the density remained unexplained at all spatial scales. The large data sets, large-scale of conducted sampling, and the use of rigorous statistical techniques at four different scales, however, add the confidence and generality for the present results.

Conclusion

Many geomorphological systems are heterogeneous at scales ranging from microstructures to continents. Drainage pattern is no exception in this regard, and its spatial pattern is controlled by different factors across multiple spatial scales. In general, most of the variation in drainage density can be explained by soil and vegetation variables in the subarctic landscape of northern Finland. This result highlights the importance of soil erodibility in drainage density, whereas the effect of relief and bedrock geology was relatively weak.

Acknowledgments

I am grateful to J. D. Vitek and G. Tucker for comments and constructive criticism made on the manuscript. Kevo Subarctic Research Institute provided high-quality facilities during the field surveys.