Soil respiration 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 (the coefficient of variation) and mean-variance power-law 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 20 through remote sensing. Combining such an 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 the 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 the large carbon stocks in permafrost (Tarnocai et al., 2009) and ongoing permafrost thaw (Grosse et al., 2011;Schuur et al., 2013Schuur 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 SR measurement to Published by Copernicus Publications on behalf of the European Geosciences Union.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 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., 2013Schuur et al., , 2015)), vast carbon 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/above ground data -across a 75 × 75 m 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 less than 1 m to greater than 40 m (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.In addition, among-study variation in the length scale of spatial autocorrelation may be partially due to differences in the spatial scale of sampling.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.Instead, we test the following qualitative hypotheses: (i) 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 hypothesize that SR will be predominantly influenced by carbon inputs such that variability in SR throughout the spatial domain will be best explained by tree-stand variables (e.g., basal area); (ii) within the drier, permafrost-free domain we hypothesize that SR will again be best explained by tree-stand variables, but within the wetter permafrost-associated domain SR will be decoupled from carbon inputs and will therefore be poorly explained by treestand variables; and (iii) cold temperatures will fundamentally constrain SR such that we hypothesize a smaller coefficient of variation and less well defined spatial gradients (i.e., weaker spatial structure) in SR under colder temperatures.
Evaluation of the above-summarized hypotheses helps elucidate how a mixture of control processes leads to aggre-gate SR rates.Such knowledge is important for improving model predictions using multi-scale approaches.

Materials and methods
The field site was within the Caribou-Poker Creeks Research Watershed (CPCRW), which comprises a relatively pristine, ∼ 100 km 2 watershed ∼ 45 km north of Fairbanks, Alaska (Fig. 1a).The CPCRW is part of the Bonanza Creek Long-Term Ecological Research site.In July 2014, six parallel transects -each running east to west -were established (Fig. 1b); each was 72 m long and covered a permafrostassociated domain (eastern half of the sampling region), a permafrost-free domain (western half of the sampling region), and a transitional zone between these conditions (thin strip in the center of the sampling region, running south to north) (Fig. 2).In broad terms, the permafrost-associated domain was characterized by black spruce (Picea mariana) and understory mosses, while the permafrost-free domain was characterized by paper birch (Betula papyrifera) and leaf litter (see Fig. 1c for overview of vegetation transition).
To enable SR measurements, 10 cm diameter PVC soil collars were installed to ∼ 5 cm depth along each transect in two phases.First, 12 collars were installed on each transect following a cyclic sampling design aimed at efficient spatial sampling (Burrows et al., 2002): spatial positions were 0,4,12,20,24,32,40,44,52,60,64, and 72 m from the eastern terminus.This arrangement deviates from regular spacing of sampling locations to enable variogram analysis by providing good sampling across all spatial lags (i.e., the distance between sampling points).The design was determined by simulating a broad range of cyclic sampling designs to find a configuration that provided a relatively uniform distribution of sample numbers across all spatial lags.Repeating this sampling design across the 6 transects resulted in 72 soil collars.SR rate measurements were taken from each collar using an EGM-4 (PP Systems, http://ppsystems.com/) with a soil respiration chamber that is designed to provide a well-mixed headspace within the chamber.These initial measurements were used to identify a transition zone across which SR increased rapidly moving from east to west.This transition zone also coincided with the transition from permafrost to permafrost-free conditions.Higher density spatial sampling was performed across this transition zone whereby collars were installed every 1 m from positions 15 to 60 m (moving east to west) on each transect.The full design resulted in 51 collars per transect and a total of 306 collars across the field site (Fig. 1d).
The full set of 306 soil collars was sampled twice from 3 to 13 August 2014 between ∼ 09:30 LT (local time) and 22:00 LT; this sampling period is referred to here as "Summer" because that time of year is usually associated with some of the warmest temperatures of the year.The minimum time between soil collar installation and SR measurement was 2 days.The strong correlation between SR measurements taken from the same soil collars in Summer and Fall (Fig. S1) indicates that the length of time between installation and measurement did not fundamentally change observed SR patterns.Co-located air and soil (5 and 15 cm depths) temperatures were collected at the time of each measurement.For temperature measurements we used an analog thermometer at the time of sampling for each soil collar.For air temperature, we shaded the thermometer to avoid elevating the temperature due to solar inputs.For soil temperature, we placed the thermometer's stem into the ground to the specified depth.Two operators each used a separate EGM-4 and respiration chamber assembly; each operator 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-33 m (moving east to west) along the second transect (from south), 31-39 m along the third transect, 37-45 m along the fourth transect, and 21-29 m along the sixth transect were removed (Fig. 1d).No collars were removed from the first and fifth transects.SR was measured twice again from 10 to 24 September 2014 between ∼ 08:30 LT and 17:30 LT 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 section).
For all SR measurements with the EGM-4 and respiration chamber assembly, the maximum change in CO 2 concentration was set to 100 ppm and the maximum time of data collection was 2 min.The manufacturer (PP Systems) set the measurement interval, which was every 4.8 s.Raw CO 2 concentration and timestamp data were collected and analyzed using custom R scripts.Plotting CO 2 concentration against time commonly showed a small degree of non-linearity.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 regres-sion 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, 2003Muggeo, , 2008)).For each measurement, the rate was taken as either 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 and had a maximum of 24 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 permafrostfree 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 used to calculate a mean and variance.In the Fall some positions did not have 6 soil collars due to soil coring (Fig. 1d); Fall-time 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 20 m of the spatial domain, which represented 54 soil collars (Fig. 1d).Preliminary analyses indicated that the number of soil collars was sufficient for a reliable variance estimate.After quantifying SR variance within the eastern-most 20 m, we added the next most western set of soil collars and recalculated 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 the following aboveground variables that may influence SR.A 150 cm soil probe was used to estimate active layer depth (ALD) to the nearest cm in September 2014; ALD was measured every 2.5 m 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 72 m from the eastern terminus) along each transect.Percent cover of feather moss was visually estimated at each soil collar within a 30 × 30 cm quadrat, with the collar at the center.Tree plots of 5 m radius were established every 10 m along each transect, starting with a plot centered at 0 m 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) greater than 1 cm were identified and their DBHs were measured.The tree plot data were used to estimate the basal area and stem density.
Relating the variance of SR to its mean suggested clustering of SR in space and/or time (see Results section).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.This package has been found to be robust via extensive use across scientific and engineering disciplines, currently with greater than 60 citations in the peer reviewed literature.To enable quantitative comparisons among SR variogram models, we used AutoKrige to fit Matern models (Pardo-Iguzquiza and Chica-Olmo, 2008) and, in turn, generate interpolations.We selected the Matern function because it has been shown to be a robust and flexible model, in general and for soil systems specifically (Stein, 1999;Minasny and McBratney, 2005).For ALD, slope, feather moss cover, tree-stem density, tree basal area, and Summer and Fall soil temperatures the AutoKrige function was used to select the best model among Matern, Gaussian, exponential and spherical structures.We allowed different model structures for these additional variables because we used the fitted variograms only to visualize the spatial structure of these variables (Fig. 2), as opposed to making quantitative comparisons of variogram parameter estimates.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).The Matern model is a special case and also includes the kappa parameter, which describes the smoothness of the fitted variogram (Minasny and McBratney, 2005).
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) re- gression and compared models -using residual 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 permafrostassociated domain, and within the permafrost-free domain.
The SR analyses were repeated separately within the Summer (August) and Fall (September).The SR data from both seasons were also pooled and a 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 use kriging due to some variables having poorly constrained variograms.Prior to GLS analysis the data were log-transformed to improve normality; feather moss percent coverage ranged from 0 to 100 such that a value of 1 was added prior to log-transformation.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 section).Doing so revealed a strong threshold at spatial positions near ∼ 40-45 m (moving east to west).At positions beyond this threshold SR variance increased rapidly and then stabilized at ∼ 55-60 m (Fig. 4a and b).ALD at the western boundary of the sampled domain also showed threshold behavior, increasing rapidly at spatial positions near ∼ 30-35 m and then reached its limit (150 cm) at about 50 m (Fig. 4a and 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 and b; Table 1).Also in both seasons, SR variance increased rapidly beyond an ALD threshold of ∼ 140 cm (Fig. 4c and d).
The variogram fits to SR revealed consistently larger nuggets in the Summer relative to the Fall 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 an especially large range (133 m) and sill (0.94) in the permafrost-free domain during Summer (Table 1).Using the fitted variograms to visualize the 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 the SR spatial structure within the permafrost-associated and permafrost-free domains -using variogram fits 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 was significantly related only to basal area when considering the full domain (Table 2).In the permafrost-free and permafrostassociated domains, multiple variables were significantly related to SR, with feather moss 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 do- 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 h −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 decline in soil respiration near 36 m, moving from west to east.
Table 2. Generalized least squares model fits for Summer soil respiration.Spherical spatial structure was used for all model fits.Models are sorted by their residual standard error (RSE), used as a metric of model fit.Model significance is reported as a p value and regression coefficients associated with p < 0.05 are in bold.Models were fit after normalizing data as standard normal deviates whereby regression coefficients can be directly compared.Entries of "NA" indicate that no model could be fit.main (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 were best explained by soil temperature, but mul-tiple variables had significant relationships with SR for all three domains (Table 4).

Spatial domain
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 and environmental correlates of SR.

The coefficient of variation (CV) 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 to 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 %).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 heterogeneous -increases 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 the 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 and 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 on 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 permafrostfree 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 carbon-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 this 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).

The Taylor power law of soil respiration
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 greater than 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.
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 processbased 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 140 cm.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 140 cm.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).
It should be recognized, however, that we were not able to characterize moisture patterns through space and time such that our interpretations of moisture impacts are speculative.In addition, field observations alone are insufficient to resolve underlying mechanisms leading to the observed thresholding behavior, especially because multiple environmental variables co-varied in space.In particular, multiple variables changed across the permafrost-to-permafrost-free transition such that multiple factors likely contributed simultaneously to the observed thresholding behavior.We strongly encourage future manipulative experiments designed to resolve the governing mechanisms.

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 permafrostassociated spatial domains.There are a number of mechanisms potentially contributing to these observed patterns.For example in the Fall, colder temperatures may place an upper constraint on microbial and root respiration (Lloyd and Taylor, 1994) and senescence of deciduous leaves -which occurred during our Fall sampling -may indicate decreases in root respiration in the deciduous-dominated permafrost-free spatial domain (Lee et al., 2003;Tomotsune et al., 2013;Noh et al., 2016).Both mechanisms could lead to weaker spatial structure and we look forward to future studies that parse their relative contributions, potentially using root-excluding soil collars.The spatial structure of SR in the permafrostassociated domain may be further influenced by high soil moisture placing an upper constraint on microbial respiration in soils with a thin ALD (Elberling et al., 2013).Climate change driven increases in temperature, shifts in forest phenology, 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 a 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 ocwww.biogeosciences.net/14/4341/2017/Biogeosciences, 14, 4341-4354, 2017 curs at ∼ 140 cm (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 feather moss cover were found along the east-to-west SR and ALD boundary (Fig. 2).Understanding how these biotic and abiotic features feed back 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 thresholds 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.
Data availability.Data are available from the corresponding author, and will be associated with the corresponding author's name via the Bonanza Creek LTER website (http://www.lter.uaf.edu/).
Author contributions.All authors contributed to components of the study design.JCS, CGA, BBL, and ARC collected field data.JCS and CGA analyzed the data.JCS drafted the manuscript and all authors contributed to the writing.
Competing interests.The authors declare that they have no conflict of interest.

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 shows the location of field sampling.(b) Satellite image of the field site with orange lines indicating transect locations.(c) Picture of the field site looking west (up-slope); gradients in tree size and the transition from black spruce (green leaves) to paper birch (yellow leaves) can be seen.(d) Spatial layout of soil collars with red points denoting collars that were removed after Summer sampling due to soil coring.

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 150 cm whereby permafrost occurring deeper than 150 cm could not be resolved.

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 section).The dashed line is the linear regression and statistics are provided, including the slope of the line and p value for model significance.

Figure 4 .
Figure 4. Relationships among soil respiration variance (CO 2 flux variance), active layer depth, and space.(a, b) Each soil respiration variance estimate (red plot, 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 plot, right-side axis) is the mean active layer at the western-most spatial position.Lines are spline-fits.(c, d) Soil respiration variance as a function of mean active layer depth -both from panels (a) and (b).The 140 cm threshold is indicated by the vertical dashed green line.The black lines are spline-fits.

Table 1 .
Soil respiration summary statistics, Matern variogram parameters, and sample sizes (n) for each season across the full, permafrostfree, and permafrost-associated spatial domains.

Table 3 .
Generalized least squares model fits for Fall soil respiration.All details as in Table2.

Table 4 .
Generalized least squares model fits for soil respiration across both seasons.The same soil collars were sampled in both seasons so spatial structure could not be introduced into the models.All other details as in Table2.