Permafrost and organic layer interactions over a climate gradient in a discontinuous permafrost zone

Permafrost is tightly coupled to the organic soil layer, an interaction that mediates permafrost degradation in response to regional warming. We analyzed changes in permafrost occurrence and organic layer thickness (OLT) using more than 3000 soil pedons across a mean annual temperature (MAT) gradient. Cause and effect relationships between permafrost probability (PF), OLT, and other topographic factors were investigated using structural equation modeling in a multi-group analysis. Groups were defined by slope, soil texture type, and shallow (<28 cm) versus deep organic (≥28 cm) layers. The probability of observing permafrost sharply increased by 0.32 for every 10-cm OLT increase in shallow OLT soils (OLTs) due to an insulation effect, but PF decreased in deep OLT soils (OLTd) by 0.06 for every 10-cm increase. Across the MAT gradient, PF in sandy soils varied little, but PF in loamy and silty soils decreased substantially from cooler to warmer temperatures. The change in OLT was more heterogeneous across soil texture types—in some there was no change while in others OLTs soils thinned and/or OLTd soils thickened at warmer locations. Furthermore, when soil organic carbon was estimated using a relationship with thickness, the average increase in carbon in OLTd soils was almost four times greater compared to the average decrease in carbon in OLTs soils across all soil types. If soils follow a trajectory of warming that mimics the spatial gradients found today, then heterogeneities of permafrost degradation and organic layer thinning and thickening should be considered in the regional carbon balance.


Introduction
The vulnerability of permafrost to climate change is a pressing question because the impacts of warmer temperatures may be greater in polar regions than in other regions (Brohan et al 2006, Trenberth et al 2007. Substantial amounts of carbon is currently stored in frozen soils which if released could affect the global carbon balance and climate change (Schuur et al 2008. The organic soil layer in particular is a large carbon pool in permafrost regions. For example, Histels and Histosols account for perhaps 30% of the total carbon pool in the top 1-m of soil (Tarnocai et al 2009). Although organic layer thickness (OLT) and permafrost distribution in the discontinuous zone are known to be tightly coupled (Harden et al 2006, Ping et al 2006, there is limited information about the linkage between permafrost degradation and the organic layer carbon pool in response to warming temperatures at regional scales (Yuan et al 2012). This is partly due to the difficulty in modeling interactions of the organic layer with both slow (e.g. top-down thaw from warming, 'press events') and rapid (e.g. fire and thermokarst, 'pulse events') disturbances (Grosse et al 2011).
Besides OLT and air temperature, the presence and persistence of permafrost is controlled in large part by soil texture, slope, and microtopography that influence site and landscape level drainage conditions (Swanson 1996a, Schuur et al 2008. When permafrost thaws, water is redistributed and localized drainage conditions are altered to form uneven patches of high and low areas (i.e. thermokarst) (Jorgenson and Osterkamp 2005). At some locations permafrost is able to persist under thick OLT because the low thermal conductivity of the organic layer acts as an insulator of soil temperature from air temperature (Hinzman et al 1991, Harden et al 2006. However, this effect may not be significant at locations with very deep organic layers such as unfrozen peats where water saturation controls the thermal properties of the soil. In the latter case, as permafrost degrades under warmer conditions, greater water saturation and lower decomposition rates may result, leading to an increase in OLT (Turetsky et al 2000). Other factors controlling permafrost include snow depth, vegetation, slope aspect and fire disturbance (Jorgenson and Osterkamp 2005, Hinzman et al 2006, Ping et al 2006. The challenge remains to quantify the impact of warming air temperatures in a way that accounts for these interactions (Schuur et al 2008).
One approach to understanding how the organic layer interacts with temperature to influence permafrost distribution is to observe soils in their natural state over a temperature gradient, i.e. a 'space for climate' analysis. An advantage to this approach is that interacting and complex soil-forming processes over environmental gradients are potentially captured by the observation data, assuming their conditions are at a near steady state or endpoint. Interpretations of relationships are strengthened when they are supported by soil-forming theories (Jenny 1941) and by exploring the data for potential thresholds that may simplify nonlinear relationships into linear relationships.
Interior Alaska lends itself to a study of permafrost response across environmental gradients because it is situated mostly within the discontinuous permafrost zone wherein there exists a variety of soil temperatures, topography, and parent material types. Our objective was to estimate changes in permafrost probability (PF), and OLT within Interior Alaska across a mean annual temperature (MAT) gradient. Permafrost probability in this context is the proportion of soil pedons having permafrost for the population or subpopulation of sampled locations. We apply structural equation modeling (SEM, i.e. path analysis) using observed variables to describe the complexity of the OLT and PF relationship for different environmental conditions (Grace 2006). Over 3000 mostly near-surface (1-m depth) observations were analyzed, emphasizing interacting effects of soil texture, landscape topography, microtopography, MAT, OLT, and PF. Finally, we discuss the implications of the results with regards to the regional carbon balance.

Field observations data sets
There are two datasets on which we base our analysis: (1) an Alaska soil carbon database (about 500 observations; Johnson et al 2011), and (2) the National Resources Conservation Service (NRCS) pedon dataset for Alaska (about 2800 additional observations; NASIS 2010) ( figure 1(a)). The latter database supports a number of georeferenced profiles that were not sampled for laboratory analysis but have recorded horizon designations and soil texture type. All Oi, Oe, and Oa horizons (excluding buried O horizons) were summed to determine the total organic layer thickness (OLT) in cm. The presence or absence of permafrost was determined by if a pedon had an 'f' or 'ff' (frozen) suffix assigned to one or more master horizons within at least the top 1-m of soil. The 'f' and 'ff' suffixes designate the presence of high moisture content and low moisture content permafrost, respectively (Soil Survey Staff 2010). The field observer relied on late summer or fall observations to make adequate estimates of the top of the permafrost. Specifically, in soils with high moisture content, the upper limit of the permafrost table was typically indicated by the top of a zone of ice lenses and masses, while low moisture content soils were judged more subjectively. Some profiles with a total depth less than 1-m were still included if the bottom horizon was bounded by bedrock ('R' or 'Cr'). If it was uncertain if a frozen horizon existed in the top 1-m (e.g. because the profile was not excavated to a 1-m depth) it was left out of the analysis. Organic layers greater than 1-m were included in the analysis. Shallow OLT (OLTs) was eventually defined as OLT <28 cm and deep OLT (OLTd) as ≥28 cm per the segmentation model below. Only profiles that were located within interior Alaska (i.e. 'Intermontane Boreal' ecoregion as defined by Nowacki et al 2002) were included in the final analysis. The sampling design of the overall dataset may be described as clustered and opportunistic, noting that individual field sampling campaigns typically aim to sample a wide variety of environmental conditions. There is poor representation of areas with MAT < −4 • C (light gray bars, figure 1(b)) which causes problems for multi-group analyses. Therefore, unless otherwise indicated, only observations with MAT ≥ −4 • C were analyzed, for a total of 3023 observations. The total area within this range of temperatures accounts for 52% of the entire interior region, or roughly the southern half (figure 1(a)). Soil texture ('particle size') was noted in the field for 87% of the pedons using field estimates of the percentage of sand, silt, and clay as well as the volume of rock fragments greater than 2 mm. Soil texture observations were grouped into more generalized particle size classes within a particle size control section, usually extending from 25 to 100 cm (Soil Survey Staff 1999). Particle size was further generalized and harmonized across datasets into four classes-silty, loamy, sandy, and variable (table S1 available at stacks.iop.org/ ERL/8/035028/mmedia). Pedons classified as 'variable' were predominantly of the 'loamy-skeletal' class in the NRCS database, which typically form on steep slopes and ridges (Ping et al 2006). 'Sandy' pedons were predominantly some variant of 'sandy-skeletal', with relatively few purely 'sandy' pedons (table S1). For the remaining pedons that had no soil texture information, the information was extracted from the parent material layer of the Permafrost Characteristics of Alaska map (Jorgenson et al 2008).
The depth of the top of the charcoal layer (n = 849) and soil temperature at 50-cm depth at the time of excavation (n = 1846) were available for only some of the pedons, limiting their application to fewer analyses. Charcoal includes large fragments and disseminated dark zones that can be identified as charcoal by visual examination. The top depth at which the frozen 'f' horizon occurred was used as the active layer depth (n = 1217). Each profile was visited and observed only one time; the visits occurred from June to September and between the years 1998 and 2008.

Spatial data sets
The raster values of several spatial data sets were extracted to the point data set of field observations. Mean annual temperature (MAT) was extracted from the 2-km Parameter-Elevation Regression on Independent Slopes Model for the years 1961-1990 (PRISM) (Daly et al 2000) ( figure 1(a)). Seasonal variants of temperature and precipitation and growing season length were also accessed from the Scenarios Network for Alaska Planning (SNAP; www.snap.uaf.edu), though they were not included in the final analysis (supplemental material available at stacks.iop.org/ERL/8/ 035028/mmedia). Slope and heat load index (HLI) at 100-m were derived from a 60-m USGS digital elevation model (Gesch 2007). Slopes were ranked and binned and then classified into 'flat' (0-0.69%), 'low slope' (0.7-3.89%), and 'high slope' (≥3.9%) according the findings of the current study. Heat load index (HLI) was calculated as: sometimes referred to as 'northness', where higher values indicate more northerly slopes (e.g. Balice et al 2000).

Statistical analysis
After compiling all the field observations data with the spatial data, statistical analyses were carried out in two parts. The first part was exploratory wherein we sought to identify predictors of PF, and directional changes of the relationships between predictors and PF ('switch points' of nonlinear relationships). We fit a segmented binomial regression, with unknown switch point, using PF as the response and OLT thickness as the single covariate. Partitioning the dataset according to the identified switch point then allowed for linear analyses of the specific effects of OLT and MAT on PF or, alternatively, MAT and PF on OLT. The observations were also partitioned into groups at the inflection points along a slope gradient. We compared mean PF and mean OLT for different classes of slope (flat, low and high), soil texture (sandy, variable, loamy and silty), and OLT depth (>28 cm, ≥28 cm) for univariate analyses. Hereafter we refer to slope and soil texture combinations as 'landforms' for simplicity. Where statistical differences of means of PF or OLT between landforms are reported, Fisher's Exact Test was used for binomial responses and Tukey's Honestly Significant Difference test for continuous responses. In the case of unequal variances for continuous responses, the non-parametric Wilcoxon Rank Sum test was used.
In the second part of the analyses, we investigated the effects of MAT, slope, and HLI on PF and OLT with structural equation modeling (SEM) for each landform. The important characteristic of SEM for this study is the path coefficients account for both direct and indirect effects of the predictors on OLT and PF (Grace 2006). The path coefficients are similar to regression coefficients (i.e. 'slopes') in that they describe the relationships between variables. Another advantage of SEM is that the direction of the effects can be modified so that for some conceptual models OLT is the predictor and PF the response, but in other cases OLT can be made the response (or result) of PF. Of the 24 landform and OLT class combinations, 18 had sufficient data (>50 observations) for SEM. For these 18 combinations, we investigated the direction and magnitude of MAT path coefficients in a multi-group analysis. The decision of model structure for each combination was based on relationships identified in literature and that made sense ecologically to the investigators (supplementary material available at stacks. iop.org/ERL/8/035028/mmedia). Therefore, this part of the study was strictly confirmatory where the goal was to test that the model structure was supported by significant correlations. In this type of SEM analysis, non-significant paths should not necessarily be deemed unimportant because the data in hand may not adequately show the relationship; therefore non-significant paths were still included in the final model but indicated with a dashed line. More details about SEM can be found in Grace (2006) and the references therein.
Nominal logistic regressions with PF make no assumptions about the normality of the response variable. A maximum likelihood estimator of standard errors and chi-square test statistic that are robust to non-normality were used to account for other non-normalities in the regression variables (Muthén andMuthén 1998-2012). This method uses a sandwich estimator to compute the standard errors, resulting in higher standard errors than those computed with the conventional method. The distribution of OLT was not normal, being skewed to the right especially at depths >28 cm. Therefore, OLT was natural log(ln) transformed to further address non-normality in the confirmatory SEM analysis. Log transformations of other model terms risk a loss of interpretability and did not greatly improve data normality. Although the effects of spatial autocorrelation were reduced by including climate and topographical variables, a separate analysis was performed to understand this effect on path coefficients (supplementary material). Non-recursive models (i.e. feedbacks) were not explored in the current analysis as they are difficult to represent in statistical models (Grace 2006). All SEM was performed with MPlus software version 6.1 (Muthén andMuthén 1998-2012) and all other statistics were performed in JMP version 9.0 (SAS Institute Inc. 2010).

Exploratory analysis-univariate relationships and interactions
3.1.1. Relationship of OLT with PF and other soil measurements. The relationship between PF and OLT was nonlinear, with a change point of 28 cm as identified by the segmented binomial model (figure 2). There was a positive correlation until a depth of 28 cm (R 2 = 0.94, RMSE = 0.07), and thereafter became negative at deeper OLT (R 2 = 0.71, RMSE = 0.17) (figure 2). For every 10-cm increase in OLT within OLTs soils, there was a 0.32 increase in PF. For OLTd soils, there was a more gradual change in the opposite direction; for every 10-cm increase in OLT there was a 0.06 decrease in PF. The 28-cm switch point was used to partition observations into classes of shallow OLT and deep OLT for further analyses. Note that OLT was log transformed in the confirmatory analysis using SEM for simplicity and lack of data at deeper depths (see section 4).
There were significant correlations of OLT with other soil variables using a subset of the OLTs pedon dataset (figure 3). Soil temperature at 50-cm was significantly negatively correlated with OLT (R 2 = 0.89, RMSE = 0.66) as was active layer depth (R 2 = 0.41, RMSE = 6.67). The depth of charcoal was significantly positively correlated with OLT (R 2 = 0.76, RMSE = 2.64). Soil temperature and charcoal depth were unequally distributed among the landform classes presented below and so were not analyzed further.
3.1.2. Relationship of PF and soil texture with slope.
Permafrost and soil texture varied nonlinearly as slope increased from very flat to low slopes to high slopes (figures 4(a) and (b)). In flat areas (0-0.7% slope) there was little or no trend in PF with increasing slope, although these areas tended to have lower PF than low slopes. From 0.7 to 3.9% there was a gradual increase of mean PF with increasing slope with sandy soils being replaced by silty and variable soils. On high slopes (>3.9%) there was a rapid decrease in PF with increasing slope and with silty soils being replaced by loamy and variable soils.
3.1.3. Relationship of PF and the interaction of slope, soil texture, and OLT depth. There were distinct patterns in both OLT and PF across landforms characterized by soil texture and topographic position (figure 5). Mean OLT was consistently lower on high slopes for all soil texture . The relationship between permafrost probability and slope (a) and soil texture type (b). The observations are grouped into bins that were ranked by slope. Every bin has at least 80 observations and no more than 130. Also depicted are the slope classes with regards to PF response over the gradient. types ( figure 5(a)). In flat and low slopes, silty soils were significantly higher in OLT than sandy or loamy soils (P < 0.0001, Wilcoxon rank sum test). Mean PF in OLTd soils was consistently higher than in OLTs soils for all situations except variable low slope and silty flat (P < 0.05) ( figure 5(b)). Furthermore, the difference in mean PF between OLTs and OLTd soils was generally smaller for flat and low slopes than for high slopes. The OLTs sandy soils had the lowest PF, regardless of slope class, while silty soils had higher PF. The PF in other OLTs soil types all decreased from low to high slopes (P < 0.05, all comparisons). In contrast, mean PF generally increased from flat to low slopes to high slopes in OLTd soils, although significant differences were observed only between flat and high slopes for loamy and silty soils (P < 0.05).
The number of observations across landforms indicated the relative frequency that OLTs and OLTd soils occurred by different landforms, ranging from 6 in sandy high slopes to 383 in silty high slopes ( figure 5(c)). For example, the occurrence of OLTd soils on sandy high slopes was uncommon. There were more OLTs observations on high slopes for all soil texture types, except sandy, than on flat and low slopes. In contrast, there were more OLTd observations in flat areas than on high slopes. The six combinations with fewer than 50 observations probably rarely occur in interior Alaska and are less important to the regional analysis, thus they were omitted from the confirmatory SEM analysis.

Confirmatory analysis-SEM
The SEM model constructs were modified for specific conditions to allow representation of the positive effect of OLT on PF in OLTs soils, and the negative relationship of OLT and PF in OLTd soils according to the above analysis and theory (supplementary material available at stacks.iop. org/ERL/8/035028/mmedia). The direction of the relationship was also modified so that in OLTs soils OLT was the predictor of PF, but in OLTd soils PF was the predictor of OLT. The paths chosen for SEM analysis, were generally (but not always) significant.

OLTs soils.
As expected, lnOLT was significantly positively correlated with PF for every soil texture and slope class ( figure 6(a)). The path coefficients between lnOLT and PF were consistently larger on high slopes than in flat areas, by a factor of 1.6-1.8, suggesting a stronger effect of OLT on PF on high slopes. Also as expected, the correlation between MAT and PF was negative and significant for almost all the landform types (figures 6(a) and 7(a)). The exceptions were variable low slopes (which had a relatively low sample number) and sandy flat and low slopes (where PF was already low) (figure 6(a)). Similar to the lnOLT and PF relationship, the path coefficients were consistently higher on high slopes than in flat areas regardless of soil texture type although the confidence intervals of the MAT on PF regressions were large ( figure 7(a)).
In none of the models was there a significant direct effect of HLI on PF (figure 6(a)). However, HLI indirectly affected PF as more northerly slopes (higher HLI) were positively correlated with higher lnOLT, which in turn led to higher PF. The exception to this was in loamy soils, although the same effect may not be detectable with the current data. The models did not confirm that the slope parameter was significant in almost all of the landforms, except in silty high slopes ( figure 6(a)). Finally, MAT was negatively correlated with lnOLT in at least six of the landforms (figures 6(a) and 7(b)). However, the models did not confirm that MAT was correlated to lnOLT in sandy and variable soils ( figure 6(a)).

OLTd soils.
As expected, there was a significant negative effect of MAT on PF for all of the landforms (figures 6(b) and 7(a)). Similar to OLTs, the confidence intervals associated with the PF on MAT relationship were large across landforms (figure 7(a)) and the paths from slope to lnOLT and PF were always non-significant. Also as expected, there was a significant negative relationship between PF and lnOLT in all the landforms that included this effect ( figure 6(b)). The effect of MAT was greater on low slopes than in flat areas, where the temperature gradient from −4 to −1 • C was 2 and 5 times higher for loamy and sandy soils, respectively. Finally, there was an indirect effect of MAT on lnOLT, because lnOLT increased simultaneously as PF decreased over the gradient ( figure 7(b)). The effects of PF and MAT on lnOLT were not included in the high slope landforms because we assumed that these were de-coupled due to a lack of thermokarst on high slopes (supplementary material available at stacks.iop.org/ERL/8/035028/mmedia).

Changes in SOC.
Structural equation results (figures 6(a), (b), 7(a) and (b)) were transformed with (2). In OLTs soils, the average decrease was −0.15 gC cm −2 as MAT increases from −4 to −1 • C (table 1). In contrast, OLTd soils increased in SOC by an average of +0.58 gC cm −2 . The average increase in carbon in OLTd soils was almost four times greater in OLTd soils compared to the average rate decrease in carbon in OLTs soils across all soil types. The largest discrepancy between the two soil types for a given landform occurred on loamy low slopes, where the change in SOC in OLTd was almost eight times higher than the change in OLTs (table 1).

Controls on permafrost distributions
Organic layer thickness and PF are tightly coupled and reflect changing conditions of microtopography, pedogenesis, vegetation, climate, and disturbance. Permafrost is highly sensitive to changes in OLTs soils, whether due to an increase (e.g. vegetation succession) or a decrease (e.g. fire) in OLT. Clearly, fire disturbance is affecting soil thermal properties as illustrated by the relationships between OLT and soil temperature, charcoal depth, and the active layer depth (figure 3). Fire typically burns the organic layer, leaving char deposits and reducing its depth. As time passes after the fire, the active layer depth decreases as the organic layer depth and char depth both increase and permafrost returns (Fenton et al 2005, O'Donnell et al 2011, Harden et al 2012. The contrast of overlying organic soil layer properties and underlying mineral soil layers directly affects the thermal impact of the organic layer (Lachenbruch 1959). Harden et al (2006) found that soil temperature at a 5-cm depth decreased by about 0.5 • C for every 1-cm increase in OLT during summer months. In another analysis Kasischke and Johnstone (2005) found about a 0.35 • C decrease in the mineral soil for every 1-cm increase in OLT. Yoshikawa et al (2002) commented on the insulation effect of OLT in terms of the persistence of permafrost after fire, noting that if the OLT remaining after fire was more than 7-12 cm then the active layer thickness would not deepen. In the current study, the decrease was 0.24 • C at the 50-cm depth ( figure 3(a)) and the point at which this thermal impact, or insulation effect, begins to lose effect is between about 25 and 30 cm OLT (figure 2). If organic layer thickness decreases in OLTs soils due to the effects of increased fire activity or change in C balance in a warmer and drier climate (Goetz et al 2005, Turetsky et al 2010 then figure 2 depicts increased vulnerability to thaw. On the other hand, OLTd soils are typically wetter and less affected by the insulation properties of the organic layer and permafrost commonly occurred when OLT was between 17 and 85-cm. Thus, vulnerability to thawing in these conditions may be less influenced by fire (Swanson 1996b, Harden et al 2001 and more directly influenced by warming temperatures (Grosse et al 2011). Furthermore, there were apparently no landforms (with their associated drainage conditions) that mediated the effect of temperature over the gradient (figure 7(a)).
Permafrost distributions were strongly controlled by localized drainage conditions and soil texture and slope that interacted to promote or limit permafrost development. Permafrost had higher occurrence in silty and loamy soils than in sandy soils for the same slope classes, owing to the lower water storage capacity and thermal conductance of sandy soils (Swanson 1996a, Osterkamp et al 2000, Jorgenson et al 2010. The insulation effect of the organic layer was evident in all landforms except sand. For example, PF was limited by a thin organic layer on higher slopes in OLTs soils due to excessive drainage and poor soil development (Ping et al 2005) ( figure 5(b)). Interestingly, the opposite was true for OLTd soils, where PF tended to be higher on higher slopes ( figure 5(b)). One interpretation of the positive effect of slope in the relatively wet OLTd soils is that better drainage on high slopes allowed for higher permafrost formation across the gradient. However, we note that most of the OLTd soils within the 'high slope' class occurred on slopes between 4 and 10% and their occurrence on more extreme slopes was far less common. As slope class and soil texture were static (i.e. not likely to change in response to warming) and seemed to capture important gradients, they were used as constraints in the multi-group analysis of warming effects on permafrost.

Temperature effects on spatial distributions of permafrost and organic layer thickness
In addition to OLT, mean air temperature was a dominant control on permafrost distribution as locations with warmer temperatures had lower PF for most landforms and regardless of OLT class (figure 7(a)). However, the differences in PF across the temperature gradient were not equal for all landforms, where the more permafrost there was to lose the more was lost. Permafrost probability decreased rapidly in silty, loamy and variable landforms (from about 0.9 to the 0.5) over only a 3 • C gradient, yet sandy landforms generally did not have permafrost and thus showed little sensitivity to MAT (figure 7). The variability in response over the temperature gradient was even more evident in the organic layer as it both thinned and thickened (figure 7). We attribute thinner organic layers at warmer temperatures to higher fire activity, higher decomposition, and/or improved soil drainage after permafrost degradation (Kasischke and Turetsky 2006, Wickland and Neff 2008, Yang et al 2010, Grosse et al 2011. Thicker organic layers in warmer temperatures are consistent with many areas where permafrost degradation creates anaerobic conditions, low decomposition rates, and low combustion losses (Osterkamp et al 2000, Harden et al 2001, Ping et al 2006. More specifically, high peat accumulation in post-permafrost thaw sites has been found in Canada (Turetsky et al 2000, Camill et al 2001 and Alaska (Jones et al 2013). OLTd soils are probably high in ice and/or water content typical of Histels, Histosols, and generally wetter conditions (Jorgenson and Osterkamp 2005). In a warmer world, sites associated with high water tables are likely to support thick peats regardless of permafrost.
Although there was indication that some soils were frozen even up to 150-cm OLT (cf O'Donnell et al 2009), permafrost was relatively uncommon in soils with OLT deeper than 85-cm. Related to this point is that we have log-transformed OLT for SEM analysis, and there may be a nonlinear relationship between OLT and PF in OLTd soils (e.g. figure S1 available at stacks.iop.org/ERL/8/035028/mmedia). However we are wary to conclude that it is truly nonlinear for lack data; there were only 21 observations available at depths ≥150-cm compared to 656 observations at depths >28 and <150 cm represented in figure 2.
In our analysis we have relied on detecting consistent trends within OLTs and OLTd soils over multiple conditions with a large dataset to make inferences about warming effects. The 28-cm value we used to simplify our statistical analysis appears to reasonably distinguish relatively drained and aerated systems (OLTs) from wetter systems (OLTd) that occur adjacently in both uplands and lowlands. It should be acknowledged, however, that the cutoff should not be considered a fixed value. For example, the same analyses were performed using a 37-cm cutoff (determined visually) with very similar resulting trends (not shown). Additionally, there are surely situations when soils nearer the cutoff correlate with PF and MAT in a way that is opposite the trends depicted in figure 7. For example, when we included both OLTd with the OLTs soils in an OLT on MAT regression for high slopes, there was still a negative relationship (not shown), suggesting that both classes respond the same. Similarly, whereas we have modeled PF and OLT to decrease with increasing MAT in OLTs, in some situations the shift may actually be toward deeper OLT. The latter complication will in particular affect flat and low slope areas that have mixed occurrences of OLTs and OLTd.

Carbon balance implications at larger scales.
A significant amount of the total soil carbon pool of interior Alaska is in the organic layer component (e.g. 17-71% of the top 1-m, Johnson et al 2011). Thus, permafrost degradation and associations with organic layer changes have obvious consequences for the carbon cycle if soils follow a trajectory of warming that mimics the spatial gradients found today (table 1). Moreover, the vulnerabilities of carbon to fire, decomposition, and hydrologic shifts are likely to vary according to changes in permafrost (Harden et al 2012). A complete regional accounting of change in the soil carbon pool requires a process-based approach. Nonetheless, a calculation of soil carbon stocks from OLT (supplementary material available at stacks.iop.org/ERL/8/ 035028/mmedia) and extrapolation to all of interior Alaska (with areas derived from Jorgenson et al 2008) provides a useful estimate of regional carbon change. Based on current data and modern climate gradients (table 1), OLTs soils decreased by −0.28 PgC over a MAT gradient of −4 to −1 • C, but OLTd soils increased +0.43 PgC, for a net gain of 0.15 PgC. For perspective, this modest sink is 0.7% of the 22 PgC total estimated soil carbon pool for the discontinuous permafrost zone in Alaska (Mishra and Riley 2012). However, despite the high number of observations, the unsystematic sample may lead to a biased estimate at the regional level.
The coldest areas of the interior (−8 to −4 • C) may not undergo permafrost degradation as dramatically as those at warmer MAT, even after an increase of several degrees. At a MAT of −4 • C, permafrost probability was already high in some landforms (about 0.8-0.9), suggesting less potential for large changes at lower temperatures ( figure 7). Furthermore, we have not investigated changes in OLT that may occur at locations that were already warm (about 0.5 • C MAT) with sporadic permafrost. In peatlands, some have suggested a rapid carbon loss in a feedback of warmer/drier climate that lowers OLT and the water table, exposing more organic material to decomposition centuries after permafrost thaw (Ise et al 2008). Others have noted the instability of permafrost in general in interior Alaska with regards to climate change (Osterkamp et al 2000). Comparison of process-based model estimates of permafrost distribution for contemporary climate with the empirical results reported in this study are important to better understand the timing and impacts of permafrost degradation on the regional ecosystem carbon balance.

Future study
Accurate mapping of OLT would allow for more accurate mapping of permafrost in interior Alaska. Given the great uncertainty of the impact of permafrost thaw on global climate, devising systems of monitoring should be considered. For example, there may exist opportunities to incorporate spatially explicit estimations of organic layer thickness from remotely sensed data, including fire activity (Pastick et al 2013b). Emerging technologies such as electromagnetic imaging also offer intriguing possibilities (Pastick et al 2013a). Finally, in our opinion, besides mapping OLT, the greatest gap remaining for spatially explicit monitoring of permafrost is the careful mapping of soil parent material and soil drainage, which are now only coarsely mapped.

Conclusions
In this study we have 'let nature perform the experiment for us' by observing how permafrost and organic layer vary over a climate gradient in their natural state. The use of multi-group analysis and SEM helped disentangle individual effects, improve insights into cause and effect relationships, and account for both direct and indirect effects. The direction of change in the organic layer depends greatly on the interaction of air temperature and drainage. The consistency in trends shown for OLT and PF suggest that, with some exceptions, as temperature increased along the gradient: (1) OLTs soils generally decreased in thickness while OLTd soils generally increased in thickness, (2) the extent of PF generally decreased, and (3) landscape level drainage conditions (indicated by landform type) further influenced variations in PF and OLT. The accumulation of soil carbon in the organic layer over temperature gradients, or lack thereof at some locations, should be considered in the spatial accounting of the regional carbon balance. Turetsky M R, Wieder R K, Williams C J and Vitt D H 2000 Organic matter accumulation, peat chemistry, and permafrost melting in peatlands of boreal Alberta Ecoscience 7 379-92 Wickland K P and Neff J C 2008 Decomposition of soil organic matter from boreal black spruce forest: environmental and chemical controls Biogeochemistry 87 29-47 Yang Z-P, Ou Y H, Xu X-L, Zhao L, Song M-H and Zhou C-P 2010 Effects of permafrost degradation on ecosystems Acta Ecol. Sin. 30 33-9