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

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

Zoom In Zoom Out Reset image size

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.

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:

$$\begin{eqnarray} &&\mathrm{HLI}={\mathrm{e}}^{\cos \left (\mathrm{aspect}\ast \frac{\pi }{1 8 0}\right )} \end{eqnarray} \tag{ 1 }$$

sometimes referred to as 'northness', where higher values indicate more northerly slopes (e.g. Balice et al 2000).

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 and Muthé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 and Muthén 1998–2012) and all other statistics were performed in JMP version 9.0 (SAS Institute Inc. 2010).

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² = 0.94,RMSE = 0.07), and thereafter became negative at deeper OLT (R² = 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).

Zoom In Zoom Out Reset image size

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² = 0.89,RMSE = 0.66) as was active layer depth (R² = 0.41,RMSE = 6.67). The depth of charcoal was significantly positively correlated with OLT (R² = 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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.

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.

3.2.1. 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)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

3.2.2. OLTd soils.

3.2.3. Changes in SOC.

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, McGuire et al 2009, 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.

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, 2007, 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.

4.2.1. Carbon balance implications at larger scales.

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.