Quantification of upland thermokarst features with high resolution remote sensing

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⁻¹ 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 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 and Osterkamp 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, Belshe et al 2012, 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 (Morgenstern et al 2011, Jones et al 2011), drained lake basins (Hinkel et al 2003, Wang et al 2012), retrogressive thaw slumps (Lantuit and Pollard 2008, Lantuit et al 2012), and thermo-erosion gullies (Godin and Fortier 2012). 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.

2.1. Site description

2.2. 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.

Zoom In Zoom Out Reset image size

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.

Table 1.  Characterization of land cover classes and the number of field training sites used for the supervised classification.

Land cover class	Training sites	Description
Hilltop shrub (HS)	10	Hill top communities dominated (85%) by medium stature (60–120 cm) shrubs (Betula nana, Betula glandulosa). Similar to DLMP land cover class closed low shrub bircha. Similar to the CAVM land cover class S2 low shrub tundrab.
Shrub tundra transition (ST)	10	Intermediate transition zone between hilltops and valleys with intermediate stature (25–55 cm) vegetation that is a mixture of shrubs (63%; Betula nana, Rubus chamaemorus, Rhododendron groenlandicum, Empetrum nigrum, Vaccinium vitis-idaea, Vaccinium uliginosum), graminoids (16%; Eriophorum vaginatum, Carex bigelowii) and non-vascular (17%; Dicranum spp., Sphagnum spp., Feather-moss, Lichen) plants. Similar to the CAVM land cover class S2 low shrub tundra and G4 tussock-sedge, dwarf-shrub, moss tundrab.
Tussock tundra (TT)	10	Mixture of short stature (<25 cm) graminoid (35%; Eriophorum vaginatum, Carex bigelowii), shrub (38%; Betula nana, Rubus chamaemorus, Rhododendron groenlandicum, Empetrum nigrum, Vaccinium vitis-idaea, Vaccinium uliginosum) and non-vascular (25%; Dicranum spp., Sphagnum spp., Feather-moss, Lichen) vegetation occurring in gently sloping valleys. Similar to DLMP land cover class low shrub sedgea. Similar to the CAVM land cover class G4 tussock-sedge, dwarf-shrub, moss tundrab.
Thermokarst pool (TP)	4	Pools of standing water created by permafrost degradation in areas of low slope. Pools are generally around 5 m wide and less than 50 cm deep.
Thermokarst (TK)	22	Thermokarst gullies that have three types of vegetation depending on the local slope and drainage. Relatively low slope areas are dominated by either non-vascular species (60%; Sphagnum spp.) or graminoid species (67%; Eriophorum spp, Carex bigelowii), while areas with greater slope and drainage are dominated by a mixture of intermediate stature (20–63 cm) shrub species (61%; Betula nana, Rubus chamaemorus, Rhododendron groenlandicum, Empetrum nigrum, Vaccinium vitis-idaea, Vaccinium uliginosum). Similar to the CAVM land cover class S2 low shrub tundra in high slope areas and W4 sedge, moss, low shrub wetland in low slope areasb.
Mixed riparian (MR)	18	Mixed vegetation near the inlet and outlet of the lake, and alongside the road. Vegetation is a matrix of patches of needle-leaf trees (Picea mariana) and large shrubs (200–600 cm; either Salix spp. or Alnus spp.). Similar to DLMP land cover classes Stunted Spruce, Alder, and Willowa. Similar to the CAVM land cover class S2 low shrub tundrab.
Road (RD)	4	Stampede Road. Similar to DLMP land cover class bare grounda.
Lake (LA)	NA	Eight Mile Lake. Similar to DLMP land cover class clear watera.

^aDenali Landcover Mapping Project (DLMP); Stephens et al (2001). ^bCircumpolar Arctic Vegetation Map (CAVM); Walker et al (2005).

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.

2.3. Data preparation

2.4. Supervised classification

2.5. Classification accuracy assessment

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², which is 12.1% of the 10.4 km² study area (table 2). Tussock tundra (TT) covers 2.4 km² (22.9% of the total area), hilltop shrub (HS) covers 1.7 km² (15.8%), shrub tundra transition (ST) covers 2.3 km² (22%), and the mixed riparian (MR) covers 2.2 km² (21%; table 2). Assuming that the thermokarst areas originally were all covered with tussock tundra, we estimate that the unthawed tussock tundra originally covered 3.7 km² and by 2008 approximately 1/3 (34.6%) had undergone permafrost thaw.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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 that are spectrally similar to the medium stature shrubs occurring on the hilltops. Additionally, the classification was not well resolved in the southern portion of the watershed and overestimates the homogeneity and dominance of the MR class at the southern tip. This incongruity was caused by the low number of training sites in this steeply sloped region. However, we achieved an overall acceptable accuracy, with thermokarst classes having comparable accuracy to other classes in the landscape.

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). The hills have gravelly soils consisting of glacial till covered by a thin aeolian silt cap with low amounts of ice within the permafrost. Alternatively, the lower slopes and valleys have fine-grained soils associated with colluvial deposits and ice-rich horizons containing up to 50% visible ice content (Osterkamp et al 2009). These ice layers and lenses formed over time as peat accumulated in the valleys and water was trapped by upward freezing of the permafrost surface (Osterkamp et al 2009). Although the peat accumulations initially acted to build and protect the permafrost, the concurrent increase in ice ultimately made permafrost in the valleys vulnerable to thermokarst formation. Through continued warming, insulation of soil by snow accumulating in subsided areas, and thermal erosion (from water being redistributed from higher to lower microtopographical areas) permafrost thaw has increased and thermokarst gullies and pools have expanded (Kane et al 2001, Osterkamp et al 2009). 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² 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⁻², which is on the upper end of the range (16–94 kg C m⁻²) 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 (Jorgenson and Osterkamp 2005, Osterkamp 2005, Kane et al 2001, Jorgenson et al 2001, Osterkamp et al 2009, Yoshikawa and Hinzman 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 (Vogel et al 2009, Belshe et al 2012, 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). On the other hand, thermokarst formation can lead to increased carbon uptake if decomposition increases plant available nitrogen (Hobbie et al 2000, Chapin and Shaver 1996, Mack et al 2004), or the plant community shifts to more productive species (Osterkamp et al 2009, 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⁻² year⁻¹ over the past few years (2008–2010; Belshe et al 2012). Although both growing season net carbon uptake and winter carbon emissions have increased in conjunction with ground subsidence at EML (Vogel et al 2009, Trucco et al 2012), growing season net sink activity is negated by winter carbon emissions (Belshe et al 2012). 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²) 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).

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.

We would like to thank Larissa Yocum from Denali National Park and Preserve for helping us locate spatial data products. We would also like to thank UNAVCO for providing GPS equipment, training, and support, which is supported by the National Science Foundation and National Aeronautics and Space Administration under NSF Cooperative Agreement No. EAG-0735156. This study was supported by grants to E Schuur: NSF grants and the DOE NICCR program. The 2009 summer field campaign was supported by a Discover Denali Fellowship awarded to E Fay Belshe, who was also supported by the Department of Energy, Graduate Research Environmental Fellowship.