Extrapolating active layer thickness measurements across Arctic polygonal terrain using LiDAR and NDVI data sets.

Landscape attributes that vary with microtopography, such as active layer thickness (ALT), are labor intensive and difficult to document effectively through in situ methods at kilometer spatial extents, thus rendering remotely sensed methods desirable. Spatially explicit estimates of ALT can provide critically needed data for parameterization, initialization, and evaluation of Arctic terrestrial models. In this work, we demonstrate a new approach using high-resolution remotely sensed data for estimating centimeter-scale ALT in a 5 km2 area of ice-wedge polygon terrain in Barrow, Alaska. We use a simple regression-based, machine learning data-fusion algorithm that uses topographic and spectral metrics derived from multisensor data (LiDAR and WorldView-2) to estimate ALT (2 m spatial resolution) across the study area. Comparison of the ALT estimates with ground-based measurements, indicates the accuracy (r2 = 0.76, RMSE ±4.4 cm) of the approach. While it is generally accepted that broad climatic variability associated with increasing air temperature will govern the regional averages of ALT, consistent with prior studies, our findings using high-resolution LiDAR and WorldView-2 data, show that smaller-scale variability in ALT is controlled by local eco-hydro-geomorphic factors. This work demonstrates a path forward for mapping ALT at high spatial resolution and across sufficiently large regions for improved understanding and predictions of coupled dynamics among permafrost, hydrology, and land-surface processes from readily available remote sensing data.

First effort to map the ALT using fine resolution remotely sensed data

A blended methodology incorporating RS data and statistical manipulation

Smaller-scale ALT is controlled by eco-hydro-geo variables

1. Introduction

Most general circulation models predict that warming observed in the Arctic over the past decades [Hassol et al., 2004] will continue and is likely to accelerate through the 21st century [IPCC, 2013]. This warming will be accompanied by a thawing of permafrost and an increase in the thickness of the active layer [Koven et al., 2013; Slater and Lawrence, 2013; Streletskiy et al., 2012; Lawrence et al., 2008], which is the layer of soil that experiences seasonal freeze and thaw. Increased active layer thickness (ALT) will be accompanied by changes in hydrology (inundation patterns and soil moisture) [Hinzman et al., 2005], vegetation [Chapin et al., 1988], surface energy balance [Porter et al., 2011], and geomorphology (thermal erosion, landslides, ground subsidence, and possibly channel development) [Gooseff et al., 2009; Rowland et al., 2010; Gangodagamage et al., 2007]. Most geomorphic, hydrological, soil microbial, and biogeochemical activities are limited to the active layer [Brown et al., 1980; Hinzman and Kane, 1991], which may vary from 30 to substantially higher than 70 cm, depending on geographic location [Hinkel and Nelson, 2003; Zhang et al., 2005; Shiklomanov et al., 2010]. The potential release of carbon into the atmosphere as the active layer thickens is of great concern; it is estimated that 1672 Pg of carbon are currently stored in permafrost, with the vast majority in the upper 3 m [Tarnocai et al., 2009]. However, in current climate-model predictions of ALT response to future warming, the permafrost landscape is represented at a relatively coarse resolution, and microtopographic variations are often ignored [Lu and Zhuang, 2012]. Methods of quantifying ALTs at fine spatial resolution across large spatial scales, therefore, have the potential to refine assessments of the rate and timing of pan-Arctic soil carbon releases.

ALT depends on many climatological parameters, such as air temperature, surface temperature, and snow thickness [Smith, 1975], as well as on land-surface properties such as surface morphology (slope, aspect) [Nelson et al., 1997], vegetation cover, and the physical, hydrological, and thermal properties of the surface cover and subsurface [Brown et al., 2000; Zhang et al., 2005]. ALT is usually directly measured by inserting a metal probe through the seasonally thawed soil to resistance at the top of permafrost, or by collecting ground temperatures from vertical arrays of temperature sensors [Brown et al., 2000]. Thaw tubes are another common method used in Canada to measure ALT [Mackay, 1973]; they primarily measure maximum annual thaw penetration and maximum annual heave and ground subsidence [Rickard and Brown, 1972; Mackay, 1973; Nixon et al., 1995]. Thaw tubes track the maximum annual thaw depth as the maximum depth penetration of a colored marker along a tube of water inserted into the ground. During the thaw season, the marker sinks through water to the depth at the top of the permafrost table where the water in tube remains frozen [Nixon et al., 1995].

Thaw depths may also be estimated using analytical procedures such as a specific form of the Stefan solution [Harlan and Nixon, 1978]. Using this method the seasonal thaw depths (X) are linked to the accumulated surface thawing-degree days (ATDD), which is the daily average of the surface temperatures from the beginning of thaw onset to the date of active layer measurement: where X is thaw depth (m), k is the unfrozen thermal conductivity of soil (Wm−1 K−1), ATDD is surface thawing-degree day (°C), and L is the volumetric latent heat of fusion (Jm−3). Here, thaw depth is a function of the accumulated seasonal thawing temperature, along with substrate moisture content and thermal properties. ATDD at a site is usually derived from accumulated air thawing-degree day (ATDD) and the thawing “n factor” (n), which are usually related to ATDD by an empirical equation given as ATDD [Lanardini, 1978]. This analytical technique is extensively used to cross validate the ALTs measured from probing or thaw tubes [Shiklomanov et al., 2010; Smith et al., 2009]. Here, n is determined using field data [Taylor, 2000]. A strong correlation between measured landscape-specific ALT and √(ATDD) was observed for the CALM (Circumpolar Active Layer Monitoring network) site at Barrow, Alaska [Shiklomanov et al., 2010].

While these methods used to measure ALT are accurate, they are labor intensive and limited in spatial coverage. At regional scales, attempts have been made to model ALT by extrapolating ground measurements using empirical and statistical relationships with air temperature, ground temperature, elevation, and surface vegetation at coarser (30–300 m) spatial resolution [Nelson et al., 1997, 1999; Zhang et al., 1997; Shiklomanov and Nelson, 1999]. Working at the same site as this study, Hubbard et al. [2013] showed that surface ground-penetrating-radar data could be used to provide very high resolution and accurate estimates of ALT. They also found that ALT, as well as other geochemical, geophysical, hydrogeological and thermal properties, varied spatially as a function of polygon type, or “geomorphic zonation.” Most recently, Liu et al. [2012] estimated ALTs across an 80 by 100 km2 area of the North Slope of Alaska, at a 100 m resolution using InSAR. They related InSAR-derived surface subsidence-to-volumetric changes associated with freeze and thaw of a saturated soil column. This method assumed the active layer was close to fully saturated, and it required some knowledge of soil and land-cover characteristics for estimating water content of the soils.

While current methods are very good at quantifying and predicting ALT at either local or regional scales, no method has been presented that is able to capture both the fine scale variation in ALT related to the microtopography of rims, troughs, and centers in ice wedge polygon landscapes, and the influence of microtopography on ALT over large spatial domains. Since ALT can vary by tens of centimeters over small spatial scales, improved representation of ALT variation has the potential to improve regional scale predictions of hydrologic response and carbon fluxes in ice wedge polygon landscapes [Minke et al., 2009; Zona et al., 2011; Hubbard et al., 2013].

Here, we present a new approach for estimating ALT at high spatial resolution (2 m), using remotely sensed data sets and ground-truth data (ALT by probing). We employed a multisensor data-fusion approach, using statistical learning from a simple, multivariate, piecewise regression technique [Berk, 2008], by combining in situ ALT measurements with (a) data sets from LiDAR-derived high-resolution surface elevation (0.25 m) to obtain local slope and landscape curvature, and (b) high resolution (2 m) Worldview-2 (DigitalGlobe, Inc.) data to obtain the Normalized Difference Vegetation Index (NDVI). We used a statistical machine-learning method to relate the surface data to the direct ALT measurements. The proposed method offers the potential to provide ALT estimates at high spatial resolution and across sufficiently large regions for improved understanding and predictions of coupled dynamics among permafrost, hydrology, and land-surface processes.

2. Study Area and Data Collection

The study site is located within the Barrow Environmental Observatory, BEO, near the city of Barrow, Alaska (71.3°N, 156.5°W), which is approximately 6 km southwest of the northernmost point of Alaska and is on the Alaskan arctic coastal plain (Figure 1). In this paper, we will refer to the study site as “Barrow” even though the BEO is located to the southeast of the city of Barrow, AK. At Barrow, the mean annual adjusted precipitation is about 173 mm yr−1, and mean annual air temperature is −12°C, [Liljedahl et al., 2011], with a minimum mean monthly temperature of −26.6°C in February and maximum mean monthly temperature of 4.7°C in July [NCDC, 2002]. Permafrost underlies the entire area and extends to depths of 300 m or greater [Brewer, 1958], with ALT varying from 20 to 70 cm [Shiklomanov et al., 2010]. The depth at which the zero annual amplitude of soil temperature (permafrost temperature) in Barrow is approximately 15 m [Romanovsky et al., 2002]. Permafrost temperature has increased by 1°C in Barrow from 1950 to 2001. However, no associated trend in ALT has been observed in Barrow, Alaska [Romanovsky et al., 2002] until 2005. The most recent active layer data collected at two CALM sites at Barrow, Alaska (see the data for site code U1 at http://www.gwu.edu/∼calm/data/webforms/u1_f.htm and site code U2 at http://www.gwu.edu/∼calm/data/webforms/u2_f.htm) show progressive increase in ALT since 2005.

Figure 1

Study area enclosed within the Barrow Environmental Observatory (BEO) (white), the LiDAR data extent (black), the geophysical transect described by Hubbard et al. [2013] (yellow), and the locations of NGEE extensive field experimental site B (blue) and C (red).

Vegetation patterns at Barrow vary with ice-wedge polygon topography and surface hydrology [Britton, 1957; Bliss and Peterson, 1992]. A variety of ice-wedge polygon morphologies exists in the area [Brown et al., 1980]. Here we focused on a 5 km2 area encompassing low centered (LC) and high centered (HC) polygons (Figures 2 and 3); this is the same region described by Hubbard et al. [2013]. Within our study area, there is a dramatic contrast in vegetation type across polygon rims, troughs, and centers. Mosses, lichen, and dwarf shrubs mainly cover elevated areas of HC polygon centers and LC polygon rims, while mosses and sedges with standing water are found in depressions of troughs and LC polygon centers [Webber, 1978; Villarreal et al., 2012]. These land-cover types can alter the surface energy balance and thermal properties in microtopographic regions, through changes in albedo, surface roughness, and evapotranspiration [Langer et al., 2010; Villarreal et al., 2012]. Low landscape relief (<7 m; Figure 2) and low hydraulic gradient characterize this landscape.

Figure 2

A schematic of how directed distance (DD) is used as a scale parameter for a high centered and a low centered polygon. DD in (a) high centered (HC) polygon and (b) the corresponding general elevation profile for a high centered (HC) polygon. DD in (c) a low centered (LC) polygon and (d) the corresponding general elevation profile for a low centered (LC) polygon.

Figure 3

(a) A shaded relief map of high centered (HC) and low centered (LC) polygons. Arrow lines indicate an example for a typical high centered polygon and a low centered polygon. (b) A map of Zone 1, Zone 2, and Zone 3 are shown. White color polygons show the delineated polygon boundaries from the LiDAR data.

The soil properties of the coastal tundra at Barrow site are products of both climate factors and the specific geological history of the area. Data gathered from nine sites during the international tundra biomes effort have shown that the soil at Barrow is similar to other tundra areas [French, 1974]. Barrow soil types are generally classified as Gelisols, which primarily includes an organic-rich soil layer underlined by a cryoturbated mineral soil (silty clay to silty loam-textured) layer, and a frozen organic-rich mineral layer [Bockheim et al., 1999; Drew, 1957].

Soil formation at Barrow is associated with the seasonal freeze and thaw processes in ice-wedge polygons at low annual temperatures, high relative moisture conditions, and low mean annual precipitations under impeded drainage conditions [Gersper et al., 1980]. Based on climate and soil factors, Rosswall and Heal [1975] documented that coastal tundra soil types are highly correlated with the microtopographic units. Their classification indicated the occurrence of uniform soil classes within meadows, polygon troughs, rims, center areas of low centered polygons, and centers of high centered polygons. The texture of the mineral soil is relatively uniform across all microtopographic zonations [Hubbard et al., 2013].

A unique characteristic of the tundra soils is the high amount of organic matter, which influences the fluxes of soil, heat, gases (like oxygen, CO2, and methane), and cation exchange [Babb and Whitfield, 1977; Chapin and Van Cleve, 1978]. One-third of this organic matter is in the upper 10 cm of soil [Gersper et al., 1980]. Variability in soil organic carbon (SOC) has also been documented in soils across ice-wedge polygon landscapes at Barrow. SOC varied from 24 to 28 kg m−3 in troughs (both in HC and LC polygons) whereas these values in polygon rims and centers have been documented at 32 and 84 kg m−3, respectively [Bockheim et al., 1999].

For this study, in situ thaw depth data were collected at three areas using the mechanical probing method (inserting a metal probe through the active layer to resistance at the top of the permafrost) in conjunction with a detailed laser level and differential GPS survey as part of the Department of Energy Next-Generation Ecosystem Experiments Arctic (DOE NGEE-Arctic) project. The survey area was (1) a 475 m long transect with one thaw-depth measurement every 3 m, three detailed short transects measurements at 0.5 m interval along the same transect (Figure 1), and two 475 m long 10 m offset transect measurements every 25 m (a total of 300 points, uniform sampling; 24 September to 1 October 2011); and (2) two areas (Site B and Site C) with multiple shorter transects (<10 m) each having one measurement every 0.5 m (108 points; Figure 1). The direction (NW-SE) of the 475 m transect was intended to capture a variety of polygon types, with the southern portion encompassing HC polygons and the northern portion LC polygons. A more spatially intensive sampling effort along the 475 m transect was performed at sites DS1, DS2, and DS3 (Figure 1; 0.5 m intervals, i.e., a total of 99 point measurements). The thaw-depth measurements were considered equivalent to ALT due to the timing of early fall field campaigns.

High-resolution LiDAR data (0.25 m spatial resolution) were collected on 4 October 2005. The LiDAR data's horizontal and vertical accuracy were approximately 0.30 and 0.15 m, respectively. This LiDAR data were quality-checked in late August 2006 by Siebert and Associates (http://seibert-ak.com/gps); 284 data points were used for this purpose. This survey indicated a systematic error of −0.143 m, which was then used to correct the LiDAR point cloud data.

3. Methodology

3.1. Creating a Spatial Database for Modeling

Four explanatory variables derived from remotely sensed data were used to predict ALT: (1) directed distance (DD, m), (2) local slope (S, m/m), (3) curvature derived from wavelet transformation (wavelet curvature; WC, 1/m), and (4) NDVI. The DD, S, and WC were used to parameterize the geomorphic variability associated with ALT. NDVI was used to characterize the vegetation cover, standing dead vegetation, and soil moisture regimes of the polygonal landscape.

Gangodagamage et al. [2011] first introduced DD in the analysis of river basins, where they related the ensemble average-basin attributes (e.g., slope, curvature) to the distance of any given pixel from the nearest drainage divide. In areas of HC polygons, DD was calculated from the elevated polygon center to the trough along the downstream surface flow paths (Figures 2a and 2b). The DD associated with LC polygons was calculated starting at the elevated rims along the downward surface flow path, i.e., from the polygon rims to the surrounding trough directly or through the polygon centers (Figures 2c and 2d). We used DD to define the topographic position of the ALT field measurements. Further, DD is used to separate topographic attributes into categories defined by values of slope, wavelet curvature, and NDVI, which are then analyzed for trends and relationships that each variable makes with measured values of ALT.

For all ALT measurements, we computed the DD along the microtopographic flow paths. We then computed the ensemble mean ALT (<ALT>) as a function of DD (Figure 2). Linear relationships between ALT field measurements and the morphological characteristics of the polygons defined three microtopographic zones (Figure 4). In LC polygons, Zone 1 (0 ≤ DD ≤ 7) encompassed the polygon rim, slope inflection point, and lower slope of the rim within the interior of the polygons (Figure 3b red area). For HC polygons, Zone 1 included the polygon center (flow characteristics are similar to polygon rim points of LC polygons), slope inflection points, and lower slope areas at the polygon-trough boundary (Figures 3a and 3b). Zone 2 (7 < DD ≤ 15) in LC polygons included the low interior portion of the polygon and typically extended to the lowest points at the polygon boundary along the downstream DD flow paths (Figures 3a and 3b). HC polygon morphology was typically different from LC polygons, and Zone 2 was not represented in HC polygons. For all polygon types, Zone 3 (DD > 15) represented the troughs.

Figure 4

The relationship between ensemble averaged ALT () and directed distance measured along the flow paths for both low centered and high centered polygons. ALT values were averaged as a function of directed distance measured along the high centered and low centered flow paths.

The relationship between ensemble averaged ALT (<ALT>) and directed distance measured along the flow paths for both low centered and high centered polygons. ALT values were averaged as a function of directed distance measured along the high centered and low centered flow paths.

The local slope was computed by considering a region composed of nine (0.75 × 0.75 m2) grid cells, whereas the curvature was computed by convolving a second derivative of the Gaussian wavelet (also called Mexican Hat Wavelet) with the elevation data set [Lashermes et al., 2007]. Imagery used for the calculation of NDVI was from the WorldView-2 (DigitalGlobe, Inc.) satellite with eight multispectral bands (coastal, Band 1, 400–450 nm; blue, Band 2, 450–510 nm; green, Band 3, 510–580 nm; yellow, Band 4, 585–625 nm; red, Band 5, 630–690 nm; red edge, Band 6, 705–745 nm; near infrared 1 (NIR 1), Band 7, 770–895 nm; near infrared 2 (NIR 2), Band 8, 860–1040 nm). We used ENVI version 4.8 (Exelis Visual Information Solutions, Boulder, Colorado) to orthorectify the WorldView-2 images using RPC (rational polynomial coefficients), LiDAR elevation data, and ground control points (GPC are collected by a GPS survey conducted by Aerometrics Inc.) to authorectify the WorldView-2 data. The pixel digital numbers (DN) of the data were converted to spectral radiances for better quantitative representation of spectral responses. The DN values were initially converted to top-of-atmosphere spectral radiance (Wm−2 sr−1 μm−1), using absolute radiometric calibration factors (Wm−2 sr−1 count−1) given for each band to get a band integrated radiance (Wm−2 sr−1), which were then normalized by effective bandwidth to get the spectral radiance [Updike and Comp, 2010]. We used the red band (Band 5, 630–690 nm) and near-infrared band (Band 7, 770–895 nm) to compute the NDVI (spatial resolution 2 m) values from an image acquired on 21 July 2010.

Thaw depths during early fall were collected using a T-handled tile probe labeled in centimeter gradations from the bottom of the probe to the top. The probe used was 1.5 m long and 1 cm in diameter with a tapered tip. The probe was strong enough to resist the bending stress generated by the resistance of the substrate material in the active layer (for ALT < 1 m). No rocks or other hard materials were found in the substrate, and the substrate soils were so wet as it is rich in ice content in this Barrow site. The probe was pushed into the soil until the downward movement ceased when the tip reached the permafrost table. Then, the tip of the probe was raised a few centimeters and pushed downward to confirm the position of the permafrost table. The depth was recorded at this position to the nearest centimeter. A TOPCON GRS-1 (Geodetic Rover System) real-time kinematic global positioning system (RTK GPS) was used to measure the latitude, longitude, and elevation at all ALT measurements locations. All positional data collection showed subcentimeter accuracy in both vertical and horizontal directions [Hubbard et al., 2013].

3.2. Relative Accuracy of LiDAR and LiDAR-Derived Products

The increasing availability of high-resolution (1 m) topography data and enhanced computational processing power presents new opportunities to study landscape processes quantitatively, in detail not possible before. Given this new power, it is now of utmost importance to understand how these LiDAR and LiDAR-derived metrics (such as slope, curvature), as well as high-resolution satellite-data-derived products (such as normalized difference vegetation index—NDVI) influence high-resolution model predictions. In the following section, we review the absolute and relative accuracies of the LiDAR data and investigate how it influences the accuracies in LiDAR-derived products such as DD, NDD, slope, and curvature.

3.2.1. Absolute Accuracy

Accuracy of a DEM (Digital Elevation Model) derived from LiDAR depends on many factors: LiDAR sensor (pulse length and rate, wavelength, and scan angle), aircraft (altitude and speed), navigation system (GPS system), LiDAR point processing and labeling, and geographic location (seasonality, land-cover, and terrain slope) [ASPRS, 2004; Hodgson et al., 2005; Shrestha et al., 1999]. Vertical accuracy is the principal criterion in determining the quality of the elevation data, with the accuracy requirement primarily depending upon the respective user applications [Hodgson and Bresnahan, 2004]. At the preliminary data-calibration stage, errors associated with yaw, pitch, and role of the aircraft are assessed and corrected by the LiDAR data provider. The vertical absolute accuracy of a data set is then derived by computing the root-mean-square error (RMSE) of the elevation data:

Here, n is the number of points used for the accuracy assessment and RMSE is a summation of LiDAR system error, error in DPGS survey, error associated with different topographic slopes, and error due to possible frozen state of landscapes, especially in tundra (e.g., frozen water in troughs). Horizontal accuracy is also an important characteristic of LiDAR elevation data [Hodgson and Bresnahan, 2004]. Horizontal inaccuracy is usually larger than vertical inaccuracy (∼two times that of vertical accuracy) and highly controlled by the vertical accuracy, in which very high vertical accuracies essentially have to be accompanied by very high horizontal accuracy. This is primarily because horizontal inaccuracies contribute significantly to computations of vertical accuracies [ASPRS, 2004; Hodgson and Bresnahan, 2004].

3.2.2. Relative Error

Because the surface morphometry represented by LiDAR data is internally consistent (spatially autocorrelated) [Gangodagamage et al., 2008], the DEM-derived products such as slope, curvature, and aspect are much more accurate than the absolute accuracies represented by the LiDAR elevation data. For example, the topographic slope derived from LiDAR data is partially independent from the absolute elevations, because the derivative of the elevation is dependent on the internal consistency of the data points. The point-to-point or relative vertical error is more important for derivative products (e.g., slope—first-order derivatives; curvature—second-order derivative) than the absolute vertical accuracy, and it only utilizes local differences among the adjacent elevation values of LiDAR elevation data [ASPRS, 2004].

We used the 284 DGPS elevation points and their corresponding LiDAR elevation to assess the relative error in elevation values of LiDAR DEM. We then used this vertical relative error to compute relative errors in LiDAR-derived products such as slopes, curvature, DD, and NDD. We computed the point-wise relative differences in elevations for the DGPS points and point-wise relative differences for LiDAR data separately. First, we considered LiDAR elevation points (284 points). The first-order differences from each LiDAR point (Δh(i)) to all other LiDAR points were computed, and then the same approach was repeated to compute the first-order differences in DGPS points (Δh(i)) as well. This procedure was carried out for all 284 points (i = 1:284), and their differences were computed as follows:

Here, Δh(i = 1:n) are the relative elevation differences in LiDAR and DGPS points where i varies from point 1 to 284, h(i) is the LiDAR elevation at point i, and h(j) is the LiDAR elevation at point j, with j varying from 1 to 284, and j always equal to or larger than i. Here, h(i) is the DGPS elevation at point i, h(j) is the DGPS elevation at point j, and n (n = 284) is the total number of points used for the analysis. Then, we computed the average relative error in elevation (δ) by taking the mean value of the elevation difference associated with each point as:

In relative error computation, we considered only vertical elevation differences between points; no horizontal distances between points were considered. We assumed that the relative horizontal error (δ) in LiDAR data is two times that of the vertical error, δ = 2 × δ, while the absolute horizontal errors of LiDAR are approximately equal or less than two times of the vertical errors [ASPRS, 2004]. We assumed that the same error ratio exists between relative vertical errors and relative horizontal errors. We considered slope (s) as the ratio of rise (height in vertical axis, a) over run (distance in horizontal axis, b). The accuracy of the slope is determined by both vertical accuracies and horizontal accuracies. The relative error in vertical elevations is denoted as δ; the relative error in horizontal distances δ. The error in slope is considered as δ. Estimated uncertainty in the vertical direction is a ± δ, estimated uncertainty in the horizontal direction is b ± δ, and estimated uncertainty in slope is s ± δ. Error in the uncertainty of slope (δ) can be expressed as (see Appendix A for details, equation (A4)):

The same derivation can be used to derive the relative error in curvature, by considering curvature as the rate of change in slope for each pixel. Change in slope for a pixel is considered as c, and relative error associated with curvature is assumed as δ, which can be expressed as:

Finally, it can be shown that the relative error of the DD, δ, is equivalent to the relative horizontal error, δ. The relative error of the NDD, δ, is equivalent to 2δ (Appendix A).

3.3. Regression Model Development

We developed a simple piecewise multivariate-linear-regression model to predict the response variable ALT using normalized DD (NDD), S, WC, and NDVI as where are the piecewise regression coefficients, and is the error term obtained for i = 1:3 cases that were derived from the three zones defined by directed distance. NDD is defined as where DD is the flow path distance of a pixel at Zone i and i represents i = 1, 2, or 3; min(DD) is the minimum DD value in Zone i; and max(DD) is the maximum value of the DD in Zone i.

We used a uniform sampling approach to obtain 407 ALT measurements (Table1), employing a 25 m sampling interval to collect ALT points (42 measurements) along two outside transects (Figure 1) and a 3 m sampling interval for the ALT measurements along the center transect (158 measurements; Figure 1). A 0.5 m sampling interval was used within detailed sampling sites (i.e., DS1, DS2, and DS3; 99 measurements) and Site B and Site C (108 measurements; Table1).

Table 1

The Distribution of ALT Point Measurements Over Different Microtopographic Zones Is Given

	Total Points	Main Transect	Side Transects	Main Transect: Detailed Sites DS1, DS2, and DS3	Site B and Site C
Sampling interval (m)		3	25	0.5	0.5
Total points	407	158	42	99	108
Zone 1 (Rims (LC) + elevated centers (HC))	145	48	14	36	47
Zone 2 (Centers (LC))	98	42	10	27	19
Zone 3 (Trough LC/HC)	165	68	18	37	42

A total of 407 ALT measurements were used for model development (Figure 1): 262 measurements were used in model training and the remaining 145 measurements were used for model validation. Among these ALT measurements, 89, 66, and 108 points were used for model training and 56, 32, and 57 points were used for model validation in Zone 1, Zone 2, and Zone 3, respectively. The distribution of microtopographic types in each zone are given in Table2.

Table 2

Distribution of ALT Points Chosen from Three Different Microtopography Classes for Training and Validation in Piecewise Regression Model Is Given

	Training				Validation
Descriptions	Total Number of ALT	Main Transect	Main Transect (DS1, DS2, and DS3)	Site B and Site C	Total Number of ALT	Side Transects	Main Transect (DS1, DS2, and DS3)	Site B and Site C
Total points	262	158	50	54	145	42	49	54
Zone 1 (Rims (LC) + elevated centers (HC))	89	48	18	23	56	14	18	24
Zone 2 (Centers (LC))	66	42	14	10	32	10	13	9
Zone 3 (Trough LC/HC)	108	68	19	21	57	18	18	21

Model uncertainty was evaluated using a coefficient of determination (r2), and model error was evaluated using analyses of residuals for Mean Average Error (MAE), Mean Bias Error (MBE), Root Mean Squared Error (RMSE), and Willmott's index of agreement (d) [Willmott, ,; Willmott and Matsuura, 2005]. Systematic and unsystematic RMSE (i.e., RMSEs and RMSEus) were also evaluated, both for training and validation data sets separately.

Relative uncertainties of predicted ALTs at each pixel location were also computed as ratio between uncertainties in ALT and absolute predicted ALT value at that pixel location. Uncertainty therein was estimated as deviation from the corresponding mean values considering their microtographic zonations: where ALT is the value of the predicted ALT in Zone i, Here i = 1, 2, or 3 and <ALT> is the mean ALT value in Zone i.

We developed three regression models representing three microtopographic classes: Zone 1, Zone 2, and Zone 3. For each multiple regression, ALT samples were chosen to represent Zone 1, Zone 2, and Zone 3, such that sample sizes satisfy the likelihood (probability) of obtaining a significant (p < 0.05) r2. To satisfy the above statistical power criteria in the regression, the minimum sample size was 50 [Cohen and Cohen, 1983; Cohen, 1988; Cohen et al., 2003; Thompson, 2012]. The sample size we used in each of three regression models was larger than 66 (Table2), which was larger than the preferred sample size as defined in the above statistical power relationship [Cohen and Cohen, 1983; Cohen, 1988]. The sample sizes used in our study also were cross-checked with the conventional tools generally used for quantifying the adequacy of sample size in regression model development. The rule proposed by Stevens [2002] and Park and Dudycha [1974] suggested using 15 measurements for each explanatory variable (ratio of 15:1). We used four variables and the preferred sample size derived from this method (60 samples). The minimum sample size in this study (i.e., 66) was larger than the preferred sample size (i.e., 60) obtained from the traditional method.

4. Results

4.1. Relationships Between ALT and Surface Parameters

A shaded relief map of HC and LC polygons is shown in Figure 3a, and a classified LiDAR elevation map for three micropotopographic regions (i.e., Zone 1, Zone 2, and Zone 3) is shown in Figure 3b [Gangodagamage and Wullschleger 2014aa, 2014bb]. Two peaks in ALT were observed when ALT values were plotted against directed distance (Figure 4). The first peak occurred at a DD of approximately 3 m, associated with the locations of troughs. The second and largest ALT peak was located at a DD of approximately 7 m, which corresponds to the distance between the rim and the LC polygon center (Figure 2). The three zones we used for the regression model development are shown in Figure 4.

We observed a systematic shift in the relationship between DD and ALT in the three zones (Figure 5). Accordingly, regression relationships were developed for each zone. In Zone 1, which represents the elevated polygon centers (HC polygons) and rims (LC polygons), ALT has a positive correlation with DD with the thaw depths increasing as one moves away from the elevated mounds (center or rim; Figure 5a). A decreasing trend in ALT was observed in Zone 2, i.e., a reduced ALT as one moves from the LC polygon center toward the trough (Figure 5b). The ALT increases with DD in Zone 3, suggesting a deepening of the active layer in the troughs compared to the edges of the HC polygon mound and the elevated LC polygon rims. However, the rate was an order of magnitude lower in Zone 3 than observed in Zone 1 (Figure 5c).

Figure 5

Ensemble averaged ALT (±1 std; ) distribution as a function of DD in (a) Zone 1 (0 ≤ DD ≤ 7), (b) Zone 2 (7 < DD ≤ 15), and (c) Zone 3 (DD > 15) regions.

Ensemble averaged ALT (±1 std; <ALT>) distribution as a function of DD in (a) Zone 1 (0 ≤ DD ≤ 7), (b) Zone 2 (7 < DD ≤ 15), and (c) Zone 3 (DD > 15) regions.

Figures 6a–6c show the ALT as a function of local slope. Across all three zones, there was a negative correlation between slope and ALT, and this trend increased from Zone 1 (220 cm per m/m) to Zone 3 (360 cm per m/m). All three zones showed a statistically significant (p < 0.05) relationship between ALT and slope. The highest correlation between slope and ALT was observed in Zone 3 (Figure 6c), whereas the lowest correlation was observed in Zone 1 (Figure 6a). All zones showed an increase in ALT as curvature increased (Figures 7a–7c). Across the study area, we observed a positive correlation between NDVI and ALT (Figure 8), with the NDVI values ranging from 0.48 to 0.51. In Zone 1, NDVI values were representative of the vegetation conditions, whereas in Zone 2 and Zone 3, NDVI values were more representative of soil moisture and inundation conditions.

Figure 6

Ensemble averaged ALT (±1 std; ) distribution as a function of slope is given in (a) Zone 1 (0 ≤ DD ≤ 7), (b) Zone 2 (7 < DD ≤ 15), and (c) Zone 3 (DD > 15) regions.

Figure 7

Ensemble averaged ALT (±1 std; ) distribution as a function of curvature is given in (a) Zone 1 (0 ≤ DD ≤ 7), (b) Zone 2 (7 < DD ≤ 15), and (c) Zone 3 (DD > 15) regions.

Figure 8

Ensemble averaged ALT (±1 std; ) distribution as a function of NDVI is given in for all Zones (i.e., Zone 1 (0 ≤ DD ≤ 7), Zone 2 (7 < DD ≤ 15), and Zone 3 (DD > 15)). Low NDVI values are associated with elevated centers and rims in HC and low centered polygons.

Ensemble averaged ALT (±1 std; <ALT>) distribution as a function of slope is given in (a) Zone 1 (0 ≤ DD ≤ 7), (b) Zone 2 (7 < DD ≤ 15), and (c) Zone 3 (DD > 15) regions.

Ensemble averaged ALT (±1 std; <ALT>) distribution as a function of curvature is given in (a) Zone 1 (0 ≤ DD ≤ 7), (b) Zone 2 (7 < DD ≤ 15), and (c) Zone 3 (DD > 15) regions.

Ensemble averaged ALT (±1 std; <ALT>) distribution as a function of NDVI is given in for all Zones (i.e., Zone 1 (0 ≤ DD ≤ 7), Zone 2 (7 < DD ≤ 15), and Zone 3 (DD > 15)). Low NDVI values are associated with elevated centers and rims in HC and low centered polygons.

The most important contribution to error in our analysis is the relative vertical error (point-to-point accuracy), because the topographic variables we used in this work are the first derivatives of the elevation (slope) and the second derivatives of the elevation (curvature). The relative vertical error of the LiDAR elevation data set was 0.0014 m (equations (3)-(6)). The same approach was used to compute the relative vertical errors in three microtopographic classes used in this study. Zone 1 showed 0.001 m relative vertical error, whereas Zone 2 and Zone 3 showed 0.0016 and 0.0020 m relative vertical errors, respectively.

We also tested how the accuracy of LiDAR elevation data set propagate into LiDAR-derived slope, curvature, DD, and NDD. The effect of the elevation error on slope was quantified using equation (5). The computed relative accuracy in slope was 0.004 m/m (∼0.004°) for the maximum slope value (i.e., 0.07 ± 0.004 m/m; p < 0.05). We also computed the relative error in slope of Zone 1, Zone 2, and Zone 3 separately: Zone 1 was 0.005 m/m (0.005°) whereas slopes in Zone 2 and Zone 3 were 0.003 m/m (0.003°) and 0.004 m/m (0.004°; p < 0.05). We used the same approach to compute the relative error in curvature by considering the curvature as the rate of change in slope for each pixel of the elevation data. This approach resulted in a relative error of curvature of 0.01 1/m for the maximum curvature of 0.04 (0.04 ± 0.01 1/m). The relative error of DD was 0.0028 m, and the relative error of NDD was 0.0056 m.

4.2. ALT Prediction

Based on the methodology outlined in section 2 above, we computed the coefficients in regression (equation (7)) for the each of three zones. Using the piecewise regression model, we predicted ALT values across a 5 km2 area (shown only for a 0.15 km2 subset from the 5 km2 for better microtopographic feature representation) for which the LiDAR elevation data were available (Figure 10b). The r2 of the regression model was 0.76 (p value < 0.001, RMSE 4.2 cm, and MAE 3.7 cm). The RMSEs was 2.24 cm, and RMSEus was 3.63 cm. The percentage of RMSEs was 27%, and the RMSEun was 73% of the overall RMSE. The MAE and MBE were 3.54 and −0.87, respectively, for the training data set.

Figure 10

(a) ALT for a ∼300 × 500 m2 of the study area, predicted using a data fusion approach adapted for use with data from remotely sensed measurements, showing the geophysical transect location (yellow line), and the location of additional ALT measurements on either side of the transect. (b) The relative uncertainty of the ALT measurement computed separately considering their microtopographic zonations (i.e., Zone 1, Zone 2, and Zone 3).

The model predicted ALT values between 20 and 70 cm, while the respective field measurements ranged from 22 to 58 cm. Lower ALT values were predicted for the HC polygons and rims of LC polygons, while higher ALT values were predicted for troughs and the center areas of the LC polygons (Figure 10b).

The model successfully predicted both the magnitude and variability in a spatially coherent way, without major systematic bias (Figure 9a). Overall r2 between training and predicted values along the transect, site-B, and site-C was 0.81 (p value 0.003; Figure 9b2). MAE and MBE were 3.54 and −0.87. RMSE was 3.99 cm, and RMSEs and RMSEus were 2.17 cm and was 3.35 cm, which were 24% and 76% from the overall RMSE, respectively. No major systematic spatial error in the model was observed in either the training or validation data sets, and errors were mainly associated with random noise. The highest correlated variable with ALT was slope; the lowest correlated variable was NDVI. We also compared the predictive power of all the parameters in ALT predictions and found that the partial covariance of slope was 158% higher than DD. Further, the correlation between DD and ALT is 144% and 121% higher than that of the NDVI and curvature with the ALT.

Figure 9

(a) Predicted ALT plotted along the measurements obtained along a geophysical transect. (b1) ALT predicted values at the calibration points are shown (green) and the measured values are shown (blue). (b2) Measured thaw-depths values at three validation sites (transect, Zone B, and Zone C) in relation to predicted ALT values. (b3) Error associated with measured and predicted ALT values at the validation sites.

The highest relative uncertainty (mean uncertainty 1.3%, standard deviation 12%, maximum <12%) in predicted ALT was shown by Zone 3 (polygons centers), whereas the lowest uncertainty (mean uncertainty 1%, standard deviation 10%, maximum <44%) was shown by Zone 1 (rims and elevated ridges; Figure 10).

5. Discussion

ALT is fundamentally a function of the warming experienced by the soil and subsurface. At regional scales, broad patterns in ALT correspond to mean annual air temperatures, with warmer regions having greater ALTs [Romanovsky et al., 2002; Jorgenson et al., 2001]. Vegetation and snow cover, organic and mineral composition of soils, and soil moisture content all influence the transfer of heat from the atmosphere into the ground [Romanovsky and Osterkamp, 1995]. These attributes vary both in time and space, and the thermal conductivity of the ground depends both on the total water content and the phase of the water [Burn and Smith, 1988; Smith and Riseborough, 2002]. Because these properties, as well as topography, vary at small scales (1–100 m), ALT may show significant local deviations from regional averages. For example, in Barrow, Alaska, a continuous time series of annual ALT records from 1994 to 2009 at the nearby CALM site (700 m east and 3500 m north of our study site) demonstrated significant variation in ALT among different landcape types, including beach ridges, drained lake basins, and polygonal uplands [Shiklomanov et al., 2010]. They found that ALT in beach ridges was 20 cm greater than that of polygonal uplands. Drained lake basins also showed on average 5 cm deeper ALTs than polygonal uplands. Thinner ALT values were documented in areas south and southwest of Barrow, on the Arctic coastal plain [Nelson et al., 1998].

Because the transfer of heat from the atmosphere into the subsurface is so strongly controlled by local surface and subsurface properties, no relationship between a deepening of ALT and observed climatic warming has been universally documented. In northern Mackenzie basin, NWT, ALT has shown a clear increase from 42 to 50 cm from 1998 to 2008 [Burn and Kokelj, 2009]. Burn and Kokelj [2009] also observed an increase in permafrost temperature of approximately 0.1°C per year from 1990 to 2000 at a depth of 28 m. In contrast, permafrost temperatures at the same location in Barrow warmed by 1°C at 15 m depth between 9 October 1950 and 9 October 2001. However, despite this warming, no associated trends in ALT were observed [Romanovsky et al., 2002] until 2005. The lack of global trends between ALT and permafrost temperature suggest local eco-hydro-geomorphology may play a significant role in ALT spatial distributions [Shiklomanov et al., 2010]. In this context, it is highly important to understand the complex interactions of microtopography, vegetation, soil moisture, etc., on ALT distributions at submeter scales.

In this study, we developed a piecewise regression model representing three microtopographic zones for ALT estimations. These zones (Zones 1, 2, and 3) were delineated using the microtopographic positions of polygonal features (quantified by directed distance) and associated ALT measurements (Figure 4). We hypothesized that these observed zones reflect the geo-hydro-thermal regimes controlled by geomorphic attributes in combination with soil moisture, water storage, and vegetation conditions. Previous studies in polygonal landscapes have documented that the zonations we observed here have unique geo-hydro-thermal process controls [Gamon et al., 2012; Liljedahl et al., 2012; Zona et al., 2011; Bockheim et al., 1999; Hinzman and Kane, 1991; Pavlov, 1975; Nelson et al., 1997; Hinzman et al., 1998; Shiklomanov and Nelson, 1999], but few studies have attempted to isolate geomorphic influences on ALT at finer scales [Minke et al., 2009]. We used LiDAR-derived topographic slope, curvature, and NDD as well as NDVI (a surrogate variable for vegetation and soil moisture variability) to capture the influence of local geomorphology, vegetation, water storage, and soil moisture on ALT. We specifically applied the regression model for three zones separately, and estimated model parameters unique to each zone. Low NDD values in Zone 1 were associated with the elevated centers and rims where low ALT values were observed. Because lower NDD values in Zone 1 were associated with low flow accumulations, high slope conditions, and high runoff conditions, this zone is well drained and has low water storage, deeper water tables, and low soil moisture conditions [Liljedahl et al., 2012; Hinzman and Kane, 1991; Pavlov, 1975]. Zone 1 is primarily a divergent landscape that warms quickly during the day and cools rapidly at night; as such, a lower heat transfer into deeper soil layers can be expected. Furthermore, because the vegetation in Zone 1 has higher albedo [Gamon et al., 2012], less heat is available to transfer into deeper soil layers. When low heat fluxes are accompanied by low moisture conditions, active layer development is limited to shallow depths [Minke et al., 2009], resulting in lower ALT values in Zone 1.

In Zones 2 and 3, areas of lower slopes and positive (concave) curvatures lead to greater water storage. Areas of positive curvature in these zones showed deeper ALT values than in areas of negative curvature in Zone 1 (Figures 7a–7c). The concave areas are associated with greater ground heat flux, resulting from high water and moisture availability [Shiklomanov et al., 2010; Zona et al., 2011]. Furthermore, concave areas in Zone 2 and Zone 3 tend to pond water and retain soil moisture, as opposed to convex areas in Zone 1 that more readily drain. Ponded waters and higher soil moisture in Zone 2 and Zone 3 likely lead to enhanced conductive transfer of heat into subsurface layers, resulting in higher ALTs (Figures 6b, 6c, 7b, and 7c). Furthermore, lateral flows of organic matter transferred through wind and surface-water runoff [Zona et al., 2011] into low-slope, concave portions of the landscape may help promote the thicker organic layers observed in Zone 2 and Zone 3 [Bockheim et al., 1999]. These thicker organic layers help retain soil moisture and high moisture conditions can increase the heat conductance into deeper layers [Hinzman and Kane, 1991; Pavlov, 1975; Williams and Smith, 1989], promoting deeper ALTs in Zone 2 and Zone 3.

Analysis of NDVI data suggest that vegetation caused NDVI values to vary from 0.32 to 0.48 across different polygon morphologies. Lower NDVI values were associated with Zone 1 and higher NDVI values were associated with Zone 2 and Zone 3. Furthermore, NDVI values inversely correspond with the ALT values (thinner ALT for lower NDVI values; Figure 8). It is also possible that NDVI values in this polygon landscape represent not only vegetation condition, but also soil moisture conditions and inundations. For example, in the regression model for Zone 1, NDVI represents vegetation conditions, whereas in Zone 2 and Zone 3, NDVI potentially represents moisture conditions as well. Because each zone has a unique regression-model representation, our approach can yield a different parameterization for the same explanatory variable in different zones.

The highest uncertainty in the ALT predictions was associated with Zone 2 (center areas), whereas the lowest uncertainty was associated with Zone 1 (rims and elevated center areas). Our results suggest that Zone 1, representing topographic highs, has a lower ALT compared to Zones 2 and 3 (centers areas and troughs of LC polygons) representing wet topographic low areas. This is primarily because of the differences in albedo, soil moisture, vegetation, and polygon morphology in ice wedge polygon dominated tundra environments. If this field campaign were conducted in the first half of the the thaw season, it might be possible to see deeper thaw depths in the HC polygons, because thawing could be much faster in HC polygons compared to LC polygons. The greater energy in HC polygons would transfer more directly as conductive heat transfer than latent heat of fusion.

LiDAR data maintain higher levels of relative vertical accuracies than absolute accuracies, and these accuracies can change with different microtopographic types. In this work, we also introduced a simple approach to compute the relative vertical errors of LiDAR elevation data, using which we showed that how these errors propagate into LiDAR-derived topographic attributes like slope, curvature, and directed distance. These relative errors can be useful in understanding model uncertainties that arise when input data contain LiDAR-derived attributes. In future studies, we recommend the collection of more ALT measurements of center areas (LC polygons) and troughs than of rims (LC polygons) and elevated center areas (HC polygons) of the ice wedge polygon features.

The ALT prediction method we investigated in this work can be applied to other parts of the tundra, where ALT is primary controlled by microtopography, vegetation, and soil moisture, and where LiDAR and WorldView-2 data available. WorldView-2 data and LiDAR data are becoming increasingly available for this region (e.g., USGS has acquired LiDAR data along the Arctic coastal belt—personal communication with Ann Gibbs USGS Pacific Coastal and Marine Center, Santa Cruz, CA). In other Arctic environments, where polygons are not present, it is essential to calibrate the model parameters with proper site-specific ALT measurements. For example, the spatial distribution of ALT in Barrow is highly dependent on the microtopographic variability introduced by the ice-wedge formation. For the landscapes where ice wedges are not the primary control, ALT depends on topographic position [Nelson et al., 1997]. In hilly environments, particularly in Arctic foothills, ALT depends on topographic positions from the top of a hillslope to the valley and channel floodplains. The topographic attributes (slope, curvature, and directed distance) used in this study can be effectively used to characterize the topographic position in hillslopes [Gangodagamage et al., 2011]. The models proposed in this work could work well for environments where microtopography primarily consists of features such as found in our Zone 1, Zone 2, and Zone 3. If this model is used in landscapes where the microtopographic features, surface, atmospheric, and subsurface dynamics are different, the model needs to be redeveloped using region specific data sets, which will possibly provide better generalizations for ALT variability for that area.

We used the regression equation to predict the ALT using three regression models representing three microtopographic classes (i.e., Zone 1, Zone 2, and Zone 3). The sample sizes we used satisfy both the likelihood (probability) of obtaining a significant (p < 0.05) and r2. Traditional rules like statistical significance do not provide insights about statistical power [Asher, 1993]. Statistical significance itself may be obtained from a large enough sample size, even though there is very little relationship between the dependent and explanatory variables [Asher, 1993]. It is important when developing predictive regression models, the sample sizes to be chosen so as to provide both statistical significance and statistical power, which then can be expected to yield better generalization of the model beyond the samples used to develop the model.

6. Conclusions

Current methods are very good at quantifying and predicting ALT at either local or regional scales. No method has been presented that is able to capture both the fine scale variation in ALT related to the microtopography of rims, troughs, and centers in ice wedge polygon landscapes, and the influence of microtopography on ALT over large spatial domains. Our results indicate that it is possible to predict ALT with a high level of accuracy at high spatial resolution in a noninvasive manner, using a combination of microtopographic and land-surface characteristics over large spatial domains.

It is generally accepted that ALT increases in response to climate warming [Lawrence et al., 2008], but how this increase is related to local topographic variability within and between polygons is still poorly quantified. The response of ALT to climate warming is likely to be complex and dependent on many parameters. As has been documented, the primary control of ALT at larger scales (∼1 km or larger) is air temperature [Zhang et al., 1997; Nelson et al., 1997], which potentially controls the average ALT values at landscape scale. Jorgenson et al. [2010] have found correlations between ALT and hillslope position, Zona et al. [2011] showed that ALT varies with polygon microtopography, and Minke et al. [2009] have observed that shallower thaw depths were associated with higher microrelief regions of LC polygons. In this study, we found that local variables, including the heterogeneity of soil moisture, water ponding, and vegetation conditions (all quantified by NDVI)—as well as geomorphic attributes quantified by slope, curvature, and NDD—play a significant role in small-scale (∼2 m) ALT variations. Consistent with previous studies, our study suggests that while broad climatic variability associated with warming will govern the regional averages of ALT, the smaller-scale variability of ALT is controlled by the aforementioned local geomorphic, hydrological, and vegetation variables.

Furthermore, this study suggests the importance of developing techniques to better landscape characterization, and the importance of identifying the process-based microtopographic zones that mainly govern the spatial distribution of ALTs. Our approaches enabled spatially classified process-based zones to capture complex geohydrothermal interactions using simple linear representations, among model variables controlling the spatial variability of ALTs. While our model does require local data for the calibration of the metrics, it should be widely applicable to polygonal terrains outside the Barrow, AK. This applicability arises from the hydrological, thermal, and geomorphic controls ALT captured via the landscape metric used in our model.