ON GENERATING DIGITAL ELEVATION MODELS FROM LIDAR DATA – RESOLUTION VERSUS ACCURACY AND TOPOGRAPHIC WETNESS INDEX INDICES IN NORTHERN PEATLANDS

Global change and GHG emission modelling are dependent on accurate wetness estimations for predictions of e.g. methane emissions. This study aims to quantify how the slope, drainage area and the TWI vary with the resolution of DEMs for a flat peatland area. Six DEMs with spatial resolutions from 0.5 to 90 m were interpolated with four different search radiuses. The relationship between accuracy of the DEM and the slope was tested. The LiDAR elevation data was divided into two data sets. The number of data points facilitated an evaluation dataset with data points not more than 10 mm away from the cell centre points in the interpolation dataset. The DEM was evaluated using a quantile-quantile test and the normalized median absolute deviation. It showed independence of the resolution when using the same search radius. The accuracy of the estimated elevation for different slopes was tested using the 0.5 meter DEM and it showed a higher deviation from evaluation data for steep areas. The slope estimations between resolutions showed differences with values that exceeded 50%. Drainage areas were tested for three resolutions, with coinciding evaluation points. The model ability to generate drainage area at each resolution was tested by pair wise comparison of three data subsets and showed differences of more than 50% in 25% of the evaluated points. The results show that consideration of DEM resolution is a necessity for the use of slope, drainage area and TWI data in large scale modelling.


Introduction
The most recent scientific assessment of climate change by the Intergovernmental Panel on Climate Change states that the world is becoming warmer, and that the extent of warming will be greater at higher latitudes (Denman et al. 2007). Northern peatlands store about 25% of the world's terrestrial carbon and emit between 10 and 20% of its natural methane sources (Gorham 1991(Gorham , 1995Moore et al. 1998). Estimating the distribution of wetness in northern peatlands is thus essential in order to determine their influence on the emission of greenhouse gases (GHG). Due to such changes in wetness, e.g. by permafrost thawing and a following change in biogeochemical function, changes in methane emissions have been reported to increase by factor of 10 to 50 (Bubier et al. 1999;Christensen et al. 2004).
At the macro scale there have been attempts to model continental and regional scale distribution of peatlands using Topographic Wetness Indices (TWIs) (Curmi et al. 1998;Gedney, Cox 2003;Kirkby et al. 1995). The success or failure of these approaches appears to be scale dependent. Wetlands within small and medium size catchments have been modelled successfully (Creed et al. 2003). However, for biogeochemical changes it is the intra-and inter-peatland complexity that is important (Baird et al. 2009). Also Baird et al. (2009) defined various scales of heterogeneity critical to peatland wetness: microtopography (1 m: S1 scale), mesotopography (10 m: S2 scale), general morphological forms (100 to 1000 m: S3 scale) which as adjacent and inter-connected systems form peatland complexes (S4 scale), and the catchment to continental scale presence or absence (S5 scale). For biogeochemistry the focus is the S1 to S3 scale, while the modelling focus to-date has been at the S5 scale.
A method has to be developed and evaluated that relates the elevation differences within peatlands (S1 to S3 scales) to characterize general wetness and the associated vegetation community structures that correlate with biogeochemical processes. Estimating the distributed wetness is one part of the effort to develop and explore innovative methods that can help us understand the relations between the physical attributes of peatland surface topography, permafrost, and hydrology, and the corresponding plant communities and their carbon storage and emissions. One approach is the use of high resolution DEMs and the spatial variability of a simulated TWI, which has been the basis of a suite of hydrological models (Quinn et al. 1995) and multiple flow algorithms (Pilesjö 2008;Pilesjö et al. 2006;Schmidt, Persson 2003).
A topographic wetness index (TWI) can be calculated using estimates of the slope and drainage area. Wetness is a geographical phenomenon that varies continuously over space. An important characteristic of peatlands is that the change in elevation between neighbouring points is relatively small, while the difference in wetness is relatively large. The use of high resolution data to estimate the wetness at many points over an area may thus provide a tool for the calculation of gas emissions in peatlands.
In order to initiate a study on the hydrology of peatlands, and at a later stage detect the important changes in wetness, it is necessary to create a digital elevation model (DEM) for a specific peatland area. In this study, high resolution LiDAR data (Fowler 2001) were used to generate a DEM for use in the estimation of wetness in a peatland area in northern Sweden.

Research questions and hypotheses
Based on the introduction presented above, the following main research question was formulated: 1. To what extent do the estimates of the slope and drainage area, and thus also the topographic wetness index, vary with the resolution of the digital elevation model for a peat land area? Our hypothesis is that estimates of the slope as well as drainage area are influenced by resolution, especially on relatively flat ground, like a peat land. Smaller terrain forms are supposed to be 'filtered' when a lower resolution is used, resulting in lower slopes. Regarding drainage area, the actual size of a cell set the minimum area, and this in combination with the fact that smaller drainage basins (depressions, often called sinks) are 'filtered' when a lower resolution is used, will result in estimations of relatively smaller drainage areas when using a higher resolution. This will together result in relatively low estimations of TWI when using high resolution data (i.e. small drainage areas divided by steep slopes) and relatively high estimations of TWI when using a DEM with a lower resolution (i.e. large drainage areas divided by less steep slopes).
When creating a DEM choice of interpolation algorithm, search radius, cell size etc. can be discussed (see e.g. Hengl 2006). These parameters, as well as the type of terrain, might influence the quality of the interpolated surface (see e.g. Hengl 2006;Kienzle 2004). This resulted in the following two, secondary, research questions: 2. To what extent do search radius and cell size influence the accuracy of a DEM created using a standard interpolation method for a peat land area?
3. Does the accuracy of the DEM change with the slope of the terrain and, if it does, to what extent? Based on studies by e.g. Burrough and McDonnel (1998); Haining (1990), we hypothesize that the accuracy of the DEM is sensitive to search radius, decreasing the quality of the DEM when the distance of significant spatial influence between points (i.e. spatial autocorrelation) is exceeded. However, cell size should not influence the accuracy since interpolation of a DEM consists of a number of point interpolations, where the cell value equals the interpolated value of the cell centre. This value is not dependent on cell size.
Regarding the relation between accuracy and slope no relationship is supposed to be found if the distribution of data point for the interpolation is homogeneous and the complexity is low. However, if the steeper terrain has a high frequency of 'hill tops' and 'valleys' , and the input data points are not evenly distributed around the point to be interpolated, we expect higher accuracy on flat terrain than in steeper areas.

Objectives and aims
The main objective of this study was to investigate possible variation in estimated slope and drainage area, and when combined also topographical wetness, using digital elevation models with different model resolutions for a northern peat land area. High resolution LiDAR data were used for the generation of the DEMs. The secondary objective was to investigate the role of search radius and cell size when using a standard interpolation algorithm in the generation of the DEM. The tertiary objective was to determine the accuracy of the DEM for different terrain with different slopes.
In order to achieve the objectives a number of specific aims were defined: 1. To create DEMs for the peat land area, using a standard interpolation algorithm with different search radiuses and spatial resolutions (cell sizes). 2. To evaluate the accuracy in the prediction of the elevation for the different DEMs using LiDAR data points as 'ground truth' . 3. To calculate and compare the accuracy in the DEM estimates of elevation for different slope intervals. 4. To study how estimated slope varies with DEM resolution. 5. To study how estimated drainage area varies with DEM resolution. 6. To conclude how estimated TWI might vary with different DEM resolution.

Elevation data and study site
This study is based on earth surface elevation data measured at the Stordalen mire and its catchment area. Stordalen is a peatland area in the Arctic region 10 km west of Abisko (68º 20' N, 19º 03' E) in northern Sweden. The hydrology and soil moisture conditions of the Stordalen mire have been reported previously . Apart from studies associated with the International Biological Program (IBP) (see e.g. Sonesson et al. (1980), the Abisko area has been included in many research programs, and its climatological records extend from 1913 to the present date (Andersson et al. 1996). This site is thus suitable for investigating methods of estimating changes in carbon storage, providing the possibility of validating tools for the prediction of changes in peatlands with past and future changes in permafrost. The mean annual temperature for the period 1913-2003 is reported to be -0.7 °C (Johansson et al. 2006) for Abisko. A regional rain shadow affects the precipitation and makes it among the lowest in Scandinavia with a mean annual precipitation of 304 mm for the period 1913-2003 (Johansson et al. 2006). An area of approximately 18 km 2 , containing the Stordalen mire, was selected as the study area in this project. An airborne LiDAR device has been used to measure the surface elevation. LiDAR is an acronym for 'Light Detection and Ranging' , and is a laser-based, remote sensing system used to collect various kinds of environmental data, including topographic data (Fowler 2001). Over the area defined above the total number of measured elevation data points (the raw data) is 76 940 341. This results in a high resolution data set with an average spatial distribution of approximately 13 points/m 2 . The accuracy of any LiDAR data points is related to the accuracy of the LiDAR device components (sensors). With recently developed LiDAR components, GPS and an Inertial Measurement Unit (IMU), range precision can reach 2-3 cm (Lemmens 2007). Airborne GPS accuracy is within 5 cm horizontally and 10 cm vertically, while the accuracy of IMU is less than a couple of centimetres. For LiDAR data in general, the root mean square error (RMSE) can get 15 cm vertically and 20 cm horizontally (BC-CARMS 2006).
The LiDAR data in the present study were retrieved with a TOPEYE S/N 425 system mounted on Helicopter SE-HJC. The altitude when sampling was 500 m. The LiDAR data have been post processed and adjusted against 54 known points connected to the national geodetic network. The mean vertical error after postprocessing corrections is +0.004 m and the average magnitude of errors is 0.022 m. The RMSE is 0.031 m and the standard deviation is 0.031 m.

General process
The LiDAR elevation data were used to generate a number of DEMs. Figure 1 illustrates the general process applied for the generation of the DEMs, as well as the evaluation of relationships between the accuracy of a DEM and different resolutions and search radiuses. The general methodology used to investigate the influence of slope and the possible relation between resolution and estimates of the slope and drainage area is also presented.

Selection of evaluation data points
The most common techniques used for the generation of evaluation points are the ʻleave one techniqueʼ within cross validation, the split-sample technique and the independent set of sample (Declercq 1996;Erdogan 2009;Smith et al. 2005). For our high density data we decided to use the split-sample technique. In this method, a part of the raw LiDAR data is omitted before performing the interpolation. Then the differences are calculated between the predicted and measured (previously omitted) values (Declercq 1996;Smith et al. 2005).
The criterion for selecting the evaluation points was that the distance between a cell centre in the DEM to be constructed and the selected point should be less than, or equal to, 10 mm. This enables us to validate the estimated elevations at the cell centres in the DEMs using data values that were measured at almost the same location (maximum 10 mm away from the point of interest).  A MATLAB (MathWorks 2008) program has been created to perform the selection process. All 76 940 341 data points were processed. The program calculates four distances from each data point to the nearest four cell centres in the DEM to be created. All data points less than or 10 mm from the nearest cell centre are then selected to be evaluation data points. Obviously, the calculation of the nearest centre points, and thus the selected evaluation data points, is dependent on the resolution (cell size) of the DEM to be created. The following six resolutions were used: 0.5, 1.0, 5.0, 10, 30 and 90 meters cell size. This implies dividing the raw data into 12 subsets, six for interpolation and six for evaluation. The number of points to evaluate the DEMs with the six different resolutions is presented in Table 1. The minimum number of evaluation points was set at 60 according e.g. to the American Society for Photogrammetry and Remote Sensing (ASPRS 2005). As can be seen from Table 1, the number of points less than or 10 mm from the nearest cell centres for the resolutions 30 m and 90 m is less than 60. To obtain more points for evaluation at these resolutions the distance from the cell centres was increased to 30 mm (30 m resolution) and 50 mm (90 m resolution). We are aware that this modification may constitute a weakness in the methodology, and thus decided to use both the original distance (10 mm) and the extended distances (30 and 50 mm) in the evaluation.

Generation of the digital elevation model
The spatial autocorrelation of the LiDAR data was tested by the generation and interpretation of semivariograms for three 225 m 2 , selected areas: one with steep terrain, one with an intermediate slope and one with flat terrain. The terrain was divided into these categories using a DEM with 0.5 m resolution. The data within the area were corrected for the influence of the major slope in the area, i.e. the general trend was removed by subtracting the average slope of the mountainside. Semivariograms of each area were then created and interpreted. None of the areas showed significant spatial influence on a point further away than 3.7 meters.
The use of different interpolation algorithms for DEM creation has been discussed by several authors (Anderson et al. 2005;Erdogan 2009;Lee 2003;Liu 2008;Myers 1994). Erdogan (2009) also investigates and reviews the role of interpolation parameters (like search radius) for the results. It is obvious that, different algorithms give different results, and that some techniques are more suitable depending on terrain and data.
In this study we have decided to use the inverse distance weighted (IDW) interpolation (Shepard 1968). This method is based on the assumption that an interpolated point is influenced more by nearby data points than points further away (Burrough, McDonnell 1998). It can be discussed if this is the most appropriate one when working in relatively flat areas, and having access to a large number of dense data points for the interpolation. Childs (2004) and Liu et al. (2007) pointed out that LiDAR data have high sampling density and even for complex terrain, the IDW approach is suitable for DEM generation, and justifies the choice. In addition, the aim of this study is not to compare different algorithms, but rather to test a commonly used one, and to look at the influence of primarily cell size on the result.
IDW is as we know by far the most common 'standard' interpolation algorithm, included in most geospatial processing software. It is likely to believe that it also in the future will be used by many researchers in the domain. Also, even if there are differences between different interpolation algorithms, we expect possible differences due to cell size in the result to be comparable. A limited number of data points (if a local algorithm is used) are used for the interpolation of each cell, and the relation between individual interpolated points (in this case cell centres) should not be heavily influenced by the interpolation algorithm.
Six different DEMs were created. The resolutions used were 0.5, 1, 5, 10, 30 and 90 meters. The value of the search radius when interpolating was varied between four different values (1, 2, 5 and 10 meters) for each resolution (cell size). The larger the search radius, the more data are obtained, but the risk of including outliers in the interpolation is increased. Too large a search radius will also include data that have no influence on the topography, as shown in the semivariogram analysis. The number of interpolations was thus 24 (six resolutions x four search radiuses).
A program was developed in MATLAB to control the interpolation process, including the problem of processing such a large number of elevation data points. In order to conduct the interpolation process for each cell it is necessary to search inside a relatively large database. Such a process can be very slow since it has a computational complexity of n 2 (in our case n = 76 940 3412 minus the number of evaluation data points), and may take several weeks. In order to speed up the process, a spatial index key called a Morton value (Orenstein, Manola 1988) was calculated for each data point to be interpolated. Adding the spatial index decreases the computational complexity to nlogn, making it possible to deal with a large number of data points within a reasonable period of time.

DEM evaluation
The purpose of evaluating DEMs with different resolutions is to detect possible differences, and identify the resolution that represents the evaluation data most accurately. In order to accomplish this, we calculated the deviations between the measured evaluation data points and the interpolated values for the overlapping cell centre for all combinations of resolution and search radiuses.
The technique most appropriate for measuring the accuracy of the DEM depends on the kind of error distribution. A normal distribution of the errors is rare in a DEM derived from data collected by LiDAR, due to e.g. filtering and interpolation errors (Höhle, J., Höhle, M. 2009). The quantile-quantile plot test (the q-q test) (Thode 2002) was used to establish the degree of deviation of the data from a normal distribution. The quantiles of the errors (Δh) are plotted against the theoretical quantiles of a normal distribution. If the data distribution is normal, the q-q plot should yield a straight line. Figure 2 shows the result of the q-q plot for the 0.5 m resolution DEM. It demonstrates that there is a strong deviation from a straight line, indicating that the data are not normally distributed at this resolution. This is also valid for all other tested resolutions. An accuracy estimation of a non-normal error distribution suggested by (Höhle, J., Höhle, M. 2009) was therefore used. This requires the calculation of four parameters, namely the median, the normalized median absolute deviation (NMAD), and two sample quantiles.
The distribution of the differences in elevation (Δh) and the absolute differences in elevation (|Δh|) were used as measures of the accuracy of the DEM. Absolute errors were used because we are interested in the magnitude of the errors, and not in their sign. Using absolute values also excludes the assumption that the distribution is symmetric.
The measurement accuracy is determined by calculating sample quantiles of the absolute differences (|Δh|) (Höhle, J., Höhle, M. 2009). The sample quantiles is the order of the sample [x(1),…, x(n)], where x(1) denotes the minimum and x(n) the maximum value in the dataset. For example, the 95% sample quantile of |Δh| means that 95% of the errors have a magnitude within the interval [0; Q |Δ| (95)]. In another way 5% of the dataset have an error larger than the 95% quantile of |Δh|. The 50% quantile is denoted as the median. The median of the error is a robust measure, which provides an estimate of a systematic shift in the DEM. Moreover, the median is less sensitive to outliers in a dataset.
The NMAD is used as a measure of the standard deviation, in comparison with the standard deviation is more resilient to outliers in the dataset (see Equation 1).
where Δh j denotes the individual errors j = 1,……, n and, m Δh denotes the median of the differences in elevation.

Accuracy of the DEM for terrain with different slopes
Many different authors have discussed alternative algorithms in estimation of slope from digital elevation models (See e.g. Grimaldi et al. 2007;Santini et al. 2009;Pilesjö et al. 2006;Skidmore 1989;Tang, Pilesjö 2011). As a result, several slope calculation algorithms employed on DEMs have been used in GIS (Geographical Information System) software (e.g. ARC/INFO and ER-DAS IMAGINE). Accounting for the importance of gradient/slope estimations in many applications, Tang and Pilesjö (2011) tested possible differences between the characteristics of eight frequently used algorithms and investigated how they behave in different terrain. They concluded that differences exist, but are in most cases not significant. Since the scope of this paper is not to test and compare slope algorithms, a standard slope estimation algorithm was chosen. At every point in a DEM the slope can be defined as a function of gradients in the X and Y direction: The key in slope estimation is the computation of the perpendicular gradients fx and fy. Different algorithms, using different techniques to calculate fx and fy yield the diversity in estimated slope. For a gridded DEM, the common approach when estimating fx and fy is by using a moving 3×3 window to derive the finite differential or local polynomial surface fit for the calculation (Florinsky 1998;Zhou, Liu 2004). In this study, we have used a polynomial surface approximation. A second-order trend surface (TS), (see e.g. Pilesjö et al. 1998), based on a least-square approximation, was applied on each 3×3 window as follows:

TS x y a a x a y a x
where i = 1, ... , 9 corresponds to the numbering of the centre cell and its eight neighbours; a 1 , .... , a 5 are the constants for the second-order trend surface; x i , y i are cell co-ordinates (mid points) in a local system. In order to estimate the gradient (slope) of the trend surface in the middle of the centre cell both the x co-ordinate and the y co-ordinate are set to zero: In order to study the accuracy of the DEM in relation to the slope of the terrain, the evaluation points were divided into six subsets, corresponding to six slope intervals. Slopes with gradients from 0 to 50 degrees were divided into five equal intervals, while the sixth interval consisted of slopes steeper than 50 degrees. A DEM with a resolution of 0.5 m created with a search radius of 1 m was used to estimate the slope. The 0.5 m resolution was chosen since this has the largest number of evaluation points (see Table 1). The evaluation points are divided into six equivalent datasets, and the accuracy of the DEM was calculated for all these datasets.

Relationship between the slope and drainage area for different DEM resolutions
A large number of different methods exist for estimation of drainage area, flow and flow accumulation (see e.g. Pilesjö 2008). Some of these estimate contributing area (up-slope) as well as dispersal area (down-slope). A few of the commonly used methods are briefly introduced below.
The single flow D8 algorithm was described by O'Callaghan and Mark (1984). It assumes that flow follows only the steepest downhill slope. Using a raster DEM, implementation of this method resulted in that hydrological flow at a point only follows one of the eight possible directions corresponding to the eight neighbouring grid cells (Band 1986;ESRI 1991;Mark 1984;O'Callaghan, Mark 1984). Here we call this approach a 'single flow' algorithm. However, for the quantitative measurement of the flow distribution, this over-simplified assumption must be considered as illogical and would obviously create significant artefacts in the results, as stated by e.g. (Freeman 1991;Holmgren 1994;Pilesjö, Zhou 1996;Wolock, McCabe 1995). More complex terrain is supposed to yield more complicated drainage patterns. The difference between the D8 algorithm and the commonly used Rho8 algorithm, presented by (Fairfield, Leymarie 1991), is that the Rho8 also includes a stochastic variable.
Attempts to solve the problem connected to the 'single flow' algorithms have led to several proposed 'multiple flow direction' algorithms (see e.g. Freeman 1991;Holmgren 1994;Pilesjö, Zhou 1996;Quinn et al. 1991;Zhou et al. 2011). These algorithms estimate the flow distribution values proportionally to the slope gradient, or risen slope gradient, in each direction. The DEMON algorithm was presented by (Costa-Cabral, Burges 1994). In order to eliminate the problem with the one-dimensional flow, present in the other algorithms, DEMON uses two-dimensional flow tubes in order to trace flow up-streams and down-streams.
In this study we have estimated drainage area using new triangular multiple flow distribution algorithm (see Pilesjö et al. 2012), partly based on the 'form-based' algorithm presented by Pilesjö et al. (1998). Given the limitations and problems of the algorithms presented above this 'multiple flow direction' approach, based on analysis of the form of individual 3×3-cell surfaces was proposed. It was assumed that flow diverge over convex surfaces and converge over concave surfaces. There is no absolute way to determine convexity and concavity of the centre cell in a 3×3-cell surface. Pilesjö et al. (2012) propose a facet based solution to do this. Around the midpoint of the centre cell in question, eight planar triangular facets are constructed, each with two corners in two adjacent cells. With the aid of these eight triangular facets, our current grid cell (centre cell) is divided to eight triangular facets. The slope direction (aspect) of each of these triangular facets can then be calculated. The area of each facet is equals to 1/8 of the cell size, directly proportional to the flow contributed by that facet. The flow on each and every triangular facet is routed towards other, neighbouring, facets or, if we have a sink, stays in the same triangular facet. If routed, depending on the aspect of a specific facet, the flow is added to a neighbour facet or split between two neighbouring facets. The sum of the flow in the eight facets covering a cell equals the estimated drainage area.
Three different DEM resolutions were used to investigate the relationship between the estimates of the slope and drainage area on the one hand, and the DEM resolution on the other. Resolutions of 10, 30 and 90 m were used, due to the overlapping evaluation centre cell points for these resolutions as illustrated in Figure 3. (i.e. the centre of a 90 meter cell denoted point P in Figure 3) is also the centre point of a 30 and 10 meter cell. The slope and drainage area evaluation cells selected from the 10 and 30 m resolution DEMs are thus the ones that have the same cell centre position as the location of the 90 m resolution evaluation cells. This results in three subsets of data that have the same number of evaluation points with the same locations, but contain slope and drainage area values estimated using different DEM resolutions. The results obtained from the three subsets are then compared to identify/reveal possible differences in the estimated slopes and drainage areas. The differences between pairs of different resolutions were calculated.

Evaluation of the DEM accuracy
Twenty four different DEMs, interpolated using six different resolutions (cell sizes) and four different interpolation search radiuses, for the same area, were to be evaluated. The aim was to determine differences between different DEMs, but also to investigate which combination of resolution and search radius gives the best accuracy. Using the robust accuracy measures appropriate for non-normal error distributions, we calculated the median, the NMAD and two quantiles (68.3% and 95%) for each combination of resolution and search radius.
The medians were all zero except for the median at 90 m resolution, which varied between -10 and -20 mm depending on the search radius. The results of the NMAD calculations are given in Table 2. The values indicate that the accuracy of the DEM is the same for different resolutions when using the same interpolation search radius. The accuracy is generally higher the shorter the interpolation search radius. The results for the two quantiles are presented in Table 3, and confirm the NMAD results. This is an expected result, as when increasing the interpolation search radius it can be expected to increase the errors in the created DEM. According to the measures of accuracy, the six most accurate combinations of resolution and search radius are those with the 1 m interpolation search radius. For these six cases, the maximum errors in the elevations are around 40 mm within the 68.3% quantile of the data (see Table 3). Moreover, the maximum errors in the DEM elevations are around 100 mm within the 95% quantile of the data.

Accuracy of the DEM for different slope intervals
The results of evaluating the relationships between the six different slope intervals and the errors represented by the NMAD are shown in Figure 4. When visually analysing the shape of the slope error curve it is obvious that there are larger errors in elevation when the terrain is steep than when it is flat. The first point in Figure 4A shows that for slopes between 0 and 9.99 degrees the error in elevation is around 0.03 m. Figure 4B illustrates two different quantiles of errors, also confirming that errors are larger in areas with steeper slopes.

Slope estimation using different DEM resolutions
Evaluation of the slopes estimated with the DEM at different resolutions shows that the medians of the differences between the slopes have negative signs, i.e. lower resolution (larger cell size) generates lower values of the slope. The frequencies of the differences in the slopes are shown in Figure 5, where the negative skewness of the distributions can also be seen. The NMAD and the quantile results also confirm the relationship between the values of the slope and the resolution. Results for the median, NMAD and the two quantiles are given in Table 4.  It should be noted that differences in estimated slope sometimes exceed 10 degrees. Also the relative differences, presented in Table 4 indicate significant differences between slope values estimated from DEMs with different resolutions. For example, the median difference in slope between the 90 and 30 meter DEM is -0.613 degrees, corresponding to a relative difference of 9%. Investigating differences for individual points (overlapping cells) we can note that more than 25% of the evaluation points differ more than 50% in estimated slope when comparing the 90 metre DEM and the 30 metre DEM. Lower resolution yields lower (less steep) slope estimations.

Estimation of the drainage area using different DEM resolutions
The medians of the differences between drainage areas estimated using different resolutions have positive signs, i.e. the lower the resolution (the larger the cell size) the higher the values of the drainage area. The frequencies of the differences in drainage area are illustrated in Figure 6, where the shift of the median towards positive values is clear. The NMAD and the quantile results also confirm the relationship between the values of the drainage areas and resolution. Results for the medians, NMAD and the two quantiles are given in Table 5. Also for the differences in estimated drainage area it should be noted that these sometimes exceed 50% (see Figure 6 and Table 5). The numbers presented in Table 5 indicate significant differences between drainage area values estimated from DEMs with different resolutions. For example, the median difference in drainage area between the 90 and 30 meter DEM is 21 480 m 2 metres, corresponding to a relative difference of 83%. Investigating differences for individual points (overlapping cells) we can note that more than 75% of the evaluation points differ more than 50% in estimated drainage area when comparing the 90 metre DEM and the 30 metre DEM. Lower resolution yields larger values in estimated drainage area.

Discussion
A number of choices influence the results of this study, some of which are discussed below. One limitation is that this is a case study, applied on one peatland area. Even if sample sizes are relatively high, it would of course strengthen the conclusions if other areas and in other regions could be included.
Regarding the interpolation algorithm, we chose the IDW algorithm, mainly because it is one of the most commonly used ones for interpolating scattered points. In order to determine whether there were large differences in the results when using other interpolation techniques or not, we created DEMs using the nearestneighbours interpolation and the bilinear interpolation algorithm. However, the results showed the same trend in errors as the IDW. This confirms that there are no significant differences in the results between IDW and other algorithms when using a high density LiDAR data, also reported by e.g. Chaplot et al. (2006);Podobnikar (2005). We also tested the possible spatial autocorrelation of the dataset in order to rule out the need for a geostatistical (Kriging) approach (Cressie 1993).
Regarding the evaluation data, ground truth field data are normally used for the evaluation of interpolated DEMs. In this study, we excluded certain data points from the raw data, and then used them for evaluation purposes, as suggested by e.g. Erdogan (2009). The excluded evaluation data were taken at a maximum distance of 10 mm from the centres of the cells, and were not used in the interpolation process. In order to increase the number of evaluation points we increased the maximum distances for the 30 and 90 m cells to 30 and 50 mm, respectively, giving 230 and 68 points, instead of 57 points and 6 points, respectively. This is a weakness of the methodology. However, the results obtained using the extended evaluation dataset show the same trend in errors as evaluation of the limited number, confirming that the use of the modified selection criteria did not affect the results significantly. The extended selection was justified for statistical reasons. Also the use of LiDAR data for the accuracy assessment can be discussed. In this study the relative errors between different DEMs, and not the absolute errors, were in focus. This justifies the use of the LiDAR data points as ground truth, even if the errors in elevation can be up to 10 cm (see e.g. Lemmens 2007). Since the same data points were used for the different DEMs, the quality of the points does not influence the relative comparisons.
When choosing the different resolutions of the DEM, we logically assumed that a high resolution should reflect reality better than a poor resolution. However, it is still interesting to create and test DEMs with different resolutions since, in most cases, DEMs with poor resolutions are more commonly available, and thus more frequently used. Evaluation of the accuracy of different resolutions showed approximately the same results regarding elevation. This was expected, since the same interpolation data set was used for all resolutions, and the evaluation points are all located close to, or very close to, the cell centres. The interpolation algorithm can then be expected to work equally well for all resolutions. However, the uncertainties between the evaluation points, in this case the centres of the cells; will be much higher at lower resolutions, where the distances between the centres of the cells are greater. This uncertainty, increasing with lower resolution, is influencing e.g. estimation of slope and drainage area.
Employing different search radiuses changes the number of data points included in the interpolation for each cell. The reason for varying the search radius in this study was that available DEMs with the same resolution are often created from different sources, having different numbers of known data points. We expected that the use of different search radiuses would have considerable effects on the results, such that the estimate of a point close to a large number of known data points would be better than including more data points further away from the point to be interpolated. This was shown by the autocorrelation study and the search radius study. The increase in search radius, from 1 to 2, 5 and 10 metres, did not give a better result since the spatial influence is limited to 3.7 meters. The reason for the 1 metre search radius to give the best results is explained by the high density of data points. Even if there was a spatial influence up to 3.7 metres, the number of closer points (<1 metre) was enough to yield the best estimates.
We also examined the influence of the slope of the terrain by dividing the data into different slope intervals (from flat to steep). We found that the errors in elevation were higher for steeper slopes, confirming the results reported by e.g. Erdogan (2009). The reason for this, which is logical, is that an accurate interpolation demands equal distribution of data points around the point to be interpolated, as well as a linear surface (if a linear interpolation algorithm like IDW is used). Even if we had equal distribution around data points, the second condition (linearity) was probably not fulfilled in our case. Worth to be noted is also that the evaluation points are located at a specified maximum distance from the cell centre, resulting in a greater deviation in steeper terrain. These slope-related differences in accuracy within DEMs are often not noted, but should not be neglected as they may be significant. The accuracy of a DEM should perhaps be associated with the slope of the terrain, especially in high accuracy modelling on a detailed scale. Further studies are needed to develop more accurate DEMs, or at least to document the weaknesses in DEMs when steep slopes are involved.
Based on the results presented in Figures 5 and 6, and Tables 4 and 5, there is a strong indication that the median of the differences in drainage area is positive, while the median of the differences in slope is negative. This indicates that high resolution DEMs will estimate lower values of the drainage area than low resolution DEMs. As mentioned in the hypothesis this is logical, since smaller drainage basins, often referred to as sinks or pits, are 'filtered' when resolution is decreased. This naturally results in larger areas. Also the resolution itself influences the estimation of drainage area. The cell size determines the minimum area, resulting in smaller drainage areas if a high resolution DEM is used. However, since there are very few 'one cell areas' the latter reason can probably almost be neglected. Our findings on slope are supported by Chang and Tsai (1991) who have reported similar results for low resolution data. Zhang and Montgomery (1994) tested the grid size impact on the TWI calculations and found that a higher resolution yields better results in a hydrological modelling. The use of data with the lower resolutions has not decreased over time. The scope of modelling larger systems, all the way to global system models, shows that this is still a highly significant issue in the matter of modelling moisture in wetness indices as well as in physical models.
Moreover, the slopes in high resolution DEMs seem to be overestimated compared to low resolution DEMs. Also this is logical, since smaller terrain forms, yielding larger slope estimations, will be 'filtered' when a lower resolution is used. These effects on estimated drainage area and slope, related to resolution, will be even more pronounced when calculating wetness indices, as these are normally based on the ratio between the slope and drainage area (see e.g. Sorenson et al. 2006). Thus, higher wetness indices will be predicted using low resolution DEMs, and relatively low values will be estimated when using high resolution DEMs. In the present study, when changing the resolution from e.g. 10 to 30 metres, we have many examples (see Figure 5 and 6 as well as Table 4 and 5) of an increase in estimated drainage area of 50% and a decrease of estimated slope of 50%, resulting in a change in estimated wetness (drainage area divided by slope) of a factor three. These differences, linked to the results presented by e.g. Bubier et al. (1999) and Christensen et al. (2004), reporting that wetness influences methane emissions in peatlands by a factor 10 to 50, strongly indicate the relevance of this study. However, more research, confirming the influence of wetness for emissions of greenhouse gases is needed.

Conclusions
Addressing the research questions presented in Section 2 we can conclude that 1. The estimates of slopes and drainage areas, and thus also the topographic wetness index, differ significantly with the resolution of the digital elevation model for a peat land area in northern Sweden. 2. The search radius, but not cell size, significantly influences the accuracy of a DEM created using a standard interpolation method (IDW) for the peat land area, and that. 3. The accuracy of the DEM differs significantly with the slope of the terrain. The first conclusion, indicating that slope values become lower and drainage areas values higher when the resolution decreases (i.e. cell size is increasing) is of critical relevance to the modelling of greenhouse gases and climate change because it alters the predicted patterns of ecosystem wetness (Zhuang et al. 2007). For example, we showed that a 10 to 30 m change in resolution of a DEM can triple the estimated topographical wetness index in certain areas. This in turn, will lead to an increase in the estimate of emissions of water dependent greenhouse gases such as methane (Bubier et al. 1999;Christensen et al. 2004).
The last two out of the three major conclusions presented above confirm and strengthen results reported by other authors (see Anderson et al. 2005;Smith et al. 2005;Podobnikar 2005;Erdogan 2009). These, as well as the first one relating to estimates of slope, drainage area and topographical wetness index, are highly relevant in hydrological as well as carbon modelling (see e.g. Walker, Willgoose 1999; Sorensen et al. 2006;Sommer et al. 2004).
Further studies using high precision field data are recommended in order to clarify the relationships between DEM resolution and estimates of topographical/ hydrological parameters such as wetness. It is obvious that high resolution field measurements, such as LiDAR data, have great potential in the development of DEMs suitable for relatively accurate estimates of hydrological processes in the landscape. In the future, one possibility is to use regional and/or global LiDAR datasets to construct DEMs with acceptable accuracy to model the wetness in peatland areas, which in turn will facilitate studies in global change.

Acknowledgments
The funding of the LIDAR survey and DEMs created from it came from a number of agencies through a partnership of researchers. We acknowledge the contributions of the Natural Sciences and Engineering Research Council of Canada, Discovery Grant to Nigel Roulet (McGill University); The Abisko scientific Research Station supported at the time by the Royal Swedish Academy of Sciences, KVA; Patrick Crill (Stockholm University) by the Swedish Research Council, VR; Torben R. Christensen (Lund University) by the Swedish Research Council, VR; Håkan Olsson (Swedish University of Agri-cultural Sciences) by the Swedish Environmental Protection Agency; and Andreas Persson's and Petter Pilesjö's research grants at the Lund University GIS Centre.