Soil carbon stock increases in the organic layer of boreal middle-aged stands

Changes in the soil carbon stock can potentially have a large influence on global carbon balance between terrestrial ecosystems and atmosphere. Since carbon sequestration of forest soils is influenced by human activities, reporting of the soil carbon pool is a compulsory part of the national greenhouse gas (GHG) inventories. Various soil carbon models are applied in GHG inventories, however, the verification of model-based estimates is lacking. In general, the soil carbon models predict accumulation of soil carbon in the middle-aged stands, which is in good agreement with chronosequence studies and flux measurements of eddy sites, but they have not been widely tested with repeated measurements of permanent plots. The objective of this study was to evaluate soil carbon changes in the organic layer of boreal middle-aged forest stands. Soil carbon changes on remeasured sites were analyzed by using soil survey data that was based on composite samples as a first measurement and by taking into account spatial variation on the basis of the second measurement. By utilizing earlier soil surveys, a long sampling interval, which helps detection of slow changes, could be readily available. The range of measured change in the soil organic layer varied from−260 to 1260 g m−2 over the study period of 16–19 years and 23±2 g m−2 per year, on average. The increase was significant in 6 out of the 38 plots from which data were available. Although the soil carbon change was difficult to detect at the plot scale, the overall increase measured across the middle-aged stands agrees with predictions of the commonly applied soil models. Further verification of the soil models is needed with larger datasets that cover wider geographical area and represent all age classes, especially young stands with potentially large soil carbon source. Correspondence to: R. Mäkipää (raisa.makipaa@metla.fi)


Introduction
Changes in the soil carbon stock can potentially have a large influence on global carbon balance between terrestrial ecosystems and atmosphere.Globally soil contains three times more carbon than atmosphere and in boreal forests soil carbon stock is three times larger than that of vegetation (Schimel, 1995;Goodale et al., 2002).Soil carbon stock, and especially its topmost organic layer, in managed boreal forests is directly (by timber harvesting and soil scarification) and indirectly (e.g. by climate change) affected by human activities.Due to the potentially large human induced changes in soil carbon balance, reporting of the changes in the soil carbon stock is an essential part of the national greenhouse gas (GHG) inventories.Currently, the majority of the countries that are able to report soil carbon apply model based approaches and only a few countries can rely on repeated soil measurements.In general, soil carbon models that can also be used in the GHG reporting (Peltoniemi et al., 2007) predict loss of carbon in regenerated young stands and accumulation of soil carbon in the middle-aged stands (Mäkipää et al., 1999;Peltoniemi et al., 2004;Palosuo et al., 2008).Such a modeled pattern is in good agreement with the chronosequence studies (e.g.Covington, 1981;Federer, 1984;Peltoniemi et al., 2004) but not confirmed with repeated measuments of permanent study sites (Yanai et al., 2000).
Verification of the modeled soil carbon dynamics with empirical data is essential for the development of reliable inventory methods.Measuring changes in soil carbon stocks is, however, challenging due to the fact soils are heterogeneous (Järvinen et al., 1993;Liski, 1995) and the rate of change is relatively small compared to the size of the stock (Yanai et al., 2003a;Peltoniemi et al., 2004).Due to the slow changes, long sampling interval may be necessary in order to Published by Copernicus Publications on behalf of the European Geosciences Union.
allow measurable changes to take place before re-sampling.At the moment, the use of earlier soil data may be the only effective way to test the soil carbon change hypothesis.Analyzing the change on the basis of repeated measurements of the earlier established sample plots is often challenging due to dissimilarities between the sampling designs applied in the first and second inventory.In general, previous soil surveys contain information on mean carbon stocks but they lack information on within-site spatial variation.Current sampling can be designed to provide representative spatial information, but following the sampling design of the first one would improve the power of statistical testing due to correlated covariances of the first and second sampling.
In addition to soil heterogeneity, measuring soil changes is also challenging due to spatial autocorrelation of soil properties.Spatial autocorrelation should be accounted for when estimating the significance of a change in single study plots.Positive autocorrelation enlarges the variance of the mean of single plots but, on the other hand, it may also decrease the estimation variances of kriging estimations.Significant small-scale autocorrelation in soil properties has been detected when the properties of mineral soil sites have been investigated using variogram analysis and kriging (Järvinen et al., 1993;Arrouays et al., 1997;Bruckner et al., 1999).Small-scale spatial autocorrelation in the amount of carbon has also been identified using variogram analysis (Liski, 1995;Möttönen et al., 1999;Schöning et al., 2006;Muukkonen et al., 2009).Spatially autocorrelated soil data have also been studied using cross-variograms and co-kriging in many geostatistical papers (Papritz and Flühler, 1994;Papritz and Webster, 1995;Lark, 2002).However, these methods are not applicable if the samples from the first inventory are composite ones with unknown variances.Since spatial variation in soil properties is widely acknowledged, the influence of spatial pattern on estimated mean values has commonly been reduced by taking numerous subsamples, which are analyzed as one composite sample (e.g.Ellert et al., 2000;Smith, 2000).
The objective of this study was to evaluate soil carbon changes in the organic layer of boreal middle-aged forest stands.In addition, the aim was to develop methods that facilitate effective use of earlier soil inventory data that was based on composite samples together with new data that carry information on the spatial variation of soil properties.The focus was on the soil organic layer, a clearly distinguished soil horizon of podzols, because it is the most dynamic part of forest soil and the most likely changes first take place in this topmost layer.This study was restricted to middle-aged stands for which stand development phase soil models predict consistently a trend of increasing soil carbon stock.However, scrutinized tests of such trends are lacking.

Soil sampling and carbon analysis
This study was based on soil data collected from a subset of 38 sample plots from a nation-wide network of forest monitoring plots.A systematic network of 3000 permanent sample plots was established by the Finnish National Forest Inventory for the monitoring of forest ecosystems (Mäkipää and Heikkinen, 2003), and the first soil sampling on mineral soil sites was performed on a sub-sample of 486 plots from the nation-wide network during 1986-1989 (Tamminen and Starr, 1990).In 2005, soil sampling was repeated on a subsample of n = 38 plots where the stand age varied between 22 and 65 years at the time of the first sampling.The 38 plots were located in Southern Finland, excluding the coastal region (Fig. 1).The tree stand was dominated by Scots pine (Pinus sylvestris L.) on 24 of the plots and by Norway spruce (Picea abies (L.) H. Karst.) on 14 plots (Table 1).The sample plots divided into two or more forest patches (according to fertility level, stand age, or management history) were not included.The tree and stand parameters were measured on circular plots of 300 m 2 (radius 9.77 m).The fertility class of the plots ranged from herb-rich to xeric heath forests (Table 1).The soil type was podzol and the organic layer was mor or moder (plots with signs of peat formation were excluded).The sample plots used in this study were a random sample that represented intermediate age classes of coniferous forest stands on mineral soils in the southern Finland.
The organic layer (excluding the litter layer) sampling in the first soil survey was based on composite samples, comprising m 1 = 30 sub-samples (Tamminen and Starr, 1990).The sub-samples were combined which resulted in only one mean value per plot and no information on variance.The second sampling was designed to provide information on the variation of the soil carbon stock and all sub-samples were analyzed separately.The organic layer was sampled with a cylinder (d = 58 mm).Above-ground parts of the living plants as well as the litter layer were excluded from the samples.Mineral soil horizon was separated out according to visual difference between the structure of organic and mineral soil layers of the podzolic soil.The instructions and the soil sampling equipment used in the first and in the second sampling were kept as similar as possible.Five different field teams participated in the first sampling (Table 1), and one person was responsible for taking the samples in the second sampling.The means of the amount of carbon in the samples taken by the different groups on the first sampling occasion were tested with simple mean tests.The results of these tests confirmed that sampling by different groups did not differ significantly from each other.
The sampling design is presented in Fig. 2. The organic layer samples were taken at points lying on a circle, r = 11 m, centered on the inventory plot.Although the first and the second round of sampling and carbon analyses were as similar as possible, sampling could not be performed at exactly the same points due to destruction of some of the first sampling points by removal of soil samples.As a result, some of the second sampling points had to be shifted.Furthermore, the second sampling of the organic layer was designed to provide information also on the spatial within-site variation of the soil carbon stock.In the second sampling, m 2 = 40 subsamples were taken instead of 30 in order to also include shorter distances between the sampling points.The organic layer consist both partially decomposed matter whose origin can be spotted on sight and well-decomposed organic matter, the origin of which is not readily visible.All the 40 subsamples were analyzed separately.
The samples were dried at a temperature of 35-44 • C, weighed, milled and sieved to pass through a 2 mm bottom sieve.The moisture content of the air-dried samples was determined on a TGA analyzer.The total carbon concentration was analyzed using a Leco CHN analyzer (Leco, St Joseph, Fig. 2. Sampling design.In the first sampling, three samples were collected at each of the 10 locations points in the square).In the second sampling, samples were collected at the same locations unless t sampling was destructive, in which case the second sampling points were shifted counterclockwise degrees.In the second sampling 4 samples were collected at each of the 10 locations.The fourth point is located consecutively at one of the three grey points shown in the square, either side by side black one or 0.2 m or 0.4 m away from it.16 Fig.2. Sampling design.In the first sampling, three samples were collected at each of the 10 locations (black points in the square).In the second sampling, samples were collected at the same locations unless the first sampling was destructive, in which case the second sampling points were shifted counterclockwise by 18 degrees.In the second sampling 4 samples were collected at each of the 10 locations.The fourth sample point is located consecutively at one of the three grey points shown in the square, either side by side with a black one or 0.2 m or 0.4 m away from it.MI, USA) in the Central Laboratory of the Finnish Forest Research Institute, which is an accredited test laboratory (in accordance with the standard SFS-EN ISO/IEC17025).

Variogram analysis
Standard geostatistical methods were used in analyzing the spatial autocorrelation and amounts of carbon on the plots (see Webster and Oliver, 2001).They were performed in the R environment (R Development Core Team, 2006) using the libraries geoR (Ribeiro Jr. and Diggle, 2001) and gstat (Pebesma, 2004).Spatial autocorrelation was studied using variograms.x tki shall denote the location of the i-th sub-sample of plot k on the t-th sampling occasion, t = 1,2, and Z t (x) shall denote the carbon concentration at the location x at the time of the t'th sampling.Plot-specific empirical variograms were estimated from the spatially explicit observations  a Fertility level of the site according to the Finnish site type classification (Cajander, 1949;Hotanen et al., 2008) second sampling using the equation where N (h) is a set of pairs (i,j ) ∈ 1,2,...,m 2 for which is the number of the pairs in the set.
A spherical model was fitted to the empirical variograms.The spherical model is defined as where c 0 is the nugget variance parameter, c the sill variance parameter, and a the range of spatial correlation.The spherical model was used because it has a well-defined range a and it exhibits linear behavior near the origin, thus making it suitable for representing properties that have high short-range variability.The experimental variograms were fitted by the restricted maximum likelihood criterion (REML) by weights The spatial autocorrelation within a plot of 300 m 2 is strong if the nugget parameter is much smaller than the sill.On the other hand, if the nugget parameter is approximately same as or bigger than the sill, there is not much actual spatial autocorrelation as the nugget explains most of the spatial variability.

Plot-specific variances for the first sampling
Because the within-plot variances of the first measurements were unknown, they were estimated on the basis of the variograms fitted to the second measurements under the assumption that the variation at the time of the first measurement was similar to that of the second one.The observations of the first sampling occasion are plot-level averages and in the plot level analysis they are considered as predictions of the unknown means 1 over circles U k with radii r and origins at the centre points of the sample plots.The variances of the prediction errors can be expressed as where The covariances in Eqs.(6-8) were obtained from the plotspecific variogram models fitted to the second measurement, and the integrals were approximated by appropriately scaled sums over dense grids discretizing the circles U k .

Change in the carbon stocks of the organic layer
For single plots the mean carbon stocks of the second measurement of the organic layer were calculated by ordinary block kriging.The ordinary block kriging estimate over a block U k is a weighted average of the data, Model-unbiasedness of the estimator 9 is ensured by the restriction N i=1 λ i = 1.The estimation variance is where and This variance was computed by replacing γ k in (10-12) with the variogram model fitted to the second measurement of the k-th plot.
The confidence intervals of the first and the second measurement were calculated from the equations z1k ± s 0.95 Var(z 1k ) and ẑ2k ± s 0.95 Var(ẑ 2k ), where s 0.95 is the statistic from Student's t distribution at the 95 % confidence level.The confidence intervals of the first and the second means of each plot were compared with each other.If they did not intersect, the change was considered statistically significant.

Variogram analysis
In most of the studied plots, spatial autocorrelation in soil carbon stock of organic layer was found on the basis of the variograms (Fig. 3), which indicates that soil spatial variation needs to be taken into account when plot-wise mean carbon stocks are estimated.The estimated range parameter was constant in two plots and in one plot (number 496702), where the estimated range was 0.70 m, the spatial variation seemed to be very small-scale -this plot could also be considered to have a constant variogram.In these cases, spatial pattern does not influence on mean estimates of plot-wise carbon stocks of organic layer.In 17 cases out of 35 where spatial autocorrelation in soil carbon of organic layer was found, it seemed to disappear at distances shorter than 7 m (range parameter < 7 m).On the other hand, the estimated range parameter was larger than the radius of the soil sampling plots  (11 m) in 7 of the 38 plots, and in 2 plots the range was larger than the diameter of the plot indicating that scale of spatial variation may be larger than we were able to evaluate with this data.Due to observed within-site spatial variation, mean carbon stocks of the organic layer in this study were estimated by kriging, which yields more realistic estimates than calculation of simple averages.

Soil carbon stock and stock change in the organic layer
We calculated simple mean stock estimates of the first and second measurements and the difference between them based on the empirical variances of the plot-specific stock estimates (Table 1).The mean carbon stock of the organic layer in the middle-aged stands (40-84 yr) of boreal forests during the second measurement was 1852 ± 67 g m −2 .The mean carbon stock of the organic layer of all the 38 plots in the first measurement, 16-19 yr earlier, was 1444 ± 95 g m −2 .The mean change of all the 38 plots was 412 ± 44 g m −2 .The amount of carbon was increased on 35 of the 38 plots (Table 1).However, the increase was statistically significant on only 6 plots (Fig. 4 and Table 1).The measured change was larger in younger stands (Fig. 5).

Discussion
The average rate of increase in the carbon stock of the organic layer (23 ± 2 g C m −2 yr −1 ) measured with repeated sampling in these managed boreal forests of intermediate age classes was higher than that reported earlier on the basis of chronosequency studies.Soil carbon accumulation of 8 g C m −2 yr −1 in the organic layer was reported in the chronosequence of windthrow pits in Alaska (Bormann et al., 1995), and long-term accumulation of the organic layer without fire resulted in an increase of 5 g C m −2 in Sweden (Wardle et al., 2003).Peltoniemi et al. (2004) measured and simulated 64 sites in boreal coniferous stands in Finland and obtained a 4.7 ± 1.4 g C m −2 yr −1 increase in carbon in the simulations, and a 4.2 ± 1.2 g C m −2 yr −1 increase in the   organic layer in measured chronosequency data.Their study sites represented a wide range of age classes, while this study was restricted to stands of intermediate age classes where a significant increase in the soil carbon stock was expected on the basis of the earlier simulation studies.Increase in the soil carbon stock in over 20 yr-old stands has been consistently predicted with several models (e.g.Mäkipää et al., 1999;Chertov et al., 2001;Yanai et al., 2003b;Peltoniemi et al., 2004).In addition, we observed that the rate of soil carbon change tends to decrease with stand age (Fig. 5), which agrees with the understanding of soil carbon dynamics based on modeling where the rate of soil change is driven by annual biomass and litter production that decrease after a fast growth period and canopy closure.
The measured soil carbon sequestration in the middleaged stands is indirectly supported by the CO 2 flux data of eddy covariance measurement sites (e.g.Valentini et al., 2000;Kolari et al., 2004).After early development of a stand and closure of the canopy they show only minor variation in the gross primary production (GPP), however, the total ecosystem respiration tends to decrease with stand age, which is an indication of soil carbon accumulation (Kolari et al., 2004).The measured rate of soil carbon sequestration in the organic layer (23 ± 2 g C m −2 yr −1 ) is relatively slow in comparison to stand scale carbon sequestration of 192 and 323 g C m −2 yr −1 measured on 40-and 75-yr-old stands, respectively (Kolari et al., 2004).Thus, in the middle-aged stands soil organic layer may contribute to less than 10 % of the forest carbon sink, but in the old-growth stands, where carbon sink is measured to be 240 g C m −2 yr −1 , soil carbon sequestration continues and living trees have only a minor role (some 40 g C m −2 yr −1 ) (Luyssaert et al., 2008).According to the earlier model predictions supported by our current results, one may interpret that a change in the age class distribution from a predominance of younger age classes to middle-aged and mature forests may result in an increase in the soil carbon stock of organic layer.The mean carbon stock of the organic layer measured in the stands of 40-80 yr (1850 ± 70 g m −2 ) is consistent with the large data set of 1248 sample plots in Southern Finland (Tamminen, 1991), where the mean carbon stocks of stand development classes 2, 3, and 4 (from young to mature) were 1450, 1630, and 1740 g C m −2 , respectively (assuming the inverse of the van Bemmel factor, 1.72, for converting the organic matter content to the carbon concentration).In this study, the average carbon stock of organic layer carbon increased by almost 30 % from 1444 to 1852 g m −2 in less than 20 yr.This accumulation of the soil organic layer takes place after a remarkable decline in the soil carbon stock after a regeneration of the stands and a release of carbon during early years of succession; Kolari et al. (2004) measured a net carbon source of 400 g C m −2 yr −1 on clear-cut site where the size of the organic layer carbon stock was similar to this study.
At the plot scale, the block kriging estimates and variances of the carbon amounts showed a significant change only on 6 of the 35 sites where the amount of carbon had increased.In general, changes smaller than one third of the stock were not significant and many large changes were also not significant due to large within-site variation (Fig. 4).This result is consistent with earlier findings which indicate that the detection of a change in soil carbon by re-measuring a single plot is challenging or impossible due to soil heterogeneity (large spatial variation) and the small temporal changes relative to the large total amount of carbon in forest soil (Yanai et al., 2000;Conen et al., 2003;Smith, 2004;Mäkipää et al., 2008).The amount of carbon appeared to have decreased on 3 plots (Table 1, Fig. 5), two of which were relatively fertile Norway spruce stands thinned some 5 years before the first sampling, and one had a management history with a lower basal area of the trees (harvesting of seedling trees), which evidently resulted in a decreased input of litter from the trees.
Possible measurement errors in this study are related to mislocated sample plots or sampling points within a plot in the second measurement, differences caused by the sampling practices of measurement groups, or inaccuracy in taking soil samples or in the carbon analysis.The exact location of the sample plots was known at both measurement times due to permanent marks of the sample plot centers, detailed descriptions of the location of the sampling points, very exact instructions and equipment.Similarly, the measurement error between the sample coordinates of the sampling points and the exact locations where they were taken was small, within 10 cm.The field personnel should not have significantly affected the results as they used the same precise instructions for both measurements.Furthermore, the different field groups were tested to ensure that they were consistent with each other in terms of their sampling practices.The carbon concentrations and moisture of the soil samples were measured in an accredited laboratory and are considered very reliable.

Conclusions
We found spatial autocorrelation in the carbon stock of the organic layer, that has already been demonstrated in many other studies for a range of soil properties (Järvinen et al., 1993;Hokkanen et al., 1995;Liski, 1995;Bruckner et al., 1999;Möttönen et al., 1999;Muukkonen et al., 2009).Therefore, the spatial structure of the data was considered in the analysis of the soil data.The spatial variation of the soil carbon stock was taken into account on the basis of spatial sampling performed during the second measurement (mean and variance of carbon stocks of second sampling were estimated by kriging) with the assumption that the within-plot variation in soil carbon was the same in the first sampling.This approach provided a realistic basis for the plot-scale analysis of the soil carbon changes that have to build on soil inventories where composite samples are used in the first measurements.With spatial sampling in the second measurement, we were able to assess the total variation and to gain a more reliable assessment of the significance of the change than with only composite samples.Since the carbon stock changes are small in comparison to the size of the soil stock, the time period needed to detect a change could be tens of years.Therefore, in the analysis of soil carbon changes, we must be able to use the data from earlier soil surveys.Regional soil surveys provide a considerable amount of data based on composite sampling (e.g.Arrouays et al., 2001;Coomes et al., 2002;Tremblay et al., 2002;Jones et al., 2005;Lettens et al., 2004), and this could be utilized in assessing changes after re-measuring the same plots.In addition, it is common in forest soil surveys that soil samples cannot be collected at every planned point because of natural or other obstacles.The methods applied in this study can also take into account sampling with an unequal number of soil samples per plot.
We measured carbon sink in the soil organic layer, which may represent some 10 % of the overall carbon sink in the middle-aged boreal forest stands.Soil carbon sequestration in the boreal forests of this age class has earlier been reported on the basis of the various soil carbon models (e.g.Mäkipää et al., 1999;Peltoniemi et al., 2004;Palosuo et al., 2008) that are also applied in the large scale forest carbon inventories including national GHG reporting under the UNFCCC.Our analysis of the empirical data from the southern boreal zone indicates that the model predicted soil carbon trend in the middle-aged stands complies with measurements.However, further verification of the models is needed with larger datasets that represent all age classes and a wider geographical area.

Fig. 1 .Fig. 1 .
Fig. 1.Sample plots were located in the southern boreal and in the central boreal vegetation z : 2 = herb-rich heath forest, 3 = mesic heath forest, 4 = sub-xeric heath forest, and 5 = xeric heath forest.b The measurement year of the first sampling.The second sampling was conducted in 2005 of all plots.c At the time of the second sampling in 2005.d The basal areas of Scots pine, Norway spruce and deciduous trees on a plot measured during the first sampling.e The field teams that collected the first samples.f The symbol * indicates a significant change (95 % confidence intervals did not intersect).

Fig. 3 .
Fig. 3. Plot-specific empirical variograms of carbon concentration (dots), numbers of pairs of observations contributing to each estimated value (attached to the dots), spherical models fitted to the empirical variograms (solid lines), and the parameter values of the models (a = range, c 0 = nugget, c = sill).

Fig. 3 .
Fig. 3.The amount of carbon at the time of the first and the second sampling with 95 % confidence intervals, connected with a line which describes the magnitude of the change.

Fig. 4 .
Fig. 4. The rate of soil carbon change in the organic layer in relation to stand age.

Fig. 4 .
Fig. 4. The amount of carbon at the time of the first and the second sampling with 95 % confidence intervals, connected with a line describes the magnitude of the change.

Fig. 3 .
Fig. 3.The amount of carbon at the time of the first and the second sampling with 95 % confidence inter connected with a line which describes the magnitude of the change.

Fig. 4 .
Fig. 4. The rate of soil carbon change in the organic layer in relation to stand age. 17

Fig. 5 .
Fig. 5.The rate of soil carbon change in the organic layer in relation to stand age.

Table 1 .
Mean carbon stocks of the organic layer according to the first (z 1k ) and the second measurement (z 2k ) and their standard deviations ( Var(z 1k ) and Var(z 2k )).