The suitability of using ASTER GDEM2 for terrain-based extraction of stream channel networks in a lowland Arctic permafrost catchment

Seasonally inundated areas and water-saturated soils are common features of lowland Arctic and sub-Arctic permafrost environments. With the onset of snow melt, and water percolation down through the snowpack, a principal factor controlling stream channel flow, aside from active layer depth, is topography. This paper investigates stream channel networks derived from the advanced spaceborne thermal emission and reflection radiometer (ASTER) global digital elevation model (GDEM) version 2 in a static terrain-based GIS-model. The suitability of using the ASTER GDEM2 for modelling the drainage network over a low-relief terrain is assessed. The aim is to use GDEM2 for the analysis of the stream channel network and to establish the network’s connectivity to previously observed spring flood patterns over the Yamal peninsula. As such, there are two parts to this study: 1) DEM validation and 2) stream channel network analysis. The results of the DEM validation study show that the root mean square error (RMSE) of the GDEM2 and reference data is approx. 10 m when compared to both reference data sets (RMSE = 12.17 m, N = 86 and RMSE = 9.64, N = 506,877), implying that the GDEM2 is sufficiently accurate for terrain-based modelling. The low connectivity between the stream channel network and seasonal inundation suggests that topographic controls play a less important role compared to the possible overbanking of lakes and basin overflow. However, drainage densities for investigated drainage basins were significantly lower than those expected from typical Arctic basins. Both more sophisticated modelling techniques as well as higher spatial resolution DEMs are needed to extract the stream channel network more accurately and hence establish a more comprehensive link between the drainage network and seasonally inundated areas.


Introduction
High latitude regions contain a combination of different water storages. Water is stored, in its liquid and solid state, in the ground in permafrost and the active layer and on the surface in lakes, the snowpack and glaciers (Barry & Gan 2011). Although not strictly surface hydrology, soil moisture, par-ticularly near-surface soil moisture (SSM) detectable with microwave remote sensing techniques (Wagner et al. 1999), is considered a crucial component of land surface hydrology due to its importance as a control parameter for ponding, runoff and infiltration (Wagner et al. 2007).
Hydrological processes at high latitudes are extremely active over a short time span in the spring The suitability of using ASTER GDEM2 for terrain-based extraction and summer season. In its liquid state, water is found in surface water bodies as well as in soils. Both surface water body extent and the degree of water saturation in soils undergo rapid changes from the onset of snowmelt until the autumn freeze-up (Vincent & Laybourn-Parry 2008). In order to establish comprehensive water balance studies that also account for rapid changes after snowmelt (Arp et al. 2011), seasonal inundation needs to be monitored. Capturing the inter-annual variability of seasonal inundation is particularly important, as it is not determined by hydrometric measurements (Zakharova et al. 2011).
To investigate the inter-annual variability of surface water extent on the landscape scale satellite data are required. In a remote sensing context spring floods may, however, result in false negatives/positives -the spurious disappearance or appearance of surface water bodies -if seasonal dynamics are not taken into consideration in change detection analyses of inter-annual water body change (Trofaier et al. 2013). In general, floodmonitoring methods have been developed for active microwave satellite sensors that can provide a means for frequently monitoring changes in spatial extent (e.g. Schumann et al. 2009;Trofaier et al. 2012;Trofaier et al. 2013). These monitoring efforts may be strengthened by an investigation of surface water flow through modelling stream channel networks. This procedure requires data on topography to describe the water pathways through the landscape (Tarboton 1997). Topography is represented by elevation data that nowadays are made available in digital formats. Regularly spaced digital elevation data in raster format are known as Digital Elevation Models (DEMs). DEMs are either digitised from compiled topographic maps (e.g. GTOPO30 -further information can be found at http: //webgis.wr.usgs.gov/globalgis/gtopo30/ gtopo30.htm) or derived from satellite sensors (e.g. SRTM and ASTER) (Rees 2012). Less common, due to the high costs (financial, logistical and time intensive) that are involved in such research, are ground (e.g. dGPS) or airborne (e.g. LiDAR) surveys with higher spatial and vertical resolutions (Hengl & Reuter 2009).
In the Arctic tundra biome, periglacial processes lead to specific landscape features associated with permafrost. In particular, surface water flow in a permafrost landscape leads to thermo-erosion and the development of gully systems (McNamara et al. 1999). In addition, thermokarst (ground subsidence after ground-ice ablation) leads to ponding and lake formation since thaw waters accumulate in thermokarst depressions (Grosse et al. 2013). Hydrological modelling may be used to identify surface water flow into these depressions as well as river systems. Morgenstern et al. (2013) investigated the modern characteristics of a drained thaw lake basin (DTLB) in the North of East Siberia by analysing morphological metrics and surface properties established from a high resolution DEM that they derived from stereo image pairs taken by the Advanced Land Observing Satellite (ALOS) Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM) sensor. The authors established different morphostratigraphic levels within the DTLB that they argue are due to climate-induced landscape evolution. These levels include the local micro-relief (ice-wedge polygons and pingos), which also plays an important role for hydrological surface flow. Liljedahl et al. (2012) examined the sensitivity of a hydrological model to topography by assimilating three artificially constructed DEMs, each representing unique tundra landscape features (low and high-centred polygons, as well as DTLBs). They found that the water balance in a polygonal permafrost landscape has significant effects on the partitioning between infiltration and runoff and that indeed permafrost micro-topography may have larger implications on the water balance that are not considered at regional scales.
Essential to all of these approaches is the accuracy of the elevation data, from which slope and curvature parameters may be deduced and, in the case of hydrological surface routing, stream channel networks may be inferred (Li & Wong 2010). Although it has become more common to produce 'home-made' DEMs using for example photogrammetric methods either from high-resolution stereo satellite imagery (e.g. Morgenstern et al. 2013) or from repeat aerial photography (e.g. Ryan, et al. 2014), open-community DEMs such as the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global DEM Version 2 (GDEM2) and the European Space Agency (ESA) Data User Element (DUE) Permafrost project DEM (D-DEM) (Santoro & Strozzi 2012) are still useful resources for geomorphological and hydrological studies dependant on relief.
Rees (2012) assessed GDEM2's predecessor, the ASTER GDEM, for applications in the Arctic by comparing it to reference data made available through airborne LiDAR campaigns and digitisation of topographic maps. The author found that the GDEM may be used tentatively using the stack-Anna Maria Trofaier and W Gareth Rees ing numbers -the number of stereo pairs -as a quality control. This present paper evaluates the applicability of the ASTER GDEM2 for hydrologic surface routing at high latitudes. The validation study takes advantage of a Soviet topographic map, as well as the aforementioned D-DEM, which are both used as reference data. This study sets out to determine whether the quality and vertical resolution of the chosen GDEM2 are good enough to continue with further topographic analysis, in the form of static terrain-based stream channel network extraction.

Regional setting
The Yamal peninsula is a flat, hummocky terrain located in the North of West Siberia in the continuous permafrost zone (Fig. 1). It is a gently rolling tableland that is divided by a number of alluviallacustrine-marine plains (Walker et al. 2011). Terraces strewn with lakes at different elevations make this permafrost area a unique but complex hydrological system.
The lower terraces are thought to be alluvial, originating from the Mordy-Yakha and Se-Yakha rivers. The III rd terrace is of alluvial marine/lacus-trine origin and is at an elevation of up to 26 m a.s.l. The IV th and V th terraces are of coastal-marine origin and lie at about 40-45 m and 58 m a.s.l., respectively. The mean elevation of the study area is approximately 16 m. A Russian research camp, known as Vaskiny Dachi, is situated in the catchment of both the Se-Yakha and Mordy-Yakha rivers (~70° 17' N, 68° 54' E). Active layer measurements and permafrost borehole temperatures at this site are continuous and records date back to the 1990s (Leibman et al. 2015).
The Vaskiny Dachi area is a landscape dissected by river valleys and their tributaries, strewn with lakes on both uplands and in DTLBs (sometimes single but often coalesced). In this paper, we refer to both drained and partially drained depressions as DTLBs. These DTLBs are known to flood and rapidly drain in the spring runoff season (Trofaier et al. 2013). Topographic controls may play a considerable role in runoff, spring flood and drainage mechanisms.

Previous research on seasonal inundation
The European Space Agency's Envisat Advanced Synthetic Aperture Radar (ASAR) instrument operating in Wide Swath Mode (WSM) was used to periodically monitor surface inundation in the Vaskiny Dachi area over the years , 2008(Trofaier et al. 2012. ASAR is a C-band in-  Brown et al. (1997). Russian administrative regions are outlined in grey (Stolbovoi & McCallum 2002) and the Yamalo-Nenets Autonomous District is highlighted by a line-fill. The suitability of using ASTER GDEM2 for terrain-based extraction strument (centre frequency 5.331 GHz); in WSM it has a spatial resolution of 120 m (Closa et al. 2003) and at high latitudes its temporal resolution is 2-3 days (using acquisitions from varying incidence angles) (Reschke et al. 2012). Through density-slicing techniques, surface water bodies were clearly distinguishable, enabling the dynamics of seasonal inundation to be monitored. This research led to the identification of recurring spring flood patterns and was the incentive for investigating topographic influences on runoff in this area.

Data and pre-processing
This study takes advantage of the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model Version 2 (GDEM2), which is available and downloadable free of charge in seamless 1° x 1° tiles via the United States Geological Survey's Earth Explorer (http:// earthexplorer.usgs.gov/). ASTER GDEM2 was released in October 2011 and retains its predecessor's nominal spatial resolution of 90 m. GDEM2 has an absolute vertical accuracy of 17 m at a 95% confidence level (ASTER GDEM Validation Team 2011). An additional 260,000 scenes and improved water body masking have strengthened GDEM2 with respect to version 1 globally; however, the AS-TER GDEM Validation Team (2011) also noted that at high latitudes stereo coverage is decreased compared to the GDEM. In addition, the GDEM2 water mask is derived from the Shuttle Radar Topography Mission (SRTM) which does not cover latitudes above 60° N (ASTER GDEM Validation Team 2011). The poor stereo coverage in association with the issue of the SRTM water mask footprint resulted in spike-anomalies in the GDEM2. These spike-anomalies are artefacts of erroneously increased elevation compared to their low surroundings and are found to be associated with areas of low stereo coverage over water bodies (Fig. 2). Quality assessment is done by inspection of ASTER stacking numbers (Rees 2012). For the Vaskiny Dachi area stacking numbers lay between 1 and 12. Masking pixels with stacking numbers ≤4 eliminated spike-anomalies over water bodies.
In addition to the ASTER GDEM2, this study also made use of a Soviet topographic map (publication date 1987, scale 1:200,000) of this area. This map was acquired from the Cambridge University Library, scanned and georeferenced in ArcMap, transformed from the original datum of Pulkovo 1942 to the more widely used WGS84. Furthermore, a second DEM was used to validate the GDEM2. This DEM was compiled under the ESA Data User Element (DUE) Permafrost project from digitised Russian Topographic Maps at 1:200,000 scale, henceforth referred to as D-DEM. Santoro and Strozzi (2012) have made this data set available on Pangaea, an open access database hosted by the Alfred Wegner Institute for Polar and Marine Research (AWI) (http://dx.doi.org/10.1594/PAN-GAEA.779748). The data set covers latitudes >55° N and is provided in Plate Carrée (WGS84 Lat/ Lon) projection. Like the GDEM2, it is provided in 1° × 1° tiles. The pixel posting is 3 arcsec (i.e. 1/1200°), which is equivalent to approximately 33 × 93 m at a latitude of 69° N.

DEM statistical validation method
The ASTER GDEM2 data were reprojected to UTM zone 42N and resampled to a 18.15 m × 18.15 m grid, representative of grid markings on the paper map, using the Nearest Neighbour algorithm in the ArcMap TM software environment. This was done to minimise aliasing artefacts and better relate its elevation distribution to the scanned paper map. Randomly generated spot-height data were extracted for both data sets. The original sample frame consisted of 100 random points. However, the sampling design was selected to only include unambiguous points which a) could be manually attributed an elevation value from the map and b) did not correspond to no data values in the DEM. 14 of the original 100 random spot-height values were hence excluded according to this sampling design. Therefore, the sample size was reduced to N = 86.
Probability density functions for the sample data (ASTER GDEM2 and Soviet map elevation values (N = 86 in both cases)) were produced in order to determine any similarities in the distribution of the elevation data. To further investigate whether the sample data are taken from populations with identical distributions a statistical test was carried out. To quantify the differences in the elevation distributions, the null hypothesis of this Kolmogorov-Smirnov test (K-S test) is that both data sets follow the same distribution. The K-S test is a non-parametric test that compares the empirical cumulative distribution function (ECDF) of two sample data sets. If the K-S statistic (D) is equal to zero, the Anna Maria Trofaier and W Gareth Rees ECDFs are indistinguishable. If D is equal to one, then the ECDFs are completely distinguishable.
The second approach was to undertake a direct comparison of the two spot-height data sets in the form of the root mean square error (RMSE) in height. The RMSE does not distinguish between random and systematic errors but indicates the overall general error (Nikolakopoulos et al. 2006). Neither can it provide any information on any spatial relationship of the error, although it has previously been found that errors in a DEM are indeed spatially autocorrelated and are coupled to the properties of the terrain and its complexity (Gao 1997;Wise 2011). Nonetheless, the RMSE is an established method in DEM error analysis that only requires the knowledge of sample spot-heights. It provides some information on the degree of correlation between two data sets. The RMSE was calculated for the sample elevation values of the GDEM2, h GDEM2 (x), and the map spot-heights, h map (x). In order to appropriately account for any biases, we further investigated the extreme difference values, by first undertaking linear regression analysis for the GDEM2 and the map spot-height data. Extreme difference values were categorised applying the conditional statement that any ASTER GDEM2 values that were greater or less than the map spot-heights ±6σ, were to be classed as extreme. As mentioned above, 86 sample points could be unambiguously used for linear regression analysis. Given this rather small sample data set, a t-test was carried out (the null hypothesis being that the sample means are equal) to further explore the statistical significance of the correlation in these two variables. Finally, the mean error in the sample points, e(x), was calculated as the difference in the elevation between the data sets The GDEM2 elevations were then analysed as a function of this sample error, including a linear regression analysis. A one-tailed t-test was undertaken to examine whether the slope of the regression fit was significantly different to zero.
In addition to the spot-height data validation technique, a DEM cross-comparison was also carried out between the ASTER GDEM2 and the D-DEM of the ESA DUE Permafrost project. The D-DEM data were reprojected into the same projection as the ASTER GDEM2 data (UTM zone 42N, WGS84) and the GDEM2 data set was resampled to the D-DEM's pixel size of 31.18 m × 92.97 m using the nearest neighbour resampling algorithm. Both data sets were resized to an entirely overlapping area, excluding any regions containing no data values. As before, a two-sided K-S test was carried out (N = 648,900). Furthermore, the RMSE The suitability of using ASTER GDEM2 for terrain-based extraction for these data were then calculated, at first without considering the effects of the water body mask that had been applied to the GDEM2 to eliminate the spike anomalies, and later accounting for this issue. This was done in order to explore to what extent these spike anomalies have an impact on the RMSE. Furthermore, as a means to visualise elevation differences, two transects were chosen for which elevation profiles were created.

GIS-based surface flow routing method
The ArcMap TM Hydrology module tools were used to extract the stream channel network. The core steps of any GIS hydrology analysis are to first amend the DEM with a pit-removal algorithm. Pits are DEM sink-cells; in a raster they are represented by pixels that indicate a depression, being surrounded by higher elevation and hence cells of higher pixel values. Pits trap flowlines, and therefore allow water to accumulate. Some of these pits are natural as they may be the location of lakes. However, most pits are gridding and resampling artefacts, which may be related to interpolation errors (Huggett 2002). It is therefore customary to run pit-removal algorithms first, ensuring the continuous flow of water downhill.
To create a local drainage direction (LDD) map, which shows the flow of water across and the connectivity between each pixel cell (Conolly & Lake 2006), the Flow Direction tool in ArcMap TM was run. This tool uses the D8-algorithm first presented by Jenson and Domingue (1988). The LDD map is a flow direction matrix and is the basis for calculating further derivatives. The LDD matrix was used to establish the catchment boundaries of the Vaskiny Dachi area. The software establishes catchments for connected matrix cells by automatically selecting pour points (the outlet points of the drainage basin) at the edge of the grid and then identifying the cells that contribute to each of these grid drains. This standard procedure was used to establish the boundaries of two adjacent catchment areas associated with the Mordy-Yakha river in the South and the Se-Yakha river in the North.
The stream channel network was extracted, by first determining the flow accumulation matrix (a LDD derivative associated with the total number of upstream cells that contribute to flow into each downslope cell). Since high accumulation values are regions of concentrated flow, stream channels may thus be identified. The crucial input parameter is the number of upstream cells that defines each watercourse. The threshold value enables catchment partitioning; if the number of upstream cells that contribute to each watercourse is high, then the stream network will be less dense, and there will be fewer sub-catchments (Huggett 2002). Given the rather coarse spatial resolutions of the GDEM2 a relatively high threshold of 500 upstream cells (approx. 0.16 km 2 ) was adopted. Stream channel density values in a 10 km neighbourhood-radius of all stream channel grid cells were calculated. Further, mean stream channel densities were calculated within areas known to be seasonally inundated (Trofaier et al. 2012); the intention being to investigate whether areas of seasonal inundation are related to a high density of stream channels and, therefore, whether seasonal floods are associated with hillslope runoff rather than lake basin overflow.
The stream channel network is a planar dendritic-system, characterised by its interconnectivity; stream channels flow into larger streams, intersecting at junctions that are defined by their surroundings. The stream channel network is defined by each of its branches, known as segment links. These links can be hierarchically organised. The ordering-system used for this analysis was the Strahler system, which orders upstream network termination links as a first-order segment. Secondorder segments are created by the interconnection of two first-order segments, third-order segments arise at the junction of two second-order segments and so on (Huggett 2002). One of the most important morphometric properties is the drainage density, Dd, of a catchment with area A. This was defined by Horton (1932Horton ( , 1945 to be: ( 2) where Li is the i th stream segment link and n is the highest order of stream segment links.
The drainage density is, therefore, the sum of the lengths of each channel segment within the drainage basin divided by the area of this basin (Melosh 2011). Drainage density is thought to be closely linked to climatically driven erosional processes (Tucker et al. 2001) as it is immediately associated with the infiltration capacity of the ground. As such, drainage density is an important geomorphic parameter that provides key information on surface flow and channel concentration, and its associated implications on sediment and mass transport (Bishop et al. 2012). Drainage den- sities per channel link for both the Mordy-and Se-Yakha catchments were calculated. Furthermore, the mean drainage density was investigated as a function of spring flood extent according to the July 2007 water body classification determined by Trofaier et al. (2012).

GDEM2 validation
The probability density functions of the two data sets indicate non-normal behaviour (Fig. 3). At first glance the ECDFs (Fig. 4) of both data sets follow a close but non-identical distribution. A normal kernel function was applied to the GDEM2 sample data, producing a kernel-estimated CDF given in black. As would be expected, the empirical and kernel-estimated CDFs of the GDEM2 data are identical within a 95% confidence interval. The ECDF of the map spot-height data fall beyond this significance level, departing substantially from the kernel estimated CDF at lower elevation values. Similar behaviour is exhibited by the ECDF of the D-DEM. To evaluate the degree of difference in the ECDFs of the GDEM2 and the map spot-heights as well as the D-DEM data sets, the two-sided K-S statistic was computed. The K-S statistic between the GDEM2 and spot-height ECDFs was found to be 0.279 with a p-value of 0.00246. The K-S statistic for the GDEM2 and the D-DEM was found to be 0.2671 with a p-value less than 2.2× 10 -16 . Given the extremely small p-value the null hypothesis has to be rejected. The two data sets appear to be statistically distinct. Conversely, the K-S statistic for the spot-height data set and the D-DEM was found to be 0.064 with p = 0.8779, showing excellent correlation between these two ECDFs.
The straightforward RMSE calculation between the ASTER GDEM2 and spot-height sample data yielded a value of 12.17 m. It has been claimed that the RMSE in the z-direction (i.e. the elevation) should be around 7 m for high vertical accuracy, and that an RMSE greater than 15 m is insufficiently accurate (Huggett 2002). According to this statement, the GDEM2 is sufficiently accurate for the Vaskiny Dachi area -albeit almost borderline. Nonetheless, given the localised application of the GDEM2 an RMSE of approx. 12 m is an acceptable result. The linear regression of the GDEM2 as a function of the spot-height data set resulted in a very low R 2 value of 0.15, with p-value < 0.0001 (Fig. 5). The red error bar indicates the confidence level that is derived from the standard deviation (σ) of the fit. A t-test was carried out, setting the degrees of freedom to 163 and only accepting a significant result if p ≥ 0.05 resulted in a critical t of 1.974. The t-test of the GDEM2 and map spotheight data determined a t of 1.116, which is smaller than the critical t. Hence the null hypothesis that the sample means are equal cannot be rejected. A boxplot (Fig. 6) clearly indicates the outliers. Figure 7 shows a scatter plot of this error versus the GDEM2 elevation.
The mean error, <|e(x)|>, in these two data sets is found to be 16 m. However, a bias towards a greater error may be found due to the anomalous sample point at h GDEM2 (x) = 62, where e(x) = 58 m. The anomaly is distinctly identifiable by the outliers of the boxplot (Fig. 6). In fact, linear regression analysis of the GDEM2 sample data versus the sample error gave a R 2 value of 0.1521, the slope Fig. 3. The probability density functions for both ASTER GDEM2 (a) and the Soviet map spot-height (b) data sets. The suitability of using ASTER GDEM2 for terrain-based extraction being 0.326. The one-tailed t-test resulted in a pvalue of 0.0001, which confirms that the slope of the regression line is significantly different from zero. This further illustrates the bias that the anomalous sample point introduces into this error analysis. The removal of the anomaly results in a mean error <|e(x)|>= 8.06 m and an RMSE of 10.505 m (N = 85). This procedure, therefore, reduces the RMSE to a more acceptable value, well within the limits suggested by Huggett (2002), indicative of fair accuracy of the ASTER GDEM2.
The RMSE of the ASTER GDEM2 with respect to the D-DEM was found to be 9.75 m (N = 648,900), excluding the water body mask. However, taking the water body mask into consideration, the RMSE is reduced to 9.64 m (N = 506,877). The DEMs of the study area are shown in Figure 8 for b) D-DEM, c) GDEM2 and d) the GDEM2-D-DEM difference. The sample statistics are given in Table 1. The standard deviation is greater for the D-DEM, but not unsurprisingly so given its lower spatial resolution. These sample statistics already give a good indication of the comparability of the data sets, which is further affirmed by the RMSE.
Elevation profiles are given in Figure 9, where the GDEM2 is given in grey and the D-DEM in black. The discrepancies between the two are clearly visible. In addition, the GDEM2 exhibits saw-tooth behaviour. This behaviour originates from the water body mask -which results in elevation values dropping to zero. For illustration purposes, the water body mask is not considered for the D-DEM in Figure 9, where a more natural, smooth elevation profile is found. A kernel estimated CDF of the GDEM2 sample is given in black with a 95% confidence interval (grey) and the maximum distance between the ECDFs is indicated by red lines with end-nodes for both graphs. Comparison between the ECDFs of D-DEM (dark green) and the spot-height values (green) (c).

Stream channel network
The borders of the Mordy-Yakha and Se-Yakha drainage basins were determined from the LDD matrix. According to these calculations, the two basins cover an area of 506.2 km 2 and 1087.9 km 2 , respectively. The total length of the stream segments in the stream channel network was found to be 0.9623 km long for the Mordy-Yakha drainage basin and 2.1603 km long for the Se-Yakha drainage basin. According to equation 2, this resulted in a drainage density of 1.901 × 10 -3 km/km 2 and 1.986 × 10 -3 km/km 2 for the Mordy-Yakha and Se-Yakha basins, respectively. Applying the Strahler system to the derived stream channel network resulted in seven stream orders. A topographic map of the two basins overlain by the stream channel network is shown in Figure 10. The gently sloping terrain of this area explains the low drainage density values. The drainage densities of each stream order for both the Se-Yakha and Mordy-Yakha basins are given in Tables 2 and 3, respectively. As is to be expected, the higher the stream order, the lower the drainage density. Slope gradient relationships between tributaries are equally as expected, with the lowest stream order associated with the flattest terrain and gentle slope gradients. Drainage density also tends to decrease with an increasing slope gradient, since drainage density per stream channel segment is inversely related to Strahler order. A total number of 4085 and 10792 stream channel segments for the Mordy and Se-Yakha basins were identified. 11.5% of the stream segments of the former basin and 10.1% of the stream segments of the latter basin were associated with first order streams and hence headwaters which according to the modelled network are initiated on the steeper hillslopes.
With respect to the hydrologic function of DTLBs and the connectivity of lakes in these DTLBs with the stream channel network, Trofaier et al. (2012) discussed seasonal water body flood and drainage   6. Boxplot of the ASTER GDEM2 and Soviet map elevation sample data. The means are not significantly different, although on average the GDEM2 sample data are slightly higher than the map sample data. The outlier affecting the bias is also clearly visible as the individual point. Fig. 7. The difference in GDEM2 and the map spot-heights (i.e. the sample error, e(x)) versus the elevation taken from GDEM2. Extreme difference values given in green. Linear regression fit given by dashed blue line. Low R 2 value was inspected by carrying out a one-tailed t-test showing that the slope of the regression line is significantly different from zero.  patterns. These patterns were determined from ESA's Envisat Advanced Synthetic Aperture Radar (ASAR) operating in Wide Swath (WS) mode, in the Vaskiny Dachi area. Floods are observed in the snow-off period and the collection of snowmelt waters in DTLBs was established but the output (drainage) mechanisms remain poorly understood.
In the present study, we analysed drainage density as a function of spring flood extent. This calculation showed a weak, non-linear correlation between the drainage network and seasonally inundated areas for both the Mordy-Yakha (R 2 = 0.0167) and the Se-Yakha (R 2 = 0.0068) basins.

DEMs for terrain-based stream channel extraction
The validation studies undertaken here imply sufficient correlation between the ASTER GDEM2 and the reference data sets. The K-S test showed rather poor correlation between the GDEM2 and both the Soviet map spot-height values and the D-DEM. It has been argued that the K-S test does not perform well for discrete data (Greenwell & Finch 2004); however, given that the second K-S test, performed on both DEMs over the entire study area (N = 649,800), gave an almost identical result to the K-S test of the 86-point sample data, the poor correlation cannot be attributed to the sample size and discretisation of the sample data. In fact, the spot-height values follow a similar ECDF to the D-DEM (Fig. 4c) and the K-S test between the spot-height sample and the D-DEM confirms that the populations were sampled from the same distribution. Given that the D-DEM was compiled from a set of Russian topographic maps (Santoro & Strozzi 2012), which may indeed have included the map used to produce the spot-height data set, the closeness of these distributions confirms the representativeness of the sample spot-height data, but it provides scope for debate as to whether the D-DEM is more representative of topography than the GDEM2. It needs to be acknowledged that we are comparing two data sets derived from different methods. The D-DEM was produced from the interpolation of rasterised topographic maps, which in turn are based on ground survey data, while the GDEM2 originates from automated stereo-correlation of the nadir and aft-telescopes (Hirano et al. 2003). Interestingly, it seems to be at low elevations where the greatest discrepancies between GDEM2 and D-DEM occur. Given that the ASTER GDEM Validation Team (2011) state the absolute vertical accuracy of the GDEM2 to be 17 m, it is likely that the D-DEM data are more accurate at low elevations. The greatest discrepancy between the GDEM2 and the D-DEM seems to occur for elevations under 16 m (Fig. 4b). The difference between the RMSEs of the GDEM2/ spot-heights (RMSE = 12.17 m) and GDEM2/D-DEM (9.75 m) was only 2.42 m. Both RMSE results are within the range suggested by Huggett (2002) for fair DEM accuracy. Considering the stacking numbers of the GDEM2, it was found that the spike anomalies of the DEM were not only related to the water bodies, but also had low stacking numbers. Hence, pixels associated with stacking numbers n ≤ 4 were excluded from the analysis. Overall, once spike anomalies have been removed, the ASTER GDEM2 produces a sufficiently accurate DEM that may be used for further topographic analysis.
As mentioned above, Trofaier et al. (2012) showed that a clear spring flood pattern could be determined (over the years 2007-2009) using the Envisat ASAR WS sensor. The spatial resolution of the instrument (120 m for each sub-swath (Closa et al. 2003)) is, however, such that melt waters in the stream channel network cannot be deduced. Nor can any hydrologic connectivity between water storages (such as flooded DTLBs and lakes) and the stream channel network be established. ASAR WS does not have the capabilities to spatially resolve stream channels. The current study was devised to model the stream channel network in order to analyse any potential interactions between the floods and the drainage system. The hypothesis was that the detected spring flood patterns were related to hillslope runoff, promoted by the limited infiltration capacity of the continuous permafrost environment (Granger et al. 1984). This would be in line with Romanovskiy (1961), who explains how the runoff of spring melt waters occurs through thermoerosion (the erosion of the ground surface through thermal and mechanical processes, e.g.   flowing water in ice-rich terrain) and thermokarst depressions (subsidence features related to ground-ice melt). Thermokarst depressions in the form of DTLBs do correspond to areas flooded after snow melt, however, there is only a weak relationship between flooded areas and the drainage network. We find that a low drainage density was associated with seasonal inundation. None of the calculations showed particularly high values over the flooded areas, but rather correspond to the mean range of drainage density values that were found over the entire Se-Yakha and Mordy-Yakha basins. Four possible explanations for this outcome may be that 1) both the GDEM2 vertical and horizontal resolutions are too crude for the de-tailed modelling of the stream channel network, 2) the D8-algorithm used to model the network is too simplified and more advanced techniques are needed to model the stream channel network in detail, 3) snow accumulated in DTLBs is not easily redistributed, melt waters that percolate through the snowpack eventually flood the DTLBs and hillslope runoff and saturated overland flow are less important contributors to these floods, and 4) possible subsurface flow through the thawed active layer is not to be neglected. Regarding the modelled stream channel network, drainage densities per Strahler order as a function of slope, follow the typical behaviour, increasing drainage density with increasing slope. However, the total drainage density of the two catchment areas is less than 0.002 km/km 2 per catchment. This is very low. According to Arp et al. (2012), a typical Arctic basin has a drainage density between 0.28-2.07 km/km 2 . The drainage densities of the Mordy and Se-Yakha basins are two magnitudes lower than the typical value. A low drainage density is a reflection of a hummocky, rolling landscape. However, the extremely low values imply the modelled drainage network represents too few stream channels. It is likely that the resolution of the DEM at horizontal scales is too coarse for detailed channel extraction. Hence small channels are not taken into consideration for the calculation of drainage density. Furthermore, as discussed above, the ECDFs of the GDEM2 and D-DEM departed considerably at low elevations and it is possible that the GDEM2 is not sufficiently accurate in order to perform the channel extraction procedure for these low elevations. In addition, Jasiewicz and Metz (2011) have pointed out the limitations of the standard D8-algorithm in low-relief terrain, which result in an unnatural representation of the stream channel network. Arp et al. (2012) further emphasise these limitations for flat, lake-abundant terrain. The automated channel extraction approach did not either allow an accurate investigation of connectivity between DTLBs and the drainage network. High spatial resolution optical data would provide a means of manually inspecting any input and output channels which could account for the drainage patterns observed by Trofaier et al. (2012).

Potential subsurface flow
As an aside, it should be noted that runoff follows a number of different pathways, which of course Fig. 10. Derived stream channel network with associated Strahler segment link order in the Se-Yakha and Mordy-Yakha catchments (northern and southern black outlines, respectively). The ASTER GDEM2 is used as the background layer. Anomalously high pixel values were not removed. 0.07% of the pixels were higher than 60 m a.s.l., resulting in the depicted colour stretch. may include the routing through (but is not solely restricted to) the stream channel network. Apart from surface flow originating from melt waters in rills adjacent to streams (Howard 1971), subsurface flow within the thawed active layer needs also be considered.
According to a study by Wang et al. (2009), active layer thawing processes were found to have a considerable impact on the spring runoff regime. They discuss how higher topsoil temperatures in shallow active layers increase surface runoff but this relationship was not identified for active layer depths greater than 60 cm. This would imply that increases in active layer depth allow water infiltration and subsurface flow rather than surface runoff. Active layer depth measurements are ongoing at the Vaskiny Dachi research camp, with 4 established grids. Measurements, however, are only taken in the summer. In the year 2007 the average active layer depths at these sites ranged from 72 to 112 cm (Leibman et al. 2015). More information on active layer depths throughout the thawing season is needed to establish whether subsurface flow pathways play a dominant role in the Vaskiny Dachi area.
Modifications in the runoff pathways due to a thawing active layer mirror an increased coupling of the hydrological system over the course of the summer season (Carey & Woo 2001). In particular, diffuse subsurface flow within organic rich active layers is a common phenomenon, the extent of or-ganic terrain significantly influencing the runoff regime (Quinton 2003). It is well known that snowmelt infiltration has a warming effect on the active layer as water flow results in the advection of heat (Hinkel & Outcalt 1994). This in turn facilitates lateral water movements through the active layer according to the hydraulic gradient reflected by the terrain's slope (Scherler et al. 2010). The present study's predicted channel network cannot account for any subsurface water movements due to the nature of the terrain-based model. However, Leibman and Streletskaya (1997) and Streletskiy et al. (2003) discuss active layer runoff as a mechanism for washing out saline deposits and ion migration in the Vaskiny Dachi area. A combination of both surface and subsurface routing of snow meltwaters is therefore plausible within the Vaskiny Dachi area. Unfortunately, due to the lack of appropriate field observations and measurements, no further discussion on lateral flow through the active layer can be made. With respect to overland water flow, it should be noted that high surface runoff is directly linked to the erosional potential of hillslopes (Cherkauer et al. 2013). Indeed, the Vaskiny Dachi area is well known to be landslideprone terrain (Leibman & Streletskaya 1997) with active-layer detachment slides being a common phenomenon (Fig. 11). Leibman and Kizyakov (2007) discuss how this cryogenic landslide activity stands in direct relation to water and sediment Fig. 11. Active layer detachment slide (a) and landslide body (b) in the Vaskiny Dachi study area. Willow shrubs are good indicators of landslide activity, as tall thickets of these willows are found on old landslide bottoms due to greater nutrient availability (Walker et al. 2011). Approximate landslide body and willow thicket outlined with dashed and solid black lines, respectively. Photographs taken by AM Trofaier on 4 July 2012. runoff, the development and evolution of stream channels that in turn dissect the landscape, creating valleys. The importance of these landslides on the stream channel network and hillslope runoff is not to be neglected, providing scope for further research and field investigations.

Conclusions
This paper has elaborated on the suitability of using the ASTER GDEM2 for analysing the drainage network in a low-relief environment by comparing the GDEM2 to reference data sets. Static terrainbased methods were then employed to delineate two catchments, the Mordy-Yakha and Se-Yakha basins. We derive the stream channel network using standard automated GIS procedures and explore the linkages between hillslope runoff and seasonal inundation in the Vaskiny Dachi area on the Yamal peninsula.
In the first instance, the validation study suggests that the ASTER GDEM2 is sufficiently accurate (after the removal of spike anomalies over water bodies). However, after extracting the stream channel network and exploring the drainage densities of the two basins, it was found that the values were significantly lower than for a typical Arctic basin. It is likely that the length of the stream channel network (and therefore the drainage densities of both basins) has been underestimated by the automated procedure, owing to a number of contributing factors. Firstly, not only vertical but also horizontal spatial resolution of the GDEM2 may in fact not be sufficient to accurately capture topographic variability at smallscales and therefore to extract a detailed stream channel network. Secondly, the D8-algorithm may be producing an unrealistic drainage network in the low-relief terrain. Investigating the connectivity between the stream channel network and seasonally flooded DTLBs by analysing drainage density as a function of spring flood extent only gave inconclusive results. There is only a weak relationship between the stream channel network and previously observed flood and drainage patterns. These floods however occur at low elevations, which in turn further highlights possible D8-algorithm related issues.
The reasons for the recurring spring flood pattern are not explored although it is likely that the greatest contributing factor of these floods comes from the overbanking of lakes and snow accumu-lated in DTLBs. It is suggested that drainage of the floods happens via small outflow channels. Stream channel extraction may have to be done manually in order to establish accurate connections between drainage network and water storages in the form of lakes and DTLBs. In general, we recommend the ASTER GDEM2 only to be used with caution in lowland, lake-abundant terrain. If possible, any surface hydrologic analysis should be supplemented with small-scale gridded elevation data.