Soil CO 2 flux across a permafrost transition zone : spatial structure and environmental correlates

Soil respiration is a key ecosystem function whereby shifts in respiration rates can shift systems from carbon sinks 10 to sources. Soil respiration in permafrost-associated systems is particularly important given climate change driven permafrost thaw that leads to significant uncertainty in resulting ecosystem carbon dynamics. Here we characterize the spatial structure and environmental drivers of soil respiration across a permafrost transition zone. We find that soil respiration is characterized by a non-linear threshold that occurs at active layer depths greater than 140cm. We also find that within each season tree basal area is a dominant driver of soil respiration regardless of spatial scale, but only in spatial 15 domains with significant spatial variability in basal area. Our analyses further show that spatial variation and scaling of soil respiration—in our boreal system—are consistent with previous work in other ecosystems (e.g., tropical forests) and in population ecology, respectively. Comparing our results to those in other ecosystems suggests that temporally-stable features such as tree stand structure are often primary drivers of spatial variation in soil respiration. If so, this provides an opportunity to better estimate the magnitude and spatial variation in soil respiration through remote sensing. Combining such an 20 approach with broader knowledge of thresholding behavior—here related to active layer depth—would provide empirical constraints on models aimed at predicting ecosystem responses to ongoing permafrost thaw.


Introduction
Given its central role in global carbon cycling (Raich and Potter, 1995), the flux of CO 2 from soil to the atmospherecollectively referred to as soil respiration (SR)-is heavily studied across a broad range of systems and spatiotemporal scales using methods ranging from small footprint flux chambers to large footprint flux towers (Bond-Lamberty and Thomson, 2010;Wolf et al., 2016).Our understanding of SR patterns is inherently linked to the scales at which measurements are made and we often lack knowledge of how variation in SR changes as we move across scales.To rigorously and mechanistically link SR measurements across scales, it is essential to understand spatial heterogeneity in SR, how spatial heterogeneity changes across scales, and the environmental features that drive those SR patterns.
Knowledge of SR spatial heterogeneity-and underlying drivers of that heterogeneity-in permafrost associated ecosystems is particularly important given large carbon stocks in permafrost (Tarnocai et al., 2009) and ongoing permafrost thaw (Grosse et al., 2011;Schuur et al., 2013;Schuur et al., 2015).Numerous studies have investigated temporal dynamics of SR and most find strong influences of temperature and moisture (e.g., Davidson et al., 1998;Hanson et al., 1993), but much less is known about spatial variation in SR, where temperature and moisture appear to have weaker influences (Yim et al., 2003;Dore et al., 2014;Ohashi et al., 2015;Song et al., 2013).
Studies that examine spatial variation in SR often find rates to be associated with tree-stand variables such as basal area, species composition, and distance from flux measurement to tree stems (Ohashi et al., 2015;Katayama et al., 2009;Soe and Buchmann, 2005;Khomik et al., 2006), though such relationships are not ubiquitous (Song et al., 2013;Saiz et al., 2006;Sotta et al., 2004).Relatively little work on spatially resolved SR has occurred, however, in permafrost-associated systems, though higher SR in permafrost-free domains has been observed (Vogel et al., 2005).Importantly, no study has provided spatially resolved SR estimates moving continuously across a permafrost transition, where the depth of soil that is above permafrost (referred to as active layer depth (ALD)) increases along a spatial gradient to the point that permafrost is lost.It is critical that we fill this knowledge gap due to a combination of ongoing permafrost thaw (Grosse et al., 2011;Schuur et al., 2013;Schuur et al., 2015), vast C stocks currently locked away in permafrost (Tarnocai et al., 2009), and the heavy contribution of SR to ecosystem carbon cycling (Raich and Potter, 1995;Longdoz et al., 2000).
Here we explore how the spatial structure and environmental correlates of SR change as we move across spatial and temporal scales.To do so we generated spatially resolved SR rates-and associated below/aboveground data-across a 75m x 75m spatial domain within boreal Alaska.This domain was selected to capture a spatially abrupt transition from a permafrost-associated system to a permafrost-free system, and at ~6 Landsat pixels, a scale relevant to potential future work using remote sensing.
Previous work in forests has revealed that SR spatial autocorrelation occurs across a broad range of length scales (referred to as the 'range' in variogram models), from <1m to >40m (Foti et al., 2014;Russell and Voroney, 1998;Song et al., 2013;Singh et al., 2008;Rayment and Jarvis, 2000), but none of these estimates come from boreal forests.It is therefore difficult to generate quantitative a priori expectations for the parameters of SR variograms (functions that describe spatial continuity or variability) in our study system.In turn, we test qualitative hypotheses whereby we expect that the length scales of SR spatial autocorrelation (i.e., the variogram range) to be more clearly defined within permafrost-free and permafrostassociated soils, and more difficult to define when simultaneously considering both conditions due to mixing abrupt and continuous environmental transitions.We also hypothesize that under warmer (Summer) conditions spatially structured environmental variables will strongly influence SR, while low temperatures experienced in the Fall season will depress SR across the whole spatial domain, leading to weaker spatial structure.
Given previous work showing that spatial variability in SR is only weakly related to temperature and soil moisture (e.g., Song et al., 2013;Yim et al., 2003), we expect that temperature and soil moisture will explain little variation in SR.We hypothesize instead that SR will be predominantly influenced by carbon inputs such that variability in SR throughout the measured each soil collar once and the order of sampling was randomized.Following the SR measurements a subset of collars were removed to enable soil coring for other purposes: positions 25-33m (moving East to West) along the second transect (moving South to North), 31-39m along the third transect, 37-45m along the fourth transect, and 21-29m along the sixth transect were removed (Fig. 1d).No collars were removed from the first and fifth transects.SR was measured twice again from September 10-24, 2014 between ~8:30am and 5:30pm in the remaining collars, along with air and soil temperatures; this sampling period is referred to here as 'Fall' because that time of year is usually associated with rapidly declining temperatures and leaf-drop from deciduous trees.Indeed, air temperatures measured during the Fall sampling period were significantly lower than those measured during the Summer sampling period.
It is important to note that analyses that directly compare results between Summer and Fall were repeated on the SR data using only locations that were sampled in both seasons.The purpose of repeating the analyses was to evaluate whether observed seasonal differences were due to differences (between seasons) in the sampling scheme; this showed that seasonal differences in SR patterns were not due to differences in the sampling scheme (see Results).
For all SR measurements with the EGM-4/Respiration chamber assembly the maximum change in CO 2 concentration was set to 100ppm and the maximum time of data collection was 2 minutes.The manufacturer (PP systems) set the measurement interval, which was every 4.8 seconds.Raw CO 2 concentration and timestamp data were collected and analyzed using custom R scripts.Plotting CO 2 concentration through time commonly showed a small degree of nonlinearity.The most pronounced non-linearity was usually between the first and second data points such that removing the first data point provided linear fits in nearly all cases.Other data excursions were noted, such as an abrupt increase in CO 2 concentration with linear trajectories on either side.To generalize the data analyses we (1) dropped the first data point from each sample, (2) fit both a linear regression and a segmented regression to the remaining data to estimate the SR rate.
Segmented regressions were estimated automatically using the function 'segmented' in R package 'segmented' (version 0.5-1.4)(Muggeo, 2003(Muggeo, , 2008)).For each measurement, the rate was taken as either as the slope of the linear regression or the slope of the first segment in the segmented regression, whichever had a higher R 2 .The only exception was when the first segment contained fewer than 4 data points; in those cases the slope from the full linear regression was used whereby all regressions contained at least 4 data points.The slope from the first segment was used instead of the second segment due to noticeable decreases in the SR rate with increasing sampling time, which was likely due to increasing CO 2 concentrations in the recirculating system.This approach significantly improved model fits, relative to using rates recorded directly from the EGM-4 instrument.For each collar in each season the duplicate SR estimates were averaged.
To facilitate comparisons to other studies investigating spatial structure in SR we quantified the coefficient of variation (CV) within the whole spatial domain, within the permafrost-associated domain, and within the permafrost-free domain; the CV was estimated in these spatial domains within each season.As a complementary approach, we estimated the slope of the power function-on log axes-relating variance in SR to mean SR.This is analogous to the classic 'Taylor power law' that relates variance in population density-of a given biological species-to mean population density (Taylor, 1961;Taylor and Taylor, 1977;Eisler et al., 2008).Within each season, the 6 soil collars at a given East/West position were Biogeosciences Discuss., doi:10.5194/bg-2016Discuss., doi:10.5194/bg- -481, 2016 Manuscript under review for journal Biogeosciences Published: 29 November 2016 c Author(s) 2016.CC-BY 3.0 License.used to calculate a mean and variance.In the Fall some positions did not have 6 soil collars due to soil coring (Fig. 1d); Falltime SR from these East/West positions were not used.To cover a broad range of mean SR conditions, the two seasons were examined together.Spatial and temporal structures were therefore examined simultaneously.
To examine how SR spatial variation might respond to the spatial advance of permafrost thaw, we iteratively quantified SR variance across spatial domains that had an increasing contribution of permafrost-free soils.Within each season, we began by quantifying SR variance using all soil collars within the Eastern-most 20m of the spatial domain, which represented 54 soil collars (Fig. 1d).Preliminary analyses indicated that number of soil collars was sufficient for a reliable variance estimate.After quantifying SR variance within the Eastern-most 20m, we added the next most Western set of soil collars and re-calculated SR variance.This procedure was repeated until reaching the far Western edge of the field domain.
At each step we also estimated mean ALD at the soil collars found at the Western edge of the domain used for variance estimation.This allowed us to connect changes in SR spatial variation to changes in ALD.
In addition to SR, we collected data on aboveground variables that may influence SR, as follows.A 150cm soil probe was used to estimate active layer depth (ALD) to the nearest cm in September 2014; ALD was measured every 2.5m along all 6 transects.A clinometer was used to estimate % slope at sampling locations (0,4,12,20,24,32,40,44,52,60,64, and 72m from the Eastern terminus) along each transect.Percent cover of feathermoss was visually estimated at each soil collar within a 30cm x 30cm quadrat, with the collar at the center.Tree plots of 5m diameter were established every 10m along each transect, starting with a plot centered at 0m at the Eastern terminus of each transect.The tree plots were hemispherical with sampling occurring to the South of each transect.This was done to minimize disturbance to soil collars, which were placed on the North side of each transect.Within each tree plot all woody stems with diameter at breast height (DBH) > 1cm were identified and their DBHs were measured.The tree plot data were used to estimate basal area and stem density.
Relating the variance of SR to its mean suggested clustering of SR in space and/or time (see Results).In turn, we conducted formal spatial analyses using the 'AutoKrige' function in R package 'automap' (version 1.1-14) (Hiemstra et al., 2009), which uses model selection to automatically fit a variogram model (Matheron, 1963) and, in turn, generates a spatial interpolation.The algorithm chooses the best fit among Gaussian, exponential, spherical, and Matern models (Pardo-Iguzquiza and Chica-Olmo, 2008).This approach was used to spatially interpolate and visualize each data type; SR, ALD, slope, feathermoss cover, tree stem density, tree basal area, and Summer and Fall soil temperatures were all interpolated (see Fig. 2 for interpolations of all but SR).We also examined parameter values of the models fit to the SR variograms.These parameters included the nugget, range, and sill.The variogram nugget represents processes that occur at spatial scales smaller than the minimum distance between samples as well as sampling error.The variogram range is an estimate of the maximum spatial distance across which there is spatial autocorrelation.The variogram sill is the variance at which the variogram model levels off (Western et al., 1998).
We further aimed to explain variation in SR using the other measured variables.Given significant spatial autocorrelation in all data types, we used generalized least squares (GLS) regression and compared models-using residual The SR analyses were repeated separately within the Summer (August) and Fall (September).The SR data from both seasons were also pooled and GLS regression was used without spatial structure to find the variable that best explained this spatiotemporal set of SR rates; spatial structure could not be used on the pooled set because measurements were made in both seasons at the same spatial locations.
Running GLS analyses required an estimate of each explanatory variable at the spatial location of each soil collar.Some data types-ALD, slope, stem density, and tree basal area-were collected at a coarser spatial grain than the soil collars.These variables were therefore linearly interpolated to provide an estimate of each at each soil collar; we did not using kriging due to some variables having poorly constrained variograms.Prior to GLS analysis the data were transformed to improve normality; feathermoss percent coverage ranged from 0 to 100 such that a value of 1 was added prior to logtransformation.To allow direct comparisons of regression coefficients across explanatory variables, log-transformed data were further normalized to have a mean of 0 and variance of 1.Initial analyses indicated that a spherical spatial structure generally lead to lower root square errors in fitted models, relative to Gaussian and exponential spatial structures.GLS models were thus fit using spherical spatial structure and models with different explanatory variables were competed against each other by comparing their root square errors.All analyses were performed with R version 3.2.1 (R Development Core Team 2016).

Results
SR measured in a soil collar in Summer was strongly correlated to SR measured in the same soil collar in the Fall (Fig. S1).
There was, however, a clear decline in both the mean and variance of SR from Summer to Fall (Table 1; Fig. S1).Across the full domain the two seasons had similar CVs (~47-50%), but across the two seasons the permafrost-associated domain CV was ~10-18% higher than the permafrost-free domain (Table 1).These results were maintained when Summer data were analyzed using the same sampling scheme used in the Fall (Table S1).The slope of the power function relating SR variance to mean SR was 1.57, with an R 2 of 0.51(Fig.3).
We plotted SR variance within a given spatial domain against the position of the Western boundary of that sampled domain (see Methods).Doing so revealed a strong threshold at ~40-45m (moving East to West), beyond which SR variance increased rapidly and then stabilized at ~55-60m (Fig. 4a,b).ALD at the Western boundary of the sampled domain also showed threshold behavior, increasing rapidly at ~30-35m and then reached its maximum value (150cm) at about 50m (Fig. 4a,b).These patterns in SR variation and ALD were found in both seasons even though the SR variance was much lower in the Fall (cf.Fig. 4a,b; Table 1).Also in both seasons, SR variance increased rapidly beyond a threshold of ~140cm ALD (Fig. 4c,d).
Variograms fit to SR revealed consistently larger nuggets in the Summer relative to the Fall-consistent with the CV results above-and modest differences in the nuggets between permafrost-associated and permafrost-free domains (Table 1, Table S1).The range and sill were both larger in the Summer, with especially large range (133m) and sill (0.94) in the permafrost-free domain during Summer (Table 1).Using the fitted variograms to visualize SR spatial structure showed elevated SR in the Summer in the permafrost-free domain, and a similar-but less pronounced-pattern in the Fall (Fig. 5).
Visualization of SR spatial structure within the permafrost-associated and permafrost-free domains-using variograms fit within those domains (Table 1)-confirmed a peak in SR near the center of the permafrost-free domain, an increase in SR from South to North in the permafrost-associated domain, and weak spatial structure across both domains in the Fall (Fig. 5).
These patterns were maintained when analyzing Summer data using the reduced sampling scheme used for the Fall (Fig. S2).
GLS regression revealed that Summer SR significantly related only to basal area when considering the full domain (Table 2).In the permafrost-free and permafrost-associated domains, multiple variables were significantly related to SR, with feathermoss percent cover and basal area having the strongest relationships, respectively (Table 2).Similar patterns were observed when analyzing Summer data using the reduced sampling scheme used for the Fall (Table S1).Fall SR was also significantly related only to basal area when considering the full domain, and was again best explained by basal area in the permafrost-associated domain (Table 3).In contrast to the Summer, Fall SR in the permafrost-free domain was best explained by active layer depth (Table 3).When considering both seasons together, SR in the full, permafrost-associated, and permafrost-free spatial domains was best explained by soil temperature, but multiple variables had significant relationships with SR for all three domains (Table 4).

Discussion
Through intensive field sampling and geospatial analyses we have revealed the spatial structure of SR across a boreal forest permafrost transition, shifts in that spatial structure between seasons and spatial scales, and the environmental variables that potentially drive spatial variation in SR within and across permafrost conditions.To the best of our knowledge this is the first direct characterization of SR spatial structure across a permafrost transition.A number of studies have characterized SR spatial structure in other ecosystems, however, and patterns revealed here show some commonality with previous work in terms of the level of variation in and environmental correlates of SR.

The coefficient of variation of soil respiration
A common approach for characterizing spatial heterogeneity in SR is through estimation of the coefficient of variation (CV) (e.g., Shi et al., 2016;Song et al., 2013;Dore et al., 2014).Our CV analyses suggest that at larger spatial scales-that traverse shifts in tree species and permafrost conditions-SR spatial heterogeneity is governed by temporally-stable environmental variables.The CV of the full (i.e., larger scale) domain was similar between seasons (~47-50%) despite declines in the mean and variance of SR from Summer to Fall (Table 1).Stability of the CV observed here is consistent with a recent study in a monsoon-impacted mixed forest in China, where the CV was 55% and 51% in the Summer and Fall, respectively (Shi et al., 2016).In the same study, the CV in the Spring was lower (42%), however, and their analyses pointed to important influences of both dynamic (soil moisture) and temporally stable (soil bulk density and tree size) environmental features.Here, we do not have SR during Spring conditions, which in our study system would correspond to snow melt and the initial stages of active layer thaw.Characterizing the SR spatial heterogeneity under such conditions may reveal more temporal dynamics in the CV.Nonetheless, our data indicate that once the active layer has thawed, larger-scale SR spatial heterogeneity is governed by environmental conditions that have consistent spatial structure between seasons (e.g., tree stand structure).
When looking across spatial and temporal scales, we found significant shifts in the CV, and the range of CV values was similar to that found in previous work from a broad range of ecosystems; across temperate forests, tropical forests, and agricultural systems the SR CV generally ranges from ~25-60% (Shi et al., 2016;Song et al., 2013;Dore et al., 2014).In both seasons the permafrost-associated domain had a higher CV (>45%) than the permafrost-free domain (~37%), and there was a modest increase in the permafrost-associated domain CV from Summer (46%) to Fall (56%).We therefore infer that temporally-stable conditions are the primary control over SR spatial heterogeneity in both the permafrost-free and permafrost-dominated domains.This is consistent with work in other systems showing important influences of tree stand structure, and in particular the distance from large trees (Ohashi et al., 2015;Katayama et al., 2009;Soe and Buchmann, 2005;Khomik et al., 2006).In the permafrost-dominated domain there may, however, be a small influence of environmental variables characterized by temporally-shifting spatial structure.One possibility is small-but spatially heterogeneousincreases in ALD from Summer to Fall.
Our results further indicate that the 10-20% higher SR CV in the permafrost-dominated domain, relative to the permafrost-free domain, was due to higher spatial variability in tree basal area in the permafrost-dominated domain.This inference is supported by two observations.First, SR in the permafrost-dominated domain was best explained by tree basal area in both seasons, while variables explaining SR in the permafrost-free domain changed across seasons and were more diffuse (i.e., there were more significant variables) (Tables 2,3).Second, the CV of basal area was much higher in the permafrost-dominated domain; examining the linearly interpolated basal area estimates that align with SR measurements revealed 45% and 80% CV for basal area in the permafrost-free and permafrost-dominated domains, respectively.These results suggest that in the permafrost-dominated domain, basal area in the vicinity of a collar-scale measurement impacts the SR rate and spatial variability of basal area at the plot-scale is a primary driver of the SR CV.These results-combined with previous work in other systems (e.g., Ohashi et al., 2015)-further support an important role of tree stand structure for SR spatial variability across biomes.It would, however, be useful to study SR CV in permafrost-free domains characterized by higher basal area CV than observed here.Doing so would help determine if there is a general threshold of basal area spatial variation above which that variation becomes an important driver of the SR CV.If a general threshold exists, it could be incorporated into ecosystem simulation models aimed at predicting C-cycling responses to environmental change.
Furthermore, if spatial variation in SR is broadly influenced by tree-stand structure, it suggests an opportunity to use remote sensing techniques to characterize stand structure (van Leeuwen and Nieuwenhuis, 2010) and, in turn, the spatial structure of SR in boreal and other high-latitude ecosystems (Kushida et al., 2004).Taylor power law analysis has been broadly used in population ecology to relate variance to mean, and extending this approach to spatial and temporal variation in SR revealed consistency between the scaling of SR and the scaling of population density.More specifically, we found a slope of the power function relating SR variance to mean SR of 1.57, which is consistent with slopes found in population ecology (Taylor, 1961;Reed and Hobbs, 2004).Slopes >1 have been interpreted as indicating aggregation or clustering (Taylor, 1984), whereby a slope of 1.57 in our dataset indicates that SR was clustered in space and/or time.

Taylor power law of soil respiration
The R 2 of the power function was only 0.51, lower than many population ecology studies.This was likely due to the small number of soil collars (6) used in each variance/mean estimate.Using more soil collars for each variance/mean estimate would generate more accurate estimates, but the number of independent estimates would drop considerably.The approach we implemented strikes a balance between accuracy and number of independent estimates (Rodeghiero and Cescatti, 2008), and our results should be interpreted as a preliminary characterization of variance-mean scaling of SR.
Nonetheless, empirically linking mean and variance may provide an opportunity to predict mean SR through process-based models and, in turn, generate appropriate spatial projections of variation in SR.For ecosystem-scale models that encapsulate non-linear or threshold-based processes, appropriately capturing spatial variation in SR could significantly improve model predictions under conditions of environmental change.

Variance in SR related to active layer depth (ALD)
Relating SR variance to ALD revealed strong thresholding behavior whereby the variance of SR increased sharply above an ALD of 140cm.This threshold was observed in both seasons suggesting underlying mechanisms that are stable through time.One potential explanation is that SR associated with thinner ALDs is constrained by relatively high soil moisture-likely due to facilitation of anaerobic conditions (Drew and Lynch, 1980)-associated with thin ALD; a small number of soil cores in the permafrost-free and permafrost-dominated domains showed higher soil moisture in the permafrost-dominated domain (~0.58 vs 0.82 gravimetric moisture content on a wet mass basis, respectively) (see also Bonan and Shugart, 1989).This is consistent with our observation that SR variance (not CV) was lower in the permafrost-dominated region (Table 1).As ALD increases and soil moisture decreases, the spatially homogeneous constraint placed on SR by high moisture should diminish, thereby allowing other-more spatially heterogeneous-mechanisms to govern SR rates.As a consequence, SR spatial variance would increase.If true, it appears that in our system this release from high-moisture-limitation occurred abruptly at an ALD of 140cm.This threshold behavior and conceptual interpretation is consistent with the representation of moisture impacts on SR in ecosystem models such as the Community Land Model (Lawrence et al., 2011), and highlights the critical need of accurately modeling changes in system hydrology as permafrost thaws (Elberling et al., 2013).

SR variograms and spatial interpolations
Variogram analyses of SR revealed shifts in the degree of small-scale heterogeneity and in the range of spatial autocorrelation across seasons and spatial scales.The variogram nugget reflects measurement error and spatial variation at distances shorter than sampled.Assuming consistent measurement error-the same instruments, operators, and protocols were used throughout-the decrease in the nugget from Summer to Fall indicates a decline in small-scale heterogeneity.A potential explanation is that lower temperatures in the Fall broadly constrained SR (Bond-Lamberty and Thomson, 2010;Davidson et al., 1998) and, in turn, limited the impact of small-scale variation in other influential environmental variables.The scale of spatial autocorrelation-as measured by the variogram range-was also much lower in the Fall and in the permafrost-associated domain during the Summer.These patterns indicate that SR has much weaker spatial structure during colder periods and across permafrost-associated spatial domains.We suggest this is due to colder temperature placing an upper constraint on microbial and root respiration during the Fall, and higher soil moisture placing an upper constraint on microbial respiration in soils with a thin ALD.Climate change driven increases in temperature and decreases in permafrost may therefore increase both small-scale heterogeneity and larger-scale spatial autocorrelation in SR.
Visualizing the spatial structure of SR-based on fitted variograms-revealed a sharp boundary in the middle of the field domain across which SR increased non-linearly moving from East to West.This pattern was most apparent in the Summer (Fig. 5), consistent with our interpretation of weaker spatial structure in Fall.This East-to-West boundary in SR was co-located with rapid thickening of the active layer, which further supports our inference that non-linear thresholds in SR should be expected as ALD thickens due to ongoing permafrost thaw.In our system this threshold occurs at ~140cm (see above); knowledge of among-system variation in the presence and magnitude of such a threshold would help constrain ecosystem models that attempt to predict elemental cycling in permafrost-associated systems under future climate scenarios.
In addition, improved understanding of the mechanisms leading to threshold behavior is needed for improved process representation in such models.For example, the highest tree stem density, soil temperatures, and feathermoss cover were found along the East-to-West SR and ALD boundary (Fig. 2).Understanding how these biotic and abiotic features feedback on each other is critical for robust predictions of ecosystem responses to permafrost thaw.

Conclusions
Characterizing the spatial structure of SR-not just spatial variation-provides a number of unique opportunities to improve our basic understanding of processes that control the carbon balance in ecosystems and to reveal key threshold in ecosystem function that are likely to arise in response to environmental change.While such efforts are relatively rare due (presumably) to logistic challenges, they provide insights that other approaches cannot.Here, for example, we revealed thresholding behavior both in terms of average SR rates (i.e., spatial discontinuities shown in Fig. 5) and spatial variation in SR (i.e., the steep increases shown in Fig. 4).These patterns provide powerful constraints on simulation models that attempt to predict carbon cycling fluxes across spatial environmental gradients and/or through time as environmental conditions change.A future challenge is therefore to couple highly-resolved spatially explicit SR measurements-like those shown here-with dynamic ecosystem models.Doing so has significant potential to reveal important model limitations that, if overcome, may facilitate improved predictions of ecosystem function under environmental change.
Biogeosciences Discuss., doi:10.5194/bg-2016-481,2016 Manuscript under review for journal Biogeosciences Published: 29 November 2016 c Author(s) 2016.CC-BY 3.0 License.standard errors-with Gaussian, exponential, or spherical spatial structure.Variogram and GLS analyses were used to study SR data across the entire spatial domain, within the permafrost-associated domain, and within the permafrost-free domain.

Figure 1 :
Figure 1: Overview of field site and spatial sampling design.The background map shows the distribution of permafrost across Alaska and indicates the approximate location of the field site.(a) Satellite image showing CPCRW domain outlined in yellow; the orange box 5

Figure 2 :
Figure 2: Interpolated environmental variables.Each panel has its own color ramp, and units are noted.The vertical and horizontal axes in all panels indicate the South-to-North and West-to-East dimensions, respectively.Interpolations are based on fitted variograms.Active layer depth of 150+ indicates that the probe had a maximum depth of 150cm whereby permafrost occurring deeper than 150cm 5

Figure 3 .
Figure 3. Mean-variance scaling of soil respiration.Each data point is derived from six soil collars used to estimate the variance and mean of soil respiration (see Methods).The dashed line is the linear regression; statistics are provided, including the slope of the line.

Figure 4 .
Figure 4. Relationships among soil respiration variance, active layer depth, and space.(a,b) Each soil respiration variance estimate (red; left-side axis) is based on data from all soil collars within a spatial domain that extends from the field site's Eastern boundary to the indicated Western-most spatial position (horizontal axis).Each active layer depth estimate (blue; right-side axis) is the mean active layer at 5

Figure 5 .
Figure 5. Interpolated soil respiration across seasons and spatial scales.Each panel has its own color ramp, and in all panels the values are CO 2 flux rates (g m -2 hr -1 ).The vertical and horizontal axes in all panels indicate the South-to-North and West-to-East dimensions, respectively.The top panels show interpolations across the full spatial domain and the bottom panels show interpolations within the permafrost-free (P-F) and permafrost-associated (P-A) domains.Stronger spatial structure in Summer is evident as is the non-linear 5

Table 1 . Soil respiration summary statistics and best-fit variogram model with associated parameters, for each season across the full, permafrost-free, and permafrost-associated spatial domains.
Biogeosciences Discuss., doi:10.5194/bg-2016-481,2016 Manuscript under review for journal Biogeosciences Published: 29 November 2016 c Author(s) 2016.CC-BY 3.0 License.