Evaluating digital terrain indices for soil wetness mapping – a Swedish case study

Introduction Conclusions References


Introduction
It is well established that forestry, agriculture, transportation corridors (roads, trails), and other land-use practices can af-fect water quality (Buttle, 2011;Ahtiainen, 1992;Laudon et al., 2009;Schelker et al., 2012).Major threats for surface waters are soil disturbances such as rutting and compaction, which subsequently lead to soil erosion and increased sediment loads.This, in turn, increases water turbidity and sediment cover of gravelly stream beds (Lisle, 1989), thereby decreasing the reproductive success of fresh-water fish (Burkhead and Jelks, 2001;Soulsby et al., 2001) and macroinvertebrates (Lemly, 1982).In forestry, primary sediment sources are logging roads, skidder trails (Sidle et al., 2006), road crossings (Kreutzweiser and Capell, 2001), and related ditching activities (Prevost et al., 1999).In their review, Moore et al. (2005) found that heavy forestry machinery traffic during forest harvesting changes water flow paths; furthermore, it increases soil wetness and soil and stream temperatures.Surface run-off from ruts along slopes and wet soils also affect water quality and aquatic habitat.For example, Munthe and Hultberg (2004) reported that rut formation and damming of a stream increased the local downstream concentration of methylmercury (MeHg) by 600 % over a period of 3 years.Subsequently, Bishop et al. (2009) estimated that 9-23 % of Hg in fish in Sweden is associated with increased Hg outputs from clear-cutting.Kronberg (2014) showed that MeHg loads in clear-cut draining streams increased by 14 %.Forest harvesting, rut formation and site preparation create conditions favourable for net MeHg production because of (i) higher soil temperature on sun exposed soils, (ii) fresh organic litter from the slash that provides an energy rich carbon source for sulfur reducing bacteria, and (iii) anaerobic conditions, due to compaction of the soil following harvesting in A. M. Ågren et al.: Evaluating digital terrain indices for soil wetness mapping combination with increased water levels and standing water pooling in the tracks (Eklöf et al., 2014).
To mitigate against soil disturbances, forest operations traffic through wet and moist areas and across flow channels should be avoided.Doing so would greatly reduce environmentally and economically costly forest operation "surprises".Among these are, e.g. the increasing costs associated with non-anticipated culvert requirements for stream crossings, inappropriate delineations of machine-free zones, increases in machine downtime, loss of wood (quality and quantity) because of poor wood-landing locations, inefficient silvicultural investments (e.g.failed plantations), errors in summer versus winter cutting allocations, and accelerated costs regarding harvest block access (Arp, 2009).Until now, areas that are sensitive to soil disturbances have not yet been mapped at resolutions sufficient to be included in forestry planning operations.But, with new and reliable high-resolution flow-channel and wet-area mapping and follow-up field inspections, best-management practices can be enhanced with added financial and economic benefits, to guide machine traffic away from wet areas through on-board navigation.In addition, this could be done in compliance with 1.The new policy from the Swedish forest industry suggesting that "driving on forest soils should be planned according to soil conditions, surface waters, and cultural heritage".In this policy, rutting is acceptable or unacceptable depending on the environmental implications for each site.Any rutting in contact with or near streams and lakes is unacceptable (Berg et al., 2010).
2. The EU Water Framework Directive, intended to establish a common framework for the sustainable and integrated management of all waters.
3. Increased needs for systematic soil and water conservation planning and related biodiversity impact-mitigation efforts (Sass et al., 2012).
The aim of this study, therefore, is to further advance and test the accuracy of flow-channel and soil wetness mapping as it is currently aided by the increasing availability of lidar (Light Detection and Ranging) data for generating bare-earth digital elevation models (DEMs).In turn, these models can be used to map flow direction, flow accumulation, and flow channel networks with increasing reliability at high-metre to sub-metre resolution.Several algorithms to do this are now available, notably D8 (Jenson and Domingue, 1988;O'Callaghan and Mark, 1984), Rho8 (Fairfield and Leymarie, 1991), DEMON (Costa-Cabral and Burges, 1994), FD8 (Quinn et al., 1991;Freeman, 1991), D∞ (Tarboton, 1997), ADRA (Lindsay, 2003), andMD∞-algorithm (Seibert andMcGlynn, 2007).Each of these algorithms have their advantages and disadvantages (Pike et al., 2009).For example, D8 produces converging flow patterns only, while D∞ deals with flow convergence as well as flow divergence.Some of the more complex multi-directional flow algorithms minimise divergence where unrealistic, and give a more natural representation of flows along ridges, pits, and flats.When used in combination with DEM-derived slope layers, any of the above algorithms can be used to develop soil wetness indices such as the topographic wetness index (T WI ; Tarboton, 1997), the topographic position index (TPI; Weiss, 2001), and the cartographic depth-to-water (D TW , Murphy et al., 2007), with T WI and D TW proving useful for mapping soil (type, drainage, chemical, and physical properties), soil trafficability, and species-or community-based vegetation distributions (Murphy et al., 2011;Zinko et al., 2006;Sørensen et al., 2006;Kuglerova et al., 2014).In this regard, some of the above flow-accumulation algorithms perform better than others.For example, Kopecky and Cizkova (2010) preferred FD8 for vegetation mapping, while Sørensen et al. (2006) concluded that using D∞ improved soil wetness mapping.Murphy et al. (2009Murphy et al. ( , 2011) ) found that T WI -based wetness maps improved from D8 to D∞, and these improvements were scale dependent, while D TW maps showed better and fairly scale-independent correlations for various soil properties, including soil and vegetation type and drainage class.
Soil wetness maps intended to guide forest planning and related operational decision making need to be as reliable as possible at metre-by-metre resolution across large areas.For that reason, this study analyses and compares lidar-derived T WI , TPI and D TW maps, and seven other DEM-derived soil wetness indicators (topographic landform, flatness, puddles, toe slopes, aspect, profile, and plan curvature) in terms of their applicability, accuracy, and conformance in emulating actual soil wetness along streams and lakes for Swedish conditions.The areas selected for this article are located within the well-studied Krycklan catchment (Fig. 1).The lidarderived mapping results for this area are analysed in terms of in-the-field wetness determinations along transects that straddle wetlands and soil drainage classes across forested valleys, floodplains, moraines, and upland tills.

Soil wetness transects
For the soil wetness survey, we did not measure soil wetness but mapped indicators along line transects on three areas within the well-studied Krycklan catchment (Laudon et al., 2013) (Fig. 1).The field survey was conducted 10-14 October 2011.During that period, discharge measured at site C7 (Laudon et al., 2013) was 0.84 (Standard deviation, SD = 0.13) mm day −1 , which matched the long-term average of 0.84 (SD = 1.53) mm day −1 for the period 1981-2013.In Area 1, eight 800-850 m transects were placed perpendicular to a number of ridges.The landscape is glaciated and till is dominating the soils, apart from a flat wetland located to the northeast.The direction of the ice flow from northwest to southeast can be seen on the DEM by the orientation of the crag and tails and drumlins.In Area 2, twelve transects were placed on a long ridge-to-valley hillslope.Till covers the hillslope, decreasing in thickness towards the top according to the Quaternary deposits map (1 : 100 000, Geological Survey of Sweden, Uppsala, Sweden).Area 2 also includes a mire at the bottom of the hill.In Area 3, eight 500 m transects were placed to cross the valley and floodplain of the Krycklan stream.The upper east side of this valley is dominated by a moraine.The floodplain is filled with ice-river alluvium, containing mostly sand, gravel and boulders.The stream cut through the ice-river sediments forming ravines that become deeper towards the south.
The geographical positions for each soil wetness class transition along the transect lines were determined using hand-held GPS, with an accuracy of < 10 m in 95 % of the measurements.Soil wetness was mapped according to the instructions for the Swedish Survey of Forest Soils (Anon, 2013).Specifically, temporal variations from dry to wet were to be ignored in favour of determining the underlying soil wetness regime and the related soil wetness classes, i.e. wet, moist, mesic-moist, mesic, and dry soil, for full definitions see http://www-markinfo.slu.se/eng/soildes/fukt/skfukt1.html.The process involved estimating the depth to the average water table level during the vegetation period, in reference to the elevation rise away from open water features (lakes, streams), wetlands, ditches, and wet obligatory (hydric) vegetation.The five moisture classes were defined as follows: on "wet soils" the groundwater table forms per-manent water standing on the soil surface.Here, (i) one cannot walk dry-footed in low shoes, (ii) soils are organic (often fens), (iii) conifers occur only occasionally.On "moist soils" the groundwater table is on average at less than 1 m depth.Here, (i) one can walk dry-footed in low shoes, provided one can step on tussocks in the wetter parts; (ii) wetland mosses (e.g.Sphagnum sp.) dominate local depressions (pits), and trees often show a coarse root system above ground (germination point above soil); (iii) ditches are common; and (iv) soils range from organic (generally fens) to mineral (generally humus-podsols).On "mesic-moist soils" the groundwater table is on average at less than 1 m depth.Here, (i) one can walk dry-footed in low shoes over the entire vegetation area, except after heavy rain or snowmelt; (ii) areas with wetland mosses (e.g.Sphagnum sp., Polytrichum commune, Polytrichastrum formosum, Polytrichastrum longisetum) are common; (iii) trees show a coarse root system aboveground (germination point above soil); (iv) soils podsolic (humoferric to humic podsols); and (v) the mineral soil is covered by a thick peaty mor (thicker than on mesic soils).On "mesic soils" the groundwater table is on average at 1-2 m depth.Here, (i) one can walk dry-footed in low shoes over the area even after heavy rains/snowmelt; (ii) the bottom layer consists mainly of dryland mosses (e.g.Pleurozium schreberi, Hylocomium splendens, Dicranum scoparium); (iii) ferric podsols with a thin (4-10 cm) humus layer (mor) are common; and (iv) the bleached horizon is grey-white and well delineated against the rust-yellow, rust-red, or brownish rustred B horizon (the darker the colour, the wetter the soil).On "dry soils" the groundwater table is deeper than 2 m.Here, (i) dry soils are found on eskers, hills, marked crowns and ridge crests; (ii) the soils tend to be coarse in texture and include lithosol, boulder soil, and iron podsol formations, generally covered with a thin humus blanket on a thin bleached horizon; and (iii) there can be significant bedrock exposure.According to the survey protocol, the depth to the groundwater was not measured but was estimated based on reading the terrain in reference to the nearest open water locations such as streams, pools, and ponds.

Lidar acquisition and digital elevation model (DEM)
Since 2009, Lantmäteriet, the Swedish Mapping, Cadastral and Land Registration Authority, is generating highresolution elevation scans using lidar technology (Light Detection and Ranging) for all of Sweden, with a point density of 0.5-1 points per m 2 , an average xy point error of 0.4 m (SWEREF 99 TM), and a vertical accuracy of 0.1 m (RH 2000).The scanning of the study area was conducted during optimal conditions: after leaf fall and before snow cover, 11-14 October 2010.A 2 m × 2 m bare-ground digital elevation model (or 2 m DEM for short), with an average elevation error of 0.5 m, was generated from the ground elevation returns of the lidar signals.This was done through triangulated irregular network (TIN) interpolation.The resulting DEM was hydrographically corrected by automatically breaching roadside impoundments and by removing DEM-wide depression artefacts.

DEM processing
All DEM processing was done with ArcGIS 9.3 modelling tools and TauDEM 5.0.The 2 m DEM was used to derive the following terrain attributes in raster format: flow direction, aspect, curvature, plan curvature, cartographic depth-towater (D TW ), flat areas, landform, puddles, toeslope, topographic position index (TPI), and topographic wetness index (T WI ).These indices were evaluated at resolutions varying from 2 to 100 m.Aspect was calculated on a 2, 4, 8, 16, and 32 m resampled DEM, using bilinear interpolation.Since the aspect is given in degrees with both 0 and 360 degrees facing north, aspect was computed in radians and then sine transformed to range from −1 to 1. Curvature was derived from the 2 m DEM in the direction of slope gradient (profile curvature) and perpendicular to the gradient (plan curvature).Profile curvature affects the acceleration and deceleration of flow, while the plan curvature affects the convergence and divergence of the flow.Both curvature types were derived using windows spanning 3, 7, and 9 cells.A flatness index was derived from the 2 m DEM using the zonal statistics function to determine the standard deviation of elevations within a radius of 10, 20, and 30 m.A low standard deviation indicates a flat area.Puddles within the DEM were identified by subtracting the 2 m DEM from a smoothed DEM, with smoothing refer-ring to the mean elevations within a rectangle of 3 × 3, 5 × 5, 7 × 7, and 9 × 9 cells.Negative differences locate local puddles.Toe slopes were DEM-derived by creating a 0, 1 raster, with toe-slope cells marked as 1 and all other cells marked as 0. This was done twice by smoothing the 2 m × 2 m DEM across 3 × 3 cells and 9 × 9 cells, and selecting those cells with a slope change of 11-20 degrees.
The topographic position index (TPI) compares the elevation of a cell to the mean elevation to the surrounding cells in a specified area.Positive values represent ridges and negative TPI values represent valleys, while flat areas have a value near zero.TPI is scale dependent and was determined from the 2 m DEM using a cell moving window average of 17, 30, and 50 cells.Topographic landform classes (TLF) were generated by classifying the TPI as follows (Weiss, 2001): 1 -ridge (SD > 1); 2 -upper slope (0.5 < SD ≤ 1 ); 3 -middle slope (−0.5 < SD < 0.5, slope > 5 • ); 4 -flat slope (−0.5 ≤ SD ≤ 0.5, slope ≤ 5 • ); 5lower slopes (−1.0 ≤ SD < −0.5); 6 -valley (SD < −1).The numerically higher landform classes refer to lower slopes or valleys and would therefore be wetter than the numerically lower landform classes (ridges, upper slopes).The topographic wetness index (T WI ) (Beven and Kirkby, 1979) was calculated using TauDEM 5.0.6.T WI was defined as where a is the D∞ specific catchment area (contributing area per unit contour length) and β is the D∞ slope, in radians (Tarboton, 1997).Flow in flat areas was calculated according to Garbrecht and Martz (1997).That means that a high T WI indicates areas were much water accumulates and the slope is low.In contrast, steep slopes drain water and are therefore drier as indicated by low T WI .Since estimating slope and therefore T WI is strongly scale dependent (Blöschl and Sivapalan, 1995), it was necessary to repeat the T WI derivation by smoothing the 2 m DEM using moving windows with 2, 4, 6, 10, 14, 24, 50, and 100 m diameters.Doing so generated 2, 4, 6, 10, 14, 24, 50, and 100 m spaced T WI grids, which were then interpolated back to 2 m resolution by way of bilinear interpolation.
The cartographic depth-to-water (D TW ) index refers to the least-cost depth or elevation difference (in metres) to the nearest open water locations such as the DEM-derived streams, lakes, pools, ponds, or shoreline where D TW is set to be 0.This index (Murphy et al., 2009(Murphy et al., , 2011) ) was derived from the 2 m DEM as follows: first the DEM was filled and the flow direction and FA data layers were generated using the D8 method (Jenson and Domingue, 1988;O'Callaghan and Mark, 1984).The resulting FA raster was used to derive the topographically defined flow channel networks with 0.5, 1, 2, 2.5, 4, 5, 8, 10, and 16 ha flow initiation thresholds.D TW was then determined for each of the resulting flow networks by determining the least elevational differences between each DEM cell and its nearest stream cell according to Hydrol.Earth Syst.Sci., 18, 3623-3634, 2014 www.hydrol-earth-syst-sci.net/18/3623/2014/ the least-elevation path between these cells.Mathematically, where dz/dx is the slope of a cell along the least-elevation path, i is a cell along the path, a is 1 when the path crosses the cell parallel to the cell boundaries and 1.414214 when it crosses diagonally; x c represents the grid cell size (m).

Data preprocessing
Each variable was tested for normality using the Kolmogorov-Smirnov test (IBM SPSS Statistics 19).
As a result, D TW , T WI , TPI, and flatness indices were log transformed to fit normality and to reduce the heteroscedasticity of model residuals.All index variables were scaled and centred using z scores (Eriksson et al., 2006a).Landform types were entered into the analyses as dummy variables.

Orthogonal projections to latent structures (OPLS)
The transect data for soil wetness were used to validate the digital terrain indices through direct point-by-point comparisons (of each 2 m × 2 m cell).This was done using the multivariate statistical program SIMCA-P+ 12.0.1,Umetrics, Umeå.The statistical tests were done using the recently developed orthogonal projections to latent structures method (OPLS; Eriksson et al., 2006a, b).This method, which is similar to principal component analysis, separates the variations of the predictors (the DEM-derived soil wetness predictors) into two parts: one part that is predictive of the fielddetermined soil wetness estimates and is plotted against the horizontal x axis (denoted "pq[1]"), and one part that is not and is plotted along the vertical y axis (the non-predictive axis, denoted "p o s o [1]").In the resulting plot, the variables with pq[1] loadings that score high or low on the predictive axis are highly positively or negatively correlated to soil wetness (Fig. 2).SIMCA-P+ also calculates the influence of each X variable in the model, called the variable importance in projection (VIP).Variables with VIP > 1 are the most relevant for predicting the variable of interest (Y ).In Fig. 2, variables with VIP > 1 appear in black text, and grey otherwise.
Variables with dots that remain closely clustered across all three study areas and also fall along the far sides of the x axis are strongly correlated with soil wetness across scales and landforms.If the dots for the same variable do not cluster, then that variable is sensitive to scale, and only those dots in the closer neighbourhood of the far sides of the x axis would achieve the VIP > 1 status.If the dot locations vary strongly from Area 1 to Area 2 and 3, those locations will be strongly influenced by landform type.

Confusion matrix
The overall conformance of D TW ≤ 1 m relative to the wet and moist soils within Areas 1, 2, and 3 was further tested by way of a confusion matrix.This was done using four D TW performance groups: (i) true positive (T P ) when the D TW ≤ 1 m correctly identified a wet area; (ii) true negative (T N ) when D TW > 1 m correctly identified a dry area; (iii) false positive (F P ) or Type I error, i.e. when D TW ≤ 1 m predicts wet soils when the soils are actually dry; and (iv) false negative (F N ), or type II error, i.e.D TW > 1 m predicts dry soils when the soils are actually wet.These tests were applied to the D TW determinations as these vary by the D TW -defining flow networks, using the flow-initiation thresholds from 0.5 to 16 ha.The accuracy (A CC ), or efficiency, of each of the D TW ≤ 1 or > 1 m locations was assessed by way of .
(3) Also used was the Matthews correlation coefficient (M CC ) for which a value of 1 indicates a perfect fit, a 0 yields a results that is no better than random prediction, and −1 indicates a perfect negative correlation.M CC was calculated as Equations ( 2) and (3) were also used to determine the A CC and M CC values for the currently used 1 : 12 500 property map for Sweden in reference to the field determined soil wetness values.This map contains all officially recognized surface water and wetland features, including mires.

Results
The OPLS model quantified the influence of each variable on rendering reliable map predictions about soil wetness by scale and across Areas 1, 2, and 3 as shown in Fig. 2.Among the variables, Fig. 2 indicates that D TW has a consistently strong influence on predicting soil wetness.In detail, all the D TW dots for Areas 1, 2, and 3 cluster closely within the same neighbourhood along the negative portion of the predictive soil wetness axis (pq[1]).In contrast, the T WI dots do not form a tight cluster but cut across the positive side of the horizontal axis, with prediction performance for T WI calculated at 24 m resolution (T WI 24) for Areas 1 and 3, and T WI 50 for Area 2. In addition, T WI calculated at 2 m resolution (i.e.T WI 2) scored low on the predictive axis and high on the orthogonal axis (p o s o [1]), thereby indicating that highresolution DEMs are not suitable for T WI -based soil wetness determinations.This is also illustrated in Fig. 3.The TPI and TLF variables both showed their best soil wetness prediction performance across the three areas when derived from the 50 m DEM (Fig. 2), but these indices were not as good soil wetness predictors as D TW or the best T WI predictors.The flatness index was a good predictor for wet areas in Area 1, but did not correspond with the soil wetness determinations in Area 2 and 3. Toe slopes, aspect, puddles, curvature, and plan curvature all ranked low along the predictive horizontal axis but high on the orthogonal axis in Fig. 2, indicating that these variables are also not good soil wetness predictors.
Figure 3 illustrates the effect of scale on using T WI as DEM-based soil wetness T WI predictor for Area 1 at 2, 24, and 100 m resolution.In detail, cells with high T WI values have large contributing drainage areas and low slopes.These areas are wetlands and are displayed in blue on the map.Cells with low T WI values are well-drained dry areas and are displayed as brown on the map.The overlay of T WI on the streams and associated wetlands previously mapped at 1 : 12 500 m scale (red crosshatch) shows that T WI -when derived from the 2 m DEM -does not correspond with the wet-area distributions.However, T WI when derived from the ≥ 24 m DEMs produced good agreements for each of the three areas (Figs. 2 and 3).This is further demonstrated in Fig. 4 by the correspondence between the mapped T WI values and the field-determined soil wetness along the transect lines.The corresponding OPLS-derived optimal flowinitiation thresholds for D TW -predicted soil wetness varied from 0.5 to 2 and 5 ha for Areas 1, 2, and 3, respectively (Fig. 2a-c).
Applying the A CC and M CC accuracy metrics to the D TWsuggested wet soil locations across the combined study areas produced best-overall values of 87.1 % and 0.52 with flowinitiation thresholds of 2 and 1 ha, respectively (Table 1).By study area, the best-attained A CC and M CC values varied from 87 to 92 % and from 0.46 to 0.70.These ranges widened from 72 to 92 % and from 0.15 to 0.72 by varying the D TW -determining flow-initiation thresholds from 0.5 to 16 ha.Across this threshold range, A CC and M CC values were more affected by study area than by resolution, with A CC and M CC generally decreasing from the 1 to the 16 ha flowinitiation threshold, while the opposite occurred for Area 3.This was also reflected by the across-area standard deviations of the A CC and M CC estimates by flow-initiation threshold: least for 1 to 2 ha for Areas 1 and 2, and least for 8 to 16 ha for Area 3. The somewhat decreasing A CC and M CC performance for D TW using the flow-initiation threshold of 0.5 ha is likely due to two reasons:  0.5 ha 1 ha 2 ha 2.5 ha 4 ha 5 ha 8 ha 10 ha 16 ha Average SD Map (1 : 12 500) A CC (%) network to smaller and smaller reaches is directly limited by DEM resolution; and (ii) with decreasing flow initiation, flow channels become drier for longer periods during each year.Figure 5 illustrates differences between the wet areas of the Swedish properties map and the D TW map: many small previously unmapped wet to moist areas along the transects conformed to the latter.Specifically, the D TW map improved M CC for Area 1 (Table 1).For Area 2, both maps produced similar A CC and M CC values, but the D TW < 1 m criterion did not fully capture the extent of wetlands below the long hillslope.This can, however, be remediated by extending the D TW -based wetland delineations towards D TW > 1 m.For the Area 3 transects across the valley, D TW improved A CC slightly, but improved M CC strongly.Generally, M CC is a better measure of model performance than A CC (Girard and Cohn, 2011).

Discussion
This study showed that the wet areas can be identified and mapped by way of DEM -derived data layers.The generally close agreement between the field-determined locations of wet soils and their corresponding T WI and D TW values confirm the underlying assumptions that water movements across the study areas are driven by gravity and that topography controls the resulting water pathways.For the boreal forest landscape, these assumptions are consistent with delineating the enduring soil drainage and wetland distributions (Rodhe and Seibert, 1999).
For the T WI determinations, several methods of calculating flow accumulation exist, from the simple D8 algorithm (O'Callaghan and Mark, 1984) to the more complex MD∞ (Seibert and McGlynn, 2007).Here, we used Tarboton's D∞ method (Tarboton, 1997) since Sørensen et al. (2006 ) showed that this method gave the best results for predicting soil wetness.Which particular T WI numbers indicate wet soils, however, vary by landscape, climate, and scale (Zinko et al., 2005;Western et al., 1999;Grabs et al., 2009;Güntner et al., 2004).The DEM scale is particularly problematic for the T WI -affecting slope derivation: with coarse DEM grids, flow channels and local depressions are not properly delineated; with fine DEM grids, local T WI variations are too strong to separate wetlands from uplands (Sørensen and Seibert, 2007).Figures 2 and 3 demonstrate that the scale as to which the T WI slope should be calculated depends on the landscape: in Area 1 and 3, the terrain was undulating and a 24 m DEM gave the best results.For Area 2, the 50 m DEM T WI derivation gave the best results along the hill slope.Since the best T WI derivations require DEM smoothing for all three areas (Figs. 2 and 3), T WI is not the best method for mapping small-scale variations of wet areas along wetlands, streams, and lakes.
In contrast to T WI , D TW -based soil wetness mapping does not require DEM smoothing, and the amount of detail so revealed is mainly limited by DEM accuracy and resolution (Murphy et al., 2007(Murphy et al., , 2009(Murphy et al., , 2011)).Figure 2 demonstrates that the D TW -based wet-area delineations remain close to the horizontal axis of the OPLS projections, and are therefore not only less sensitive to DEM scale but also remain fairly robust by stream flow initiation.Hence, by varying the threshold for stream flow initiation, spatial and temporal variability of the stream network and adjacent wet soils can also be modelled, with setting lower and larger threshold values for wet and dry seasons, respectively.For example, a 4 ha flow-initiation threshold tends to reflect end-of-summer flow and soil wetness.In comparison, D TW maps based on 1 m DEMs and setting 1-2 ha flow initiation thresholds are useful (i) for planning or locating road-stream channel crossings except for sandy landforms such as, e.g.floodplains and glacial outwash (Campbell et al., 2013); and (ii) for estimating the distribution of wet-area obligatory species (Hiltz et al., 2012;White et al., 2012).Lowering the flow-initiation threshold further to 0.5 and 0.25 ha would emulate soil wetness and soil trafficability during wet summer weather and the snowmelt season (Murphy et al., 2011).Going from flat to montane areas may also require a downward change in the flow initiation threshold from 4 ha, as a reported by Jaeger et al. (2007).In arid regions, flow initiations and related depthto-water mapping may increase to 1000 ha or more during dry and drought seasons.
In reference to the case study area, setting an appropriate flow-initiation threshold for the DEM-derived flow accumulation pattern depends, in part, (i) on the surface expressions of the landscape or landform as these vary from hummocky to flat, and (ii) on substrate permeabilities as these vary from fast to slow (Jutras andArp, 2011, 2013).In terms of the undulating terrain covered by compacted glacial till about 5-10 m deep in Areas 1 and 2 (Seibert et al., 2009), saturated hydraulic conductivities would decrease rapidly with depth (Nyberg et al., 2001), as is mostly the case on till deposits across Scandinavia and elsewhere (Nyberg, 1995;Beven and Germann, 2013).As a result, (i) 90 % of lateral water movements on compacted tills would occur within the top 50 cm (Bishop et al., 2011), and (ii) the soil portions of these till deposits would saturate quickly during wet weather and wet seasons.Hence, subsurface flow would remain shallow and converge quickly into down-slope flow channels, each with low flow accumulation requirements.In Area 3, the terrain has high hydraulic conductivities (Koch et al., 2011), since it cuts across ice-river alluvium, which is mostly composed of sand, gravel, and boulders, and the adjacent slopes are dissected by steep and deep ravines.Here, (i) water would drain quickly and deeply (Aneblom and Persson, 1979), (ii) more flow accumulation drainage area would be required for water to re-appear within the surface channels, and (iii) only the areas immediately next to the water-filled channels would remain wet.As a result, D TW maps with flow-initiation thresholds from 8 to 16 ha were the most accurate for this area (Table 1).
Many of the transect-determined wet areas along streams and more diffuse areas along valley bottoms of this case study were not marked on Sweden's current property map.But, with the D TW methodology, these areas could be mapped (Fig. 5) with an overall mapping improvement for Areas 1 and 3 (Table 1).For Area 2, the D TW < 1 m criterion would need to be extended to D TW > 1 m to capture all of the wetland areas.In some cases, using a flatness index (e.g.White et al., 2012: mean elevational standard deviation < 0.1 m within 30 m neighbourhood of each grid cell) can also be used to extend the D TW -located wetland areas towards the actual wetland borders.For the Swedish context, however, large wetlands are already well mapped, but using these delineations in combination with D TW , as in Fig. 5, improves the high-resolution delineation of all the smaller wet areas next to streams and lakes.
It is suggested that the D TW -derived soil wetness maps can be used to reduce environmentally and economically costly off-road traffic surprises such as unacceptable rutting.For example, a recent D TW advance dealt with mapping potential and actual soil disturbance impacts for the purpose of off-road soil trafficability assessment (Campbell et al., 2013).Additional forestry benefits refer to improving harvest scheduling (summer versus fall versus winter), in-field harvest navigation, selecting tree seedlings by species for within-block planting dry versus moist and wet sites, optimising block access routes, and guiding in-block harvesting and wood-forwarding trails (Arp, 2009).Further operational benefits refer to knowing how to lay-out the harvest and skidding trails depending on current or forecasted weather conditions by changing the D TW index by flow-initiation threshold (Vega-Nieva et al., 2009), i.e. when and where to harvest and/or to scarify depending on the prevailing weather and soil conditions.The D TW -determining flow-channel net-works are of interest in themselves as these determine where trails should avoid ephemeral streams, or where to use brush mats to avoid soil compaction and subsequent sub-surface flow blocking.(For practitioner accounts of D TW -derived benefits within the forestry, lay-out for oil and gas extraction, and park management, see http://watershed.for.unb.ca/media/.)Elsewhere, D TW -generated data layers already proved useful in systematic wetland border delineations and wetland classification (Murphy et al., 2007;White et al., 2012).Similarly, D TW -derived maps can be useful to forecast upslope flooding and associated soil wetness distribution once streams are blocked by, e.g.roads, trails, beaver dams, and logs falling across streams.In terms of vegetation indexing, Hiltz et al. ( 2012) was able to relate a plot-based indexing of vegetation soil-moisture regime preferences to log 10 (D TW ).To that effect, Zinko et al. (2005) and Kuglerova et al. (2014) found a strong relationship between plant species richness, groundwater levels, and local groundwater discharge areas.
The highest resolution bare-earth DEMs produce the best D TW results across natural landscapes.However, where roads cross streams, it is necessary to ensure natural flow connections across roads.This can be done by breaching the DEM (i) at all known culvert and bridge locations, and/or (ii) at potential road blocks and related ditch diversions.Where there are flow-path uncertainties, or where the streams braid, flow accumulation algorithms other than D8 can be used to generate the D TW -determining flow channel network.Where the land is drained by way or subterranean infrastructure, it becomes important to "burn" this infrastructure into the bare-earth DEM to prescribed depths.Seepage locations and springs that occur outside upslope topographic control (e.g.artesian wells; seepage emerging from permeable bedrock) can also be incorporated into the D TW -generating algorithm by adding these locations as additional D TW = 0 defining locations.

Concluding remarks
D TW and T WI are both useful soil wetness predictors.However, T WI soil wetness delineations are sensitive to scale and landscape variations, and are limited in providing soil wetness detail at less than their optimal resolution of 24 m.In contrast, D TW is fairly scale-independent in predicting wet areas and produces a resolution-consistent wet-area delineation accuracy of at least A CC = 80 % and M CC = 0.40.In addition, D TW can be further optimised by choosing appropriate flow-initiation thresholds according to seasons (including climatic region) and landforms.Doing so would enable a systematic reduction of false positive and false negative wetarea determinations.In conclusion, D TW maps have the potential to form the next generation of high-resolution wetarea maps, and the process of doing that would find many applications in forestry and elsewhere.

Figure 1 .
Figure 1.Locator map for Areas 1, 2, and 3 with the Krycklan catchment with its Quaternary deposits.The black lines show the location of the study transects.

Figure 2 .
Figure2.OPLS loading plots for Area 1, 2, and 3 and their DEMderived terrain indices regarding soil wetness prediction.Variables in black/grey text have a higher/lower influence on soil wetness prediction, respectively.Variables that cluster closely within the same neighbourhood along the far sides of the x axis are the more robust soil wetness predictors across DEM scales and landforms.For reference, the OPLS projection for soil wetness moisture class (black dot) is also shown.

Figure 4 .
Figure 4. Soil wetness transects (coloured lines) for Areas 1, 2, and 3 on top of the 24 m DEM-derived T WI maps (left panel).The hill-shaded 2 m DEM is overlain by the cartographic depth-to-water (D TW ) classes ranging from 0 (dark blue) to 1 m (light blue), with flow channels mapped using a 1 ha flow initiation threshold (right panel).The lower panels show close-ups for Areas 1 (1 ha flow initiation) and Area 3 (10 ha flow initiation).

Table 1 .
Accuracy (A CC , %) and Matthews correlation coefficient (M CC ) calculated for Areas 1, 2, and 3 by D TW determining flowinitiation threshold, also showing the averages and standard deviations across the areas and the flow-initiation thresholds.The best results are highlighted in bold.For comparison, numbers for the Swedish Property map (1 : 12 500) are also given.The comparison between map data and ground truth were done on a cell by cell basis (area 1 n = 3903 cells, area 2 n = 3218 cells, and area 3 n = 2268).

Figure 5 .
Figure 5. Field-mapped soil wetness superimposed on today's most high-resolution map (property map 1 : 12 500).The wet areas of the Swedish property map are marked in beige (left) and are superimposed by the D TW < 1 m map (right; 1 ha flow initiation).