Biogeosciences, 13, 4439–4459, 2016 www.biogeosciences.net/13/4439/2016/ doi:10.5194/bg-13-4439-2016 © Author(s) 2016. CC Attribution 3.0 License. Underestimation of boreal soil carbon stocks by mathematical soil carbon models linked to soil nutrient status Boris Tˇupek1, Carina A. Ortiz2, Shoji Hashimoto3, Johan Stendahl2, Jonas Dahlgren4, Erik Karltun2, and Aleksi Lehtonen1 1Natural Resources Institute Finland, P.O. Box 18, 01301 Vantaa, Finland 2Swedish University of Agricultural Sciences, P.O. Box 7014, 75007 Uppsala, Sweden 3Forestry and Forest Products Research Institute, Tsukuba, Ibaraki 305-8687, Japan 4Swedish University of Agricultural Sciences, Skogsmarksgränd, 90183 Umeå, Sweden Correspondence to: Boris Tˇupek (boris.tupek@luke.fi), Aleksi Lehtonen (aleksi.lehtonen@luke.fi) Received: 21 December 2015 – Published in Biogeosciences Discuss.: 18 January 2016 Revised: 6 July 2016 – Accepted: 7 July 2016 – Published: 10 August 2016 Abstract. Inaccurate estimate of the largest terrestrial carbon pool, soil organic carbon (SOC) stock, is the major source of uncertainty in simulating feedback of climate warming on ecosystem–atmosphere carbon dioxide exchange by process- based ecosystem and soil carbon models. Although the mod- els need to simplify complex environmental processes of soil carbon sequestration, in a large mosaic of environments a missing key driver could lead to a modeling bias in predic- tions of SOC stock change. We aimed to evaluate SOC stock estimates of process- based models (Yasso07, Q, and CENTURY soil sub-model v4) against a massive Swedish forest soil inventory data set (3230 samples) organized by a recursive partitioning method into distinct soil groups with underlying SOC stock develop- ment linked to physicochemical conditions. For two-thirds of measurements all models predicted ac- curate SOC stock levels regardless of the detail of input data, e.g., whether they ignored or included soil properties. How- ever, in fertile sites with high N deposition, high cation ex- change capacity, or moderately increased soil water content, Yasso07 and Q models underestimated SOC stocks. In com- parison to Yasso07 and Q, accounting for the site-specific soil characteristics (e. g. clay content and topsoil mineral N) by CENTURY improved SOC stock estimates for sites with high clay content, but not for sites with high N deposition. Our analysis suggested that the soils with poorly predicted SOC stocks, as characterized by the high nutrient status and well-sorted parent material, indeed have had other predomi- nant drivers of SOC stabilization lacking in the models, pre- sumably the mycorrhizal organic uptake and organo-mineral stabilization processes. Our results imply that the role of soil nutrient status as regulator of organic matter mineralization has to be re-evaluated, since correct SOC stocks are decisive for predicting future SOC change and soil CO2 efflux. 1 Introduction In spite of the historical net carbon sink of boreal soils, 500 Pg of carbon since the last ice age (Rapalee et al., 1998; DeLuca and Boisvenue, 2012; Scharlemann et al., 2014), boreal soils could become a net source of carbon dioxide to the atmosphere as a result of long-term climate warm- ing (Kirschbaum, 2000; Amundson, 2001). They have the potential to release larger quantities of carbon than all an- thropogenic carbon emissions combined (337 Pg; Boden et al., 2010). In order to preserve the soil carbon pool and to utilize the soil carbon sequestration potential to mitigate anthropogenic CO2 emissions, mitigation strategies of cli- mate forcing aim to improve soil organic matter management (Schlesinger, 1999; Smith, 2005; Wiesmeier et al., 2014). Supporting soil management decisions requires an ac- curate quantification of spatially variable soil organic car- bon (SOC) stock and SOC stock changes (Scharlemann et al., 2014). The initial level of SOC stock is essential in order to estimate SOC stock changes (Palosuo et al., 2012; Todd-Brown et al., 2014), especially when estimating carbon emissions due to land-use change, e.g., afforestation of grass- Published by Copernicus Publications on behalf of the European Geosciences Union. 4440 B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility lands (Berthrong et al., 2009). Process-oriented soil carbon models like CENTURY, Roth-C, Biome-BCG, ORCHIDEE, JSBACH, ROMUL, Yasso07, and Q are important tools for predicting SOC stock change, but there are also risks for poor predictions (Todd-Brown et al., 2013; DeLuca and Boisv- enue, 2012). The models need further validation and im- provement as they show poor spatial agreement on fine scale and moderate agreement on regional scale against SOC stock data (Todd-Brown et al., 2013; Ortiz et al., 2013). Despite the potentially quantitative importance of CO2 emissions the ex- pected change will be small in relation to the SOC stock. Therefore, the uncertainty of measurements and/or model estimates could prevent conclusions on SOC stock changes (Palosuo et al., 2012; Ortiz et al., 2013; Lehtonen and Heikki- nen, 2015) especially for the soils with the largest SOC stocks, which are the most sensitive to carbon loss. Beside large uncertainties, the poor agreement between the modeled and measured SOC stocks (Todd-Brown et al., 2013) could also indicate missing biotic or abiotic drivers of long-term carbon storage (Schmidt et al., 2011; Averill et al., 2014). For example, ignoring the essential role of soil nutrient availability in ecosystem carbon use efficiency (Fernández- Martínez et al., 2014) could lead to missing important con- trols of plant litter production and soil organic matter stabi- lization mechanisms. Soil nutrient status is linked to the mo- bility of nutrients in the water solution (Husson et al., 2013), production, quality and microbial decomposition of plant lit- ter (Orwin et al., 2011), and formation of the soil organic matter (SOM). The SOM affects soil nutrient status by re- cycling of macronutrients (Husson et al., 2013), and water retention and water availability (Rawls et al., 2003). In spite of state of the art soil carbon modeling based on the amount and quality of plant litter “recalcitrance”, affected by climate and/or soil properties as in the Yasso07, Q, and CENTURY models, these types of process-based models do not include mechanisms for SOM stabilization by (a) the or- ganic nutrient uptake by mycorrhizal fungi; (b) humic or- ganic carbon interactions with silt-clay minerals; and (c) the inaccessibility of deep soil carbon and carbon in soil aggre- gates to soil biota (Orwin et al., 211; Sollins et al., 1996; Torn et al., 1997; Six et al., 2002; Fan et al., 2008; Dun- gait et al., 2012; Clemente et al., 2011). Although the models do not contain aforementioned mechanisms and controls for changes in SOM stabilization processes, they have been pa- rameterized using a wide variety of data sets and can treat soil biotic, physicochemical, and environmental changes im- plicitly. The Yasso07 model (Tuomi et al., 2009, 2011) is an advanced forest soil carbon model and it is used for Kyoto protocol reporting of changes in soil carbon amounts for the United Nations Framework Convention on Climate Change (UNFCCC) by European countries, e.g., Austria, Finland, Norway, and Switzerland. The Q model (Ågren et al., 2007) is a mechanistic litter decomposition model developed in Sweden and used, e.g., to compare results produced with Swedish national inventory data (Stendahl et al., 2010; Ortiz et al., 2011) and also with other models at national or global scales (Ortiz et al., 2013; Yurova et al., 2010). The CEN- TURY model (Parton et al., 1987, 1994; Adair et al., 2008) is one of the most widely applied models and it is used for soil carbon reporting to the UNFCCC by Canada, Japan, and USA. Although individual parameters and functions vary, mathematical models such as Yasso07, Q, and CENTURY have similar structures. For example, these models are driven by the decomposition rates of litter input and SOM. Decom- posing litter and SOM is divided into pools based on litter quality, and its transfer from one pool to another is, apart from model functions and parameters, affected by temper- ature (Q), and/or water (Yasso07), and/or soil texture and structure (CENTURY). The Q model does not include ex- plicit moisture functions, whereas precipitation affects de- composition for the Yasso07 and CENTURY models (Tuomi et al., 2009; Adair et al., 2008). On the other hand, the mod- els do not explicitly or by default include mechanisms that reduce decomposition by excessive precipitation/moisture (Falloon et al., 2011). We hypothesized that (1) soil carbon estimates of the Yasso07, Q, and CENTURY models would deviate for soils where SOC stabilization processes not implicitly accounted by the models are predominant, (2) the Yasso07 and Q mod- els ignoring soil properties would fail on the nutrient-rich sites of the southwestern coast of Sweden and on occasion- ally paludified clay and silt soils, and (3) the CENTURY model outperforms the Yasso07 and Q models due to fact that it includes soil properties as input variables. We grouped Swedish forest soil inventory data into ho- mogenous groups with specific soil physicochemical con- ditions using a regression tree and recursive partitioning modeling methods. After that we ran the models until they reached an equilibrium with a litter input that was derived from the Swedish forest inventory. Thereafter, we compared the model estimates against data by groups that were ob- tained from the regression tree model. In discussion we ad- dress the reasons why the models deviate and indicate direc- tions of further improvements. 2 Material and methods 2.1 Measurements We analyzed data from the Swedish forest soil inventory (SFSI), which is a stratified national grid survey of vege- tation and physicochemical properties of soils (SLU, 2011; Olsson et al., 2009). The soil data distinguished between the organic, B (0–5 cm of B horizon), BC (45–55 cm below ground surface), and C (55–65 cm from the top of the mineral soil) horizons (Olsson et al., 2009). All analysis was done using R software for statistical computing and graphics (R Core Team, 2014). The soil data were identical to a data set used in Stendahl et al. (2010). We restricted our sample plots Biogeosciences, 13, 4439–4459, 2016 www.biogeosciences.net/13/4439/2016/ B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility 4441 Table 1. Description of the Swedish Forest Soil Inventory (SFSI) data reduction of soil sorting of parent material and humus types; SFSI conversion estimate of soil classes of soil moisture to numerical representation of soil water content; and SFSI conversion estimate of classes to numerical representation of soil texture (sand, silt, and clay content for sediments by Lindén (2002) and for tills by Albert Atterberg’s distribution of the different grain size fractions). Sorting parent material Humus type Moisture SFSI Reduced SFSI Reduced SFSI SFSI Numeric Bedrock Bedrock Moder No-peat Water Long-term Poorly sorted sediments Unsorted Mor 1 No-peat level (m) moisture % Tills Unsorted Mor 2 No-peat Dry < 2 10 Well-sorted sediments Sorted Mull No-peat Fresh 1–2 20 Mull-Moder Peat Fresh-moist < 1 30 Peat Peat Moist < 0.5 50 Peat-Mor Peat Texture SFSI Numeric Sediments Tills Sand % Silt % Clay % Sand % Silt % Clay % Bedrock 0 0 0 0 0 0 Boulder 0 0 0 0 0 0 Gravel 10 0 0 10 0 0 Coarse sand 40 5 0 40 5 0 Sand 80 10 0 45 10 0 Fine sand 70 25 5 55 15 0 Coarse silt 50 40 10 65 20 5 Fine silt 10 75 15 55 35 10 Clay 0 65 35 0 85 15 Peat 0 0 0 0 0 0 to minerogenic soils since the Q, Yasso07, and CENTURY models were not developed for use on peat soils, and only to plots for forest land use with Swedish forest inventory data (SFI). We also excluded samples with total SOC stock be- low 2.8 and above 470.5 (tCha−1), i.e., samples with SOC stock below 0.01 and above 99.9 percentile. Measurement data originated from 1993 to 2002, which constitute a full inventory, and from 2020 sample plots located around Swe- den, and in total it included 3230 samples. For each sample plot the weather (years 1961–2011) and N deposition (years 1999–2001) data were retrieved from the nearest stations of Swedish Meteorological and Hydrological Institute (SMHI) network (Fig. 1). The plots, which were linked by the closest distance to the given weather station had the same weather and N deposition data, and the number of soil samples per station ranged between 10 and 70. The mean total SOC stock of samples corresponding to weather stations ranged from 40 to 200 (tCha−1), and the SOC stock level decreased from southern to northern Sweden (Fig. 1). Each sample plot contained categorical data from the field survey on the sorting of soil parent material, humus type, soil texture, and soil moisture. In our analysis we reduced categorical classes by basing them on the sorting of soil par- ent material and humus type (Table 1). We determined nu- meric values for silt, clay, and sand content from soil texture categories by Albert Atterberg’s distribution of the different grain size fractions in tills and distributions for sediments by Lindén (2002) (Table 1). We also determined numeric val- ues of volumetric soil water content (SWC) from categorical field data classified according to the depth of the ground wa- ter level (WL; Table 1). As is typical for soil carbon inventories, the variation of data was large (Table 2). For example, the mean total SOC stock of all samples was 93 (tCha−1) while 1st and 99th per- centiles were 17 and 309 (Table 2). The mean SOC stock was 33.3 and 66.8 (tCha−1) for the humus horizon and the mineral soil. The mean values of cation exchange capacity (CEC) (23.9 mmolc kg−1), the base saturation (36.4 %), and the C /N ratio (16.5) indicated conditions of medium fertil- ity, although the soils were mostly acidic (mean pH was 5.2). The mean prevailing soil water content (22.3) was typical for the well-drained forest soils. The mean annual temperatures ranged from below 0 to above 8 ◦C, and annual precipita- tion varied between 392 and 1154 mm (Table 2). Total SOC stock for all the samples generally increased for peat and peat like humus forms, for well-sorted sediments, for soils with high fraction of silt and clay and with increasing soil mois- ture (Fig. S1 in the Supplement). www.biogeosciences.net/13/4439/2016/ Biogeosciences, 13, 4439–4459, 2016 4442 B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility Table 2. Descriptive characteristics (mean, confidence interval, 1st, 50th, and 99th percentile) of selected variables (n= 3230 samples). The values of the bulk density, cation exchange capacity, base saturation, C /N ratio, and pH are shown only for BC soil horizon (fixed 45–50 cm depth below the ground surface) due to the strong correlation to the total soil carbon stock. The soil was cut off at 1 m. The site productivity index (H100, m) is an approximation of the site fertility expressed as the height of trees at 100 years of age. Stand and understory biomass, and litter input are modeled values for approximated equilibrium conditions based on observations. Mean CI 1st percentile 50th percentile 99th percentile Total soil carbon stock (tCha−1) 93.24 1.95 17.02 79.68 308.68 Humus carbon stock (tCha−1) 33.29 1.17 3.89 22.82 176.66 Mineral soil carbon stock (tCha−1) 66.82 1.7 6.92 54.81 273.91 Depth of humus (cm) 10.52 0.27 1 8 36 Depth of soil (cm) 93.37 0.6 18 99 99 Stoniness (%) 39.91 0.54 3.96 42.37 65.05 Bulk density of BC (gdm−3) 1267.1 5.5 790.55 1294.9 1522.13 Cation exchange capacity of BC (mmolc kg−1) 23.94 1.28 1.53 12.33 203.25 Base saturation of BC (%) 36.44 1.02 4.33 25.73 100 C /N ratio of BC 16.5 0.35 3.33 14.98 62.45 pH of BC 5.17 0.02 4.36 5.08 7.26 Silt content (%) 19.98 0.57 0 15 85 Clay content (%) 3.16 0.25 0 0 35 Sand content (%) 51.25 0.63 0 55 80 Long-term soil moisture (%) 22.36 0.2 10 20 30 Mean air temperature (◦C) 4.63 0.09 −0.44 5.34 8.47 Total precipitation (mm) 697.87 7.13 392.54 637.11 1154.55 Nitrogen deposition (kgNha−1 y−1) 7.17 0.14 2.35 6.56 17.67 Productivity class (H100, m) 23.61 0.21 12 23 36 Total stand biomass (tCha−1) 56.02 1.39 1.34 51.14 156.52 Total understory biomass (tCha−1) 2.69 0.05 0.96 2.37 6.02 Total litterfall input (tCha−1) 3.17 0.03 1.65 3.07 5.28 2.1.1 Biomass and litterfall estimates For the biomass and litterfall estimation we adopted a stan- dard method of national greenhouse gas inventories for esti- mating soil carbon stock changes (Statistics Finland, 2013). In order to model SOC stocks of forest in equilibrium (not SOC stocks changes), we modified the method by estimating the long-term litterfall of forest in equilibrium. Forest stand biomass was estimated by allometric biomass functions for stem with bark, branch, foliage, stump, coarse roots and fine roots applied to basic tree dimensions (breast height diame- ter, total height of tree, number of trees) of SFI stands (Mark- lund, 1988; Pettersson and Ståhl, 2006; Repola, 2008; Lehto- nen et al., 2016a). In order to simulate “equilibrium” soil car- bon stock, we estimated long-term mean forest biomass, re- ferred to as “equilibrium forest” below. We adopted an observed fraction of photosynthetically ac- tive absorbed radiation (fAPAR; Fig. A1 in Appendix A) as a relative indicator of a site’s capacity to produce biomass (minimum is 0, maximum is 1) by accounting for the forest stand structure, ranging from the absent stand fAPAR = 0 to the closed canopy stand fAPAR = 1, through its major role on limiting of the potential gross primary production (Pel- toniemi et al., 2015). The fAPAR was calculated based on SFI measurements of basic tree dimensions as in Härkönen et al. (2010) and for the main tree species (pine, spruce, de- ciduous) it was well correlated with the stand basal area (Ap- pendix A). The equilibrium forest fAPAR values were assumed to be in a range between the median and the maximum fraction of the observed state forest fAPAR for a given species, latitu- dinal degree, and site productivity index (Appendix A). We selected equilibrium fAPAR as the 70th percentile (fAPAR70) out of a range from the 50th to 95th, because the modeled soil carbon distributions with a litter input from the fAPAR70 biomass agreed best with the measured soil carbon distri- butions (Fig. S2). The fAPAR70 was the estimated 70th per- centile of the observed fraction of absorbed radiation specific for a given species, latitudinal degree, and site productivity index H100 (height of trees at 100 years of age; m; Fig. B1 in Appendix B). The site index H100, that can be translated to a specific productivity (m3 ha−1 yr−1), was for Swedish for- est inventory plots determined based on height development curves and observed site properties by using the methodol- ogy of Hagglund and Lundmark (1977) (Swedish Statistical Yearbook of Forestry, 2014). Instead of modeling of equi- librium biomasses for every tree stand component separately for the species, latitude, and site productivity index, we sim- plified the biomass modeling first by estimating only equilib- Biogeosciences, 13, 4439–4459, 2016 www.biogeosciences.net/13/4439/2016/ B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility 4443 10 15 20 25 56 58 60 62 64 66 68 68 l34 and number of Meteo.station nearest soil samples (t ha−1) Soil carbon La tit ud e (°) Longitude (°) 50 100 150 Figure 1. Geographical locations of meteorological stations with corresponding number of nearest soil samples (n, size of the circle) and their mean measured soil organic carbon stock (tCha−1, color of the circle) across Sweden. rium forest stand structure for the species, latitude, and pro- ductivity (fAPAR70, Table A1 in Appendix A) and secondly by using fAPAR70 with fAPAR biomass models (Table B1 in Appendix B) to estimate the biomass components. We modeled the equilibrium biomass by applying the fit- ted exponential functions between the observed state forest biomass components (stem, branch, foliage, stump, coarse roots, fine roots, estimated by tree stand measurements and the allometric biomass functions) and the observed fraction of absorbed radiation (fAPAR; Appendix B) to the estimated fAPAR70 of the equilibrium forest. The understory vegeta- tion of the equilibrium forest was estimated by applying our ground vegetation models (Appendix C) to the modeled equi- librium forest characteristics, and plot-specific environmen- tal conditions. In order to derive the litter inputs, the annual turnover rate (TR), the fraction of living biomass that is shed onto the ground per year of biomass components, was applied to the modeled biomass components of the equilibrium for- est. The needle litter TR was a linear function of latitude for pine and spruce and a constant for deciduous species (Ågren et al., 2007). The TR of branches and roots were from Mukkonen and Lehtonen (2004) and Lehtonen et al. (2004) and the TR of stump and stem were from Viro (1955), Mälkö- nen (1974, 1977) and Liski et al. (2006). For tree fine roots, we assumed there was a difference between tree species and between southern and northern Sweden. For pine, spruce, and birch the TR fine roots were 0.811, 0.868, and 1.0, re- spectively, as reported by Maidi (2001), Kurz et al. (1996), and Liski et al. (2006). Kleja et al. (2008) and Leppälampi- Kujansuu et al. (2014) reported different fine-root TR for southern (1 and 0.83) and northern Finland (0.5). We interpo- lated TR according to the mean annual temperature gradient between TR of fine roots in the south and the north. The fine- roots TR of 0.811, 0.868, and 1.0 in the warmest southern- most soil plots were thus reduced down to 0.5 in the coldest northernmost soil plots. The understory TR was applied as in Lehtonen et al. (2016b). The major part of the litter input originated from the tree stand biomass components, which were modeled by the non- linear functions with R2 values close to 0.9 (Fig. B1, Ta- bles A1 in Appendix A, and B1 in Appendix B). The linear understory vegetation models had low R2 values (Table C1 in Appendix C). However, when the understory models (Ap- pendix C) were applied only to plots close to equilibrium for- est, as in our application, the R2 values of predicted and ob- served understory components were larger (Fig. S9). In com- parison to major understory litterfall originating from rea- sonably well-predicted dwarf shrubs and mosses (Figs. S9 and S10), the influence of poorer understory models (for herbs, grass, and lichens) was small on predictions of the un- derstory litter and marginal on predictions of the total forest litterfall (Fig. S10). The main improvement on the accuracy of total litter input was achieved by avoiding the confounding effect of management on observed forest state by modeling the biomass/litterfall estimates representing the mean long- term conditions (defined by estimated equilibrium fAPAR70) for small regions (defined by degree of latitude and produc- tivity index for dominant species; Fig. A1 in Appendix A). Thus the estimates accurately reflected the long-term spatial variability in dominant species, nutrient status and climate (Fig. S11) and lacked higher spatial and temporal precision; as attempts for high precision of the estimates applied for the period of the last few thousand years would be uncertain due to high variation of factors affecting plot history. 2.1.2 Correlation analysis Overall our data consist of 3230 soil samples and their car- bon stocks linked to soil physicochemical variables, stand and ground vegetation biomass and litterfall components, and nearest weather station environmental variables. We per- formed the Spearman’s rank correlation analysis between the total soil carbon stock and the other soil variables, site, cli- mate, and vegetation characteristics. As expected the total soil carbon stock most strongly correlated with the measured www.biogeosciences.net/13/4439/2016/ Biogeosciences, 13, 4439–4459, 2016 4444 B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility CEC.BC < 17 CEC.BC < 7.9 parent.material: bedrock,unsortd CN.BC < 16 parent.material: unsortd CEC.BC < 33 N.deposition < 9.9 Bound.H2O.C < 2.4 Humus: no.peat >= 17 >= 7.9 Sorted >= 16 Sorted >= 33 >= 9.9 >= 2.4 Peat 65 n=959 30 % 82 n=909 28 % 130 n=136 4 % 86 n=335 10 % 126 n=182 6 % 104 n=296 9 % 137 n=180 6 % 269 n=8 0 % 144 n=142 4 % 203 n=83 3 % Measured soil carbon stock (t C ha−1)(a) Carbon low medium high medium high medium high extra high extra Moisture dry−fresh fresh moist−fresh fresh moist−fresh fresh fresh fresh fresh moist−fresh Fertility low medium medium medium high low high high medium medium Soil group 1 2 3 4 5 6 7 8 9 10 (b) n=959 n=909 n=136 n=335 n=182 n=296 n=180 n=8 n=142 n=83 Figure 2. (a) Classification/regression tree for the measured soil carbon stock (tCha−1), soil physicochemical properties, and site envi- ronmental characteristics; the cation exchange capacity of BC horizon (CEC.BC, (mmolc kg−1)), the C /N ratio (CN.BC), the nitrogen deposition (N.deposition kgNha−1 y−1), the highly bound soil water of C horizon (bound.H2O.C, %), and soil class variables as type of sorted or unsorted soil parent material and humus type. Note that variables used to calculate the soil carbon stock (bulk density, carbon content, depth, and stoniness) were excluded from the regression tree analysis. The values in the leaves of the tree show for the distinct en- vironmental conditions mean soil carbon stock (tCha−1), number and percentage of samples. (b) The interpretation of 10 physicochemical soil groups of the regression tree model into the levels of carbon, soil moisture, and fertility roughly increasing from left to right. variables used for its calculation, e.g., bulk density, depth of humus and mineral soil, carbon content, and stoniness. These variables were excluded from further regression tree analysis, which aimed to group data according to the processes of soil carbon stock development. 2.1.3 Regression trees In order to organize SOC data into groups according to the physicochemical soil variables and to better understand the nature of measured data, we generated regression trees of SOC stocks by using recursive partitioning (RPART; Th- erneau and Atkinson, 1997). RPART is based on developing decision rules for predicting and cross-validation of continu- ous output of soil carbon stocks (regression tree). The clas- sification tree was built by finding a single variable, which best splits the data into two groups. Each sub-group was re- cursively separated until no improvement could be made to the soil carbon stock estimated by using the split-based re- gression model. The complex resultant regression tree model was cross-validated for a nested set of sub-trees by comput- ing the estimate of soil carbon stock to trim back the full tree. When building the regression tree models, we excluded variables such as bulk density, carbon contents of soil lay- ers, soil depth, and stoniness, since these measured variables were used for determining the total soil carbon stock. The selected variables for the RPART data mining were based on the correlations analysis (see Sect. 2.1.2), the processes of soil organic matter formation (e.g., Husson et al., 2013) and decomposition, and represented the soil categorical vari- ables (sorting of parent material, soil texture, long-term soil moisture, and humus form), soil physicochemical variables Biogeosciences, 13, 4439–4459, 2016 www.biogeosciences.net/13/4439/2016/ B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility 4445 (sand, clay, and silt content, long-term soil moisture, highly bound water, C /N ratio, pH, CEC of organic, B, BC, and C horizons), climatic variables (annual mean air temperature, annual precipitation sum), and stand and site characteristics (tree species coverage of pine, spruce and deciduous, total foliar litter input, productivity class and N deposition). Alter- natively, we also ran regression and classification analysis by excluding all measured soil variables because soil variables are often unavailable for landscape level modeling. The regression tree model separated the measured to- tal SOC stocks (tCha−1) into 10 groups. The cation ex- change capacity of the BC horizon (CEC, mmolc kg−1) di- vided all the samples into two-thirds of lower SOC stock groups (means between 65 and 130 tCha−1) and one-third of larger groups (means between 86 and 269 tCha−1; Fig. 2a). The group of the smallest SOC stock consisted of 959 sam- ples compared to eight samples of the group with the largest SOC stocks. We acknowledge that this is a small distinct group based only on eight observations. However, we did not have any reasons to exclude these data points as outliers. These observations indicated highly fertile conditions (high N deposition, the largest H100 among groups (31 m), second largest litter input, the highest temperature and precipitation on well-drained soil) (Fig. 2, Table S1 in the Supplement). Two-thirds of samples with smaller SOC stocks were sub- divided by CEC and the type sorting of soil parent material (sorted or unsorted). One-third of samples with larger SOC stocks was subdivided by the C /N ratio, CEC, N deposition among others. Roughly generalized, groups from left to right or from 1 to 10 formed a gradient in levels of SOC stock, moisture, nutrient status, and production (Fig. 2, Table S1). The alternative regression tree model was built with vari- ables other than soil properties. The regression tree with the annual mean air temperature, the annual precipitation sum and the percentage of pine trees in the stand, and the ni- trogen deposition separated measured SOC stocks (tCha−1) into five groups (Fig. S3). Colder groups with smaller SOC stocks (with means 67 and 85) had less litter input (below 3 tC ha−1) and a low site productivity index (H100 < 20 m; Table S2). 2.2 Soil carbon stock modeling The Q model (Rolff and Ågren, 1999) is a continuous mech- anistic litter decomposition model describing change of soil organic matter over time. The decomposition rate for the branch, stem, needle, fine root, and woody litter fractions is controlled by the temperature, litter quality, microbial growth, and litter invasion rate. The model has been cali- brated for seven climatic regions of Sweden in order to ac- count for Swedish temperature and precipitation gradients (Ortiz et al., 2011; Table 3). The Q model was applied in several studies of SOC stock and change estimation in Swe- den (e.g., Stendahl et al., 2010; Ortiz et al., 2013; Ågren et al., 2007). The Q model was run for seven Swedish cli- matic regions (Ortiz et al., 2011). The mean regional parame- terization from the calibration of the 2011 Q model was used for the plot simulations. Thus, the simulations in each re- gion represent variations in climate and litter input and not parameter variations. The equilibrium soil carbon stocks are estimated in the model using the equation for equilibrium soil carbon stock, which is derived from the decomposition functions with constant amounts and quality of litter input. The Yasso07 model (Tuomi et al., 2009, 2011) is one of the most widely applied SOC models. The model was calibrated based on almost 10 000 measurements of litter decomposi- tion from Europe, North and South America (Table 3). The required annual inputs of litterfall, its size and chemical com- position, temperature, and precipitation determine the de- composition and sequestration rates of soil organic matter. Yasso07 estimates SOC stock to a depth of 1 m (organic and mineral layers), change of SOC stock, and heterotrophic soil respiration. Species-specific chemical composition of differ- ent litter compartments of Yasso07 were used according to Liski et al. (2009). The initial soil organic matter of Yasso07 was zero. The simulated soil carbon stock corresponding to equilibrium between the litter input and decomposition was achieved by a Yasso07 spin-up run of 10 000 years. Yasso07 runs used litter inputs of the equilibrium forest biomasses (see Sect. 2.1.1) and climate variables (annual air tempera- ture, monthly temperature amplitude, and annual precipita- tion). The global parameter values of decomposition rates, flow rates, and other dependencies of the Yasso07 soil car- bon model were adopted from Tuomi et al. (2011) and the estimates of Yasso07 SOC stocks were used in comparison with measurements and other models. We did not use the SOC stocks simulated with the more recent Yasso07 param- eters based on the litter decomposition data from the Nordic countries (Rantakari et al., 2012), because the SOC stocks simulated with the global parameter values produced a better fit with SFSI measurements. The CENTURY mathematical model originally devel- oped for grassland systems (Parton et al., 1987, 1992) has been since modified for various ecosystems including boreal forests (Nalder and Wein, 2006). The CENTURY is also one of the most widely applied models. The soil organic mat- ter in the model consists of active, slow, and passive pools, which have different TR (Table 3). The decomposition rates are modified by temperature and moisture, and in addition the decomposition rates of the slow and passive pools rely on lignin to N and C to N ratios, while the active pool de- composition rate relies on soil texture. The model simulates soil organic matter to a depth of 20 cm. The model simu- lates plant production and pools of living biomass, while TR for biomass pools determine the litterfall inputs to soil. To compare the performance of the soil sub-model with other soil carbon dynamics models, Q and Yasso07, we only used the CENTURY soil sub-model. We used the same litterfall inputs as used by the Q and Yasso07 simulations, which were estimated by our litterfall modeling (see Sect. 2.1.1). www.biogeosciences.net/13/4439/2016/ Biogeosciences, 13, 4439–4459, 2016 4446 B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility Table 3. Description of models and data inputs relevant for this study. Model Yasso07 Q CENTURY v4.0 soil submodel Time step Year Year Month Parameterization Global Scandinavian Combined global with site specific Carbon pools Labile (acid -, water -, and ethanol- sol- uble and non-soluble), recalcitrant (hu- mus) Cohorts (foliage, stems, branches, coarse roots, fine roots, “grass”), soil organic Litter (surface structural and metabolic, belowground str. and met.), surface microbial, soil organic matter (active, slow and passive) Biomass Biomass components estimated by allometric biomass functions and provided stand data for litter input estimation Litter amount Annual or monthly fractions of biomass components (species specific, same total litter inputs for all models) Litter quality Literature-based solubilities Estimated cohorts qualities C /N ratios and lignin /N ratios Temperature air Annual mean, monthly amplitude Annual mean Max and min monthly mean Precipitation Annual total – Monthly total Soil properties – – Bulk density, sand, silt, and clay content Soil depth (m) 1 – 0.2 The litter inputs reflected N deposition and site productivity (Fig. S11). For CENTURY we adopted general parameters from the parameter file “tree.100”, parameters of site “AND H_J_ANDREWS” for conifers, and site “CWT Coweeta” for deciduous trees. The N dynamics in CENTURY sub-model included tuning site-specific parameters of topsoil mineral N relative to N deposition (Throop et al., 2004) and reduction of C /N ratio of the litterfall up to 15 % for most productive sites (Merilä et al., 2014). We also accounted for site-specific soil drainage by varying its parameter between 1 and 0.6 rel- ative to long-term soil water content ranging between 10 and 50 % (Raich et al., 2000). The CENTURY SOC stocks sim- ulation were run with equilibrium forest litter inputs, site- specific C /N ratios of litterfall, site-specific soil parameters (specific bulk density, sand, silt, and clay content, mineral N in topsoil, and drainage) and climate variables (monthly air temperature, and monthly precipitation). In order to account for the deep soil carbon (Jobbágy and Jackson, 2000), we scaled CENTURY estimates representing the topsoil horizon by adding 40 % of estimated site-specific SOC stock. The simulated equilibrium SOC stocks were estimated by a spin- up run of 5000 years. The number of years to reach equi- librium (equilibrium between the litter input and decomposi- tion) was sought empirically on 100 random sites, and differs from Yasso07 and Q models. 3 Results The distributions of Yasso07, Q, and CENTURY model es- timates of total SOC stocks (tCha−1) were in agreement for two-thirds of the measured data with lower SOC stock (Fig. 3, distributions of groups 1, 2, and 4). The remain- ing one-third of SOC data were underestimated by models. This one-third of data were separated into seven physico- chemical soil groups (means of groups ranging from 104 to exceptionally large 269 tCha−1, see Fig. 3, distributions of groups 3, and 5–10). The linear regression of mean levels of all 10 physicochemical soil groups (weighted by the number of samples in each group) between the modeled and mea- sured SOC stocks showed smaller underestimation of CEN- TURY compared to Yasso07 and Q models (Fig. 4). The weighted root mean square error (RMSE) was 27.5 (tCha−1) for CENTURY and 31.6 and 38.8 for Yasso07 and Q, respec- tively. The proportion of explained variance was larger for Q (r2 = 0.58) than for Yasso07 and CENTURY (r2 = 0.42 and 0.32; Fig. 4). The deviation of the distributions of CENTURY SOC stocks, simulated using soil bulk density, sand, silt, and clay content, were lower than those of Yasso07 and Q esti- mates for 10 physicochemical soil groups (Fig. 3). Account- ing for site-specific soil texture (clay, silt, and sand content) and structure (bulk density) by the CENTURY model im- proved SOC stock estimates for fertile sites with high clay content, but not for sites with high N deposition. Varying CENTURY parameters of site-specific topsoil mineral nitro- gen and C /N ratio of the litterfall showed that this impact on SOC stocks estimates was small in comparison to sensitivity of SOC stock estimates to litterfall (Fig. S12). The applica- tion of site-specific drainage on our mostly well-drained soils showed minor impact on estimated CENTURY SOC stocks. As expected, the models clearly showed less variation than the measurements. The shift of the mean values from the cen- ter of distribution, the width of confidence intervals of means, and the width of the tails of distributions were clearly larger for the measurements than for the modeled estimates (Fig. 3). The modeled distributions agreed for the poor–medium fer- tility soils with low and medium measured SOC stocks, low and medium CEC, unsorted parent material, low tempera- tures, and low production (groups 1, 2, and 4; Figs. 2, 3, Table S1). Disagreement between modeled and measured SOC stock distributions were formed on fertile soils with sorted parent material (groups 3 and 5), soils with higher water content (groups 3, 5, and 10), where nitrogen depo- sition was large (groups 7 and 8), and where CEC was me- dian or large (Figs. 2, 3). The largest deviation between the Biogeosciences, 13, 4439–4459, 2016 www.biogeosciences.net/13/4439/2016/ B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility 4447 0 10 0 20 0 30 0 Measurements Yasso07 CENTURY Q 0 10 0 20 0 30 0 n=959 n=909 n=136 n=335 n=182 n=296 n=180 n=8 n=142 n=83 1 2 3 4 5 6 7 8 9 10 So il ca rb o n st o ck t ha ( − 1 ) Soil groups Figure 3. Bean plot of distributions of the soil carbon (tCha−1) measurements (gray fill) and estimates for 10 physicochemical groups. The full and dashed horizontal lines represent the group means and their confidence intervals. The n is the number of samples. For description of group levels of SOC stocks, moisture, and fertility see Fig. 2 and Table S1. measured and modeled distributions was found for the rela- tively small physicochemical groups of soils (3 %) typical for highly bound water and peat humus types (groups 8 and 10; Figs. 2, 3). The distributions of measured total SOC stocks (tC ha−1) generally increased for the groups with higher nu- trient status (Figs. 3, S4). The distributions of SOC stocks in mineral soil were larger than those in humus horizon, and distributions of mineral SOC stocks increased with fertility slightly more than distributions of SOC stocks in humus hori- zon (Fig. S4). After excluding all the soil physicochemical characteris- tics from the recursive partitioning, the SOC stock distribu- tions of five group regression tree models (Fig. S3, Table S2) were in agreement between the measurements and model es- timates for three groups (77 % of samples) and deviated for two groups (23 %; Fig. S5). The modeled SOC stock dis- tributions agreed with measurements for all models on sites with low annual temperatures < 3 ◦C in northern sites (low- C.cold.pine, low-C.cold.other) and for warmer conditions in middle Sweden on sites with low nitrogen deposition and me- dian SOC stocks (Fig. S5). However, the models underesti- mated SOC stocks on sites with high (> 10 kgNha−1 y−1) N deposition (21 % of samples) and on sites with warm and dry climate (2 % of samples; Fig. S5). The variation of density functions of modeled SOC stocks for 10 physicochemical groups (Fig. 3) was similar to the variation of the total annual plant litter input (tCha−1; Fig. S6) indicating that litterfall was the main driver of SOC accumulation in the models . The mean levels of annual plant litter input and mean SOC stocks for 10 soil groups were more strongly correlated for Yasso07 and Q models (with r2 values 0.86 and 0.96, respectively) than for CENTURY (r2 = 0.52). Although, models performed reasonably well for the largest soil groups of nutrient and production levels (Figs. 3 and 4), none of the models was able to predict vari- ation of individual samples (Fig. S7). The model estimates were well correlated between Yasso07 and CENTURY with r2 ranging from 45 to 73 % for individual samples of 10 soil groups, whereas the correlations of estimates between Q and the other two models were lower (Fig. S8). 4 Discussion 4.1 SOC stock distributions linked to mechanisms of SOM stabilization It has been suggested that process-based soil carbon mod- els with the current formulation lacking major soil environ- mental and biological controls of decomposition would fail for conditions where these controls predominate (Schmidt et al., 2011; Averill et al., 2014). Even so, the effect of the soil properties on SOC stocks, e.g., soil nutrient status in the widely used models such as Yasso07, Q, and CENTURY, have not previously been quantitatively evaluated. We found that in comparison with Swedish forest soil inventory data, the models based on the amount and quality of inherent struc- tural properties of plant litter (Q, Yasso07, and CENTURY) produced accurate SOC stock estimates for two-thirds of northern boreal forest soils in Sweden. Two-thirds of the distributions of SOC stocks measurements of SFSI agreed with distributions of SOC stock estimates of the Q, Yasso07, and CENTURY soil carbon models (Fig. 3, distributions of groups 1, 2, and 4). However, the SOC stocks underestima- tion by these models for one-third of the data (Fig. 3, distri- www.biogeosciences.net/13/4439/2016/ Biogeosciences, 13, 4439–4459, 2016 4448 B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility Figure 4. Scatter plot between mean measured and mean mod- eled soil organic carbon stocks (tCha−1) for 10 physicochemical groups for Yasso07, CENTURY and Q models. Data were fitted with weighted linear regression (lines). The number of samples in each group was used as weights for fitting and also as weights for the weighted mean of squared differences between the modeled and measured values (MSE, tCha−1). The RMSE is the square root of MSE. The r2 is the proportion of explained variance. The p value is the calculated probability that the fit is significant. butions of groups 3, and 5–10) indicated that some drivers other than molecular structure, especially site nutrient status, play an important role in higher SOC stocks sequestration. Some level of deviation from measurements and poorly explained spatial variation (Fig. S7) was expected from the uncertainties of the SOC measurements, annual plant lit- ter inputs and climate variability for the model SOC stock change estimates (Ortiz et al., 2013; Lehtonen and Heikki- nen, 2015). For the long-term SOC stock development the model uncertainties are less known than for the short-term litter decomposition. Previously reported fine-scale compari- son also showed poor agreement between Earth system mod- els and the Northern Circumpolar Soil Carbon Database (Todd-Brown et al., 2013), although drivers of the devia- tion still remained open. Our results showed that if models strongly depend on the litter inputs (Fig. S6) then the spa- tial differences between measured and modeled SOC stock distributions could be linked to sites with rich nutrient status through cation exchange capacity, C /N ratio, N deposition, drainage (sorting of parent material) among other factors (Figs. 2 and 3). Additionally, when the soil properties were excluded from the regression, the estimates of SOC stocks also deviated for the fertile groups (Fig. S5). However, the rich nutrient status for these groups was linked to differences in species composition, N deposition, and climate (tempera- ture, precipitation) instead of soil properties (Fig. S3). Larger net soil carbon accumulation in nutrient-rich sites could be attributed to the relative differences in litterfall com- ponents (relatively more leaves and branches with higher N content than fine roots) and, to higher N availability and car- bon use efficiency of decomposers, reduction of respiration per unit of C uptake (Ågren et al., 2001; Manzoni et al., 2012; Fernández-Martínez et al., 2014). The largest deviation be- tween measured and modeled data in our study was found for fertile presumably N rich and fresh to fresh-moist sites. The soils with large N deposition were also highly productive and showed high to exceptionally high SOC stocks (Figs. 2, 3, soil groups 7 and 8). This was in agreement with fertiliza- tion and modeling study of Franklin et al. (2003) showing an increase in soil C accumulation with N addition. Our forest biomass and litterfall estimates were based on for- est inventory and modeling, but the site nutrient status and N deposition was only partially reflected in the amount of biomass/litterfall (Fig. S11) and its quality. The quality was only reflected through the biochemical differences between species and plant litter components. The relative differences between the biomass/litterfall components or between C /N ratios of litterfall in relation to site fertility are not accounted for by the current biomass models, but soil fertility could be considered in an attempt of SOC stock modeling (in- cluded in CENTURY but missing in Yasso07 and Q models). For example the proportion of acid-, water-, and ethanol- soluble and non-soluble litter inputs for Yasso07 could be re-evaluated by allowing it to vary depending on site fertil- ity, in addition to currently used variation specific for species and the litter components. Although CENTURY SOC stocks were sensitive to the amount of clay, the variation of topsoil mineral N and C /N ratio of litterfall did not improved SOC stock predictions for sites with high N deposition (Fig. 3 and Table S1). The litter decomposition and SOC stabilization rates in Yasso07, Q, and CENTURY based on the litter quality “re- calcitrance” originating from the litter bag mass loss mea- surements have major drawbacks. The mass loss from the litter bags is assumed to be fully mineralized, although the litterbags are subjected to non-negligible leaching (Rantakari et al., 2012; Kammer and Hagedorn, 2011). The SOC stabi- lization represented in models by the remaining litter mass is thus underestimated due to the fraction of particulate organic matter and dissolved organic carbon that is lost from the lit- terbags but later immobilized, e.g., through organo-mineral stabilization. The use of stable isotopes seems to determine the field carbon mineralization and accumulation rates from the labile (high C quality and N concentration) or recalcitrant (low C quality and N concentration) litter more accurately than litter bags (Kammer and Hagedorn, 2011). A higher amount of more recalcitrant fine roots compared to more labile leaves (Xia et al., 2015) heavily increased the soil carbon sequestration in CENTURY model simulations, Biogeosciences, 13, 4439–4459, 2016 www.biogeosciences.net/13/4439/2016/ B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility 4449 which was in line with McCormack et al. (2015). Though, the contribution of fine roots to SOC stabilization is still not settled due to the significant role of mycorrhizal fungi in SOC accumulation (Averill et al., 2014; Orwin et al., 2011). Xia et al. (2015) claimed that more recalcitrant fine roots contribute to stable SOC more than leaf litter, because fine roots degrade slower. This would be supported by the fact that the derivatives of fine roots from degradation by fungi are more stable than the derivatives of leaves from degra- dation by microbes. However, more recalcitrant plant litter has been also suggested to stabilize fewer SOC stocks (Kam- mer and Hagedorn, 2011). This is a result of recalcitrant lit- ter satisfying less of the microbial N demands promoting respiration and reducing the long-term production of micro- bial products, precursors for the organo-mineral stabilization (Cotrufo et al., 2013, Castellano et al., 2015). According to the microbial efficiency-matrix (MEM) stabilization mecha- nism (Cotrufo et al., 2013) fertile sites with relatively more labile plant litter, but with larger absolute production and larger microbial activity than poor sites, would in long-term stabilize more carbon through organo-mineral stabilization. Our results supported MEM stabilization theory by showing larger carbon stocks in mineral soil than in humus horizon, and by relatively more SOC stocks in mineral soil in fertile groups than in poor conditions (Fig. S4). Expanding on the CENTURY model structure, the MySCaN model incorporating the organic nutrient uptake by mycorrhizal fungi estimated a positive effect on SOC accumulation, relatively larger in poor than in fertile sites (Orwin et al., 2011). Therefore, not accounting for the or- ganic nutrient uptake by mycorrhizal fungi by the Yasso07, Q, and CENTURY models probably led to the underestima- tion of SOC stocks in sites with higher nutrient status. This hypothesis needs to be tested in further studies. We did not have all input data and the source code to include MySCaN into our model intercomparison. The spatial trends of N and P data of litter in Sweden that would be needed to make such a study were not available. However, adjusting biomass turnover rates, used for the litter input estimation, in depen- dence to site fertility would lead into larger inputs for fer- tile sites and increased SOC stock accumulation as a result of increasing plant productivity and inputs. It is well estab- lished that SOM increases soil fertility by improving the soil water and nutrient holding capacity; recycling of SOM in- creases CEC, humic substances and nutrient availability for plant resulting in larger biomass/litter production (Zandonadi et al., 2013). As an alternative to adjusting turnover rates with site fertility, we suggest that a feedback link in models between increasing fertility due to SOC stock accumulation (e.g., due to increased CEC relative to humus, increased ni- trogen availability), increasing litter inputs, and reduced rates of SOC decomposition per unit of litter input (e.g., through satisfying more microbial N demand with less respiration, limited oxygen in increased moisture conditions) would also increase SOC stock accumulation. Increased moisture and more frequent water saturation due to SOC accumulation limits soil oxygen availability and slows rates of microbial decomposition, which increases the rate of SOC stabilization. Our results, which were derived from mostly well-drained soils, suggest that measured high SOC stocks may be partly caused by reduction of decom- position at increased water content (Fig. 2). The CENTURY model has an optional function that represents the reduction of decomposition caused by anaerobic conditions. The func- tion becomes active when a controlling parameter, “drain”, is changed, and the value of the parameter has to be arbitrar- ily determined through parameter fitting against SOC data (e.g., Raich et al., 2000). However, this function was meant for anaerobic conditions in poorly drained soils; therefore, it was not applicable to the prevailing conditions of our sites. Accounting for drainage only on some sites slightly affected decomposition, when precipitation increased and potential evapotranspiration decreased in late spring or early autumn. Water availability affecting soil fertility and SOC formation is beside climate also affected by topography (Clarholm and Skyllberg, 2013), which was not accounted for by CEN- TURY. Detailed modeling of soil water conditions requires specific functions and many parameters, which are not in- cluded in simpler SOC models like Q and Yasso07. However, appropriate modeling of soil water conditions and reduction of decomposition in wet conditions (not necessarily at satu- ration) would potentially improve the performance of SOC models in particular for soils with high SOC stocks. 4.2 Intercomparison of models The similarities between the variations of modeled SOC stocks and litterfall inputs for the soil groups with different fertilities (Figs. 3, S6) could be expected for the Yasso07 and Q models, which ignore the soil properties. These models run organic matter decomposition and humus stabilization with litterfall, temperature, and/or precipitations input data. Litter quality as input in Yasso07 and Q implicitly includes some information on soil properties, but as we saw litter quality hardly mapped any of soil fertility. Although, the impact of soil properties on the estimates was seen in the more complex CENTURY model for sites with high clay content, the SOC stock of sites with high N deposition were underestimated. The CENTURY model depended less on the amount of litter input. In testing multiple soil carbon models with the same litter inputs, Palosuo et al. (2012) observed larger variation in modeled SOC stocks at the early stage of the litter decom- position (10 years) but later on at 100 years the variation de- creased. Although the variations of SOC stocks were similar between the models, the estimated CENTURY SOC stocks distributions were lower than the Yasso07 estimates when we did not accounted for deep soil carbon. CENTURY in its original configuration simulated SOC stock up to 20 cm soil depth (Metherell, 1993), whereas the Yasso07, Q, and mea- sured SOC stock data represented up to 100 cm of the soil www.biogeosciences.net/13/4439/2016/ Biogeosciences, 13, 4439–4459, 2016 4450 B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility (Tuomi et al., 2009; Stendahl et al., 2010). In Yasso07 model parameters were calibrated based on soil age chronosequence data of SOC stocks for soil depths up to 30 cm, which was assumed to represent 60 % of the total SOC stocks up to 100 cm soil depth (Liski et al., 1998, 2005 as cited by Tuomi et al., 2009). Therefore, when 40 % of the missing deep car- bon (Jobbágy and Jackson, 2000) were added on top of the original CENTURY estimates as was done when calibrating Yasso07, the SOC stock levels for CENTURY were larger than those for the Yasso07 and Q models. Although estimated SOC stocks of CENTURY were gen- erally larger than those of Yasso07, the correlation between CENTURY and Yasso07 estimates was stronger than for Q model compared to two other models (Fig. S8). The reason was probably similar global parameterizations of Yasso07 and CENTURY whereas Q was specifically parameterized and applied for the regions in Sweden (Ågren and Hyvö- nen, 2003; Ortiz et al., 2013). Furthermore the Q model SOC stock estimates were more sensitive to differences in species coverage e.g., to pine and spruce (Ågren and Hyvö- nen, 2003) and formed two distinct point cloud distribu- tions (one for pine and broadleaves, the other for spruce) when compared with the CENTURY or Yasso07 estimates (Fig. S8). In spite of similarities in Yasso07 and CENTURY SOC stocks estimates, Yasso07 was more sensitive to species coverage through species-specific litterfall solubility (Liski et al., 2009) than CENTURY, which treated conifers in a single group (Metherell et al., 1993). Pine and other species (spruce) coverage was shown to affect measured low and me- dian SOC stocks of colder climate if the soil properties were not considered (Fig. S5). Therefore, the pattern of increased accumulation of SOC stock on sites with larger spruce cov- erage partially observed in distribution of Yasso07 estimates, and missing in the CENTURY estimates, could be related to the slightly lower solubility/decomposability of spruce com- pared to pine litterfall. However, the CENTURY model SOC stocks were also highly sensitive to accurate estimation of fine-root litterfall (McCormack et al., 2015) typically in- creasing with colder climate and increasing the C /N ratio of the organic layer (Lehtonen et al., 2016a), which is driven by the dominant tree species (Cools et al., 2014). Large SOC stock measurements on sites with high long- term nitrogen deposition over 10kgNha−1 y−1 (Figs. 3 and S4) were underestimated by the Q, Yasso07, and CEN- TURY models. A positive correlation between nitrogen de- position and SOC stocks measurements in Sweden had been previously reported by Olsson et al. (2009), and the model- ing study by Svensson et al. (2008) indicated that Swedish soil carbon was decreasing in the north and increasing in the south mainly as a result of different nitrogen inputs. The Q and Yasso07 models do not have nitrogen processes. As for CENTURY, it is reported that large N input could enhance plant productivity and then increase SOC (Raich et al., 2000). The purpose of our study was to evaluate the performance of soil carbon models against the SOC data using the same litter input, and the feedback of nitrogen input to plant productiv- ity was primarily included in this study indirectly, through estimated equilibrium litter input based on site productivity index, which strongly correlated with N deposition (Figs. A1 in Appendix A and S11). In spite of a slight increase of SOC stock estimates when CENTURY accounted for the site- specific topsoil mineral N, the C /N ratio of litterfall (Fig. S12), in sites with large N deposition CENTURY still un- derestimated. However, as in the case of drainage discussed earlier, the CENTURY incorporates more detailed processes than the relatively simpler soil carbon models do, Q and Yasso07, and hence the CENTURY could potentially repro- duce a wider range of SOC stocks if it was parameterized with more detailed data. 5 Conclusions In this study we presented the reasons to re-evaluate the con- nection between the soil nutrient status and performance of widely applied soil carbon models (Yasso07, Q, and CEN- TURY). As previously described in detail, our simulation was based on the widely used process-based SOC models, accurate driving data including litter inputs, and massive SOC data points (Swedish inventory data, N= 3230). The models differed in the main controls and functions and their performance was expected to depend on model complexity (CENTURY outperforming Q and Yasso07). The intercom- parison of SOC stocks between Yasso07, Q, and CENTURY models and Swedish soil carbon inventory data revealed that these process-based mathematical models developed for pre- dicting short-term SOC stock changes can all in their cur- rent state predict accurate long-term SOC stocks for most soils. However, in medium–highly productive sites of south- ern Sweden for conditions where the high nutrient status pre- dominates soil carbon accumulation, the models with their current formulation (lacking nutrient status-related controls of decomposition and SOC accumulation) underestimated SOC stocks. The estimates of CENTURY fitted generally better to measurements than those of Yasso07 and the Q model. Although the Yasso07 model, which requires fewer parameters and less input data, showed similar performance than CENTURY, except for sites with high clay content. Through the intercomparison of three different widely used SOC models with massive data points, we identified that re-evaluation of the impact of nutrient status would im- prove the model development towards their accuracy. Partic- ularly, the relationship between the soil nutrient status and the mechanism of soil organo-mineral carbon stabilization needs to be re-evaluated, because larger SOC stocks were found more in the mineral than in the humus soil horizon. We suggest evaluating enhanced microbial transformation of soil organic matter and the mycorrhizal organic nutrient up- take in relation to larger plant biomass/litter production in nutrient-rich sites resulting in higher SOC stock accumula- Biogeosciences, 13, 4439–4459, 2016 www.biogeosciences.net/13/4439/2016/ B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility 4451 tion in deeper soil layers. In addition to the organo-mineral carbon stabilization, we also suggest further model develop- ment accounting for the soil nutrient status through evaluat- ing the effect of topography on sorting of the parent material, and its silt and clay complexes. Our study is very useful for developing accurate soil car- bon and Earth system models. Furthermore, developing ac- curate models that would account for the soil nutrient sta- tus as one of the key controls affecting the soil organic mat- ter production and SOC stabilization improves estimation of feedback of global warming on SOC stock temperature sen- sitivity and soil CO2 efflux, national reporting of soil carbon stock changes for UNFCCC, and implications of decisions mitigating the climate change effects on soil carbon stocks. 6 Data availability The source codes of the Yasso07, Q, and CENTURY models used in this paper are available through the Supplement. Data used in this study can be available directly by contacting the authors. www.biogeosciences.net/13/4439/2016/ Biogeosciences, 13, 4439–4459, 2016 4452 B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility Appendix A: Models of fraction of absorbed radiation for observed and equilibrium forest The fraction of photosynthetically active absorbed radiation (fAPAR) for the observed state forest was calculated based on measurements of Swedish forest inventory as in Härkönen et al. (2010). For the main tree species fAPAR was also well correlated with the stand basal area (r2 was 0.85, 0.86, and 0.88 for pine, spruce, and deciduous stands, respectively, co- efficients of regressions in Table A1 in Appendix A). The ob- served state forest fAPAR varied between 0 and a maximum close to 1 (Fig. A1 in Appendix A). The equilibrium forest fAPAR values were assumed to be ranging between the median and the maximum fraction of observed state forest fAPAR for given species, latitudinal de- gree, and site productivity index (indicated by the height of largest tress at 100 years of stands age). The equilibrium for- est fAPAR values were set to 70th percentile of maximum (fAPAR70) for given species, latitudinal degree, and site pro- ductivity index. We selected 70th percentile from the range between 50th and 95th, because the modeled soil carbon dis- tributions with the litter input from biomass of fAPAR70 best agreed with measured soil carbon distributions (Fig. S2). Table A1. Parameter estimates and their standard errors of the fAPAR regressions with the stand basal area (BA, m2 ha−1), and the fAPAR70LAT and fAPAR70H100 regressions with the latitude (LAT, ◦) and with the site productivity index (H100, m) for Scots pine, Norway spruce, and deciduous stands. fAPAR = a×BA/(b+BA) a±SE b±SE c±SE adj.R2 Pine 1.00± 0.03 11.75± 0.81 0.85 Spruce 1.17± 0.03 10.67± 0.87 0.86 Deciduous 1.13± 0.06 7.41± 1.15 0.88 fAPAR70LAT = LAT/(a+ b×LAT)+ c Pine −9976± 3691a 143± 54b 0.72± 0.02 0.92 Spruce −2689± 3507c 35± 50d 0.97± 0.09 0.74 fAPAR70LAT = a+ b×LAT Deciduous 1.36± 0.28 −0.01± 0.01e 0.26 fAPAR70H100 = a× e(b/H100) Pine 0.86± 0.02 −5.22± 0.41 0.89 Spruce 0.97± 0.01 −2.85± 0.22 0.86 Deciduous 0.94± 0.02 −2.63± 0.50 0.51 p < 0.001 for all parameters except for a 0.023, b 0.024, c 0.461, d 0.498, and e 0.076. The fAPAR70 values specific for pine, spruce, and decidu- ous stands were modeled with latitude and site productiv- ity index (H100) in two steps. First, the fAPAR70LAT and the fAPAR70H100 values were modeled separately by regression models with latitude and with site productivity index (Ta- ble A1 in Appendix A). Second, the fAPAR70LAT was reduced by the difference between the fAPAR70H100 and the maximum fAPAR70H100 (fAPAR70 = fAPAR70LAT+ fAPAR70H100− max- imum fAPAR70H100). The fAPAR70 equaled the fAPAR70LAT only for the maximum site productivity index, otherwise it was reduced. Biogeosciences, 13, 4439–4459, 2016 www.biogeosciences.net/13/4439/2016/ B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility 4453 Figure A1. The fraction of photosynthetically active absorbed radiation (fAPAR; estimated as in Härkönen et al., 2010) observed fAPAR and equilibrium fAPAR (fAPAR70, set to 70th percentile of maximum fAPAR for given species, latitudinal degree, and site productivity index). Panels (a), (b), and (c) show relation between fAPAR and latitude (◦) for forest stands dominated by Scots pine, Norway spruce, and deciduous species, whereas panels (d), (e), and (f) show relation between fAPAR and site index H100 (height of dominant trees at 100 years in meters). www.biogeosciences.net/13/4439/2016/ Biogeosciences, 13, 4439–4459, 2016 4454 B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility Appendix B: Models of forest dry weight biomass with fAPAR We fitted species-specific exponential regression models between the biomass components (stem, branch, foliage, stump, coarse roots, fine roots, all in kgha−1) of observed state forest and the observed fraction of absorbed radiation (fAPAR) (statistics of the regression models in Table B1 in Appendix B). The biomass components derived with allo- metric models (measured) and those derived with fAPAR models (modeled) showed strong correlations (Fig. B1 in Appendix B). In order to model the long-term mean forest biomass “equilibrium forest biomass” we applied the fAPAR biomass models to the modeled fAPAR70 values. Table B1. Parameter estimates and their standard errors for the co- efficients of the dry weight biomass (kgha−1) models with the frac- tion of absorbed radiation (y = abfAPAR ) for Scots pine, Norway spruce, and deciduous stands. y = abfAPAR Species a±SE b±SE adj.R2 Branch pine 610± 21 122± 6 0.92 spruce 877± 35 54± 2 0.92 deciduous 290± 26 156± 16 0.89 Fine root pine 422± 13 21± 1 0.84 spruce 317± 14 15± 1 0.80 deciduous 453± 28 14± 1 0.82 Foliage pine 361± 24 86± 8 0.71 spruce 766± 40 33± 2 0.83 deciduous 141± 28 71± 16 0.56 Root pine 703± 26 183± 10 0.92 spruce 629± 32 113± 7 0.90 deciduous 359± 33 150± 16 0.89 Stem and bark pine 1793± 84 254± 17 0.89 spruce 974± 72 229± 19 0.86 deciduous 972± 98 161± 18 0.88 Stump pine 232± 10 214± 13 0.89 spruce 171± 10 129± 9 0.88 deciduous 80± 8 216± 25 0.87 p < 0.001 for all parameters. Biogeosciences, 13, 4439–4459, 2016 www.biogeosciences.net/13/4439/2016/ B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility 4455 0 50 100 150 0 50 10 0 15 0 Stem+bark Modeled (tC ha−1) M ea su re d (tC ha − 1 ) r2 = 0.95 0 10 20 30 0 10 20 30 Branch Modeled (tC ha−1) M ea su re d (tC ha − 1 ) r2 = 0.97 0 5 10 15 0 5 10 15 Foliage Modeled (tC ha−1) M ea su re d (tC ha − 1 ) r2 = 0.94 0 5 10 15 0 5 10 15 Stump Modeled (tC ha−1) M e a su re d tC ha ( − 1 ) r2 = 0.96 0 10 30 50 0 10 20 30 40 50 Root Modeled (tC ha−1) M e a su re d tC ha ( − 1 ) r2 = 0.97 0 1 2 3 4 0 1 2 3 4 Fine root Modeled (tC ha−1) M ea su re d (tC ha − 1 ) r2 = 0.96 Figure B1. Scatter plots (n= 3698 in each panel) for the dry weight tree biomass components (tCha−1) between “modeled” (estimated based on fraction of absorbed radiation, fAPAR, and our fAPAR models) and “measured” (estimated based on basic tree dimensions and allometric biomass models). The r2 values represent the coefficient of determination indicating how close the modeled values fit the measured values. www.biogeosciences.net/13/4439/2016/ Biogeosciences, 13, 4439–4459, 2016 4456 B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility Appendix C: Models of understory vegetation We used Swedish forest inventory ground vegetation cov- erage (%) data visually monitored between 1993 and 2002 on 2440 plots around Sweden with altogether 4472 observa- tions separately for species of forest floor vegetation or their classes (Table S3). In order to derive the ground vegetation biomass and to apply the coverage/biomass conversion func- tions (Lehtonen et al., 2016), we grouped the species cov- erage observations into five functional types (dwarf shrubs, herbs, grasses, moss, and lichen; Table S3). The applied coverage/biomass conversion functions estimated separately the above- and below-ground biomass components for dwarf shrubs, herbs, and grasses, and total biomass for moss, and lichen. Except the understory coverage, the forest inventory data also contained basic tree dimensions (diameter and height of trees) and stand variables (species dominance, age, basal area, site productivity class indicated by the height of largest tress at 100 years of stands age), and also we linked the plots by their closest proximity to SMHI weather stations with weather data (air temperature, precipitation) and location at- tributes of the weather stations (latitude, longitude, altitude). Table C1. Parameter estimates and their standard errors for the coefficients of the forest understory vegetation dry weight biomass (kgha−1) models (Eq. C1) for functional types (1 – dwarf shrubs, 2 – herbs, 3 – grasses, 4 – mosses, and 5 – lichens) with intercept (a) and n – number of predictors (b1 – age (years), b2 – basal area (m2 ha−1), b3 – annual air temperature (◦C), b4 – latitude (◦), b5 – H100 (height of trees at 100 years of age, m), b6 – H100 of spruce trees (m), b7 – H100 of pine trees (m), b8 – pine dominance (0/1), and b9 – spruce dominance (0/1)). For the latin names of species included into understory functional types see Table S3. W a±E b1±SE b2±SE b3±SE b4±SE b5±SE b6±SE b7±SE b8±SE b9±SE adj.R2 Above 1 24.28± 0.32 0.13± 0.01 -0.43± 0.02 7.13± 0.33 0.29 ground 2 −82.13± 6.8 −0.1± 0.1a 1.23± 0.1 0.77± 0.03 0.12 3 4.07± 0.30 −0.16± 0.01 0.27± 0.01 −1.36± 0.15 0.21 4 32.9± 0.62 −0.78± 0.04 0.48± 0.06 3.66± 0.3 5.76± 0.29 0.22 5 19.91± 0.57 −0.13± 0.01 −0.45± 0.02 6.31± 0.29 0.25 total 43.68± 0.29 0.12± 0.01 −0.41± 0.01 6.34± 0.3 0.30 Below 1 −256.3± 3.5 0.1± 0.01 −0.35± 0.02 5.05± 0.06 8.56± 0.35 0.75 ground 2 −89.34± 7.85 −0.03± 0.1b 1.4± 0.12 0.78± 0.04 −4.97± 0.27 0.19 3 5.97± 0.37 −0.19± 0.01 0.32± 0.01 −1.78± 0.19 0.21 total −251.9± 3.3 −0.2± 0.01 5.15± 0.05 0.7 Total −222.7± 4.0 0.12± 0.01 −0.44± 0.02 4.9± 0.07 0.67 p < 0.001 for all parameters except for ap = 0.44, and bp = 0.84. We built linear models for dry weight biomass of understory vegetation (kgha−1) in a two level selection of the predictors from stand, weather and location variables. First, we selected the predictors into linear models by using R package “Mass” and its stepwise model selection by exact Akaike’s informa- tion criterion (AIC; Venables and Ripley, 2002). Second, we refined the model by using “relaimpo” R package estimating usefulness (Grömping, 2006), or relative importance for each of the predictors in the model, and by selecting only predic- tors with relative importance ≥ 0.1. The general form of the models was yi = a+ b1x1+ . . .+ bnxn+ ε, (C1) where yi is the understory dry weight biomass (kgha−1), x1 . . . xn are the predictors, a, b1 . . . bn are parameters of the ith understory functional type (Table C1 in Appendix C), and ε is the residual error. Statistics of the models are shown in Table C1 in Appendix C. Scatter plots between the measured coverage derived biomass and modeled dry weight biomass (kgha−1) of the functional types of ground vegetation for the forests in their observed state close to the estimated equilib- rium are shown on Fig. S9. Biogeosciences, 13, 4439–4459, 2016 www.biogeosciences.net/13/4439/2016/ B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility 4457 The Supplement related to this article is available online at doi:10.5194/bg-13-4439-2016-supplement. Acknowledgements. We thank the Finnish Ministry of Environ- ment and the Finnish Ministry of Agriculture and Forestry for funding this work through the Metla project 7509 “Improving soil carbon estimation of greenhouse gas inventory”, and Academy of Finland for funding the mobility projects 276300 and 276602. We would like to thank the editor and the reviewers for their valuable comments improving the manuscript. Edited by: A. V. Eliseev Reviewed by: three anonymous referees References Adair, E. C., Parton, W. J., Del Grosso, S. J., Silver, W. L., Harmon, M. E., Hall, S. A., Burke, I. C., and Hart, S. C.: Simple three-pool model accurately describes patterns of long-term litter decompo- sition in diverse climates, Global Change Biol., 14, 2636–2660, 2008. Ågren, G. I., Bosatta, E., and Magill, A. H.: Combining theory and experiment to understand effects of inorganic nitrogen on litter decomposition, Oecologia, 128, 94–98, 2001. Ågren, G. I. and Hyvönen, R.: Changes in carbon stores in Swedish forest soils due to increased biomass harvest and increased tem- peratures analysed with a semi-empirical model, Forest Ecol. Manage., 174, 25–37, 2003. Ågren, G., Hyvönen, R., and Nilsson, T.: Are Swedish forest soils sinks or sources for CO2—model analyses based on forest inven- tory data, Biogeochemistry, 82, 217–227, 2007. Amundson, R.: The carbon budget in soils, Annu. Rev. Earth Planet. Sci., 29, 535–562, 2001. Averill, C., Turner, B. L., and Finzi, A. C.: Mycorrhiza-mediated competition between plants and decomposers drives soil carbon storage, Nature, 505, 543–545, 2014. Berthrong, S. T., Jobbágy, E. G., and Jackson, R. B.: A global meta- analysis of soil exchangeable cations, pH, carbon, and nitrogen with afforestation, Ecol. Appl., 19, 2228–2241, 2009. Boden, T. A., Marland, G., and Andres, R. J.: Global, re- gional, and national fossil-fuel CO2 emissions, Carbon Diox- ide Information Analysis Center, Oak Ridge National Labo- ratory, US Department of Energy, Oak Ridge, Tenn., USA, doi:10.3334/CDIAC/00001_V2010, 2010. Castellano, M. J., Mueller, K. E., Olk, D. C., Sawyer, J. E., and Six, J.: Integrating Plant Litter Quality, Soil Organic Matter Stabiliza- tion and the Carbon Saturation Concept, Glob. Change Biol., 21, 3200–3209, doi:10.1111/gcb.12982, 2015. Clarholm, M. and Skyllberg, U.: Translocation of metals by trees and fungi regulates pH, soil organic matter turnover and nitrogen availability in acidic forest soils, Soil Biol. Biochem., 63, 142– 153, 2013. Clemente, J. S., Simpson, A. J., and Simpson, M. J.: Association of specific organic matter compounds in size fractions of soils un- der different environmental controls, Org. Geochem., 42, 1169– 1180, 2011. Cools, N., Vesterdal, L., De Vos, B., Vanguelova, E., and Hansen, K.: Tree species is the major factor explaining C : N ratios in Eu- ropean forest soils, Forest Ecol. Manage., 311, 3–16, 2014. Cotrufo, M. F., Wallenstein, M. D., Boot, C. M., Denef, K., and Paul, E.: The Microbial Efficiency-Matrix Stabilization (MEMS) framework integrates plant litter decomposition with soil organic matter stabilization: do labile plant inputs form stable soil or- ganic matter?, Glob. Change Biol., 19, 988–995, 2013. Deluca, T. H. and Boisvenue, C.: Boreal forest soil carbon: distri- bution, function and modelling, Forestry, 85, 161–184, 2012. Dungait, J. A. J., Hopkins, D. W., Gregory, A. S., and Whitmore, A. P.: Soil organic matter turnover is governed by accessibility not recalcitrance, Glob. Change Biol., 18, 1781–1796, 2012. Falloon, P., Jones, C. D., Ades, M., and Paul, K.: Direct soil mois- ture controls of future global soil carbon changes: An impor- tant source of uncertainty, Global Biogeochem. Cy., 25, GB3010, doi:10.1029/2010GB003938, 2011. Fan, Z., Neff, J. C., Harden, J. W., and Wickland, K. P.: Bo- real soil carbon dynamics under a changing climate: A model inversion approach, J. Geophys. Res.-Biogeo., 113, G04016, doi:10.1029/2008JG000723, 2008. Fernández-Martínez, M., Vicca, S., Janssens, I. A., Sardans, J., Luyssaert, S., Campioli, M., Chapin III, F. S., Ciais, P., Malhi, Y., Obersteiner, M., Papale, D., Piao, S. L., Reichstein, M., Roda, F., and Penuelas, J.: Nutrient availability as the key regulator of global forest carbon balance, Nature Climate Change, 4, 471– 476, 2014. Franklin, O., Högberg, P., Ekblad, A., and Ågren, G. I.: Pine for- est floor carbon accumulation in response to N and PK addi- tions: bomb 14C modelling and respiration studies, Ecosystems, 6, 644–658, 2003. Grömping, U.: Relative importance for linear regression in R: the package relaimpo, J. Stat. Softw., 17, 1–27, 2006. Hagglund, B. and Lundmark, J.: Site index estimation by means of site properties, Scots pine and Norway spruce in Sweden, Stud. For. Suec., 138, 38, 1977. Härkönen, S., Pulkkinen, M., Duursma, R., and Mäkelä, A.: Esti- mating annual GPP, NPP and stem growth in Finland using sum- mary models, For. Ecol. Manage., 259, 524–533, 2010. Husson, O.: Redox potential (Eh) and pH as drivers of soil/plant/microorganism systems: a transdisciplinary overview pointing to integrative opportunities for agronomy, Plant Soil, 362, 389–417, 2013. Jobbágy, E. G. and Jackson, R. B.: The vertical distribution of soil organic carbon and its relation to climate and vegetation, Ecol. Appl., 10, 423–436, 2000. Kammer, A. and Hagedorn, F.: Mineralisation, leaching and sta- bilisation of 13C-labelled leaf and twig litter in a beech for- est soil, Biogeosciences, 8, 2195–2208, doi:10.5194/bg-8-2195- 2011, 2011. Kirschbaum, M. U. F.: Will Changes in Soil Organic Carbon Act as a Positive or Negative Feedback on Global Warming?, Biogeo- chemistry, 48, 21–51, 2000. Kleja, D. B., Svensson, M., Majdi, H., Jansson, P., Langvall, O., Bergkvist, B., Johansson, M., Weslien, P., Truusb, L., and Lin- droth, A.: Pools and fluxes of carbon in three Norway spruce ecosystems along a climatic gradient in Sweden, Biogeochem- istry, 89, 7–25, 2008. www.biogeosciences.net/13/4439/2016/ Biogeosciences, 13, 4439–4459, 2016 4458 B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility Kurz, W. A., Beukema, S. J., and Apps, M. J.: Estimation of root biomass and dynamics for the carbon budget model of the C ana- dian forest sector, Can. J. Forest Res., 26, 1973–1979, 1996. Lehtonen, A. and Heikkinen, J.: Uncertainty of upland soil carbon sink estimate for Finland, Can. J. Forest Res., 45, 1–13, 2015. Lehtonen, A., Sievänen, R., Mäkelä, A., Mäkipää, R., Korhonen, K. T. and Hokkanen, T.: Potential litterfall of Scots pine branches in southern Finland, Ecol. Modell., 180, 305–315, 2004. Lehtonen, A., Linkosalo, T., Peltoniemi, M., Sievänen, R., Mäkipää, R., Tamminen, P., Salemaa, M., Nieminen, T., Tupek, B., Heikki- nen, J., and Komarov, A.: Soil carbon stock estimates in a na- tionwide inventory: evaluating performance of the ROMUL and Yasso07 models, Geosci. Model Dev. Discuss., in review, 2016a. Lehtonen, A., Palviainen, M., Ojanen, P., Kalliokoski, T., Nöjd, P., Kukkola, M., Penttilä, T., Mäkipää, R., Leppälammi-Kujansuu, J. and Helmisaari, H.: Modelling fine root biomass of boreal tree stands using site and stand variables, Forest Ecol. Manag., 359, 361–369, doi:10.1016/j.foreco.2015.06.023, 2016b. Leppälammi-Kujansuu, J., Aro, L., Salemaa, M., Hansson, K., Kleja, D. B., and Helmisaari, H.: Fine root longevity and carbon input into soil from below-and aboveground litter in climatically contrasting forests, Forest Ecol. Manage., 326, 79–90, 2014. Lindén, A., G.: Swedish Geological Survey report, pp. 10, avail- able at: http://resource.sgu.se/produkter/ae/ae118-beskrivning. pdf (last access: 3 August 2016), 2002. Liski, J., Tuomi, M., and Rasinmäki, J.: Yasso07 user-interface manual, Finnish Environment Institute, Helsinki, 2009. Liski, J., Lehtonen, A., Palosuo, T., Peltoniemi, M., Eggers, T., Muukkonen, P., and Mäkipää, R.: Carbon accumulation in Fin- land’s forests 1922–2004–an estimate obtained by combination of forest inventory data with modelling of biomass, litter and soil, Ann. Forest Sci., 63, 687–697, 2006. Majdi, H.: Changes in fine root production and longevity in relation to water and nutrient availability in a Norway spruce stand in northern Sweden, Tree Physiol., 21, 1057–1061, 2001. Malkonen, E.: Annual primary production and nutrient cycle in a birch stand, Commun. Inst. For. Fenn., 91, 1–35, 1977. Manzoni, S., Taylor, P., Richter, A., Porporato, A., and Ågren, G. I.: Environmental and stoichiometric controls on microbial carbon use efficiency in soils, New Phytol., 196, 79–91, 2012. Marklund, L.: , Biomassafunktioner för tall, gran och börk i Sverige, Biomass functions for pine, spruce and birch in Sweden, Swedish University of Agricultural Sciences, Department of Forest Sur- vey, Report 45, ISSN 0348–0496, 1988. McCormack, M. L., Crisfield, E., Raczka, B., Schnekenburger, F., Eissenstat, D. M., and Smithwick, E. A.: Sensitivity of four eco- logical models to adjustments in fine root turnover rate, Ecol. Model., 297, 107–117, 2015. Metherell, A. K.: Century: Soil Organic Matter Model Environ- ment: Technical Documentation: Agroecosystem Version 4.0, Colorado State University, 1993. Merilä, P., Mustajärvi, K., Helmisaari, H., Hilli, S., Lindroos, A., Nieminen, T. M., Nöjd, P., Rautio, P., Salemaa, M., and Ukon- maanaho, L.: Above-and below-ground N stocks in coniferous boreal forests in Finland: Implications for sustainability of more intensive biomass utilization, For. Ecol. Manage., 311, 17–28, 2014. Muukkonen, P. and Lehtonen, A.: Needle and branch biomass turnover rates of Norway spruce (Picea abies), Can. J. Forest Res., 34, 2517–2527, 2004. Nalder, I. A. and Wein, R. W.: A model for the investigation of long-term carbon dynamics in boreal forests of western Canada: I. Model development and validation, Ecol. Model., 192, 37–66, 2006. Olsson, M. T., Erlandsson, M., Lundin, L., Nilsson, T., Nilsson, Å., and Stendahl, J.: Organic carbon stocks in Swedish Podzol soils in relation to soil hydrology and other site characteristics, Silva Fennica, 43, 209–222, 2009. Ortiz, C., Karltun, E., Stendahl, J., Gärdenäs, A. I., and Ågren, G. I.: Modelling soil carbon development in Swedish coniferous for- est soils—An uncertainty analysis of parameters and model esti- mates using the GLUE method, Ecol. Model., 222, 3020–3032, 2011. Ortiz, C. A., Liski, J., Gärdenäs, A. I., Lehtonen, A., Lundblad, M., Stendahl, J., Ågren, G. I., and Karltun, E.: Soil organic carbon stock changes in Swedish forest soils – A comparison of uncer- tainties and their sources through a national inventory and two simulation models, Ecol. Model., 251, 221–231, 2013. Orwin, K. H., Kirschbaum, M. U., St John, M. G., and Dickie, I. A.: Organic nutrient uptake by mycorrhizal fungi enhances ecosys- tem carbon storage: a model-based assessment, Ecol. Lett., 14, 493–502, 2011. Palosuo, T., Foereid, B., Svensson, M., Shurpali, N., Lehtonen, A., Herbst, M., Linkosalo, T., Ortiz, C., Rampazzo Todorovic, G., Marcinkonis, S., Li, C., and Jandl, R.: A multi-model comparison of soil carbon assessment of a coniferous forest stand, Environ. Modell. Softw., 35, 38–49, 2012. Parton, W. J., Schimel, D. S., Cole, C. V., and Ojima, D. S.: Analysis of Factors Controlling Soil Organic Matter Levels in Great Plains Grasslands, 51, 1173–1179, 1987. Parton, W., Ojima, D., and Schimel, D.: Environmental change in grasslands: Assessment using models, Climate Change, 28, 111– 141, 1994. Parton, W.J., McKeown, R., Kirchner, V., and Ojima,D.: CEN- TURY Users’ Manual, Natural Resource Ecology Laboratory, Colorado State University, Ft. Collins., 1992. Peltoniemi, M., Pulkkinen, M., Aurela, M., Pumpanen, J., Kolari, P. and Mäkelä, A.: A semi-empirical model of boreal-forest gross primary production, evapotranspiration, and soil water– calibration and sensitivity analysis, Boreal Environ. Res., 20, 2015. Petersson, H. and Ståhl, G.: Functions for below-ground biomass of Pinus sylvestris, Picea abies, Betula pendula and Betula pubescens in Sweden, Scand., J. For. Res., 21, 84–93, 2006. Raich, J. W., Parton, W. J., Russell, A. E., Sanford Jr, R. L., and Vi- tousek, P. M.: Analysis of factors regulating ecosystemdevelop- ment on Mauna Loa using the Century model, Biogeochemistry, 51, 161–191, 2000. Rantakari, M., Lehtonen, A., Linkosalo, T., Tuomi, M., Tammi- nen, P., Heikkinen, J., Liski, J., Mäkipää, R., Ilvesniemi, H., and Sievänen, R.: The Yasso07 soil carbon model–Testing against re- peated soil carbon inventory, Forest Ecol. Manage., 286, 137– 147, 2012. Rapalee, G., Trumbore, S. E., Davidson, E. A., Harden, J. W., and Veldhuis, H.: Soil Carbon stocks and their rates of accumulation Biogeosciences, 13, 4439–4459, 2016 www.biogeosciences.net/13/4439/2016/ B. Tˇupek et al.: Underestimation of modeled SOC stocks linked to soil fertility 4459 and loss in a boreal forest landscape, Global Biogeochem. Cy., 12, 687–701, 1998. Repola, J.: Biomass equations for birch in Finland, Silva Fenn., 42, 605–624, 2008. Rawls, W. J., Pachepsky, Y. A., Ritchie, J. C., Sobecki, T. M., and Bloodworth, H.: Effect of soil organic carbon on soil water re- tention, Geoderma, 116, 61–76, 2003. Rolff, C. and Ågren, G. I.: Predicting effects of different harvesting intensities with a model of nitrogen limited forest growth, Ecol. Model., 118, 193–211, 1999. R Core Team R: A language and environment for statistical com- puting. R Foundation for Statistical Computing, Vienna, Aus- tria, available at: http://www.R-project.org/ (last access: 3 Au- gust 2016), 2014. Scharlemann, J. P., Tanner, E. V., Hiederer, R., and Kapos, V.: Global soil carbon: understanding and managing the largest ter- restrial carbon pool, Carbon Management, 5, 81–91, 2014. Schlesinger, W. H.: Carbon Sequestration in Soils, Science, 284, 2095–2095, 1999. Schmidt, M. W. I., Torn, M. S., Abiven, S., Dittmar, T., Guggen- berger, G., Janssens, I. A., Kleber, M., Kogel-Knabner, I., Lehmann, J., Manning, D. A. C., Nannipieri, P., Rasse, D. P., Weiner, S., and Trumbore, S. E.: Persistence of soil organic mat- ter as an ecosystem property, Nature, 478, 49–56, 2011. Six, J., Conant, R. T., Paul, E. A., and Paustian, K.: Stabilization mechanisms of soil organic matter: Implications for C-saturation of soils, Plant Soil, 241, 155–176, 2002. SLU: Markinfo, available at: http://www-markinfo.slu.se/eng/ index.html (last access: 3 August 2016), 2011. Smith, P.: An overview of the permanence of soil organic carbon stocks: influence of direct human-induced, indirect and natural effects, Eur. J. Soil Sci., 56, 673–680, 2005. Sollins, P., Homann, P., and Caldwell, B. A.: Stabilization and destabilization of soil organic matter: mechanisms and controls, Geoderma, 74, 65–105, 1996. Statistics Finland: Greenhouse gas emissions in Finland 1990– 2011, in: National Inventory Report to the UNFCCC Secretariat, Ministry of the Environment, Helsinki, Finland, 285–286, 2013. Stendahl, J., Johansson, M., Eriksson, E., Nilsson, Å., and Lang- vall, O.: Soil organic carbon in Swedish spruce and pine forests—differences in stock levels and regional patterns, Silva Fenn., 44, 5–21, 2010. Swedish Statistical Yearbook of Forestry: Official Statistics of Swe- den, 370 pp., Skogsstyrelsen, 2014. Svensson, M., Jansson, P., and Kleja, D. B.: Modelling soil C se- questration in spruce forest ecosystems along a Swedish transect based on current conditions, Biogeochemistry, 89, 95–119, 2008. Therneau, T. M. and Atkinson, E. J.: An introduction to recursive partitioning using the RPART routines, 1997. Throop, H. L., Holland, E. A., Parton, W. J., Ojima, D. S., and Keough, C. A.: Effects of nitrogen deposition and insect her- bivory on patterns of ecosystem level carbon and nitrogen dy- namics: results from the CENTURY model, Glob. Change Biol., 10, 1092–1105, 2004. Todd-Brown, K. E. O., Randerson, J. T., Hopkins, F., Arora, V., Ha- jima, T., Jones, C., Shevliakova, E., Tjiputra, J., Volodin, E., Wu, T., Zhang, Q., and Allison, S. D.: Changes in soil organic carbon storage predicted by Earth system models during the 21st cen- tury, Biogeosciences, 11, 2341–2356, doi:10.5194/bg-11-2341- 2014, 2014. Todd-Brown, K. E. O., Randerson, J. T., Post, W. M., Hoffman, F. M., Tarnocai, C., Schuur, E. A. G., and Allison, S. D.: Causes of variation in soil carbon simulations from CMIP5 Earth system models and comparison with observations, Biogeosciences, 10, 1717–1736, doi:10.5194/bg-10-1717-2013, 2013. Torn, M. S., Trumbore, S. E., Chadwick, O. A., Vitousek, P. M., and Hendricks, D. M.: Mineral control of soil organic carbon storage and turnover, Nature, 389, 170–173, 1997. Tuomi, M., Rasinmäki, J., Repo, A., Vanhala, P., and Liski, J.: Soil carbon model Yasso07 graphical user interface, Environ. Modell. Softw., 26, 1358–1362, 2011. Tuomi, M., Thum, T., Järvinen, H., Fronzek, S., Berg, B., Harmon, M., Trofymow, J. A., Sevanto, S., and Liski, J.: Leaf litter de- composition – Estimates of global variability based on Yasso07 model, Ecol. Model., 220, 3362–3371, 2009. Viro, P.: Investigations on forest litter, Valtioneuvoston Kirjap, 1955. Wiesmeier, M., Hübner, R., Spörlein, P., Geuß, U., Hangen, E., Reischl, A., Schilling, B., von Lützow, M., and Kögel-Knabner, I.: Carbon sequestration potential of soils in southeast Germany derived from stable soil organic carbon saturation, Glob. Change Biol., 20, 653–665, 2014. Xia, M., Talhelm, A. F. and Pregitzer, K. S.: Fine roots are the dominant source of recalcitrant plant litter in sugar maple- dominated northern hardwood forests, New Phytol., 208, 715– 726, doi:10.1111/nph.13494, 2015. Yurova, A. Y., Volodin, E. M., Agren, G. I., Chertov, O. G., and Komarov, A. S.: Effects of variations in simulated changes in soil carbon contents and dynamics on future climate projections, Glob. Change Biol., 16, 823–835, 2010. Venables, W. N. and Ripley, B. D.: Modern applied statistics with S-PLUS, Springer Science & Business Media, 2013. Zandonadi, D. B., Santos, M. P., Busato, J. G., Peres, L. E. P. ,and Façanha, A. R.: Plant physiology as affected by humified organic matter, Theoretical and Experimental Plant Physiology, 25, 13– 25, 2013. www.biogeosciences.net/13/4439/2016/ Biogeosciences, 13, 4439–4459, 2016