Quantification of upland thermokarst features with high resolution remote sensing

Climate-induced changes to permafrost are altering high latitude landscapes in ways that could increase the vulnerability of the vast soil carbon pools of the region. Permafrost thaw is temporally dynamic and spatially heterogeneous because, in addition to the thickening of the active layer, localized thermokarst features form when ice-rich permafrost thaws and the ground subsides. Thermokarst produces a diversity of landforms and alters the physical environment in dynamic ways. To estimate potential changes to the carbon cycle it is imperative to quantify the size and distribution of thermokarst landforms. By performing a supervised classification on a high resolution IKONOS image, we detected and mapped small, irregular thermokarst features occurring within an upland watershed in discontinuous permafrost of Interior Alaska. We found that 12% of the Eight Mile Lake (EML) watershed has undergone thermokarst, predominantly in valleys where tussock tundra resides. About 35% of the 3.7 km2 tussock tundra class has likely transitioned to thermokarst. These landscape level changes created by permafrost thaw at EML have important implications for ecosystem carbon cycling because thermokarst features are forming in carbon-rich areas and are altering the hydrology in ways that increase seasonal thawing of the soil.


Introduction
Temperature is increasing in northern high latitudes in response to the radiative forcing associated with increasing greenhouse gases and changes in albedo (Chapin et al 2005). Surface air temperature has increased 0.35 • C decade −1 since the 1970s, with dramatic increases since the 1990s (Polyakov et al 2002, Serreze and Francis 2006, Euskirchen et al 2007, Hinzman et al 2005 and temperatures are predicted to rise Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. by up to 6.5 • C in northern high latitudes by the year 2100 (IPCC 2007). This warming affects the spatial distribution of permafrost, which currently occurs within 24% of the ice-free land area in the northern hemisphere (Zhang et al 1999). It is estimated that 25-90% of areas with near-surface permafrost will transition to seasonally frozen ground by the year 2100 (Anisimov and Nelson 1996, Lawrence et al 2008, Saito et al 2007, Jafarov et al 2012, and rate and areal extent of permafrost thaw has already begun to increase (Jorgenson et al 2006, Camill 2005, Hinzman et al 2005, Osterkamp et al 2009, Romanovsky et al 2010, Schuur et al 2008. This thaw could result in a strong positive feedback to climate change if a portion of the estimated 1670 Pg of soil organic carbon in permafrost regions is released to the atmosphere (Schuur et al 2008, Tarnocai et al 2009.
Understanding the relationships among warming, permafrost thaw, and carbon cycling is a challenge because permafrost thaw is temporally dynamic and spatially heterogeneous. Various disturbances and feedbacks with hydrology, vegetation, and geomorphology additionally complicate the response of permafrost to warming (Grosse et al 2011). In addition to the predicted thickening of the active layer (seasonally thawed surface layer) (Anisimov et al 2007), thawing of ice-rich permafrost can result in thermokarst formation (Hinzman et al 2005, Jorgenson andOsterkamp 2005). Thermokarst landforms often have very different physical and biogeochemical environments compared to the surrounding un-affected permafrost landscape, with impacts on both rates of carbon uptake and emission (Lee et al 2010, Vogel et al 2009, Walter et al 2007, Jones et al 2012. Ultimately, net landscape-scale carbon exchange will depend on the cumulative response of the mosaicked landscape created by permafrost thaw. Therefore, to estimate the rate and magnitude of carbon efflux from the landscape it is imperative to quantify the size and distribution of thermokarst landforms. The wide diversity of thermokarst features depends on local permafrost and ground ice conditions, landscape position, hydrology, soil texture, and surrounding vegetation (Jorgenson and Osterkamp 2005). Sixteen primary modes of surface permafrost degradation have been identified, ranging from features with substantial thaw settlement (2-20 m) such as glacial thermokarst and drained lake basins to features with relatively low subsidence (0.2-2 m) such as thermokarst gullies and non-patterned thawed ground (Jorgenson and Osterkamp 2005). One of the major limitations for estimating the size, number, and spatial distribution of thermokarst features is attaining adequate spatial coverage because high latitude ecosystems are vast, isolated, and often inaccessible. Remote sensing allows us to observe these areas in a cost-effective manner and on the timescales needed to answer questions about the rapidly changing landscape. So far, the remote detection of thermokarst has been limited to relatively large and visually distinct features, such as thermokarst lakes . However, to capture the true magnitude of thermokarst, it is also necessary to expand observation capabilities towards detection of small, subtle thermokarst forms. With recent advancements in spectral and spatial resolution this is becoming possible.
The objective of this study is to assess the feasibility of very high resolution remote sensing for detecting small irregular thermokarst features (pits, gullies, and long water tracks) in an upland permafrost area in the discontinuous permafrost zone where such features are numerous. Similar thermokarst gullies and pits occur throughout the ∼64 000 ha Toklat Basin that extends 80 km to the west of our study site in Interior Alaska, and in many other high latitude upland areas. By performing a land cover classification targeting thermokarst features in IKONOS satellite imagery, we quantify the spatial extent of thermokarst and determine their spatial distribution and landscape position. Using this information, along with field observations of thaw depth and organic soil depth, we assess the vulnerability of soil carbon pools in a typical upland area in the discontinuous permafrost zone of Interior Alaska.

Site description
Our study site is the Eight Mile Lake (EML; 63 • 52 42 N, 149 • 15 12 W) watershed and adjacent hill slopes located in the northern foothills of the Alaska Range near Denali National Park and Preserve, Interior Alaska (Schuur et al 2007). The study areas covers 10.4 km 2 and is located at ∼700 m elevation on a gently sloping (∼5%), north facing glacial terminal moraine that dates back to the Early Pleistocene. Nenana gravel constitutes the major bedrock formation of the area and due to multiple glacial advances and retreats a landscape of rolling hills and valleys has formed (Wahrhaftig 1958). The hilltops have gravelly soils associated with glacial till and are covered by a thin aeolian silt cap, while the lower slopes and valleys have thick peat accumulations over fine-grained soils associated with colluvial deposits (Osterkamp et al 2009).
The study area is underlain by permafrost but it is also mapped within the vulnerable band of discontinuous permafrost near the point of thaw due to the combination of its elevation and geographic position (Yocum et al 2007, Romanovsky et al 2007. Permafrost at EML is classified as climate-driven, ecosystem protected (Shur and Jorgenson 2007) because the permafrost was formed under an earlier cold climate independent of ecological conditions and later protected from degradation by vegetation and soil development (Osterkamp et al 2009). The long-term mean annual air temperature ) of the area is −1.0 • C and the growing season mean temperature is 11.2 • C, with monthly averages ranging from −16 • C in December to +15 • C in July. The long-term annual precipitation is 378 mm with a growing season precipitation of 245 mm (National Climatic Data Center, NOAA). Deep permafrost temperature has been measured in a borehole at the site since 1985 and during this time mean annual ground temperature rose by 1 • C (1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998) and thermokarst terrain has developed and expanded (Osterkamp et al 2009).

Image acquisition and field sampling
A multispectral IKONOS satellite image covering the EML watershed and adjacent hillslope was collected on 25 June 2008 (figure 1). The commercial IKONOS satellite was launched in September of 1999 and has a polar, circular and sun-synchronous 681 km orbit with a swath-width of 11.3 km, resulting in a revisit rate of three to five days for off-nadir and 144 days for true-nadir. The multispectral sensor utilized in this study has four spectral bands (blue 450-520 nm; green 520-600 nm; red 625-695 nm; near infrared 760-900 nm). IKONOS Geo products have a spatial resolution of 3.28 m for multispectral bands with a position accuracy of 5 ± 3 m. To improve this georectification we collected eleven widely distributed ground control points in the study area with a differential global positioning system (DGPS). We then identified these ground control points in the IKONOS image and rectified the image to these DGPS-based coordinates using a second order polynomial transformation with nearest neighbor resampling. With this rectification we achieved an average root mean square error (RMSE) of 0.98 m.
First, an unsupervised Isodata classification (Tou and Gonzalez 1974) of the IKONOS image with 15 initial classes was conducted to understand the spectral variability of land cover in the study area. Based on this information and our general field knowledge of the study area we defined eight potentially spectrally separable land cover types. We then collected detailed field information on vegetation and soils within these eight a priori land cover types (table 1). Approximately 10 field training sites were sampled within each class, resulting in a total of 78 training sites distributed throughout the watershed (table 1, figure 1). Our sampling strategy focused on spatial coverage, therefore rapid assessment quantification metrics were used. At each site, four contiguous (2 × 2 m) plots were established within homogeneous areas that were representative of the site's vegetation. Within each plot, the per cent of ground cover (vegetation, water, bare soil) was estimated and geographic location was taken with DGPS. Vegetation was further categorized by the per cent cover of functional groups (graminoid, herbaceous, non-vascular, shrub, and tree) and per cent cover of each species. Due to identification difficulties and time limitations, Salix spp., Alnus spp., Dicranum spp. and Sphagnum spp. were only identified to genus and rare non-vascular plants were grouped into a non-vascular 'other' category. These data (per cent cover of vegetation functional groups, species composition, vegetation height) were used to further characterize the eight land cover types (table 1), which we compared to the Denali landcover mapping project (DLMP; Stephens et al 2001) and the Circumpolar Arctic vegetation map (CAVM; Walker et al 2005). At the end of the growing season (late August) active layer thickness (depth to permafrost) was measured at each site by inserting a metal rod until the impermeable permafrost layer was reached. Differences in thaw depth among land cover classes were evaluated using an ANOVA with a Tukey post-hoc test and analysis was performed in R (R Core Development Team 2011). In addition, depth of the near-surface organic horizon was measured in small soil pits dug at each site. Because deep organic horizons (40-60 cm) occurred in some classes, and due to field time limitations, pits were only dug to ∼20 cm. If the organic horizon extended below 20 cm the site was considered highly organic and grouped into a >20 cm category.
Fine-scale differences in elevation of the study area were measured with a survey-grade DGPS. One GPS unit (Trimble 5400) was placed at a nearby USGS geodetic marker (WGS84, 63 • 53 16.56 N, 149 • 14 17.92 W), which acted as the reference receiver. Using a second GPS unit (Trimble 5400) secured to a backpack, a kinematic survey was conducted by walking transects of the 10.4 km study site. Geographic position and elevation were collected at 5 s intervals, yielding a total of 68 000 data points. These data were post processed with methodology developed by UNAVCO using Trimble Geomatics Office software (Dayton, Ohio). A combination of these DGPS point measurements of elevation and USGS Hydrologic Line Graphs Maps were used to create a digital elevation model (DEM) with 5 m grid cell size of the study area with the Topo to Raster function in the 3D Analyst Tool of ArcGIS 9.3. Slope was derived from the resulting DEM after it was smoothed using mean focal statistics (10 × 10 cell rectangle) with the Spatial Analyst tool in GIS.

Data preparation
IKONOS bands 3 and 4 were used to calculate the normalized difference vegetation index (NDVI), which is a measure of photosynthetically active vegetation and can be used as an indicator for 'greenness' and productivity (e.g., Raynolds et al 2006). We expected that NDVI could potentially help separate thermokarst classes because vegetation within thermokarst at our study site has been shown to have higher productivity than surrounding tussock tundra (Schuur et al 2007). To reduce the spectral overlap between some training classes, we also performed a principal components analysis (PCA) on IKONOS bands 1-4. Based on visual inspection of spectral separation, we used the first three PCA bands in subsequent analysis.

Supervised classification
For the final classification we combined six data layers (PCA1, PCA2, PCA3, NDVI, elevation, slope) and performed a maximum likelihood classification based on spectral signatures of the training sites. The spectral separability of individual training data were visually evaluated in scatter plots combining the different data layers. Classes were merged if their spectra appeared to have insufficient separation. Originally the spruce, alder and willow classes were evaluated independently, but because they spectrally were very similar they were merged into a mixed riparian class. Post classification majority filtering was done to remove small numbers of pixels within larger contiguous classes and to achieve a smoother final classification. All classification analysis was done with ERDAS Imagine 9.1 software (Atlanta, GA).

Classification accuracy assessment
Independent validation sites representative of each land cover class were collected during the field campaign in 2009 (78 sites total; ∼10 within each class) and from high resolution true-color aerial images taken of the watershed in 2003 (203 total). Cross-tabulated error matrices of the mapped classes and validation sites were used to assess classification accuracy (Congalton and Green 1999) and estimate overall accuracy, user's and producers accuracies, and the Kappa statistic. The producer's accuracy indicates the probability of a reference pixel being correctly classified (omission error), while the user's accuracy is the probability that a pixel on the map accurately represents that category on the ground (commission error). The Kappa statistic incorporates the off-diagonal elements of the error matrices (i.e. classification errors), which represents agreement obtained after removing the proportion that could be expected by chance (Jensen 2005).

Results
The EML watershed grades in elevation from 680 m in the south to 632 m in the north, where the lake is situated. Superimposed upon this overall gradient in elevation are rolling hills and valleys that differ in elevation by roughly 15 m. The final thematic map derived from the supervised classification shows that hilled areas correspond with the Hilltop Shrub (HS) class, which transition to the tussock tundra (TT) class found in the valleys (figure 2). The two classes associated with permafrost degradation (Thermokarst, TK, and Thermokarst Pool, TP) mainly occur in the valleys of the watershed within the tussock tundra class (figure 2). In general, thermokarst features at EML are elongated (∼300-700 m), slightly incised (<1 m) gullies that drain from the valleys towards the lake (see inset of figure 2 for a close-up of typical thermokarst patterns). On average, ground subsidence within thermokarst features is less than 1 m, and while individual gullies are relatively narrow (1-5 m wide), aggregated together they can affect areas up to 75 m wide (figure 2; Belshe et al 2012). The mapped thermokarst classes combined (TK and TP) cover 1.31 km 2 , which is 12.1% of  Classes with organic depths less than 20 cm (MR, HS, ST) are shown in box-plots, with solid line denoting mean and box denoting the 25th and 75th quartile. Classes grouped into the >20 cm category (TT, TK, TP) have boxes that start at 20 cm and fade towards 40 cm to represent that previous site measurements show organic depths are greater than 20 cm (ranging from 37 to 54 cm) for these classes.
All training sites of thermokarst classes (TK, TP) and tussock tundra (TT) had soil organic layers greater than 20 cm thickness, while all other classes had organic layers shallower than 20 cm (figure 3). Our field method of categorizing all training sites with organic layers >20 cm as highly organic does not allow us to statistically differentiate between classes because we did not measure variance within sites with thick organic layers due to the challenge of sampling into frozen soil with limited field time. However, previous work in the same watershed reported that soil organic layers range from 37 ± 3.1 cm at minimally thawed tundra sites to 54 ± 4.8 cm at extensively thawed sites (Hicks Pries et al 2012). These site-specific data along with our field data across the watershed show that permafrost subsidence at EML is occurring in areas with relatively large soil organic layers.
Thermokarst classes (TK, TP) have greater maximum seasonal thaw depths than other land cover classes (figure 4). On average, thaw depth within TK and TP classes was 73 and 85 cm, respectively. The mixed riparian class (MR) did have a few sites near the lake inlet that exhibited deep thaw depths (66 and 81 cm), but on average all other land cover classes had thaw depths of 44 cm or lower. MR sites near the lake inlet most likely exhibited deeper thaw depths due to the thermal influence of a thaw bulb developing beneath the lake. Notably, the tussock tundra (TT) class had much lower and less variable thaw depths than the thermokarst classes, even though they are adjacent to each other. This confirms that soils within thermokarst are seasonally thawing to a greater extent than any other land cover class.
The overall classification accuracy of the IKONOS image was 83%, and there was adequate agreement between the final map and training data (Kappa value of 0.79). The user's and producer's accuracies of individual classes were consistently high, ranging from 69% to 100% (table 3). Thermokarst classes had a producer's accuracy of 79% and 83% and a user's accuracy of 75% and 100%, for TK and TP classes respectively. For the thermokarst classes, the highest confusion occurred between the TK and ST class, with 15% of TK training pixels being classified as ST. This misclassification most likely occurred because sloped sites within the TK class have very similar vegetation compositions as the ST class. Both are dominated by short stature shrubs and have similar species compositions (table 1). Previous work at EML has shown that the vegetation composition shifts from tussock-forming sedges to shrubs as thermokarst develop (Schuur et al 2007). This is similar to the vegetation change that occurs at the transition between tussock tundra valleys and shrubby hilltops. The other obvious misclassification occurred between the MR and HS classes, with 31% of MR training pixels classified as HS (table 3). We believe this confusion occurred because the MR class is a mosaic of vegetation types that include tall shrubs Table 3. Summary of supervised classification accuracy assessment that includes a confusion matrix of the classification and ground validation data, along with the producer's and user's accuracies (%), the overall accuracy, and the Kappa statistic.

Discussion
Twelve per cent of the landscape within the EML watershed has undergone permafrost thaw resulting in thermokarst formation (table 2). Although thermokarst features are widespread throughout the watershed, they predominately occur within the low-lying valleys dominated by tussock tundra ( figure 2). This pattern of ground subsidence is dictated by the initial amount and distribution of ice-rich permafrost, which in turn resulted from the geology and soils of the area. The EML watershed experienced multiple glacial advances and retreats that created a landscape of rolling hills and valleys (Wahrhaftig 1958 Based on our land cover classification dendritic patterns of thermokarst gullies have formed throughout the valleys of the watershed. We estimate approximately 1/3 of the original 3.7 km 2 tundra land cover class has undergone near-surface permafrost thaw that resulted in thermokarst formation. It is important to note that this is likely an underestimate of total permafrost degradation since thawing that does not result in subsidence (such as the deepening of active layers) may not be remotely detected. Thermokarst and tussock tundra classes occur in areas with relatively thick soil organic layers compared to other land cover classes within the EML watershed (figure 3). Although our soil data is semi-quantitative where we did not quantify variability >20 cm, previous intensive work within tundra at EML reports soil organic layers range from 37 to 54 cm (Hicks Pries et al 2012; sites occur within the subset of figure 2). At these same sites, carbon inventories of the top meter of soil ranged from 55 to 69 kg C m −2 , which is on the upper end of the range (16-94 kg C m −2 ) reported for similar tundra soils across Alaska (Michaelson et al 1996). They also found no difference in carbon pools at minimally and extensively subsided tundra sites (Hicks Pries et al 2012). By extrapolating this intensive plot data, we estimate that both tundra and thermokarst classes have organic soil layers that are twice as deep as other classes in the watershed. In spite of the thermal protection provided by the extensive organic soil accumulations in the valleys (Shur and Jorgenson 2007), under present climate conditions, this protection is no longer sufficient to keep the underlying near-surface permafrost from thawing and thermokarst from forming and expanding.
The soil profile of thermokarst classes is seasonally thawing to a greater extent than any other land cover class, and thaw depth in thermokarst is twice as deep as in the surrounding tussock tundra (figure 4). This increase in thaw depth results from water being redistributed into subsided areas, which increases in soil thermal conductivity as organic matter becomes saturated, and increases ground heat flux within depressions while decreasing heat flux in higher, dryer areas ( Yoshikawa andHinzman 2003, Bonan 2002). The net result is deeper thawing in subsided, saturated areas and a decrease in thaw in the high relief areas of the adjacent tundra from which water is drained. This pattern of increased thaw depth within subsided areas has been widely documented at EML , Trucco et al 2012, Lee et al 2010.
The landscape level changes created by permafrost thaw within the EML watershed have important implications for its carbon cycle because thermokarst features are forming in the valleys with soils rich in organic carbon, and are altering the physical environment in ways that increase the vulnerability of the soil carbon pool. Through associated changes in temperature and hydrology, thermokarst causes the active layer of soil to thaw to a greater extent. This exposes more previously frozen organic carbon to above-freezing temperatures and extends the time it remains unfrozen because it takes longer to refreeze the greater volume of unfrozen soil. These abiotic changes can stimulate decomposition and nitrogen mineralization, and result in increased carbon emissions, unless buffered by changes in soil hydrology (Goulden et al 1998, Schuur et al 2008, Shaver et al 1992, Davidson and Janssens 2006 Schuur et al 2007). At EML, shifts in the plant community from a dominance of tussock-forming sedges to more productive shrubs have already occurred in response to permafrost thaw (Schuur et al 2007).
Currently, tundra undergoing permafrost thaw at EML is a carbon source on an annual basis, emitting on average 54 g C m −2 year −1 over the past few years (2008-2010Belshe et al 2012). Although both growing season net carbon uptake and winter carbon emissions have increased in conjunction with ground subsidence at EML , Trucco et al 2012, growing season net sink activity is negated by winter carbon emissions . We believe this current imbalance is created because decomposition can continue to increase as permafrost thaw exposes more of the vast soil organic pool, while the magnitude of potential carbon uptake is mediated by size limitations of species currently residing in the tundra. This trajectory will continue as more soil is exposed with thaw, unless offset in the future, on the timescale of plant succession, by boreal trees advancing from the mixed riparian land cover class into the tundra. Above ground biomass is predicted to increase by 15-68% by the 2050s as Arctic vegetation composition, density and distributions change in response to future warming (Pearson et al 2013). In addition, long-term (since the early 1980s) remote sensing observations show a greening Arctic (Myneni et al 1997, Jai and Epstein 2003, Nemani et al 2003, Goetz et al 2005, Sitch et al 2007, Jai et al 2009 and increased shrub encroachment (Stow et al 2004, Strum et al 2005, both suggesting a regional response of accelerated carbon uptake. An important next step is to build upon this analysis to detect changes in thermokarst extent as well as shifts in vegetation cover with remote sensing image time series. The considerable time and work necessary to collect ground validation data within our relatively small (10.4 km 2 ) study area indicate that mapping of such thermokarst features remains challenging despite their importance for understanding and projecting tundra ecosystem transitions. Clearly, homogeneous very high resolution remote sensing data will be useful to investigate such features in other areas. But available elevation data, covering widespread areas in high latitudes, have horizontal resolution (15-30 m) and vertical accuracy (>5-10 m) that make it difficult to extend this type of analysis to larger areas. With increasing efforts and advances in technology, such as LIDAR, successful mapping of such small, but widespread upland thermokarst features over larger study regions is becoming possible (e.g., Stevens and Wolfe 2012).

Conclusions
Using a supervised maximum likelihood classification approach and a very high resolution IKONOS image, we detected and mapped small, irregular thermokarst features occurring within the EML watershed (figure 2). Overall we achieved an acceptable classification accuracy of 83% when compared to ground validation data, with thermokarst classes having comparable accuracy to other classes in the landscape (table 3). This demonstrates the feasibility of using high resolution remote sensing data to detect a diversity of thermokarst features, including small, subtle features that are more difficult to detect. Our analysis of this site, typical for uplands with discontinuous and partially ice-rich permafrost in Interior Alaska, shows that 12% of the landscape at EML has undergone subsidence, and thermokarst-induced physical changes have important implications for the carbon balance of this ecosystem.