Scaling impacts on environmental controls and spatial heterogeneity of soil organic carbon stocks

Introduction Conclusions References


Introduction
Soil organic carbon (SOC) stocks have large spatial heterogeneity globally (Todd-Brown et al., 2013) and in the northern circumpolar permafrost region (Ping, 2013;Hugelius et al., 2014).Achieving an accurate representation of the spatial heterogeneity of existing SOC stocks in Earth system models (ESMs) is a prerequisite for predicting future carbon-climate feedbacks (Todd-Brown et al., 2013), given that uncertainty resulting from the distribution of SOC stocks accounts for a large proportion of the overall uncertainty (Burke et al., 2012) in predictions of future greenhouse gas concentrations and associated climate changes.Reliable estimates of re-Published by Copernicus Publications on behalf of the European Geosciences Union.
U. Mishra and W. J. Riley: Scaling impacts on environmental controls and spatial heterogeneity gional SOC stocks and their spatial heterogeneity are also essential to gaining a better understanding of existing environmental controls of SOC stocks and their vulnerability to changing climate (Johnson et al., 2011).
At present, large differences exist between observationbased and ESM-based SOC stock estimates at high latitudes (Mishra et al., 2013).Several factors likely account for these differences (Koven et al., 2013;Todd-Brown et al., 2013); the relevant factors for the current analyses are that the natural spatial variability of some environmental controls observed to affect SOC stocks are not represented in the models.One potential approach to account for spatial heterogeneity in ESMs is to incorporate the impacts of environmental controls on SOC stocks at a spatial scale consistent with the observations.
Among environmental variables, climate, topography, organisms, parent material, time, and spatial position (Jenny, 1941;McBratney et al., 2003) have widely been used to predict the spatial variability of soil properties from plot to regional scales (Thompson and Kolka;2005;Vasques et al., 2012;Minasny et al., 2013;Mishra and Riley, 2014).At high latitudes where solar radiation and evaporation are strongly related to the geometry of the land surface (McKenzie et al., 2000), topography has been reported to play an important role in determining the level of SOC at a specific location (Jenny, 1980;Birkeland, 1984;Hobbie et al., 2000).Hill slope processes, such as erosion and deposition, are related to terrain attributes, such as the soil wetness index (SWI) and sediment transport index (STI).Both the SWI and STI are strong predictors of SOC stocks and have been used extensively in digital mapping of SOC stocks at various scales (e.g., Zhu et al., 2001;Thompson et al., 2006;Adhikari et al., 2013).
Despite several global and regional studies on SOC inventories, quantitative data on the relationship between SOC stocks and the environmental controlling factors do not exist at multiple scales (Lal, 2004;Torn et al., 2009;Chaplot et al., 2010).In the spatial prediction of SOC concentration and stocks, different primary and secondary topographic attributes are calculated from a digital elevation model (DEM) and used, along with other soil-forming factors such as climate, land cover type, and parent material, to predict SOC stocks across different scales (Hancock et al., 2010;Vasques et al., 2010;Kumar et al., 2012).Usually, fine-scale environmental data are preferred over coarser spatial scales as they are more representative of natural landscapes.However, it has also been reported that the highest-scale DEMs do not always produce highest accuracy in predicting soil properties and that knowledge of the appropriate DEM scale for a particular landscape is important (Smith et al., 2006;Roecker and Thompson, 2010).Calculated values of topographic indices of a study area differed depending on the spatial scale of the DEM used (Wilson and Gallant, 2000;Thompson et al., 2001;Kienzle, 2004).In a recent study, Vasques et al. (2012) showed that the quality (R 2 ) of predicted carbon stocks consistently decreased with increasing grid size of the environmental data used by the models.
Information collected at one scale can be used to infer properties at either smaller or larger scales (Beven, 1995), such as when using point observations to estimate regional SOC stocks, or taking areal averages of SOC stocks from ESMs and disaggregating them.The essence of successful scaling is to infer the key patterns from information collected at one scale and use those patterns to make inferences at another scale that maintain consistency in the desired metrics across the scales (Western et al., 2002).We believe that the statistical spatial structure of SOC stocks may be useful in (1) downscaling and inferring fine-scale variability from coarse-scale ESM predictions; (2) model benchmarking; and (3) developing reduced order model formulations relating biogeochemical processes across scales.By model benchmarking, we mean using observation-based scaling relationships to test the land biogeochemical representations in Earth system models.
A large body of literature exists on spatial scaling of soil moisture, which is a dominant control on SOC dynamics and stocks.Several studies have modeled the spatial variability of soil moisture patterns from relatively smaller to larger scales, and have attempted to characterize the spatial structure based on system properties (Rodriguez-Iturbe et al., 1995;Peters-Lidard et al., 2001;Western et al., 2002;Riley and Shen, 2014;Pau et al., 2014).For some systems, soil moisture spatial variance follows a power law decay as a function of spatial area (Manfreda et al., 2007), although in other systems there are clear scale breaks in this relationship (e.g., Das and Mohanty, 2008;Joshi and Mohanty, 2010;Pau et al., 2014).Gebremichael et al. (2009) reported that, in a watershed located in the southern Great Plains of the USA, soil moisture showed scale invariance, and that if the scaling parameters could be estimated from large-scale soil moisture fields, it might be feasible to transform spatial soil moisture fields between scales.Despite the progress made in modeling the scaling behavior of soil moisture, to our knowledge no study exists that has examined the statistical structure of the scaling behavior of SOC stocks.
Scaling approaches might also improve prediction of the spatial structure of SOC stocks.The variogram, which relates the variance of SOC stock differences between two points to the distance between the two points (Webster and Oliver, 2007), is one approach to representing the spatial structure of SOC stocks (Western et al., 2002).The main structural parameters of a variogram are the nugget, sill, and range (Western and Bloschl, 1999;Webster and Oliver, 2007).A nugget shows the unexplained portion of the variance and the sill is the level at which a variogram flattens and therefore represents the total variance of the variable.The nugget-to-sill ratio can be used to quantify the percentage of overall variance found at a distance smaller than the smallest lag interval, and to classify the spatial dependence of soil properties (Cambardella et al., 1994;Mora-Vallejo et al., 2008).
Current ESMs operate at spatial scales larger than 50 km and use a nested subgrid hierarchy approach to represent the land surface heterogeneity (Lawrence et al., 2012).ESM spatial scales are becoming increasingly finer to more accurately represent the localized features of Earth systems that affect energy, water, and greenhouse gas fluxes.For example, nextgeneration ESMs such as those being developed by the US Department of Energy (Accelerated Climate Modeling for Energy; ACME) will have a spatial scale for the land model of ∼ 10 km (Bader et al., 2014).In this study, we successively increased the spatial scale (from s = 50 m to 10 km) of environmental variables and used observations and geospatial approaches to predict SOC stocks.Throughout this paper, we refer to the "scale" (s) as either the distance across which properties are assumed to be homogeneous or the square root of the pixel area satisfying that criteria, and note that the terms "scale" and "resolution" are often interchangeable in this context.The idea was to model the change in predicted heterogeneity of SOC stocks resulting from changing scale of the environmental data.In this study, by the term "scaling" we mean the transfer of information about environmental controls and statistical properties of SOC stocks from one scale to another.We also investigated how environmental controls on SOC stocks and the spatial structure of SOC stocks (the correlation length (range), total variance (sill), and unstructured variability (nugget)) changed as a result of change in spatial scale of environmental variables.

Soil organic carbon observations and environmental variables
We used 556 georeferenced soil profile observations by combining data from two recently published studies of soil organic carbon stocks across Alaska (472 soil profiles from Mishra and Riley (2012) and 84 additional samples from a recently published database of Michaelson et al., 2013).This database combined the US Department of Agriculture National Resource Conservation Service's pedons with the soil pedon observations collected from the University of Alaska Fairbanks under its northern latitude soil program (Michaelson et al., 2013).This database included measured representative soil profiles from Alaska and covered all soil types at the soil suborder level (18 suborders) and all 27 major land resource areas of Alaska.The observed SOC stocks showed a wide range (0.35-296 kg m −2 ) due to varying sampling depths (0.05-4.5 m depth), various soil types with differing carbon content levels (Inceptisols to Gelisols), and the large environmental heterogeneity present in the study area.About 50 % of the observed samples were from the 0.05-1 m depth interval, 41 % of the samples were from the 1-2 m depth interval, and the remaining samples (9 %) were from the 2-4.5 m depth interval.A large number of observations had SOC stocks ≤ 50 kg m −2 , and about 49 samples had SOC stocks greater than 100 kg m −2 .A majority of pedons with high SOC stocks were located in the Arctic region of Alaska, which contains areas with soils impacted by cryoturbation, and the majority of pedons with lower SOC stocks were located in interboreal and coastal rainforests.The SOC stock data has a unimodal (kurtosis, K = 9.4) positively skewed distribution (coefficient of skewness, S = 2), and is unevenly distributed throughout Alaska (Fig. 1).The fine-scale (spatial scale, s = 50 m) spatially heterogeneous SOC stock estimates have uncertainty, as described by the prediction errors in Mishra and Riley (2012).However, that study indicated that the SOC stock estimates have reasonable fidelity with the observations withheld from the estimation procedure, and therefore we believe are sufficiently accurate to be used in the scaling analyses described here.As more data become available, the approaches described here to evaluate spatial scaling controls can be readily re-applied.
We collected environmental variables of the study area from various sources.The topographic attributes were derived from a digital elevation model (DEM) with a 50 m spatial scale (s = 50 m) obtained from the US Geological Survey (USGS) database (Gesch et al., 2009).From the DEM, 13 topographic attributes were calculated that have been reported to be useful predictors of the SOC stocks across different environmental conditions using the Spatial Analyst function of ArcGIS (ArcGIS version 10.2, Environmental Systems Research Institute, Inc., Redlands, CA, USA).
Climate data for long-term ) average annual temperature, precipitation, potential evapotranspiration (PET), and summer shortwave radiation were obtained from the Scenario Network for Alaska and Arctic Planning (SNAP, 2014) database, with a spatial scale of s = 2 km.
Land cover data of 60 m spatial scale (s = 60 m) was extracted for Alaska from the National Land Cover Database (NLCD) database (Homer et al., 2007).In this study, to reduce the number of categorical variables with potentially redundant information, we reclassified the 19 NLCD land cover types found in Alaska into 9 major categories.Thus, developed open space, low intensity, medium intensity, and  high intensity were classified as a developed category; deciduous, evergreen, and mixed forests were classified as a forest category; dwarf scrub and shrub scrub were classified as a scrub category; shrub, sedge, or moss were classified as a herbaceous category; pasture and cultivated lands were classified as a cultivated category; and woody and herbaceous wetlands were classified as a wetland category.In the reclassified map, the category with the largest land area was the scrub category (43 %), followed by forest (25 %), barren (8.5 %), herbaceous (7 %), and wetlands (7 %).The remaining surface area (9.5 %) was under open water, perennial ice, barren land, and moss vegetation.Indicator variables for the presence or absence of seven land cover types (except open water and perennial ice) were created and used in the model selection process.
The surficial geology data of Alaska was obtained from a USGS database (Karlstrom, 1964) in which the entire state of Alaska was classified into 27 different types of surficial geology.The category with the largest land area was mountain alluvium and colluvium (22.5 %), followed by coarse rubble (19 %), coastal alluvium (8.5 %), glacial moraine (7 %), and undifferentiated mosaics (6 %).The remaining land area (37 %) was placed in the remaining 22 surficial geology types.

Spatial scaling and environmental controls
We investigated changes in environmental controls on SOC stocks at spatial scales of 100, 200, and 500 m and 1, 2, 5, and 10 km.For this purpose, the DEM of a 50 m spatial scale was successively aggregated to 100, 200, and 500 m and 1, 2, 5, and 10 km using the bilinear convolution algorithm of ArcGIS 10.2 (Fig. 2).Topographic attributes were calculated at each spatial scale and included in predicting SOC stocks at specific scales.The land cover types and surficial geology data were aggregated at similar scales using the nearestneighbors algorithm of ArcGIS 10.2 and were included in predicting SOC stocks at various spatial scales.
Median values of the geographically weighted regression coefficients (β) of log-transformed (natural log) SOC stock and environmental variables (that were significant at all spatial scales) were calculated for each spatial scale.Log trans-formations of SOC stocks were applied to meet the linear regression assumptions of normality and constant variance of errors.The change in the strength of environmental controls (β values) on SOC stocks resulting from a change in spatial scale was investigated by plotting the median values against the spatial scale at which they were calculated.Several mathematical functions were used to represent the statistical nature of scaling of the environmental controls on SOC stocks.

Selection of environmental predictors, prediction of SOC stocks, and variogram modeling
We used best subset regression in SAS 9.3 (SAS Institute, 2011) to generate all possible combinations of environmental variables to predict the SOC stocks across Alaska.We used Mallows' Cp criterion to select three candidate linear models for each spatial scale for which Cp values were close to the number of predictors (Kutner et al., 2004;p. 358), from which we subsequently selected one linear model for each spatial scale with uncorrelated and statistically significant environmental predictors.The presence of multicollinearity among selected environmental predictors was tested by calculating the variance influence factors (VIFs) for each of the selected variables.The VIFs for all the variables included in models selected at different spatial resolutions were < 5. High levels of multicollinearity (VIF > 10) can falsely inflate the least squares estimates (Kutner et al., 2004;O'Brien, 2007;Gomez et al., 2013).The selected environmental predictors were then used in a geographically weighted regression approach (Mishra andRiley, 2012, 2014) to predict the whole profile SOC stocks across Alaska at different spatial scales using where C is the SOC stock at a certain location, u i v i are the geographical coordinates, X k is the environmental predictor, p is the number of independent variables, β o and β k are the geographically weighted regression coefficients, and ε is the error term.We used adaptive kernel bandwidth in this study given that the sample density varied over the study area.The optimal bandwidth was determined by minimizing the corrected Akaike Information Criterion (AICc) as described in Fotheringham et al. (2002).
We calculated several statistical parameters, including mean (µ), variance (σ 2 ), skewness (S), and kurtosis (K) from the predicted SOC stocks at each specific spatial scale.The predicted variances were plotted against spatial scale and a best-fit mathematical function was determined.Relationships between the µ and other statistical parameters were investigated by plotting the µ SOC stocks calculated at each spatial scale across Alaska versus its other statistical parameters calculated at the same scale.
To study the change in spatial structure of SOC stocks due to spatial scaling, we calculated the SOC stock variograms (Webster and Oliver, 2007) at different spatial scales: where Z(X i ) and Z(X i + h) are the measured SOC stocks at X i and X i + h, respectively; h is the lag; n is the number of paired comparisons at that lag; and γ (h) is the semivariance.By varying h in discrete steps, we obtain an ordered set of semivariances (Webster and Oliver, 1992).Suitable mathematical functions were fitted to the calculated semivariance values using weighted least squares (Webster and Oliver, 2007).The best-fit model parameters for each spatial scale are provided in Table 2.We examined the change in correlation length (range), total variance of SOC stocks (sill), and the unstructured variability (nugget) resulting from the change in the spatial scales of environmental variables.A variable is considered to have a strong spatial dependence if the nugget-to-sill ratio is less than 0.25; a moderate spatial dependence if the ratio is between 0.25 and 0.75; and otherwise a weak spatial dependence (Cambardella et al., 1994;Mora-Vallejo et al., 2008).The range, which is also referred as correlation length, is the distance at which the sill is reached and represents the maximum extent of the spatial correlation between predicted SOC stocks.The change in the spatial structure of variability of SOC stocks was evaluated using the nugget-to-sill ratio and the correlation length.

Inferred environmental controls on SOC stocks change with scale
The set of significant environmental predictors and their individual strengths varied with spatial scale.Statistically significant environmental predictors at different spatial scales and the predicted variance of SOC stocks are summarized in Table 1.From the 13 topographic attributes, 5 attributes -elevation (meters), slope (degree), aspect (degree from north), soil wetness index (SWI), and sediment transport index (STI) -were significant predictors of SOC stocks at various spatial scales.Among land cover types, barren, scrub, herbaceous, and cultivated land covers were significant predictors.Among different climatic factors, surface air temperature and potential evapotranspiration were significant predictors.Among surficial geology types, the glacial moraine, fluvial, undifferentiated mosaic, coastal delta, and volcanic mountain types were significant predictors at various spatial scales.Among all of the environmental variables investigated in this study, only elevation, temperature, potential evapotranspiration, and scrubland cover types were significant predictors of SOC stocks at all spatial scales.In general, at finer spatial scales, controls of topographic attributes were more   prominent, whereas land cover and surficial geology types were more significant predictors at coarser scales (≥ 500 m spatial scale).Fewer environmental variables controlled the variability of SOC stocks at finer spatial scales, whereas the number of predictors increased at larger spatial scales.
Our results showed that the strength of the control (median geographically weighted regression coefficient across Alaska) of elevation on SOC stocks decreased by 31 % as the spatial scale increased from 50 m to 1 km.Beyond this scale, we found no change in the control of elevation on predicted SOC stocks (Fig. 3a).The control of temperature on SOC stocks decreased with spatial scale between 50 m and 500 m and became constant at larger scales (Fig. 3b).The controls of elevation and temperature on SOC stocks across spatial scales can be accurately modeled by using exponential functions with R 2 = 0.83 and 0.94, respectively.The control of potential evapotranspiration decreased by 36 % as the spatial scale increased from 50 m to 500 m, and became constant beyond 500 m (Fig. 3c).The control of potential evapotranspiration on SOC stocks across spatial scales can be modeled by using an exponential decay function (R 2 = 0.97).The control of scrub vegetation (shrubs less than 5 m tall) on SOC stocks decreased linearly (from 0.3 to 0.13) with scale over the range 50 m to 10 km (Fig. 3d).At finer scales, scrub vegetation showed higher control on SOC stocks, which decreased by about 57 % as the scale increased.due to scaling, the relationship between scrub vegetation and SOC stocks remained positive.The control of scrub vegetation on SOC stocks can be modeled using a linear function (R 2 = 0.85).

Scaling impacts on the spatial structure of SOC stocks
The variogram summarizes the statistical structure of spatial variability and is known as a summarizing function of the investigated property.For the predicted SOC stocks, the variograms at finer scales showed significant amount of unexplained variance (nugget) (Table 2) due to measurement errors and/or microscale spatial heterogeneity not captured across the study area.Both the nugget and sill decreased as the scale increased from finer to coarser spatial scales.The nugget and sill decreased by 75 and 64 %, respectively, as s increased from 50 m to 500 m and remained constant beyond 500 m scale (Fig. 4a).The relationship between these variogram parameters and scale can be accurately described using exponential functions with R 2 of 0.98.The correlation length of SOC stocks remained relatively constant up to about 1100-1400 km across spatial scales (Fig. 4b).
The nugget-to-sill ratio of predicted SOC stocks at the 50 and 100 m spatial scales showed moderate spatial dependency (i.e., a nugget-to-sill ratio > 25 %; Cambardella et al., 1994).However, strong, and relatively similar, spatial dependency (a nugget-to-sill ratio ≤ 25 %) was predicted for spatial scales between s = 200 m and s = 10 km (18-21 %; Table 2).These results suggest that the predicted SOC stocks show different spatial structure below s = 100 m as compared to s ≥ 200 m.

Statistical behavior of spatial scaling of SOC stocks
We found different statistical properties in predicted SOC stocks as the spatial scale of environmental variables increased.The variance of predicted SOC stocks decreased by 45 % as the scale increased from 50 to 500 m and remained constant beyond 500 m (Fig. 5).An exponential function based on spatial scale accounted for 98 % of variability in the variance of predicted SOC stocks across Alaska.These results suggest that the spatial heterogeneity of SOC stocks decreases exponentially with scale until about 500 m.Beyond this scale, no change in the spatial heterogeneity of SOC stocks was predicted.The higher-order moments (skewness S and kurtosis K) of the predicted SOC stocks decreased as the spatial scale increased (Fig. 6a).S decreased from 1.4 to 0.8 as the scale increased from 50 to 200 m; and fluctuated between 1.25 and 0.74 as the scale increased further.K decreased from 5.7 to 3.6 as the scale increased from 50 to 200 m and fluctuated between 5.1 and 3.68 beyond a scale of 200 m.We evaluated the relationships of mean SOC stocks (µ) generated at each spatial scale with σ 2 , S, and K. Within the investigated range of µ SOC stocks (45.7 to 50.4 kg m −2 ), we found moderately accurate linear relationships between µ and σ 2 , S, and K of predicted SOC stocks, with R 2 values of 0.59, 0.63, and 0.55, respectively (Fig. 6b-d).Significant negative slopes of these relationships indicate that as µ increases, the higher-order moments of SOC stock decrease,  and therefore the µ of SOC stocks can be used to predict SOC spatial heterogeneity.

Discussions
Soil organic carbon lies at the interface between the atmosphere and pedosphere and can be altered at time scales rele-vant to climate change by anthropogenic and climatic factors.
As a result, the spatial heterogeneity of SOC stocks may affect the magnitude of greenhouse gas fluxes from the land surface and the terrestrial carbon cycle.Several studies have demonstrated the impact of DEM scale on calculated topographic attributes (Thompson et al., 2001;McMaster, 2002;Smith et al., 2006;Roecker and Thompson, 2010).However, very little attention has been paid to the impact of scaling on environmental controls and the predicted spatial heterogeneity of soil properties.We found that apparent significant environmental controllers of high-latitude SOC stocks changed with scale.Our findings are consistent with the results of Vasques et al. (2012), who showed inconsistent controls of various tropical system environmental factors in predicting SOC stocks as the scale of environmental variables increased from 30 to 1920 m.Due to a lack of other scaling studies examining the change in environmental controllers of SOC stocks, we compared our results with findings from two studies that reported environmental controllers of SOC stocks at two different scales.Johnson et al. (2011) reported that at the plot scale (∼ 1 m 2 ), surface air temperature, topographic attributes, and soil texture were primary controllers of SOC stocks across Alaska.Similarly, Martin et al. (2011) predicted SOC stocks at the 12 km scale over France, and reported land use and clay content to be important drivers of SOC stock spatial variability.Our results at the 50 m and 10 km scales are partially consistent with these findings, showing significant controls of temperature, land cover types, and topographic attributes on SOC stocks.Our results showed that the strength of controls of environmental factors that were significant decreased as spatial scale increased.The largest decrease was found in the control of scrubland cover type, and the smallest decrease was found in the control of temperature.The rates of these changes can be modeled using simple mathematical functions.
The scaling properties of soil moisture have been widely studied (Western and Boschl, 1999;Isham et al., 2005;Famiglietti et al., 2008;Li and Rodell, 2013).Some studies reported that the variance of soil moisture follows a power law relationship with surface area (Rodriguez-Iturbe et al   1995; Manfreda et al., 2007), while others have shown nonsimple scaling (Joshi and Mohanty, 2010;Riley and Shen, 2014;Pau et al., 2014).Beven (1995) pointed out that it is unlikely that any general scaling theory can be developed because of the dependence of hydrological systems on geological perturbations, and advocated for a disaggregation approach to develop scale-dependent models to represent hydrological heterogeneity.These studies suggested that further work focusing on the factors influencing soil moisture (topography, vegetation, soil properties, and rainfall) is expected to increase understanding of the mechanisms affecting the scaling properties of soil moisture.
Keeping these inferences and that soil moisture is an important control on SOC dynamics in mind, in this study we used the soil-forming factors that have been documented to affect SOC stocks to study the spatial scaling behavior of SOC stocks.However, in contrast to soil moisture, our results show that the variance of SOC stocks follows an exponential decay with spatial scale up to about 500 m, and then remains constant thereafter.We interpreted environmental controllers on SOC stocks in an equilibrium sense; i.e., over hundreds to thousands of years these controllers are expected to be related to observed SOC stocks.However, transient responses are not expected to follow the same scaling properties, given likely changes in thermokarst formation and other hydrological and geomorphological changes at high latitudes.
The spatial structure of SOC stocks results from the spatial structure of soil-forming factors (Jenny, 1941;McBratney et al., 2003).Besides the state factors described in the Introduction, SOC stocks of Alaska and other arctic regions are also impacted by cryopedogenic features (pingos, frost boils, hummocks) and thermal surface erosion resulting from the freezing and thawing of the polygonal land surface (Tarnocai and Bockheim, 2011;Ping et al., 2008;Ping, 2013).In this study, the spatial organization of predicted SOC stocks in-   creased with spatial scale as indicated by moderately lower nugget-to-sill ratios for spatial scales > 200 m.However, the correlation length remained relatively constant in the upscaled SOC stocks across Alaska.
One potential approach to represent the spatial heterogeneity of soil properties in ESMs is to relate their statistical properties to the mean state.As an analogy, several soil moisture studies investigated the relationships between soil moisture mean and higher-order statistics, such as σ 2 , S, and K (Famiglietti et al., 1999(Famiglietti et al., , 2008;;Ryu and Flamiglietti, 2005;Li and Rodell, 2013).The results from these studies suggest that the mean of soil moisture is often related to its σ 2 , S, and K, although with different functional forms depending on various system properties.Our results indicate a moderate but statistically significant linear relationship between the mean and higher-order moments of SOC stocks (i.e., σ 2 , S, and K decrease linearly with the mean SOC stock within the range of mean SOC stocks investigated in this study (45.5-50.4kg m −2 )).The strength of the linear relationship was highest between µ and S (R 2 = 0.63, P < 0.01), and lowest between µ and K (R 2 = 0.55, P < 0.05).However, available observations are insufficient to determine whether the same relationship holds with dynamically changing mean SOC stocks, highlighting an area for further investigation.
Current land models integrated in ESMs, such as CLM4.5 (Lawrence et al., 2012;Koven et al., 2013;Tang et al., 2013), use a nested subgrid hierarchy approach to represent the land surface heterogeneity.In this approach the grid cells (∼ 1 × 1 • ) are divided into nonspatially explicit land units (such as vegetated, lakes, urban, glacier, and crops), columns (with variability in hydrological, snow, and crop management), and plant functional types (accounting for variations in broad categories of plants and bare ground).We have shown here that this type of representation does not characterize the environmental controls and scaling properties of SOC stocks; rectifying this problem would require substantial restructuring of the model's subgrid hierarchy.One potential application of the relationships we developed in this study could be to apply them with coarse-resolution ESM results in order to generate fine-scale spatial heterogeneity parameters of SOC stocks that are more representative of the natural landscape.Currently the land models of most ESMs have a spatial scale of ≥ 50 km.In the next 5-10 years, we believe ESM land models will function much closer to the resolution we identify in our study as being representative of the SOC landscape heterogeneity (∼ 10 km; Bader et al., 2014).As the model resolution becomes finer in nextgeneration ESMs, data sets such as the one we describe in this study will be critical for model benchmarking.We note that many environmental factors that we found significant at various scales are not represented in current land models.However, representing these factors in future land model developments can improve our understanding of SOC dynamics of arctic/boreal systems.
Median values of regression coefficients of different environmental controls might change if a global or continental scale study would have been conducted.Other environmental factors which are not significant predictors of Alaskan SOC stocks might become significant in different study domains such as in temperate and tropical ecosystems.Therefore extending this kind of scaling study in larger spatial domains might produce more important information for benchmarking ESM results.We note this as a limitation to our study and recommend further studies to investigate these factors in other systems and at larger spatial extents.
Soil texture has also been reported to be related to SOC stocks.To investigate the use of soil texture in our scaling study, we collected the soil texture data currently used in CLM 4.5 (Bonan et al., 2002) that was available at ∼ 8 km spatial resolution across Alaska.However, this data has approximately one soil texture value in each Alaska ecoregion implying that the data must have been derived from a source with a much coarser spatial resolution (International Geosphere-Biosphere Programme soil data set that had 4931 soil mapping units globally; CLM 4.5 Technical notes).Because of this limitation, we were unable to include soil texture in the current scaling study.However, if soil texture information becomes available at finer spatial scales, it could be readily incorporated in future studies using the methods described here.

Summary and conclusions
Understanding the spatial scaling structure of SOC stocks and their relationships with environmental factors is important for prediction of carbon-climate feedbacks under a changing climate.Here, we estimated the spatial scaling properties of environmental factors and statistical properties of SOC stocks.We conclude that environmental controllers of SOC stocks change with scale, and the strength of environmental controls can be accurately modeled using simple mathematical functions.The variance of predicted SOC stocks decreases exponentially with scale up to about 500 m, and then remains constant thereafter.Our results showed that the mean predicted SOC content is linearly related with its variance, skewness, and kurtosis.If the objective of a study is to represent the spatial heterogeneity of SOC stocks resulting from environmental factors, the prediction spatial scale should be finer than ∼ 500 m.Similarly, the choice of environmental predictors should be based on the intended final spatial scale of predicted SOC stocks.Therefore, the current subgrid structure of many ESM-scale land models does not allow for characterization of the environmental controls and scaling properties affecting SOC stocks.The findings of this study may help to develop scale-appropriate techniques for ESM-scale land models and ultimately reduce the existing uncertainties in carbon-climate feedback predictions.

Figure 2 .
Figure 2. The spatial scaling approach used in this study.The black dot represents the location of an SOC sampling site, and the rectangles denote different pixel sizes of environmental variables (rectangles are not to scale). 4 Figure 3:

Figure 5 .
Figure 5. Variance of soil organic carbon stocks as a function of spatial scale.Each dot is predicted variance of SOC stock at each spatial scale across Alaska.The fitted function is SOC variance = 260.61+ 422.42 × exp(−0.0318× s), adjusted R 2 = 0.98, P < 0.001. 7 Figure 6:

Figure 6 .
Figure 6.Higher-order moments of predicted SOC stocks as a function of spatial scale (a), and relationships between predicted mean SOC stocks and its (b) variance, (c) skewness, and (d) kurtosis within the range of predicted mean values of SOC stocks.