The importance of tree species and soil taxonomy to modeling forest soil carbon stocks in Canada

The importance of tree species and soil taxonomy to modeling forest soil carbon stocks in Canada

Geoderma Regional 4 (2015) 114–125 Contents lists available at ScienceDirect Geoderma Regional journal homepage: The...

2MB Sizes 0 Downloads 126 Views

Geoderma Regional 4 (2015) 114–125

Contents lists available at ScienceDirect

Geoderma Regional journal homepage:

The importance of tree species and soil taxonomy to modeling forest soil carbon stocks in Canada C.H. Shaw a,⁎, K.A. Bona b, W.A. Kurz c, J.W. Fyles d a

Natural Resources Canada, Canadian Forest Service, 5320 122nd Street, Edmonton, AB T6H 3S5, Canada Macdonald Campus of McGill University, 21 111 Lakeshore Road, Ste-Anne-de-Bellevue, QC H9X 3V9, Canada Natural Resources Canada, Canadian Forest Service, 506 West Burnside Road, Victoria, BC V8Z 1M5, Canada d Macdonald Campus of McGill University, 21 111 Lakeshore Road, Ste-Anne-de-Bellevue, QC H9X 3V9, Canada b c

a r t i c l e

i n f o

Article history: Received 11 July 2014 Received in revised form 22 December 2014 Accepted 9 January 2015 Available online 13 January 2015 Keywords: Forest soil carbon Model Initialization Soil taxonomy Tree species Redundancy analysis Mollisols Alfisols Entisols Inceptisols Spodosols Gelisols Albeluvisols Solonetz Luvisols Podzols Cambisols Fluvisols Regosols Gleysols Planosols Cryosols

a b s t r a c t Accurate initialization of soil and dead organic matter carbon (C) stocks in forest ecosystem models is challenging but critical to forest C estimation, assessing current and future responses to climate change, and evaluation of management options for climate change mitigation strategies. We identified opportunities to improve the accuracy of soil C estimates from the Carbon Budget of the Canadian Forest Sector (CBM-CFS3) — a model of forest C dynamics used to support greenhouse gas emission reporting. Accuracy of soil C stocks estimated by models is very dependent on the initialization process. Here, we used redundancy analysis (RDA) and ordinations in an exploratory analysis to compare the variance structures of soil C estimates determined by model variables used in the initialization process, in two different soil C datasets; one derived from the model, the other obtained from 2391 ground plots. We also used the ground plot data to determine if soil taxonomy (information currently not used in the CBM-CFS3) could be used to explain variation in addition to that already accounted for by variables in the model. Total variance of the plot C dataset was about twice as large as the variance of the model C dataset confirming that currently the model does not represent all factors that control variation in soil C stocks. Soil C stocks in the mineral soil were highly correlated with C stocks in soil organic horizons in the model dataset but not in the plot dataset, suggesting that the variables included in our assessment controlling C stocks in the mineral soil horizons are different than in the organic soil horizons. Tree productivity (maximum yield curve volume per hectare) explained a much larger proportion of the total variation in the model dataset than in the plot dataset, whereas the leading tree species explained more variation in the plot dataset than in the model, suggesting that accuracy of initialization of soil C stocks could be improved by including leading tree species to stratify soil C modeling parameters. Leading species that are in greatest need of improved representation were identified by ordination. The results from the RDA showed that soil taxonomy explained 4 (order) to 13% (subgroup) of plot soil C variance, in addition to that explained by variables currently used in the model that determine initial soil C stocks. Soil taxonomy and leading species can compensate for one another to explain variance in soil C stocks. Our results suggest the potential of using the combination of leading tree species and soil taxonomy to improve soil C stocks initialized by forest C models, but this remains to be tested. Crown Copyright © 2015 Published by Elsevier B.V. All rights reserved.

1. Introduction Soil carbon (C) is indisputably important as a C stock and as a contributor to the global C budget (Schlesinger and Andrews, 2000; Gottschalk et al., 2012). Forest ecosystems cover approximately 3.8 billion ha globally (Pan et al., 2011) and their soil represents a significant component of forest C budgets (Lal, 2005). The contribution of

⁎ Corresponding author at: Northern Forestry Centre, 5320 122 Street, Edmonton, AB T6H 3S5, Canada. E-mail addresses: [email protected] (C.H. Shaw), [email protected] (K.A. Bona), [email protected] (W.A. Kurz), [email protected] (J.W. Fyles). 2352-0094/Crown Copyright © 2015 Published by Elsevier B.V. All rights reserved.

soil to forest C budgets is particularly high for boreal forest ecosystems (Lal, 2005; Kurz et al., 2013) — the dominant forest biome in Canada (Brandt, 2009). Soil C models in various forms (Peltoniemi et al., 2007) are integrated into forest ecosystem models like the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3; Kurz and Apps, 1999; Kurz et al., 2009), FORCARB (Harmon and Marks, 2002), BIOME-BGC (Running and Hunt, 1993), CO2FIX (Schelhaas et al., 2004) and EFIMOD 2 (Komarov et al., 2003), and these models have been designed or adapted to estimate forest C budgets or to contribute to the understanding of forest C dynamics. Important applications of these models include meeting national and international reporting requirements (Stinson et al., 2011), assessing current and predicting future responses to climate change (Boisvenue and Running, 2010;

C.H. Shaw et al. / Geoderma Regional 4 (2015) 114–125

Shanin et al., 2010), and/or evaluating forest management activities and their contributions to climate change mitigation objectives (Smyth et al., 2014). However, the uncertainty of soil C stock and stock change estimation by these models is often high (Palosuo et al., 2012; Todd-Brown et al., 2013; Shaw et al., 2014). This study compares the relative importance of factors influencing soil C stock estimation from plot data and a forest ecosystem C budget model (CBM-CFS3) with the goal of improving the accuracy of initial estimates of soil C stocks derived from the model spin-up procedures (Kurz et al., 2009). Soil C models vary in structure, number and type of pools, and the approach used to initialize the pools (McGill, 1996; Peltoniemi et al., 2007). Here we focus on variables that affect soil C stocks in the initialization process that is used to estimate stocks in modeled pools prior to running simulations or scenario analyses. Initialization of soil pools as part of the modeling process, rather than using inventory-based estimates, is necessary because the data required to estimate initial soil C stocks are typically not measured in forest inventories, soil data that is collected is often insufficient to represent areas that are large or that have diverse soils, and because model pools are often defined in such a way that their measurement in field samples is difficult or time consuming (Zimmerman et al., 2007). Most approaches to initialization depend on the validity of the steady-state assumption; that soil C inputs (primarily from plants) are equivalent to soil C outputs (primarily from heterotrophic respiration). One of the most common approaches to initialization is to run the model to a steady state in which the slower pools cease to change (Hashimoto et al., 2011). The approach to initialization of soil C in models is important because it affects the outcome of contemporary and predictive simulations. It influences soil C stock and stock change estimation (Foereid et al., 2012; Xu et al., 2011; Palosuo et al., 2012), calibration of model pools with slow turnover rates (Peltoniemi et al., 2007; Wutzler and Reichstein, 2007), and the range in soil C stocks or sequestration rates (Foereid et al., 2012; Wutzler and Reichstein, 2007). The CBM-CFS3 differs from many other models in that it includes pools for soil organic horizon C as well as soil mineral horizon C and it simulates natural (e.g., fire, insects) and anthropogenic (e.g., harvest) disturbance effects on biomass, dead organic matter (DOM: standing and downed dead wood) pools, and soil organic horizon C; in turn these can affect mineral soil C. Including disturbance effects has been an important and unique feature of the CBM-CFS3 since its inception (Kurz et al., 1992) because they are important for estimating greenhouse gas (GHG) emissions from Canada's forested area. In the CBM-


CFS3 all DOM and soil pools are initialized (starting at zero) through a spin-up process in which cycles of stand development, initiated by stand-replacing wildfire (or another user-defined stand-replacing disturbance such as insects or harvesting), are repeated for as many times as required to reach a quasi-steady state defined as the point where the sum of the slow pools at the end of two successive cycles differs by less than 0.1% (Fig. 1). Thereafter, the model simulates the last stand-replacing disturbance and grows the stand to the age in the inventory. Thus, depending on the time since last disturbance, the pools can be far from an equilibrium state (e.g., recently disturbed stands) or close to an equilibrium state (e.g., in very old stands which are rare in boreal forest ecosystems). Estimates of initial soil C stocks by the CBM-CFS3 are typically lower than estimates from the most common equilibrium approach, because it accounts for losses of soil C due to disturbances (Fig. 1). Despite these differences, the CBM-CFS3 sometimes over- or underestimates initial soil C stocks, and estimates a much narrower range (lower total variance) in soil C stocks compared with ground plot data (higher total variance) (Shaw et al., 2014). Thus, the main objectives of this exploratory study are to determine if the CBM-CFS3 makes full and appropriate use of variables used in the initialization process to express variation in soil C stocks, and if soil taxonomy can account for additional variation in predicted soil C stocks. We included soil taxonomy because currently it is not used in the model, but may be useful for improving modeled soil C estimates because it has been linked to differences in forest soil C stocks (Shaw et al., 2008) and pedological factors and processes important to soil C stabilization and accumulation (Shaw et al., 2008; Schaetzl and Anderson, 2005; Jandl et al., 2007). 2. Methods We compared the variance structure of soil C stocks estimated from 2391 plots located across Canada to the variance structure of soil C estimated by the CBM-CFS3 in national-scale simulations, as explained by variables used in the model that determine initial soil C stocks. In the CBM-CFS3 variables that influence initial soil C stock estimates (e.g., decomposition, disturbance impacts, and biomass estimation) are defined at different spatial scales (Table 1; Kurz et al., 2009). We used redundancy analysis (RDA) (ter Braak, 1987; ter Braak and Prentice, 1988) to partition total variance of the model-based and the plot-based estimates of mineral and organic horizon soil C as explained by the model's spatial structure and associated variables. Here we describe the two datasets (plot- and model-derived) and the RDA methods used for their analysis. 2.1. Model-derived dataset for the RDA

Fig. 1. Comparison of the accumulated sum of slow soil C pools during initialization of the CBM-CFS3 with disturbance from a stand-replacing fire using a 200-year historical fire return interval (solid line) and with no disturbances (dashed line).

Canada's National Forest Carbon Monitoring, Accounting and Reporting System (NFCMARS; Stinson et al., 2011) combines national scale forest inventory, growth and yield curve, disturbance and landuse change activity data with the CBM-CFS3 to estimate C stocks and fluxes for Canada's managed forest area. In this application forest stands (ranging in size from 1 to 370,000 ha) are represented by aspatial inventory data located within polygons of varying sizes depending on the source of the inventory data (Kurz et al., 2009). Nearly three million records representing stands in the managed forest area (Fig. 2) are used as input to the CBM-CFS3 when it is used in the NFCMARS. The model generates, among other information, estimates of soil C stocks. Output from the CBM-CFS3 model runs for 1990 from the 2010 National Inventory Report (2010NIR; Environment Canada, 2010) was compressed to a manageable 21,994 records by averaging results for stands with similar classifier set values. Classifier sets vary by jurisdiction and include information such as national or regional ecological classification, leading species or forest type, management units, management history, or productivity class. Classifier sets are used in the CBM-CFS3 to describe stands and to assign the appropriate yield curves. Estimates for organic


C.H. Shaw et al. / Geoderma Regional 4 (2015) 114–125

soil horizon (ORGSOIL) and mineral soil horizon (MINSOIL) C density (Mg C ha−1) as well as the explanatory variables and scales (Table 1) were extracted for each record in the compressed dataset. Modeled C stocks for ORGSOIL were calculated as the sum of C stocks estimated by the CBM-CFS3 for the Aboveground Slow and Aboveground Very Fast pools and for MINSOIL as the sum of C stocks estimated for the Belowground Slow and Belowground Very Fast pools. Soil pools in the CBM-CFS3 are described using kinetic labels reflecting their relative rates of decomposition and types of organic matter inputs, in combination with their location relative to the mineral soil surface (above- or belowground). The Aboveground Very Fast pools include dead fine roots and foliar litter and the Belowground Very Fast pool includes dead fine roots in the mineral soil. Above- and Belowground Slow pools represent humified organic matter in the soil organic and mineral layers, respectively (see Kurz et al., 2009 for detailed pool descriptions).

order (e.g., Histosols) and organic soils in the Cryosolic order; SCWG, 1998; see Table 2 for World Reference Base and US Soil Taxonomy taxonomic equivalents) because the CBM-CFS3 is primarily designed for upland forest types. Data were compiled from seven sources (Siltanen et al., 1997; (RDB) Resource Data Branch, 2002; CEMA, 2004; ELCG, 2005; Shaw et al., 2005; (CanSIS) Canadian Soil Information Service, 2007; NFI, 2008, 2010, 2011) comprising previously compiled and new data. Care was taken to avoid duplication of records. Carbon stocks (Mg ha−1) for ORGSOIL (L, F, H and O horizons; SCWG, 1998) and MINSOIL (mineral soil C to 100 cm depth) were estimated for the pedons located at the plots. Organic C stocks for each horizon were calculated using Eq. (1),

2.2. Plot-derived dataset for the RDA

 horizon thickness ð cmÞ: ð1Þ The C stocks for organic soil horizons were summed to provide an estimate for the combined organic horizons (ORGSOIL) and the estimates for the mineral soil horizons were summed to provide an estimate for the combined mineral soil horizons (MINSOIL). All profiles

A dataset (n = 2391), similar in structure to the model-derived dataset, was compiled using data from ground plots (Fig. 2). We used plot data for upland mineral soils only (excluding soils in the Organic

  −1 Organic C stocks Mg ha ¼ organic carbon ð%Þ   −3  bulk density g cm

Table 1 Description of parameters and variables affecting soil C stock estimates associated with different scales used in the CBM-CFS3, and description of additional soil variables linked to taxonomy in the Canadian system of soil classification (SCWG, 1998) that are currently not used in the model. We refer to groups of parameters and variables associated with the ecozone scale as E-variables, with the province or territory scale as P-variables, with the reconciliation unit as R-variables and with the stand scale as S-variables. Variables and parameters used in the CBM-CFS3 CBM-CFS3 scale

Parameter or variable



Base decay rate

Scope of the CBM-CFS3 for the 2010NIR is the managed forest area of Canada (Fig. 2) One base decay rate for each of the four CBM-CFS3 soil pools. Default rates are AGSlow = 0.015; AGVFast = 0.355; BGSlow = 0.0033; BGVFast = 0.5. As defined by Kull et al. (2011) The CBM-CFS3 default fire return intervals (mean stand age) for initialization (Fig. 1) were extracted from the 2010NIR database files for the model dataset. FRIs were assigned to plots for the plot dataset based on the ecozone in which they are located. Provinces or territories of Canada Merchantable hardwood and softwood proportions for stumps and tops (Appendix 2 in Kull et al., 2011). Units defined by the intersection of provinces or territories (P) with terrestrial ecozones (E) Proportions of the CBM-CFS3 AGSlow pool retained after wildfire were extracted from disturbance matrices (Kurz et al., 2009) in the 2010NIR database files for the model dataset and assigned to plots based on the reconciliation unit in which they are located. The AGSlow pool in the CBM-CFS3 comprises most of the soil C found in the organic soil horizons. Leading species (GENUS; genus only; GSPECIES; genus and species) were extracted from the 2010NIR databases for the model dataset. They were based on species with the largest proportion of merchantable volume or biomass for the plot data. Hardwood, softwood or mixedwood stand type. Extracted from the 2010NIR databases for the model dataset. For the plot dataset the composition is hardwood or softwood if the volume or biomass of hardwood or softwood species, respectively, is ≥75% of the total plot volume or biomass. If the composition was not determined to be softwood or hardwood — it was mixedwood. Individual records in the 2010NIR compressed database and individual plots for the plot dataset Maximum volume at or below the age equal to the FRI, or at or below the plot's stand age, whichever is greater, from the plot yield curve for the plot dataset. Maximum volume at or below the age equal to the FRI, or at or below the averaged stand age, whichever is greater, from the averaged yield curve for a record in the compressed 2010NIR database for the model dataset. The age (years) estimated for the stand was extracted from the 2010NIR database for the model dataset and estimated from tree core data for the plot dataset.

Ecozone Fire return interval (years)

Province or territory Merchantability limit Reconciliation unit Proportion of AGSlow retained

Leading species


Stand Maximum volume (m3)

Stand age

Mean annual temperature (°C)

Mean annual temperatures are estimated for the 2010NIR using the model of McKenney et al. (2001) and were extracted from the 2010NIR database files for the model dataset. The methods of McKenney et al. (2001) were used to estimate mean annual temperatures for the plots. Base decay rates are modified by mean annual temperature (Kurz et al., 2009).

Codea BDR E FRI (years) (75 to 300) P MLIM R FFRET (0.28 to 0.95)



S MVOL (m3 ha−1) (model 5 to 2402) (plot 1 to 2554) STANDAGE (years) (model 1 to 600) (plot 11 to 208) MAT (°C) (model −9.9 to 9.1) (plots −12.3 to 10.1)

Soil taxonomic variables not in the CBM-CFS3 Variable



Soil order

Highest level of organization for soil taxonomy in the Canadian system of soil classification (SCWG, 1998); conveys the least taxonomic information Soil taxa organized into strata separated by soil C stocks analyses (Shaw et al., 2008) Lowest level of organization for soil taxonomy in the CSSC (SCWG, 1998); conveys the most taxonomic information


Soil carbon modeling category Great group subgroup a

Range of parameter or variable values in the CBM-CFS3 are shown in parentheses where appropriate.


C.H. Shaw et al. / Geoderma Regional 4 (2015) 114–125


Fig. 2. Distribution of plots used to produce the plot-derived dataset. The gray zone is the managed forest area for reporting forest greenhouse gas emissions and removals in the NIR2010 using the CBM-CFS3. Output from the CBM-CFS3 runs used for the NIR2010 was used to produce the model-derived dataset. The area in the southern portion of central Canada without plot data is comprised of prairie ecozones that are dominated by agricultural land-uses.

had data for horizon thicknesses. We only used profiles if % organic C values were available for most of the horizons. In some cases % organic C analyses were not available for very thin or transitional horizons. In these cases we substituted values using a mean value for similar horizons in the same soil taxonomic class. It was very common that bulk density values were not commonly measured in legacy data. Where bulk density values were missing we substituted values from Table 8 in Shaw et al. (2005) that is a look-up table for bulk density values based on an analysis of a dataset (n = 3736) where soil taxonomy and horizon names were known and bulk density was measured. Mean bulk density values were estimated for horizons, or groups of pedogenically similar horizons, for different taxonomic classes. Bulk density values for thick organic horizons were estimated from a separate database (n = 6989) (Zoltai et al., 2000) based on a similar type of analysis. Carbon stock values summarized by soil order are provided in Table 3. Terrestrial ecozone (ECOZONE), province or territory (PROVINCE), reconciliation unit (RU, the intersection of Province and Ecozone), leading tree species genus (GENUS), leading tree species genus and species (GSPECIES), and an estimate of stand age (STANDAGE) (Table 1) were determined from the plot data. Maximum input volume (MVOL) (see Table 1 for detailed definition) was determined for a subset of plots (n = 785) where yield curves had been assigned to the plots for other research purposes. Mean annual temperature (MAT) (Table 1) was estimated using the methods of McKenney et al. (2001). Based on plot location in a terrestrial ecozone and/or reconciliation units, values for mean fire return interval (FRI) (Appendix 3 in Kull et al., 2011), merchantability limits (MLIM) (Appendix 2 in Kull et al., 2011), and the proportion of the soil organic horizons retained after a wildfire (FFRET; Table 1) (defined in fire disturbance matrices accessible through the model interface) were assigned to plots. Soil order, great group and subgroup (SCWG, 1998) were compiled for each plot and also used to assign each plot to a taxonomically-based stratification of soil carbon modeling categories (SCMCs; Shaw et al., 2008).

2.3. Full and partial redundancy analysis Redundancy analysis (RDA) (ter Braak, 1987; ter Braak and Prentice, 1988) was utilized to explain the variation in ORGSOIL and MINSOIL with model and soil taxonomic explanatory variables. RDA is an extension of a principal component analysis where the response variables are constrained by the linear combinations of explanatory variables, such as in a multiple regression. RDA was also chosen over other ordination techniques because we suspected monotonic responses (versus unimodal) and because RDA has been shown to be an appropriate method for analysis of forest soil data (Odeh et al., 1991). We first conducted RDAs with different combinations of variables used to explain the variance in ORGSOIL and MINSOIL, which for clarity in this paper we will call the full RDAs. We then conducted a set of partial RDAs used to partition the variance explained in ORGSOIL and MINSOIL by groups of variables associated with different spatial scales (Table 1) defined by the model's structure, and their interactions. The full RDAs were conducted using the CANOCO for Windows 4.5 software package (Biometris; Wageningen, Netherlands). Two full RDAs were conducted with the plot- and model-derived estimates of ORGSOIL and MINSOIL as the response variables. Different combinations of explanatory variables were used to explain the variation in soil C from the list of available site and soil plot characteristics: ECOZONE, PROVINCE, RU, stand type (COMP), FFRET, MVOL, MAT, GENUS, and GSPECIES (Table 1). For the plot RDAs, variables based on soil taxonomy (soil order, ORDER; soil carbon modeling category, SCMC; soil subgroup and great group, SGGG) (Table 1) were added to the analysis. The qualitative explanatory variables were converted to dummy variables (binary data), while quantitative variables were centered and standardized. For each combination of explanatory variables used to explain the response variables (ORGSOIL and MINSOIL) a correlation bi-plot was generated. RDA correlation bi-plots (e.g., results Figs. 4 and 5), can be interpreted by observing the angle between variables that reflect their


C.H. Shaw et al. / Geoderma Regional 4 (2015) 114–125

Table 2 Taxonomic correlation at the Canadian order and great group levelsa. Adapted from Table 8 in “The Canadian system of soil classification” (SCWG, 1998). Canadian system (SCWG, 1998)

U.S. soil taxonomy (Soil Survey Staff, 1992)

WRB (FAO) system (Deckers et al., 1998)



Brown Chernozem Dark Brown Chernozem Black Chernozem Dark Gray Chernozem

Luvisolic Gray Brown Luvisol

Aridic Boroll subgroups Typic Boroll subgroups Udic Boroll subgroups Boralfic Boroll subgroups, Albolls Natric great groups, Mollisols & Alfisols Boralfs & Udalfs Hapludalfs or Glossudalfs

Kastanozem, Chernozem, Greyzem, Phaeozem Kastanozem (aridic) Kastanozem (Haplic) Chernozem Greyzem

Gray Luvisol


Podzolic Humic Podzol Ferro-Humic Podzol

Eutric Brunisol

Spodosols, some Inceptisols Cryaquods, Humods Humic Cryorthods, Humic Haplorthods Cryorthods, Haplorthods Inceptisols, some Psamments Cryochrepts, Eutrochrepts, Hapludolls Cryochrepts, Eutrochrepts

Sombric Brunisol

Humbric Dystrochrepts

Dystric Brunisol Regosolic Regosol Humic Regosol Gleysolic Humic Gleysol

Dystrochrepts, Cryochrepts Entisols Entisols Entisols Aqu-suborders Aquolls, Humaquepts

Gleysol Luvic Gleysol

Aquents, Fluvents, Aquepts Argialbolls, Argiaquolls, Aqualfs Gelisols Turbels Orthels


Humo-Ferric Podzol Brunisolic Melanic Brunisol

Cryosolic Turbic Cryosol Static Cryosol a

Solonetz Luvisol Albic Luvisol, Haplic Luvisol Albic Luvisol, Gleyic Luvisol Podzol Humic Podzol Orthic Podzol Orthic Podzol Cambisol Cambisol, Eutric Cambisol Eutric Cambisol, Calcic Cambisol Dystric Cambisol, Umbric Cambisol Dystric Cambisol Fluvisol, Regosol Regosol Fluvisol, Regosol Gleysol, Planosol Mollic, Umbric, Calcic Gleysol Eutric, Dystric Gleysol Planosol Cryosol Cryosol Cryosol

Only the nearest equivalents are indicated.

correlation (90° is analogous to having a correlation coefficient of 0, while 0° or 180° is analogous to having a correlation coefficient of 1 or −1, respectively), with the length of the arrow expressing the strength of the relationship between variables (Legendre and Legendre, 1998). The partial RDAs were conducted using the open source R software and vegan package (R. Development Core Team, 2011). Partial RDAs, as in partial regressions, can be used to explain the variance in Y (ORGSOIL and MINSOIL) with X (one group of explanatory variables) while controlling for the confounding effects of W (a separate group of explanatory variables) on Y. By using different combinations of variables in the X and W matrices we can partition the variance explained Table 3 Carbon stocks (Mg ha−1) by soil order in the Canadian system of soil classification (SCWG, 1998) for organic soil horizons (ORGSOIL) and mineral soil horizons to 100 cm depth (MINSOIL). Soil order

Brunisolic Chernozemic Cryosolic Gleysolic Luvisolic Podzolic Regosolic Solonetzic




Mean Standard deviation n

Mean Standard deviation

918 2 34 216 652 452 101 16

33 16 55 66 31 39 45 26

80 116 194 127 85 124 156 112

26.4 5.6 56.6 67.4 26.6 29.7 60.2 10.2

893 2 34 215 638 452 91 16

by groups of variables, to estimate how much of the variation is explained by one variable without the confounding effects of another variable. For this study, we wanted to examine how the explained variation was partitioned among the different spatial scales and associated variables used in the model, for soil C estimates obtained from the model and from plot data. For both model and plot data we defined four groups of variables (E-variables, P-variables, R-variables, S-variables) associated with the four different scales (Ecozone, Province, Reconciliation Unit, Stand) in the model (Table 1). As in the full RDAs, qualitative data were converted to dummy variables and quantitative data were centered and standardized.

74.4 40.9 124.2 124.3 66.8 113.1 130.6 54.2

3. Results The total variance for MINSOIL and ORGSOIL estimated for the model was about half (2536 (Mg ha− 1)2) of that for the plot data (5959 (Mg ha−1)2). Results from the full RDA showed that scales (E, P, R), and model variables associated with those scales, explained a lower percent of the total variance for plots (26.2%) compared with the model (Table 4). However, this percent is similar to the percent variation explained by variables linked to these scales (24.1%) which supports the underlying premise used in the model of linking variables to scales. The percent of the model's total variance explained by scale (131.4%; Table 4), or by variables associated with scales, exceeded 100% indicating significant interactions in the model between scales or between variables associated with scales. Based on this outcome we re-analyzed the data using partial RDA to attribute variance components to variables associated with scales, and associated with interactions of scales. All of the results describe from here on are based on the partial RDA with the exception of one clear reference back to the full RDA, and the soil taxonomy analysis where we used results from the full RDA. A much larger percent of the total variance was accounted for by model variables in the model dataset (88.2%) than by the same variables in the plot dataset (24.1%) (Table 5). Seventy-eight percent of the total variance accounted for in the plot dataset was explained by variables applied at the scale of the reconciliation units (R-variables; Fig. 3a), another 5% by variables applied at the scale of the stand (S-variables), and another 12% by an interaction of R and S-variables with variables applied at the scale of the ecozone (E-variables) (Fig. 3a). The remaining 5% is distributed in small percentages among the remaining 12 partitions (Fig. 3a; Table 5). Because few significant interactions were identified in the partial RDA from the plot data we can refer back to the full RDA to conclude that of the R-variables (genus, genus and species, forest floor retained after fire, hardwood or softwood), leading tree genus, or genus and species, are the ones that explains most of the variance (94% or 98%) for R-variables (Table 4). Significant percentages of the total variance accounted for in the model dataset (Fig. 3b) were distributed among more partitions when compared with the plot dataset (Fig. 3a). R-variables explained only 10%, and partitions to which they contributed (RES, RESP) an additional 40%, of the total variance accounted for in the model dataset. Of the primary sets of variables (E-, P-, R-, S-variables), stand-scale (S) variables explained the largest percentage of variance accounted for in the model dataset (28%), and the partitions stand-scale variables contribute to (RES, ES, ESP, RESP) explain an additional 60% (Fig. 4b). Stand-scale (S) variables (primarily productivity (MVOL)), account for a small amount of variance in the plot data (65.4 (C Mg ha− 1)2) compared with the model (622.9 (C Mg ha−1)2), while the converse was true for reconciliation unit (R) variables (primarily leading tree species), where more variance was explained by R-variables in the plot data (952.2 (C Mg ha−1)2) than for the modeled data (230.2 (C Mg ha− 1)2) (Table 5). Thus, Rvariables (primarily leading species) explained most of the variance in the plot data whereas S-variables (primarily MVOL) explained most of the variance in the model data, followed by R-variables. We compared the variance structure of the R-variables for the model- and plot-based datasets using ordination diagrams (Fig. 4)

C.H. Shaw et al. / Geoderma Regional 4 (2015) 114–125 Table 4 Percent variance explained (R2) by spatial scales, variables associated with spatial scales and soil taxonomic variable based on the full RDA for soil C estimates from plot data and the model.

Spatial scale Ecozone (E) Province or Territory (P) Reconciliation Unit (R) Total Variables associated with spatial scales Ecozone Fire return interval (FRI) Province or territory Merchantability limit (MLIM) Reconciliation unit with GENUS Leading tree genus only (GENUS) Proportion of AGSlow pool retained after wildfire (FFRET) Hardwood, softwood or mixedwood stand type (COMP) GENUS, FFRET, COMP combined Reconciliation unit with GSPECIES Leading tree genus and species (GSPECIES) Proportion of AGSlow pool retained after wildfire (FFRET) Hardwood, softwood or mixedwood stand type (COMP) GSPECIES, FFRET, COMP combined Stand (S) Maximum volume of input yield curve (MVOL) Age of the stand (STANDAGE) Mean annual temperature (MAT) MVOL, STANDAGE, MAT combined

Plot R2a

Model R2b

8.9 4.3 13 26.2

51.1 21.1 59.2 131.4





10.9 0.6 0.4 11.6

33.8 1.8 0.5 33.9

15.7 0.6 0.4 16

44.7 1.8 0.5 45.1

1 1.1 3.4 4.2

68.7 9.2 10.1 78

Soil taxonomic variables Sub group-great group (SGGG) Soil carbon modeling category (SCMC) Soil order (ORDER)

18.8 9.7 7.5


15 17 24.3 20 21.5 29.7

a n = 2321 for plot data except for the stand scale where n = 477 because of missing data mainly for MVOL; p-values for plot R2 are 0.002 or 0.006; plot total variance = 5956.5 (Mg ha−1)2 where n = 2321, and 5095.9 (Mg ha−1)2 where n = 477. b n = 21,994 for the model data; p-values for model R2 are 0.002 or 0.005; model total variance = 2536.1 (Mg ha−1)2.

from the partial RDA. The 17 leading species shown in the ordination graphs are representative of the leading species in Canada's forests. The suite of species is similar in composition to the 17 most common leading species in the inventory data from the 2010NIR and in national forest inventory ground plot data that were located based on a sampling design representative for Canada's forested area (Shaw et al., 2014). Because MINSOIL and ORGSOIL are strongly correlated in the model (co-linear in Fig. 4b), but are not in the plot data (Fig. 4a), some relationships between ORGSOIL and MINSOIL with R-variables seen in the plot diagram (Fig. 4a) simply cannot occur in the model diagram (Fig. 4b). For example, in the plot diagram some variables are strongly related to ORGSOIL (e.g., black spruce [Picea mariana (Mill.) BSP], white spruce [Picea glauca (Moench) Voss], eastern white cedar [Thuja occidentalis L.], SW, FFRET) but are not strongly related (orthogonal) to MINSOIL. In the model diagram any variables strongly related to ORGSOIL (e.g., Amabalis fir [Abies amabilis [Dougl. (ex Loud.) Dougl. ex J. Forbes], western hemlock [Tsuga heterophylla (Raf.) Sarg.], Douglas-fir [Pseudosuga menzziesii (Mirb.) [Franco]), must also be strongly related to MINSOIL because ORGSOIL and MINSOIL are co-linear. Despite this outcome, some relationships in the plot diagram also appear in the model diagram. For example, a positive relationship between western hemlock and Amabalis fir with MINSOIL, a negative relationship between jack pine (Pinus banksiana [Lamb]) and lodgepole pine (Pinus contorta [Dougl. Ex Loud. var. latifolia Engelm.]), with both MINSOIL


Table 5 Total and percent variance explained in plot and model data estimated by the partial RDA for hierarchical scales defined for the CBM-CFS3 and partitions from their overlap. All RDAs are significant at p ≤ 0.01; Model n = 21,994; Plot n = 477. Partition (same as for results in Fig. 3)

Province (P) Reconciliation Unit (R) Ecozone (E) Stand (S) PS PR ES RS RE PE PES RES SPR PER PRES All Unaccounted (UN)

Percent variance explained

Total variance explained (C t ha−1)2





0.2 18.7 0.9 1.3 0.1 0.9 -0.4 0.8 −0.8 0.3 −0.1 2.8 0.5 −0.3 −0.7 24.1 75.9 100.0

0.2 9.1 0.3 24.6 −0.2 −0.1 13.2 0.6 0.7 0.5 4.4 25.7 0.2 −0.5 9.5 88.2 11.8 100.0

11.5 952.2 46.4 65.4 5.2 46.6 -20.6 39.2 −42.1 13.7 −5.5 143.1 27.2 −17.3 −37.7 1227.3 3868.7 5096.0a

6.1 230.2 8.5 622.9 −4.3 −2.5 335.5 15.4 17.7 12.3 110.3 650.9 6.1 −13.4 240.7 2236.3 299.7 2536.0

a The total variance from the partial RDA is lower that for the full RDA (5956 (C t ha−1)2) because the sample size was lower for the partial RDA due to missing values in the S partition.

and ORGSOIL (Fig. 4). This suggests that the model accurately reflects some, but not all, of the relationships between soil C stocks and leading species. MINSOIL and ORGSOIL in the model are more often positively correlated with higher productivity leading species (e.g., Douglas-fir, western hemlock, Amabalis fir) and negatively correlated with lower productivity leading species (e.g., black spruce, white spruce). These relationships are not consistently observed in the plot data and the variables in these relationships represent examples of variables available in the model, but not currently used to express variation in soil C stocks. Soil taxonomy is currently not used in the CBM-CFS3, but because of the links between soil development, soil C stocks, and soil taxonomy we used the full RDA for plot data to estimate the proportion of total variance that could be explained by soil order (ORDER), soil carbon modeling categories (SCMC) and soil great group and subgroup (SGGG) that we assume to be taxonomic classifications ranging from a lower to higher degree of co-variance with soil C stocks. We found that ORDER, SCMC and SGGG explained 7.5, 9.7 and 18.8% respectively, of total variance for MINSOIL and ORGSOIL in the plot data, confirming our assumption and demonstrating that soil taxonomy at the SGGG level explains more variation than all R-variables whether genus (11.6%), or genus and species (16%), were used to describe leading tree species (Table 4). When combined with all R-variables the increments for the percentage of total variance explained by ORDER, SCMC and SGGG were 4, 5.5 and 13.7%, respectively if both genus and species were used to identify the leading species. We compared the variance structure of SCMC and SGGG where we observed that as more taxonomic information is included (SGGG), more taxa separate out, explaining more variation in relation to soil C stocks (Fig. 5b). For example, the Humo-Ferric Podzol (HFP) SCMC is negatively correlated with ORGSOIL but when analyzed at the subgroup level (SGGG) Orthic Humo-Ferric Podzols (O.HFP) were negatively correlated, and Ortstein Humo-Ferric Podzols (OT.HFP) positively correlated, with ORGSOIL. The ortstein separates out from the Orthic, Humo-Ferric Podzols most likely because the impermeable Ortstein layer (Bockheim, 2011) creates a seasonally perched water table and restricts root penetration into the mineral soil which would favor accumulation of C in the organic horizons. Thus opportunities are available to improve model accuracy, and express more variation in modeled soil C estimates, by representing effects


C.H. Shaw et al. / Geoderma Regional 4 (2015) 114–125

Fig. 3. Comparison of results for the partial RDA of explained variance partitioned among variables associated with different scales used in the model based on the a) plot- and b) modelderived datasets for soil C. The largest proportion of explainable variance (24.1% of the total plot-based variance) in the plot data is explained by R-variables (Table 1) and interactions with R-variables. The largest proportion of explainable variance (88.2% of the total model-based variance) in the model data is explained by S-variables (Table 1) and interactions with S-variables. Interactions with no percentage value explained ≤1% of variance in the plot-derived, or model-derived, data. E, ecozone; P, province or territory; R, reconciliation unit; S, stand.

of leading species on soil C, a variable already available in the model, and to explain additional variance by using soil taxonomy, a variable currently not in the model but one that could be added in the future. 4. Discussion We undertook this exploratory study to understand the sources of variation that could help us improve the accuracy of predictions for soil C from the CBM-CFS3. Our results indicate this might be achieved by shifting attribution of soil C stock estimation to lower the proportion explained by tree productivity and raise the proportion explained by leading species effects apart from productivity. This is consistent with a recent inter-model comparison where Todd-Brown et al. (2013) determined that net primary productivity and soil temperature explained 62–93% of variation in modeled soil C stocks, but only 10% of variation in soil C data. In our analyses we calculated similar percentages with 78% of soil C variation being explained by productivity related variables and MAT in the model data, but only 4.2% in the plot data. The sizable portion of the variance (amount and proportion) explained by leading species in the plot data, aside from productivity,

may relate to variations in ecological processes that can be associated with leading species, but are not yet captured in the model. Leading species effects on soil and soil C were recently described in papers contributing to a special issue journal devoted to the topic (Prescott and Vesterdal, 2013). In that issue Vesterdal et al. (2013) reviewed tree species effects on soil C in boreal and temperate forests and concluded that there was consistent evidence that tree species affect soil C stocks but that the mechanisms are unclear, especially for the mineral soil C. Studies have shown that leading species influence forest floor (ORGSOIL) C stocks and dynamics because they differ in the quality and timing of inputs to dead organic matter pools (Finzi et al., 1998; Letang and de Groot, 2012; Hilger et al., 2012) which in turn influences the direction and intensity of pedological processes (van Cleve and Powers, 1995; Schaetzl and Anderson, 2005). The quality of dead organic matter can be a major determinant of decay rates (Flanagan and van Cleve, 1983; Trofymow et al., 2002). Although some relationships explaining variation in soil C were similar in the plot and modeled data our primary goal in comparing the plot- and model-derived ordination graphs was to look for specific relationships between leading species and soil C in the plot data that are not

C.H. Shaw et al. / Geoderma Regional 4 (2015) 114–125


Fig. 4. A comparison of ordinations based on results from the partial RDA for R-variables using the a) plot and b) model-derived datasets for soil C. Solid black arrows show that ORGSOIL and MINSOIL are highly correlated in the model (are collinear) and not correlated (orthogonal) in the plot data. Gray dashed arrows show the relationships between R-variables and soil C. The percent variation of ORGSOIL and MINSOIL explained by R-variables is 21.8% in the plot data, and 45.2% in the model data. The percent variation explained by axis 1 = 18.5% and axis 2 = 3.3% in the plot data, and 44.8% and 0.4%, respectively, in the model data. Arrows associated with MINSOIL and ORGSOIL are shorter in the plot- than model-derived ordinations because the proportion of total variance explained by R-variables in the ground plot data is less than half of that explained by R-variables in the model. Gray dashed arrows show the relationships between R-variables and soil C. In both model- and plot-derived partial RDAs leading species accounted for most of the variation and FFRET and COMP (HW or SW) combined accounted for a very small proportion of the total variance. MINSOIL, total C stocks in the mineral soil to 100 cm depth; ORGSOIL, total C stocks in the organic soil horizons.

observed in the model data, and therefore could be used to improve model accuracy. MINSOIL and ORGSOIL were consistently orthogonal in all plot-derived ordination graphs implying that the variables examined affected ORGSOIL and MINSOIL C stocks differently. This result is consistent with conclusions in a recent review by Vesterdal et al. (2013), a meta-analysis by Boca et al. (2014), and by detailed analyses of Olsson et al. (2009), Hobbie et al. (2007), Johnson et al. (2011), and Hoffmann et al. (2014) who found that trends in organic and mineral horizon C stocks or dynamics differed with pedological factors. Our analysis indicates that using information about leading species offers more opportunities to improve ORGSOIL than MINSOIL C estimates. For example, white spruce and black spruce had a strong positive correlation with ORGSOIL in the plot data, but not in the model data. The high ORGSOIL C stocks may be explained by slow decomposition rates associated with the litter and woody debris contributions to organic soil horizons for black and white spruce (Laiho and Prescott, 2004; Lang et al., 2009), and by the fact that they are also commonly associated with mosses in the boreal forest of Canada, especially when they

occur on less well or poorly drained sites (Bona et al., 2013). It is likely that accurate modeling in these, and some other boreal forests with significant moss associations, will require representation of moss contributions to ORGSOIL that are currently not included in the CBM-CFS3 (Bona et al., 2013). Like black and white spruce, Eastern white cedar was also positively correlated with ORGSOIL. Eastern white cedar is not as widespread or as well studied as the spruce species (Hofmeyer et al., 2007), however it is known to commonly occur in less well-drained sites transitional between bogs and upland hardwood forests. In poorly drained sites it is associated with woody swamp peat accumulation (Burns and Honkala, 1990), as opposed to the moss-derived peat accumulations often associated with black and white spruce. ORGSOIL had a negative correlation with trembling aspen (Populus tremuloides Michx.) in the plot data, but not in the model data. The negative correlation in the plot data is likely related to the high decay rates of trembling aspen litter and dead wood relative to other leading species found in Canada's forests (Prescott and Vesterdal, 2005; Angers et al., 2012; Braise and Drouin, 2012) and that in general, angiosperm wood


C.H. Shaw et al. / Geoderma Regional 4 (2015) 114–125

Fig. 5. Ordinations showing the relationship between soil taxonomy and soil C based on the results from the full RDA for the plot soil C data only. Two levels in detail of soil taxonomy are compared; a) the simplified soil carbon modeling categories (SCMC) and b) the full taxonomic detail to the level of the subgroup and great group (SGGG). The percent variation of ORGSOIL and MINSOIL explained by SCMCs is 9.7% and by SGGGs is 18.8%. The percent variation explained by axis 1 = 8.4% and axis 2 = 1.3% for SCMCs, and 13.1% and 5.7%, respectively, for SGGGs. Solid black arrows show the relationship between ORGSOIL and MINSOIL. The arrows associated with MINSOIL and ORGSOIL are larger for the SGGGs than for the SCMCs because providing more resolution in soil taxonomy provides more explanation of variance observed for soil C stocks. Gray dashed arrows show the relationships between soil taxonomy and soil C. MINSOIL, total C stocks in the mineral soil to 100 cm depth; ORGSOIL, total C stocks in the organic soil horizons.

decays at faster rates than gymnosperm wood (Weedon et al., 2009). Douglas-fir and ORGSOIL are negatively correlated in the plot data, but positively correlated in the modeled data. Although high decomposition rates could account for low ORGSOIL C stocks for Douglasfir the available evidence indicates that Douglas-fir litter (Prescott et al., 2004; Moore et al., 2006) and wood (Means et al., 1985) do not decompose at rates significantly higher than other leading species. It is more likely that ORGSOIL C stocks are low for this leading species because it is commonly associated with complex fire histories characterized by low frequency high intensity fires in

combination with high frequency low intensity fires (Gray and Riccius, 1999; Wetzel and Fonda, 2000). Frequent low intensity fires do not kill the overstory Douglas-fir trees but reduce ORGSOIL horizons, and dead organic matter (e.g., downed dead wood) available for input to ORGSOIL horizons. Currently, the CBM-CFS3 only simulates stand-replacing fires during the spin-up phase, but not the high frequency, low intensity fires that may be controlling accumulation of C in the ORGSOIL under Douglas-fir. Through our analysis we found that soil taxonomy could explain variance in addition to that already accounted for by productivity and

C.H. Shaw et al. / Geoderma Regional 4 (2015) 114–125

leading species, and that the magnitude of the increment in the percent variation explained increased as soil taxonomic information became more precise. Results from the ordinations indicate that there are opportunities to improve soil C stock estimation for both ORGSOIL and MINSOIL in relation to soil taxonomy. This result is consistent with site-level studies where vegetation (e.g., tree species) was held constant and soil type was allowed to vary, or in factorial experiments designed to examine interaction effects between tree species and soil type, C stocks varied depending on soil type (Kulmatiski et al., 2004; Kelly and Mays, 2005). At regional or national scales Vesterdal et al. (2013) concluded that mineral soil C was more strongly influenced by soil type or climate than by tree species. We also found that soil taxonomy could compensate for explanatory power of leading species if leading species identification was less precise (i.e., genus only, or generic hardwood or softwood) than genus and species. The ability of leading species and soil taxonomy to compensate for one another in explaining variation in soil C stocks, is potentially valuable in cases where leading species may not be precisely identified (e.g., remote locations, noncommercial forests). It would be impossible to discuss all pedological processes explaining the relationships between soil C and taxonomy here, but readers are referred to Shaw et al. (2008) for a description of forest soil C stocks in relation to soil taxon for Canada, and a special issue of the Canadian Journal of Soil Science (Pennock et al., 2011) where current understanding of the genesis, distribution and classification of soil orders in Canada are discussed in detail. However, we will briefly highlight some of the relevant major trends and observations from this study. Most of the variation in C stocks explained by soil taxonomy is associated with a gradient of increasing soil C stocks and soil taxa ranging from well-drained Luvisolic and Brunisolic soils to poorly drained soil taxa such as Gleysols. This trend has been reported by others (Little et al., 2002; Bois et al., 2009; Johnson et al., 2011) and knowledge about the relationships between geomorphology and hydrology linking soil taxa to the landscape is entrenched in the fundamentals of pedology (Chapter 13 in Schaetzl and Anderson, 2005). Modeling of hydrological effects on forest soil C dynamics at a national scale in Canada is challenged by the lack of suitable national-scale hydrologic models and data. Recent conceptual and modeling efforts (Devito et al., 2005; Fan and Miguez-Macho, 2011) show promise but have not yet been integrated with forest C models. Therefore, a useful interim solution to improve soil C models may be to stratify by soil taxa linked to a range in drainage conditions (Shaw et al., 2008). The greater separation of ORGSOIL and MINSOIL with increased soil taxonomic information (from SCMCs to SGGGs) provides further insight to processes that should be considered in soil C modeling that differ between organic and mineral soil horizons. All Gleysols and gleyed Regosols had a strong positive correlation with ORGSOIL. Gleysolic Turbic Cryosols, Podzols with ortstein horizons, and Orthic Humic Podzols were also positively correlated with ORGSOIL. These are soil taxa where classification is based on criteria in the mineral soil horizons that are indicative of sufficient long-term (Gleysols) or periodic soil saturation to influence accumulation of C in the organic soil horizons. Gleyed features were significant enough to determine the classification of other soils in this analysis (GLE.DYB; GL.SC) but they did not have a positive correlation with ORGSOIL C stocks, suggesting that the degree of saturation by water in these soils is not sufficient to influence ORGSOIL C stock accumulation. The presence of impermeable layers such as ortstein (Bockheim, 2011) or soils with poor physical properties such as the gleyed Regosols, likely played a prominent role in the relationship between gleyed subgroups and C accumulation in the ORGSOIL because they can create seasonal perched water tables and restrict root penetration into the mineral soil. Other processes that our analysis suggests be considered in modeling of MINSOIL C include podzolization (with Humo-Ferric Podzols distinguished from Ferro-Humic Podzols), accumulation of buried mineral A horizons in floodplains associated with Cumulic Regosols, and genesis of mull Ah horizons (e.g., Sombric Brunisols). Soils negatively correlated


with MINSOIL or ORGSOIL are primarily associated with lessivage and transport of dissolved or particulate organic matter out of the pedon (e.g., Luvisols; Howitt and Pawluk, 1985) or associated with high fire frequencies (Brunisols) with reduced inputs of dead organic matter to soils (Shaw et al., 2008). A large proportion of total soil C variance in the plot data still remains unexplained. Although we cannot expect to explain all the variation that occurs at the scale of a pedon, there is reason to believe that through further analyses of a larger and more comprehensive dataset, more variation could be explained. Based on other studies additional factors to consider include soil texture, mineralogy, cation chemistry and geological parent material (Giardina et al., 2001; Sanger et al., 1998; Torn et al., 1997; Hobbie et al., 2007; Heckman et al., 2009, 2011) because of their effects on C stabilization and decay dynamics.

5. Conclusions One objective of the initialization process in many soil and ecosystem models is to distinguish soils with different C density. Our exploratory analysis demonstrates that leading species, a variable already available in the CBM-CFS3, and soil taxonomy, a variable currently not in the model but that could be added in a future version, have the potential to be used to stratify the landscape into meaningful modeling categories in a way that differentiates C dynamics for organic and mineral soil horizons. Processes active over pedological time frames are a major determinant of accumulated C stocks. Even if those processes do not measurably influence stock changes over a time frame meaningful to current C assessments (decades to hundreds of years) they have created soil environments that influence current C stocks and dynamics. Once parameters are determined for the central tendency of C retention for leading tree species/soil taxonomy combinations, additional parameters can be developed to reflect contemporary processes responding to changes in climate, disturbances, management, or land-use change. Most soil C models are derivatives of models developed for grassland or agricultural systems where soils are united by a similar pedological theme (central concept of landscapes with low relief, grasses as the major source of C input, organic rich mineral Ah horizon). However, modeling of forest soil C dynamics is more complex than previously thought (Dungait et al., 2012) leading the soil modeling community to call for a change to fundamental modeling approaches because current constructs of soil models will be unable to accurately predict response to climate change as they do not directly represent soil (e.g., microbial activity, stabilization with inorganic components, and aeration) and landscape processes (e.g., permafrost thaw, thermokarst, change in insulation from peat, forest, or snow cover, and hydrology) controlling heterotrophic respiration (Allison et al., 2010; Conant et al., 2011; Schmidt et al., 2011). These processes are more common and important in the forested areas of Canada compared with other land types. It may be some time before these new constructs are developed and tested for implementation in national-scale reporting. The empirical approach described here can act as an interim solution that accounts for more variation in soil C stocks than existing approaches, until such time as soil modeling knowledge is sufficient to quantitatively integrate dynamic pedological processes into soil C models.

Acknowledgments We acknowledge Sue Mayer for graphics production; Brenda Laishley for editorial assistance, Michelle Filiatrault for map production and assistance with database management, and Dan Thompson and Paul Sanborn for feedback on an earlier version of this manuscript. Funding for this project was provided by the Natural Resources Canada under Canada's Clean Air Agenda.


C.H. Shaw et al. / Geoderma Regional 4 (2015) 114–125

References (CanSIS) Canadian Soil Information Service, 2007. Soil Pedon Data From Yukon Territory. Canada, Agriculture and Agri-Food Canada, Ottawa. (CEMA) Cumulative Environmental Management Association, 2004. Data Acquired in 2004. Contact on-line at. (ELCG) Ecological Land Classification Group, 2005. Ontario Terrestrial Assessment Program, Ontario Ministry of Natural Resources. Sault Ste, Marie, Ontario, Canada. (NFI) National Forest Inventory, 2008. Canada's National Forest Inventory — Ground Sampling Guidelines, Version 5.0 (Available from (NFI) National Forest Inventory, 2010. Canada's National Forest Inventory, National Standards for Ground Plots Compilation Procedures Version 1.7.2 [Draft]. Canadian Council of Forest Ministers, Ottawa, ON. (NFI) National Forest Inventory, 2011. Canada's National Forest Inventory National Standard For Ground Plots, Data Dictionary February 2011, Version 5.1.3. Canadian Council of Forest Ministers, Ottawa, ON (Available from: documentation/ground_plot/Gp_data_dictionary_v5.1.3.pdf (cited 12.23.14)). (RDB) Resource Data Branch, 2002. Alberta Ecological Site Data and Ecological Site Description Manual. Alberta Sustainable Resource Development (Now, Environment and Sustainable Resource Development), Edmonton, AB. (SCWG) Soil Classification Working Group, 1998. The Canadian System of Soil Classification. 3rd ed. Publication 1646, Agric. Agri-Food Can, Ottawa, ON. Allison, S.D., Wallenstein, M.D., Bradford, M.A., 2010. Soil-carbon response to warming dependent on microbial physiology. Nat. Geosci. 3 (5), 336–340. Angers, V.A., Drapeau, P., Bergeron, Y., 2012. Mineralization rates and factors influencing snag decay in four North American boreal tree species. Can. J. For. Res. 42, 157–166. Boca, A., Van Miegroet, H., Gruselle, M.-C., 2014. Forest overstory effect on soil organic carbon storage: a meta-analysis. Soil Sci. Soc. Am. J. 78, S35–S47. Bockheim, J., 2011. Distribution and genesis of ortstein and placic horizons in soils of the USA: a review. Soil Sci. Soc. Am. J. 75 (3), 994–1005. Bois, C.H., Janzen, D.T., Sanborn, P.T., Fredeen, A.L., 2009. Contrasting total carbon stocks between ecological site series in a subboreal spruce research forest in central British Columbia. Can. J. For. Res. 39, 897–907. Boisvenue, C., Running, S.W., 2010. Simulations show decreasing carbon stocks and potential for carbon emissions in Rocky Mountain forests over the next century. Ecol. Appl. 20 (5), 1302–1319. Bona, K.A., Fyles, J., Shaw, C.H., Kurz, W.A., 2013. Are mosses required to accurately predict upland black spruce forest soil carbon in national-scale forest C accounting models? Ecosystems 16, 1071–1086. Braise, S., Drouin, P., 2012. Interactions between deadwood and soil characteristics in a natural boreal trembling aspen–jack pine stand. Can. J. For. Res. 42, 1456–1466. Brandt, J.P., 2009. The extent of the North American boreal zone. Environ. Rev. 17, 101–161. (tech. coords)Burns, R.M., Honkala, B.H., 1990. Silvics of North America: 1. Conifers. Agriculture Handbook 654. U.S. Dep. Agric., For. Serv., Washington, DC (Accessed on-line April 2, 2013. occidentalis.htm). Conant, R.T., Ryan, M.G., Ågren, G.I., Birge, H.E., Davidson, E.A., Eliasson, P.E., Evans, S.E., Frey, S.D., Giardina, C.P., Hopkins, F.M., Hyvönen, R., Kirschbaum, M.U.F., Lavallee, J.M., Leifeld, J., Parton, W.J., Megan Steinweg, J., Wallenstein, M.D., Martin Wetterstedt, J.Å., Bradford, M.A., 2011. Temperature and soil organic matter decomposition rates — synthesis of current knowledge and a way forward. Glob. Chang. Biol. 17 (11), 3392–3404. Deckers, J., Spaargaren, O.C., Nachtergaele, F.O. (Eds.), 1998. World Reference Base for Soil Resources. International Society of Soil Science Working Group RB, 1st ed. . Acco. vols. 1 and 2. KULeuven Academic Press (93 and 96 pp.). Devito, K., Creed, I., Gan, T., Mendoza, C., Petrone, R., Silins, U., Smerdon, B., 2005. A framework for broad-scale classification of hydrologic response units on the Boreal Plain: is topography the last thing to consider? Hydrol. Process. 19, 1705–1714. Dungait, J.A.J., Hopkins, D.W., Gregory, A.S., Whitmore, A.P., 2012. Soil organic matter turnover is governed by accessibility not recalcitrance. Glob. Chang. Biol. 18 (6), 1781–1796. Environment Canada, 2010. National Inventory Report 1990–2008: Greenhouse Gas Sources and Sinks in Canada. Environ. Can., Ottawa, ON (Available from: (cited 04.05.13)). Fan, Y., Miguez-Macho, G., 2011. A simple hydrologic framework for simulating wetlands in climate and earth system models. Clim. Dyn. 37, 253–278. Finzi, A.C., Van Breemen, N., Canham, C.D., 1998. Canopy tree–soil interactions within temperate forests: species effects on soil carbon and nitrogen. Ecol. Appl. 8 (2), 440–446. Flanagan, P.W., van Cleve, K., 1983. Nutrient cycling in relation to decomposition and organic-matter quality in taiga ecosystems. Can. J. For. Res. 13, 795–817. Foereid, B., Bellamy, P.H., Holden, A., Kirk, G.J.D., 2012. On the initialization of soil carbon models and its effects on model predictions for England and Wales. Eur. J. Soil Sci. 63, 32–41. Giardina, C.P., Ryan, M.G., Hubbard, R.M., Binkley, D., 2001. Tree species and soil textural controls on carbon and nitrogen mineralization rates. Soil Sci. Soc. Am. J. 65, 1272–1279. Gottschalk, P., Smith, J.U., Wattenbach, M., Bellarby, J., Stehfest, E., Arnell, N., Osborn, T.J., Smith, P., 2012. How will organic carbon stocks in mineral soil evolve under future climate? Global projections using RothC for a range of climate scenarios. Biogeosci. Discuss. 9, 411–451. Gray, R., Riccius, E., 1999. Historical fire regime for the Pothole Creek interior Douglas-fir research site. B. C. Minist. For. Res. Program Working Paper 38. Victoria, BC, Canada.

Harmon, M.E., Marks, B., 2002. Effects of silvicultural practices on carbon stores in Douglas-fir-western hemlock forests in the Pacific Northwest, USA: results from a simulation model. Can. J. For. Res. 32 (5), 863–877. Hashimoto, S., Wattenbach, M., Smith, P., 2011. A new scheme for initializing process-based ecosystem models by scaling soil carbon pools. Ecol. Model. 222, 3598–3602. Heckman, K., Welty-Bernard, A., Rasmussen, C., Schwartz, E., 2009. Geologic controls of soil carbon cycling and microbial dynamics in temperate conifer forests. Chem. Geol. 267, 12–23. Heckman, K., Vazquez-Ortega, A., Gao, X., Chorover, J., Rasmussen, C., 2011. Changes in water extractable organic matter during incubation of forest floor material in the presence of quartz, goethite and gibbsite surfaces. Geochim. Cosmochim. Acta 75, 4295–4309. Hilger, A.B., Shaw, C.H., Metsaranta, J.M., Kurz, W.A., 2012. Estimation of snag carbon transfer rates by ecozone and lead species for forest in Canada. Ecol. Appl. 22 (8), 2078–2090. Hobbie, S.E., Ogdahl, M., Chorover, J., Chadwick, O.A., Oleksyn, J., Zytkowiak, R., Reich, P.B., 2007. Tree species effects on soil organic matter dynamics: the role of soil cation composition. Ecosystems 10, 999–1018. Hoffmann, U., Hoffmann, T., Johnson, E.A., Kuhn, N.J., 2014. Assessment of variability and uncertainty of soil organic carbon in a mountainous boreal forest (Canadian Rocky Mountains, Alberta). Catena 113, 107–121. Hofmeyer, P.V., Knefic, L.S., Seymour, R.S., 2007. Northern white-cedar (Thuja occidentalis L.). An annotated bibliography. Cooperative Forestry Research Unit Research Report CFRU RR 07-01. University of Maine, USA (Accessed on-line April 2, 2013. http:// Howitt, R., Pawluk, S., 1985. The genesis of a Gray Luvisol within the boreal forest region. II. Dynamic pedology. Can. J. Soil Sci. 6, 9–19. Jandl, R., Lindner, M., Vesterdal, L., Bauwens, B., Baritz, R., Hagedorn, F., Johnson, D.W., Minkkinen, K., Byrne, K., 2007. How strongly can forest management influence soil carbon sequestration? Geoderma 137 (3–4), 253–268. Johnson, K.D., Harden, J., McGuire, A.D., Bliss, N.B., Bockheim, J.G., Clar, M., NettletonHollingsworth, T., Torre Jorgenson, M., Kane, E.S., Mack, M., O'Donnell, J., Ping, C.-L., Schuur, E.A.G., Turetsky, M.R., Valentine, D.W., 2011. Soil carbon distribution in Alaska in relation to soil-forming factors. Geoderma 167–168, 71–84. Kelly, J.M., Mays, P.A., 2005. Soil carbon changes after 26 years in a Cumberland Plateau hardwood forest. Soil Sci. Soc. Am. J. 69, 691–694. Komarov, A., Chertov, O., Zudin, S., Nadporozhskaya, M., Mikhailov, A., Bykhovets, S., Zudina, E., Zoubkova, E., 2003. EFIMOD 2 — a model of growth and cycling of elements in boreal forest ecosystems. Ecol. Model. 170 (2–3), 373–392. Kull, S.J., Rampley, G.J., Morken, S., Metsaranta, J., Neilson, E.T., Kurz, W.A., 2011. Operational-scale Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3) Version 1.2: User's Guide. 2011. Nat. Resour. Can., Can. For. Serv., North. For. Cent., Edmonton, AB, Canada (370 pp.). Kulmatiski, A., Vogt, D.J., Siccama, T.G., Tilley, J.P., Kolesinskas, K., Wickwire, T.W., Larson, B.C., 2004. Landscape determinants of soil carbon and nitrogen storage in southern New England. Soil Sci. Soc. Am. J. 68, 2014–2022. Kurz, W.A., Apps, M.J., 1999. A 70-year retrospective analysis of carbon fluxes in the Canadian forest sector. Ecol. Appl. 9 (2), 526–547. Kurz, W.A., Apps, M.J., Webb, T.M., McNamee, P.J., 1992. The carbon budget of the Canadian forest sector: phase 1. For. Can., Northwest Reg., Edmonton, AB. Inf. Rep. NOR-X326. Kurz, W.A., Dymond, C.C., White, T.M., Stinson, G., Shaw, C.H., Rampley, G.J., Smyth, C., Simpson, B.N., Neilson, E.T., Trofymow, J.A., Metsaranta, J., Apps, M.J., 2009. CBMCFS3: a model of carbon-dynamics in forestry and land-use change implementing IPCC standards. Ecol. Model. 220, 480–504. Kurz, W.A., Shaw, C.H., Boisvenue, C., Stinson, G., Metsaranta, J., Leckie, D., Dyk, A., Smyth, C., Neilson, E.T., 2013. Carbon in Canada's boreal forest — a synthesis. Environ. Rev. 21, 260–292. Laiho, R., Prescott, C.E., 2004. Decay and nutrient dynamics of coarse woody debris in northern coniferous forests: a synthesis. Can. J. For. Res. 34 (4), 763–777. Lal, R., 2005. Forest soils and carbon sequestration. For. Ecol. Manag. 220, 242–258. Lang, S.I., Cornelissen, J.H.C., Klahn, T., van Logtestijn, R.S.P., Broekman, R., Schweikert, W., Aerts, R., 2009. An experimental comparison of chemical traits and litter decomposition rates in a diverse range of subarctic bryophytes, lichen and vascular plant species. J. Ecol. 97, 886–900. Legendre, P., Legendre, L., 1998. Chapter 11, Canonical Analysis. In: Numerical Ecology. 2nd ed. Elsevier Science BV, Amsterdam, pp. 579–593. Letang, D.L., de Groot, W.J., 2012. Forest floor depths and fuel loads in upland Canadian forests. Can. J. For. Res. 42, 1551–1565. Little, T.I., Pluth, D.J., Corns, I.G.W., Gilmore, D.W., 2002. Post-fire forest floor development along toposequences of white spruce-trembling aspen mixedwood communities in west-central Alberta. Can. J. For. Res. 32 (5), 892–902. McGill, W.B., 1996. Review and classification of ten soil organic matter (SOM) models. NATO ASI series I global environmental change. Heidelberg, Germany. NATO 38, 111–132. McKenney, D.W., Hutchinson, M.F., Kesteven, J.L., Veneir, L.A., 2001. Canada's plant hardiness zones revisited using modern climate interpolation techniques. Can. J. Plant Sci. 81 (1), 129–143. Means, J.E., Cromack Jr., K., MacMillan, P., 1985. Comparison of decomposition models using wood density of Douglas-fir logs. Can. J. For. Res. 15, 1092–1098. Moore, T.R., Trofymow, J.A., Prescott, C.E., Fyles, J., Titus, B.D., CIDET Working Group, 2006. Patterns of carbon, nitrogen, and phosphorus dynamics in decomposing foliar litter in Canadian forests. Ecosystems 9, 46–62. Odeh, I.O.A., Chittleborough, D.J., McBrathney, A.B., 1991. Eluciation of soil-landform interrelationships by canonical ordination analysis. Geoderma 49, 1–32.

C.H. Shaw et al. / Geoderma Regional 4 (2015) 114–125 Olsson, M.T., Erlandsson, M., Lundin, L., Nilsson, T., Nilsson, A., Stedahl, J., 2009. Organic carbon stocks in Swedish Podzol soils in relation to soil hydrology and other site characteristics. Silva Fenn. 43 (2), 209–222. Palosuo, T., Foereid, B., Svensson, M., Shurpali, N., Lehtonen, A., Herbst, M., Linkosalo, T., Ortiz, C., Rampazzo Todorovic, G., Marcinkonix, S., Li, C., Jandl, R., 2012. A multimodel comparison of soil carbon assessment of a coniferous forest stand. Environ. Model Softw. 35, 38–49. Pan, Y., Birdsey, R.A., Fang, J., Houghton, R., Kauppi, P.E., Kurz, W.A., Phillips, O.L., Shvidenko, A., Lewis, S.L., Canadell, J.G., Ciais, P., Jackson, R.B., Pacala, S.W., McGuire, A.D., Piao, S., Rautiainen, A., Sitch, S., Hayes, D., 2011. A large and persistent carbon sink in the world's forests. Science 333 (6045), 988–993. Peltoniemi, M., Thurig, E., Ogle, S., Palosuo, T., Schrumpf, M., Wutzler, T., Butterback Bahl, K., Chertov, O., Komarov, A., Mikhailov, A., Gardenas, A., Perry, C., Liski, J., Smith, P., Makipaa, R., 2007. Models in country scale carbon accounting of forest soils. Silva Fenn. 41 (3), 575–602. Pennock, D., Arocena, J.M., Scott, C.A., 2011. Soils of Canada—preface. Can. J. Soil Sci. 91 (5), 671–673. Prescott, C.E., Vesterdal, L., 2005. Effects of British Columbia tree species on forest floor chemistry. Chapter 2. In: Binkley, D., Menyailo, O. (Eds.), Tree Species Effects on Soils. Implications for Global Climate Change. Springer, Dordrecht, The Netherlands, pp. 17–29. Prescott, C.E., Vesterdal, L., 2013. Tree species effects on soil in temperate and boreal forests; emerging themes and research needs. For. Ecol. Manag. Spec. Issue 309 (1), 1–3. Prescott, C.E., Blevins, L.L., Staley, C., 2004. Litter decomposition in B.C. forests: controlling factors and influences of forestry activities. BC J. Ecosyst. Manag. 5 (2), 30–43. R. Development Core Team, 2011. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria (http://www.R-project. org/). Running, S.W., Hunt, E.R.J., 1993. Generalization of a forest ecosystem process model for other biomes, BIOME-BGC, and an application for global-scale-models. In: Ehleringer, J.R., Field, C.B. (Eds.), Scaling Physiological Processes: Leaf to Globe. Academic Press, San Diego, CA, USA, pp. 141–157. Sanger, L.J., Cox, P., Splatt, P., Whelan, M., Anderson, J.M., 1998. Variability in the quality and potential decomposability of Pinus sylvestris litter from sites with different soil characteristics: acid detergent fibre (ADF) and carbohydrate signatures. Soil Biol. Biochem. 30 (4), 455–461. Schaetzl, R., Anderson, S., 2005. Soils: Genesis and Geomorphology. Cambridge Univ, Press, New York, NY. Schelhaas, M.J., van Esch, P.W., Groen, T.A., de Jong, B.H.J., Kanninen, M., Liski, J., Masera, O., Mohren, G.M.J., Nabuurs, G.J., Palosuo, T., Pedroni, L., Vallejo, A., Vilén, T., 2004. CO2FIX V 3.1 — A Modelling Framework for Quantifying Carbon Sequestration in Forest Ecosystems. Alterra-rapport, Wageningen, Alterra, Wageningen. Schlesinger, W.H., Andrews, J.A., 2000. Soil respiration and the global carbon cycle. Biogeochemistry 48, 7–20. Schmidt, M.W.I., Torn, M.S., Abiven, S., Dittmar, T., Guggenberger, G., Janssens, I.A., Kleber, M., Kogel-Knabner, I., Lehmann, J., Manning, D.A.C., Nannipieri, P., Rasse, D.P., Weiner, S., Trumbore, S.E., 2011. Persistence of soil organic matter as an ecosystem property. Nature 478 (7367), 49–56. Shanin, V.N., Komarov, A.S., Mikhailov, A.V., Bykhovets, S.S., 2010. Modelling carbon and nitrogen dynamics in forest ecosystems of Central Russia under different climate change scenarios and forest management regimes. Ecol. Model. 222 (14), 2262–2275. Shaw, C.H., Bhatti, J.S., Sabourin, K.J., 2005. An ecosystem carbon database for Canadian forests. Nat. Resour. Can., Ca. For. Serv., North. For. Centr., Edmonton, Alberta, Inf. Rep. NOR-X-403. Shaw, C., Banfield, E., Kurz, W.A., 2008. Stratifying soils into pedogenically similar categories for modeling forest soil carbon. Can. J. Soil Sci. 88, 501–516.


Shaw, C.H., Hilger, A.B., Metsaranta, J., Kurz, W.A., Russo, G., Eichel, F.G., Stinson, G.C., Smyth, C.M., Filiatrault, M., 2014. Evaluation of simulated estimates of forest ecosystem carbon stocks using ground plot data from Canada's National Forest Inventory. Ecol. Model. 272, 323–347. Siltanen, R.M., Apps, M.J., Zoltai, S.C., Mair, R.M., Strong, W.L., 1997. A Soil Profile and Organic Carbon Data Base for Canadian Forest and Tundra Mineral Soils. Natural Resources Canada, Canadian Forest Service, Northern Forestry Centre, Edmonton, AB. Smyth, C.E., Stinson, G., Neilson, E., Lempriere, T., Hafer, M., Rampley, G.J., Kurz, W.A., 2014. Quantifying the biophysical climate change mitigation potential of Canada's forest sector. Biogeosciences 11, 3515–3529. Stinson, G., Kurz, W.A., Smyth, C.E., Neilson, E.T., Dymond, C.C., Metsaranta, J.M., Boisvenue, C., Rampley, G.J., Li, Q., White, T.M., Blain, D., 2011. An inventory-based analysis of Canada's managed forest carbon dynamics, 1990 to 2008. Glob. Chang. Biol. 17 (6), 2227–2244. Survey Staff, Soil, 1992. Keys to soil taxonomy. SMMS Technical Monograph No. 19, 5th ed. Pocahontas Press, Inc., Blacksburg, Virginia (556 pp.). ter Braak, C.J.F., 1987. Ordination. In: Jongman, R.H.G., ter Braak, C.J.F., Tongeren, O.F.R. (Eds.), Data Analysis in Community and Landscape Ecology. Center for Agricultural Publishing and Documentation, Wageningen, The Netherlands, pp. 91–173. ter Braak, C.J.F., Prentice, I.C., 1988. A theory of gradient analysis. Adv. Ecol. Res. 18, 259–267. Todd-Brown, K.E.O., Randerson, J.T., Post, W.M., Hoffman, F.M., Tarnocai, C., Schuur, E.A.G., Allison, S.D., 2013. Causes of variation in soil carbon simulations from CMIP5 Earth system models and comparison with observations. Biogeosciences 10, 1717–1736. Torn, M.S., Trumbore, S.E., Chadwick, O.A., Vitousek, P.M., Hendricks, D.M., 1997. Mineral control of soil organic carbon storage and turnover. Nature (London) 389, 170–173. Trofymow, J.A., Moore, T.R., Titus, B., Prescott, C., Morrison, I., Siltanen, M., Smith, S., Fyles, J., Wein, R., Camire, C., Duschene, L., Kozak, L., Kranabetter, M., Visser, S., 2002. Rates of litter decomposition over 6 years in Canadian forests: influence of litter quality and climate. Can. J. For. Res. 32, 789–804. van Cleve, K., Powers, R.F., 1995. Chapter 9. Soil carbon, soil formation, and ecosystem development. In: McFee, W.W., Kelly, J.M. (Eds.), Carbon Forms and Function in Forest Soils. Soil Sci. Soc. Am, Madison, WI, pp. 155–200. Vesterdal, L., Clarke, N., Sigurdsson, B.D., Gundersen, P., 2013. Do tree species influence soil carbon stocks in temperate and boreal forests? For. Ecol. Manag. 309 (1), 4–18. Weedon, J.T., Cornwell, W.K., Cornelissen, J.H.C., Zanne, A.E., Wirth, C., Coomes, D.A., 2009. Global meta-analysis of wood decomposition rates: a role for trait variation among tree species? Ecol. Lett. 12, 45–56. Wetzel, S.A., Fonda, R.W., 2000. Fire history of Douglas-fir forest in the Morse Creek drainage of Olympic National Park, Washington. Northwest Sci. 74 (4), 263–279. Wutzler, T., Reichstein, M., 2007. Soils apart from equilibrium—consequences for soil carbon balance modelling. Biogeosciences 4 (1), 125–136. Xu, X., Liu, W., Kiely, G., 2011. Modeling the change in soil organic carbon of grassland in response to climate change: effects of measured versus modelled carbon pools for initializing the Rothamsted Carbon model. Agric. Ecosyst. Environ. 140 (3–4), 372–381. Zimmerman, M., Leifeld, J., Schmidt, M.W.I., Smith, P., Fuhrer, J., 2007. Measured soil organic matter fractions can be related to pools in the RothC model. Eur. J. Soil Sci. 58, 658–667. Zoltai, S.C., Siltanen, R.M., Johnson, J.D., 2000. A wetland data base for the western boreal, subarctic, and arctic regions of Canada. Nat. Resour. Can., Can. For. Serv., North. For. Cent., Edmonton, Alberta. Inf. Rep. NOR-X-368.