Data-mining analysis of the global distribution of soil carbon in observational databases and Earth system models

Geosci. Model Dev., 10, 1321–1337, 2017 www.geosci-model-dev.net/10/1321/2017/ doi:10.5194/gmd-10-1321-2017 © Author(s) 2017. CC Attribution 3.0 License. Data-mining analysis of the global distribution of soil carbon in observational databases and Earth system models Shoji Hashimoto1, Kazuki Nanko1, Boris ˇ Tupek2, and Aleksi Lehtonen2 1Forestry and Forest Products Research Institute (FFPRI), Tsukuba, Japan 2Natural Resources Institute Finland, Latokartanonkaari 9, Helsinki, Finland Correspondence to: Shoji Hashimoto ([email protected]) Received: 27 May 2016 – Discussion started: 30 June 2016 Revised: 8 February 2017 – Accepted: 28 February 2017 – Published: 28 March 2017 Abstract. Future climate change will dramatically change the carbon balance in the soil, and this change will affect the terrestrial carbon stock and the climate itself. Earth system models (ESMs) are used to understand the current climate and to project future climate conditions, but the soil organic carbon (SOC) stock simulated by ESMs and those of observational databases are not well correlated when the two are compared at fine grid scales. However, the specific key processes and factors, as well as the relationships among these factors that govern the SOC stock, remain unclear; the inclusion of such missing information would improve the agreement between modeled and observational data. In this study, we sought to identify the influential factors that govern global SOC distribution in observational databases, as well as those simulated by ESMs. We used a data-mining (machine-learning) (boosted regression trees – BRT) scheme to identify the factors affecting the SOC stock. We applied BRT scheme to three observational databases and 15 ESM outputs from the fifth phase of the Coupled Model Intercomparison Project (CMIP5) and examined the effects of 13 variables/factors categorized into five groups (climate, soil property, topography, vegetation, and land-use history). Globally, the contributions of mean annual temperature, clay content, carbon-to-nitrogen (CN) ratio, wetland ratio, and land cover were high in observational databases, whereas the contributions of the mean annual temperature, land cover, and net primary productivity (NPP) were predominant in the SOC distribution in ESMs. A comparison of the influential factors at a global scale revealed that the most distinct differences between the SOCs from the observational databases and ESMs were the low clay content and CN ratio contributions, and the high NPP contribution in the ESMs. The results of this study will aid in identifying the causes of the current mismatches between observational SOC databases and ESM outputs and improve the modeling of terrestrial carbon dynamics in ESMs. This study also reveals how a data-mining algorithm can be used to assess model outputs. 1 Introduction Soil is the largest organic carbon stock in terrestrial ecosystems (Batjes, 1996; IPCC, 2013; Köchy et al., 2015). The soil organic carbon (SOC) stock represents a balance between carbon inputs to soil and carbon losses from soil via decomposition and dissolved organic carbon, and this influx and efflux of soil carbon is controlled directly and indirectly by environmental conditions (Carvalhais et al., 2014; Schimel et al., 1994). Future climate change will dramatically affect the global soil carbon balance (Bond-Lamberty and Thomson, 2010; Crowther et al., 2016; Friedlingstein et al., 2006; Hashimoto et al., 2011, 2015), and this change will affect terrestrial carbon and consequently the climate itself (Cox et al., 2000; Zaehle, 2013). Earth system models (ESMs) were developed to understand the current climate and provide future climate projections, and these models incorporate the terrestrial carbon cycle, including SOC (Arora et al., 2013; Friedlingstein et al., 2014). In ecosystem carbon cycle models of ESMs, SOC is calculated as the balance between carbon inputs via dead organic matter and carbon emissions via organic matter decomposition, with both processes influenced by temperature and water conditions. SOC dynamics have a critical influence on Published by Copernicus Publications on behalf of the European Geosciences Union. 1322 S. Hashimoto et al.: Data-mining analysis of the global distribution of soil carbon the land carbon sink in ESM simulations (Friedlingstein et al., 2014). Observational global soil databases are often used as benchmarks to examine whether ESMs successfully describe the global distribution of soil carbon stocks (Anav et al., 2013; Hararuk et al., 2014; Todd-Brown et al., 2013; Wieder et al., 2014). Several global soil databases have been developed over the past two decades, and several are undergoing further improvement (Scharlemann et al., 2014). Certain databases describe the global distribution of soil physiochemical properties and enable calculations of the global distribution of SOC stocks (e.g., Harmonized World Soil Database – HWSD), whereas others provide SOC stocks by default (e.g., International Geosphere Biosphere Programme’s (IGBP) Data and Information System (DIS) database). These databases incorporate observed data points with global coverage, although there are biases in the spatial distribution or density of the data points. In these databases, gridded SOC data have been generated by linking the soil properties to soil maps or by inter-extrapolating the model outputs derived from analyses of observed SOC data points. Compared with the SOC distribution derived from ESMs, SOC estimates derived from SOC observations are more data oriented; however, even the observational databases include significant uncertainty because of errors in the source data and building processes (Köchy et al., 2015; Todd-Brown et al., 2013). A recent study (Todd-Brown et al., 2013) found that although ESM results are moderately consistent at the biome level, the correlation between the distribution of soil carbon stocks simulated by ESMs and observational databases is poor when the two are compared at fine scales (e.g., a 1◦ scale). Furthermore, estimates of SOC by ESMs and terrestrial biosphere models exhibit high uncertainty (Nishina et al., 2014, 2015; Tian et al., 2015). Several studies have examined the cause of the inconsistency in data derived by observational databases and ESMs, and the high variation of SOC outputs from ESMs (Exbrayat et al., 2013; ToddBrown et al., 2013; Wieder et al., 2013). Todd-Brown et al. (2013) analyzed the soil carbon outputs from 11 ESMs from the fifth phase of the Coupled Model Intercomparison Project (CMIP5) and soil carbon data from the HWSD, and found that net primary productivity (NPP) and temperature could explain the SOC spatial variations in the ESM output but not in the HWSD output. These authors also found that the differences in SOC from the ESMs were driven by differences in the simulated NPP and the parameterization of soil heterotrophic respiration and not by differences in the soil model structure of the ESMs. The key influence of parameterizing soil heterotrophic respiration (e.g., turnover time) on SOC in the CMIP5 ESM has also been discussed by Exbrayat et al. (2013). Anav et al. (2013) examined the relationships between simulated SOC and vegetation carbon and compared them with reference data obtained from observational databases, and they found that although simulated values clustered around the reference values, the ratio of SOC to vegetation carbon differed among the models. This finding suggests that the parameterization of plant production, mortality, and decomposition vary greatly among ESMs. More realistic representations of turnover times (Koven et al., 2015) and focusing on the model treatment of the hydrological cycle on the carbon cycle (Shao et al., 2013) are suggested as future improvements. Despite these research results, the key processes and factors that govern the SOC stock and the relationships among them remain unclear. The appropriate inclusion of these processes/factors would improve the consistency between the model results and observational data. In this study, we sought to identify the key factors that govern the global SOC distribution in observational databases as well as those simulated by ESMs. We applied a data-mining (machine-learning) scheme (boosted regression trees – BRT) to identify the influential factors and explore how they relate to SOC stocks (Elith et al., 2008). The BRT method is based on regression trees and boosting. We combined the potentially influential variables from many data products and SOC data from observational databases and ESMs, and examined the factors influencing the distribution of SOC and the relationships between these factors and SOC stocks. We assessed how closely ESMs could match the influential factors and their relationships with factors obtained from observational databases. By comparing the influential factors in the observational databases with those in the ESMs, we clarified the model–data discrepancies and the areas in which ESMs can be improved. 2 Materials and methods 2.1 Observational global SOC database We used SOC data from two global databases and one northern observational database. The first global database was the HWSD (FAO/IIASA/ISRIC/ISSCAS/JRC, 2012). The HWSD is a global database of soil physiochemical properties that has been developed by the International Institute for Applied Systems Analysis (IIASA) and the Food and Agriculture Organization of the United Nations (FAO) in collaboration with the International Soil Reference and Information Centre (ISRIC) World Soil Information, the European Commission Joint Research Centre (JRC), and the Institute of Soil Science, Chinese Academy of Sciences (ISSCAS). The database was constructed by compiling the European Soil Database (ESDB), a 1 :1 million soil map of China, various regional SOTER databases (SOTWIS database), and a soil map of the world from the FAO. We used an SOC stock database obtained with HWSD from the Joint Research Centre (JRC) (Hiederer and Köchy, 2011) (Fig. 1a). The second database included global gridded surfaces of selected soil characteristics (IGBP-DIS) (Global Soil Data Task Group, Geosci. Model Dev., 10, 1321–1337, 2017 www.geosci-model-dev.net/10/1321/2017/ S. Hashimoto et al.: Data-mining analysis of the global distribution of soil carbon 1323 60˚ S 0˚ 60˚ N (a) HWSD 60˚ S 0˚ 60˚ N (b) IGBP−DIS 180˚ 90˚ W 0˚ 90˚ E 180˚ 60˚ S 0˚ 60˚ N 0 20406080100 kg C m−2 (c) NCSCD Figure 1. Soil carbon stock in the upper 100 cm (kg C m−2)from the observational databases (HWSD, IGBP-DIS, and NCSCD). 2000) (Fig. 1b), which contains gridded soil physiochemical properties. The database has been developed by the Global Soil Data Task Group of the International Geosphere Biosphere Programme’s (IGBP) Data and Information System (DIS), and the database was generated by linking the pedon records in the Global Pedon Database to the FAO/UNESCO digital soil map of the world. The third database was the Northern Circumpolar Soil Carbon Database, version 2 (NCSCD) (Hugelius et al., 2013; Tarnocai et al., 2009) (Fig. 1c). This database is a spatial database of SOC stock of the northern circumpolar permafrost region. The soil map data were obtained from different regions/countries (USA, Canada, Russia, etc.) and were harmonized. The NCSCD were based on 1778 pedon data points. We used the HWSD and IGBP-DIS to analyze the global distribution of SOC stocks; then, we extracted a database of northern circumpolar regions from the three above databases and analyzed the SOC stocks in the northern region. The relationships among the databases are shown in Fig. S1 in the Supplement. The SOC in the upper 100 cm in each database was used. 2.2 Global SOC estimated using Earth system models The global distribution of SOC stocks estimated by ESMs was obtained from CMIP5. We examined the results of 15 ESMs (Fig. 2; Table 2). When more than one result was obtained by the same model family (e.g., MIROC-ESM and MIROC-ESM-CHEM), we generated an ensemble average database for each family (e.g., average of MIROC-ESM and MIROC-ESM-CHEM): Todd-Brown et al. (2013) showed through a hierarchical cluster analysis that SOC distributions were very similar among ESMs from the same climate center. The mean values from 1980–2004 were calculated. The results of the historical and ensemble member r1i1p1 were used in this study. The notation “r1i1p1” is an identifier of the model simulation and is an ensemble member that is often used for analyses (Chang et al., 2012; Dirmeyer et al., 2013; Jiang et al., 2015; Kumar et al., 2014). The overviews of SOC submodels in the ESMs have been previously described (Exbrayat et al., 2014; Todd-Brown et al., 2013, 2014) and are also shown in Table 2. In general, each soil submodel consisted of one to nine pools and incorporated the effects of temperature and moisture. Some ESMs have litter carbon pools; these were excluded from this study. A comparison between the mean of ESMs and global observational databases in a 1◦grid is shown in Fig. S2. 2.3 Other databases We used five groups of variables/factors to examine their effects on global SOC: climate, soil property, topography, vegetation, and land-use history. Detailed data sources for the databases are described in Table 1. The mean annual temperature and annual precipitation were used as the climate variables, and the clay content, carbon-to-nitrogen (CN) ratio, and texture (Appendix Table A1) were used as the soil variables (0–30 cm). The compound topographic index, elevation, slope, and wetland ratio were used as the topographic indices. The CN ratio was calculated by dividing the carbon density by the nitrogen density. The wetland ratio was calculated by dividing the number of wetland grids at 30 s by the total grids at 1◦. The lake, reservoir, and river were not quantified as wetlands and were excluded from the total grids. The land cover type (Appendix Table A2) and NPP were adopted as vegetation indices, and the cropland ratio and human appropriation of net primary production percentage, which is a percentage of human consumption of NPP to local NPP (Imhoff and Bounoua, 2006), were used as the indices of land-use history. The average human appropriation of the NPP percentage was calculated at 1◦. Histograms of the variables are shown in Fig. S3. www.geosci-model-dev.net/10/1321/2017/ Geosci. Model Dev., 10, 1321–1337, 2017 1324 S. Hashimoto et al.: Data-mining analysis of the global distribution of soil carbon 60˚ S 0˚ 60˚ N (a) BCC−ensemble (b) BNU−ESM (c) CanESM2 60˚ S 0˚ 60˚ N (d) CCSM4 (e) CESM1−ensemble (f) CMCC−CESM 60˚ S 0˚ 60˚ N (g) GFDL−ESM2M (h) GISS−ensemble (i) HadGEM2−CC 60˚ S 0˚ 60˚ N (j) INM−CM4 (k) IPSL−ensemble (l) MIROC−ensemble 180˚ 90˚ W 0˚ 90˚ E 180˚ 60˚ S 0˚ 60˚ N (m) MPI−ensemble 180˚ 90˚ W 0˚ 90˚ E 180˚ 0 20406080100 kg C m−2 (n) MRI−ESM1 180˚ 90˚ W 0˚ 90˚ E 180˚ (o) NorESM1−ensemble Figure 2. Soil carbon stocks (kg C m−2)from Earth system models (CMIP5). The term “ensemble” indicates the result of an ensemble of family members. 2.4 Database handling All global databases, except for the databases with a spatial resolution of 1◦by default, including observational and ESM model outputs, were regridded to a spatial resolution of 1◦for the analyses. Regridding of data in the NetCDF format was performed using the Climate Data Operators (CDO) software, version 1.6.9, provided by the Max Plank Institute for Meteorology (https://code.zmaw.de/projects/cdo). A bilinear interpolation, which is one of the most widely used algorithms, was used (remapbil in CDO). 2.5 Boosted regression trees (BRT) scheme To identify the influential factors and their relationships with SOC stocks, BRT scheme was used in this study (Elith et al., 2008). This technique involves a data-mining (machinelearning) algorithm that combines the advantages of a regression tree (decision tree) algorithm and boosting. Regression trees are a classification algorithm that classify data through recursive binary splits, and boosting is a machine-learning algorithm that generates many rough models and combines them to improve their predictive capability. The main advantages of this method are that the BRT scheme can analyze different types of variables and interaction effects among variables, and is applicable to nonlinear relationships. In recent years, the BRT technique has been used to examine the distribution of soil characteristics at a regional scale (Aertsen et al., 2011; Cools et al., 2014; Martin et al., 2011). Major outputs from BRT analyses can identify the following: (1) the relative importance (percentage of influence or contribution) of predictor variables (explanatory variables) on the basis of the weighted and scaled number of times a variable is selected for splitting (Elith et al., 2008) and (2) the relationships among variables and the explained variable shown in partial dependence plots. Geosci. Model Dev., 10, 1321–1337, 2017 www.geosci-model-dev.net/10/1321/2017/ S. Hashimoto et al.: Data-mining analysis of the global distribution of soil carbon 1325 Table 1. Variables used in the analyses and their sources. Variable Abbreviation Source (database) Original resolution Reference Mean annual temperature1MAT ISLSCPII (CRU05) 1◦New et al. (2011) Mean annual precipitation1MAP ISLSCPII (CRU05) 1◦New et al. (2011) Clay content (0–30 cm) Clay ISLSCPII 1◦Scholes and Brown de Colstoun (2011) CN ratio (0–30 cm)2CN ratio ISLSCPII 1◦Scholes and Brown de Colstoun (2011) Soil texture (0–30 cm) Texture ISLSCPII 1◦Scholes and Brown de Colstoun (2011) Compound topographic index3CTI ISLSCPII 1◦Verdin (2011) Elevation3Elev. ISLSCPII 1◦Verdin (2011) Slope3Slope ISLSCPII 1◦Verdin (2011) Wetland ratio Wetland Global Lakes and Wetlands Database 30 s Lehner and Döll (2004) Land cover LandCover ISLSCPII 1◦Friedl et al. (2010) Net primary production NPP ISLSCPII 1◦Prince and Zheng (2011) Cropland ratio Cropland ISLSCPII 1◦Ramankutty and Foley (2010) Human appropriation of NPP percentage HANPPpct HANPP collection 0.25◦Imhoff et al. (2004) 1The original database provides monthly data. Annual means were calculated by the authors. 2The CN ratio was calculated by dividing the carbon density by the nitrogen density. 3The native database is hydro1k, and its resolution is 1 km. The mean value of 1 km was used in this study. We used the open-source BRT package (brt.functions.R) in R software versions 3.2.1 and 3.2.2 (R Core team, 2013) developed by Elith et al. (2008). The R code for the BRT algorithm is available in the supplementary material of Elith et al. (2008). The gbm package was used (version 2.1.1) to run the BRT package. The calculations were performed in Mac OS X (version 10.9.5 and version 10.10.5). To do so, the “windows” function in the “brt.functions.R” needed to be replaced with the “quartz” function in R. In practice, three parameters in the BRT package – the learning rate (lr), tree complexity (tc), and bag fraction (bg) – control the BRT performance. The lr determines the contribution of each tree, the tc controls the number of splits, and the bg is the proportion of data selected at each step. The number of trees was determined using the cross-validation method in the R package. The maximum number of trees was set to 15 000. The tc value was set to 5. We tested different lr (0.001, 0.005, 0.01, 0.05, 0.1) and bg values (0.5, 0.6, 0.7) and used the best parameter set for each database, but the changes in parameter values had little effect on the model performance. 2.6 Model performance The goodness of fit between the BRT model and data was assessed by using the linear relationship between the predicted and observed values, the coefficient of determination (R2), and the root mean square error (RMSE); it is shown in Tables S2 and S3 in the Supplement. For both the observational databases and ESM databases, the BRT models exhibited good performance, with high R2values in most of the databases, but the performance was relatively lower for NCSCD and CMCC (northern soils). 3 Results 3.1 Observational databases 3.1.1 Global soil The relative contributions of variables in the BRT model of global SOC stocks to the observational databases are shown in Fig. 3a and b. In HWSD, the contributions of land cover, mean annual temperature, CN ratio, and wetland ratio were high. For IGBP-DIS, the mean annual temperature, followed by clay content, CN ratio, and land cover also highly contributed. In particular, the mean annual temperature was very influential. The contribution of elevation to each HWSD and IGBP-DIS was 6 and 7 %, respectively. The NPP contributed 5 % in both databases. The relationships between the influential variables and SOC are shown in Fig. 4a–e. In general, the two databases showed similar relationships. For example, the SOC decreased with increasing mean annual temperature, particularly at sites with a mean annual temperature >0◦C (Fig. 4a), but increased with increasing clay content and CN www.geosci-model-dev.net/10/1321/2017/ Geosci. Model Dev., 10, 1321–1337, 2017 1326 S. Hashimoto et al.: Data-mining analysis of the global distribution of soil carbon 0 10 20 30 40 50 (a) Global, HWSD Contribution (%) 0 10 20 30 40 50 MAT MAP Clay CNratio Texture CTI Elev Slope Wetland LandCover NPP Cropland HANPPpct (b) Global, IGBP-DIS Contribution (%) Variables 0 10 20 30 40 50 (c) North, HWSD Contribution (%) 0 10 20 30 40 50 (d) North, IGBP-DIS Contribution (%) 0 10 20 30 40 50 MAT MAP Clay CNratio Texture CTI Elev Slope Wetland LandCover NPP Cropland HANPPpct (e) North, NCSCD Contribution (%) Variables Figure 3. Relative contribution (influence) of predictive variables for the model of soil carbon stocks in the global observational databases (left) and northern observational databases (right). ratio (Fig. 4b and c). The SOC increased rapidly with an increasing CN ratio. Relationships with the mean annual temperature were similar (Fig. 4a). The relationship with clay was steeper in IGBP-DIS than in HWSD, but the opposite was true for the CN ratio (Fig. 4b and c). With respect to land cover, evergreen needleleaf forests and permanent wetlands had higher SOC (Fig. 4e). 3.1.2 Northern soils In the northern region, the dominant contributors differed among northern soil databases and from those identified in the global database analyses described above (Fig. 3c–e). In HWSD, the CN ratio was the dominant contributor, followed by the wetland ratio, clay content, and mean annual precipitation. In IGBP-DIS, clay content, CN ratio, and elevation were the most important contributors. For NCSCD, elevation contributed the most (∼25 %), but all of the variables except for the cropland ratio and HANPPpct contributed 5–15 %. The mean annual temperature was not as influential as the global databases. The relationships between variables and SOC stock varied more among the databases for northern soils than those of global databases (Fig. 4f–k). Furthermore, because the northern regions were extracted, the ranges of variables were narrower than the global databases. In NCSCD, the SOC decreased with increasing temperature (Fig. 4f) and increased with increasing precipitation (Fig. 4g). The SOC increased with increasing clay content and CN ratio in HWSD and IGBP-DIS (Fig. 4h and i), which was consistent with the findings obtained from the global databases. The increasing trend with increasing CN ratio was also observed in NCSCD. The SOC decreased with increasing elevation in all databases but showed considerable variability at low elevations (Fig. 4j). 3.2 Earth system models 3.2.1 Global soil The contributions of some variables varied among ESMs, but the mean of the results of the ESMs showed that the mean annual temperature, land cover, and NPP clearly contributed to SOC distribution (Fig. 5a and b). Large inconsistencies between the observational databases and ESMs were found in the low contributions of clay content and the CN ratio and in the high contributions of NPP in ESMs (Fig. 5a and b). The contribution of NPP to ESMs was greater than in the observational databases. The relationships between SOC and certain variables substantially varied among the ESM databases (Fig. 6a–e), particularly in the mean annual temperature (Fig. 6a). The SOC decreased with increasing mean annual temperature (Fig. 6a) but increased with increasing precipitation (Fig. 6b) and NPP (Fig. 6e). The mean of the relationship with mean annual temperature for ESMs was highly consistent with that in the Geosci. Model Dev., 10, 1321–1337, 2017 www.geosci-model-dev.net/10/1321/2017/ S. Hashimoto et al.: Data-mining analysis of the global distribution of soil carbon 1327 -10 0 10 -20 -10 0 10 20 30 (a) Fitted function (kg C m− 2 ) MAT ( C) -10 0 10 0 10 20 30 40 50 (b) Clay (%) -20 -10 0 10 20 0 10 20 30 40 50 (c) C / N ratio -10 0 10 0 0.5 1 (d) Fitted function (kg C m− 2 ) Wetland ratio -10 0 10 0 5 10 15 (e) Land cover HWSD IGBP-DIS -10 0 10 -20 -10 0 10 20 30 (f) Fitted function (kg C m− 2 ) MAT ( C) -10 0 10 0 1000 2000 3000 4000 (g) MAP (mm) -20 -10 0 10 20 0 10 20 30 40 50 (h) Clay (%) -20 -10 0 10 20 0 10 20 30 40 50 (i) Fitted function (kg C m− 2 ) C / N ratio -10 0 10 0 1000 2000 3000 (j) Elev (m) -10 0 10 0 0.5 1 (k) Wetland ratio HWSD IGBP-DIS NCSCD o o Figure 4. Effects of the most influential variables in the model of the soil carbon stock for each global (a–e) and northern (f–k) observational databases. The fitted functions were centered by subtracting their means. See Table A2 for land cover classifications. Because of the small number of data points, the results for “permanent snow and ice” are not shown (e). The y-axis scales for clay and the CN ratio are different from those of other factors (c, h, i). HWSD and IGBP-DIS databases of the temperature range −5 to 15 ◦C (Fig. 6a). The increasing trend with increasing NPP in ESMs was consistent with that of the HWSD, particularly below approximately 500 g C m−2of NPP (Fig. 6e). Although the wetland ratio did not contribute to the ESMs (Fig. 6a) with respect to land cover, permanent wetlands had higher SOC (Fig. 6d). 3.2.2 Northern soils The mean of the ESMs showed that for northern soils, the main contributors (mean annual temperature, land cover, and NPP) were mainly the same as in the ESM global outputs (Fig. 5c and d). The contribution of the mean annual temperature was lower than that of the global results of the ESMs (mean of 14 % for the northern and 29 % for the global temperatures). The relatively large discrepancy between the observational databases and ESMs included the lower contribution of clay content, CN ratio, and elevation, and the higher contribution of the mean annual temperature, land cover, and NPP in the ESMs. The relationship between SOC and variables in ESMs as well as the results of the observational databases are shown in Fig. 6f–i. The mean of the ESMs indicated that the SOC in the northern region increased with increasing NPP, and the relationship was similar to that in HWSD (Fig. 6i), although the contribution of NPP in the ESMs differed from www.geosci-model-dev.net/10/1321/2017/ Geosci. Model Dev., 10, 1321–1337, 2017 1328 S. Hashimoto et al.: Data-mining analysis of the global distribution of soil carbon (b) (d) % % 0 10 20 30 40 50 60 70 80 MAT MAP Clay CNratio Texture CTI Elev Slope Wetland LandCover NPP Cropland HANPPpct (a) Global Contribution (%) ESM mean HWSD IGBP-DIS 0 10 20 30 40 50 60 70 80 MAT MAP Clay CNratio Texture CTI Elev Slope Wetland LandCover NPP Cropland HANPPpct (c) North Contribution (%) Variables ESM mean HWSD IGBP-DIS NCSCD Databases NorESM1−ensemble MRI−ESM1 MPI−ensemble MIROC−ensemble IPSL−ensemble INM−CM4 HadGEM2−CC GISS−ensemble GFDL−ESM2M CMCC−CESM CESM1−ensemble CCSM4 CanESM2 BNU−ESM BCC−ensemble ESM−mean IGBP−DIS HWSD MAT MAP Clay CNratio Texture CTI Elev Slope Wetland LandCover NPP Cropland HANPPpct 0 20 40 60 80 100 Variables Databases NorESM1−ensemble MRI−ESM1 MPI−ensemble MIROC−ensemble IPSL−ensemble INM−CM4 HadGEM2−CC GISS−ensemble GFDL−ESM2M CMCC−CESM CESM1−ensemble CCSM4 CanESM2 BNU−ESM BCC−ensemble ESM−mean NCSCD IGBP−DIS HWSD MAT MAP Clay CNratio Texture CTI Elev Slope Wetland LandCover NPP Cropland HANPPpct 0 20 40 60 80 100 Figure 5. Relative contribution (influence) of predictive variables for the model of the soil carbon stock from ESMs and a comparison with those of observational databases. Box plots show the results of ESMs, and the purple, green, light blue, and blue marks indicate the mean of the ESMs and results from observational databases (a: global; c: north). Mosaic plots of detailed relative contributions for each ESM (b: global; d: north) are shown. those of the observational database (Fig. 5c). The decreasing trend with elevation was not replicated in the ESMs (Fig. 6g). 4 Discussion and concluding remarks 4.1 Identified influential factors Compared with previous studies, we examined the contributions of a wider variety of factors to SOC distributions. Our analyses revealed that the most distinct differences between the observational database data and the ESM outputs were the effects of the CN ratio and clay content (Fig. 5). For both global observational databases, the CN ratio was a substantial contributor (Fig. 3a and b). The important contribution of the CN ratio was the same in the northern databases (Fig. 3c–e). The SOC in the observational databases increased with increases in the CN ratio (Fig. 4c), whereas the SOC values of the ESMs were insensitive to the CN ratio. Our results support the importance of properly incorporating the nitrogen (N) cycle into SOC models (e.g., control over decomposition, soil fertility, nutrient availability, and plant litter quality) (Berg et al., 2001; Cotrufo et al., 2013; Fernández-Martínez et al., 2014; Liski et al., 2005; Tuomi et al., 2009; ˇ Tupek et al., 2016). None of the ESMs except for the CESM1 and NorESM in CMIP5 included terrestrial nitrogen processes (Todd-Brown et al., 2013); however, including this parameter has been suggested as a key improvement for the next model intercomparison (CMIP6) (Hajima et al., 2014; Zaehle et al., 2015). The results of our analysis support the importance of including the N cycle in ESM models. Clay content is also often used as a regulator of the decomposability of organic matter in the soil (e.g., CENTURY and RothC) (Coleman and Jenkinson, 1999; Parton et al., 1987). Generally, high clay content inhibits organic matter decomposition in the soil. Furthermore, high clay contents often result in low drainage and anaerobic soil conditions, which also inhibit organic matter decomposition. For the IGBPDIS data, the contribution of the clay content was as high as that of the CN ratio. The control of decomposability by the Geosci. Model Dev., 10, 1321–1337, 2017 www.geosci-model-dev.net/10/1321/2017/ S. Hashimoto et al.: Data-mining analysis of the global distribution of soil carbon 1329 -10 0 10 -20 -10 0 10 20 30 (a) Fitted function (kg C m− 2 ) MAT ( C) -10 0 10 0 1000 2000 3000 4000 (b) MAP (mm) -10 0 10 0 1000 2000 3000 (c) Elev (m) -10 0 10 0 5 10 15 (d) Fitted function (kg C m− 2 ) Land cover -10 0 10 0 500 1000 1500 (e) NPP (g C m−2) ESM mean HWSD IGBP-DIS -10 0 10 -20 -10 0 10 20 30 (f) Fitted function (kg C m− 2 ) MAT ( C) -10 0 10 0 1000 2000 3000 (g) Elev (m) -10 0 10 0 5 10 15 (h) Land cover -10 0 10 0 500 1000 1500 (i) Fitted function (kg C m− 2 ) NPP (g C m−2) ESM mean HWSD IGBP-DIS NCSCD o o Figure 6. Effect of the most influential variables in the model for global (a–e) and northern (f–i) outputs from ESMs and a comparison with those of observational databases. Grey lines show the results of each ESM, and the purple line indicates the mean of the ESMs. The fitted functions were centered by subtracting their means. See Table A2 for land cover classifications. Because of the small number of data points, the results for “permanent snow and ice” are not shown (d, h). clay content has been previously incorporated in site-scale process-based models (Parton et al., 1987) and may be incorporated in certain ESMs because the soil carbon submodels in these ESMs are based on the CENTURY model (see the soil model history reported in Todd-Brown et al., 2014). However, regardless of whether the control of decomposability by clay is incorporated, our results suggest that the influence of clay on the carbon cycle is not well captured in most ESMs. The mean annual temperature was identified as an influential factor in the global databases (Fig. 3a and b) but not in the northern soil databases (Fig. 3c–e). Temperature is a main factor controlling both plant production (source of carbon input to soil) and soil organic matter decomposition, which are already incorporated in ESMs. Based on an analysis of the output of heterotrophic respiration, the temperature sensitivity (e.g., Q10 value) of soil organic matter decomposition in the ESMs has been reported as 1.4 to 2.2 (Todd-Brown et al., 2014). In addition, our analyses identified diverse relationships between the mean annual temperature and SOC. The lower contribution of the mean annual temperature in the northern soils likely occurred because temperature sensitivity is an exponential process, and the magnitude of observed changes under changing temperature is relatively small at a low temperature range, as shown in a comparison of temperature functions in common biogeochemical models (Sierra et al., 2015). The relationships between the SOC and temperature obtained in this study include the integration of the temperature sensitivity of both plant production and soil organic decomposition and thus do not provide the temperature sensitivity parameter of individual processes for ESMs. However, the results of this study can be used to examine www.geosci-model-dev.net/10/1321/2017/ Geosci. Model Dev., 10, 1321–1337, 2017 1336 S. Hashimoto et al.: Data-mining analysis of the global distribution of soil carbon Liski, J., Palosuo, T., Peltoniemi, M., and Sievänen, R.: Carbon and decomposition model Yasso for forest soils, Ecol. Model., 189, 168–182, doi:10.1016/j.ecolmodel.2005.03.005, 2005. Luo, Y., Ahlström, A., Allison, S. D., Batjes, N. H., Brovkin, V., Carvalhais, N., Chappell, A., Ciais, P., Davidson, E. A., Finzi, A., Georgiou, K., Guenet, B., Hararuk, O., Harden, J. W., He, Y., Hopkins, F., Jiang, L., Koven, C., Jackson, R. B., Jones, C. D., Lara, M. J., Liang, J., McGuire, A. D., Parton, W., Peng, C., Randerson, J. T., Salazar, A., Sierra, C. A., Smith, M. J., Tian, H., Todd-Brown, K. E. O., Torn, M., van Groenigen, K. J., Wang, Y. P., West, T. O., Wei, Y., Wieder, W. R., Xia, J., Xu, X., Xu, X., and Zhou, T.: Toward more realistic projections of soil carbon dynamics by Earth system models, Global Biogeochem. Cy., 30, 40–56, doi:10.1002/2015GB005239, 2016. Manzoni, S. and Porporato, A.: Soil carbon and nitrogen mineralization: Theory and models across scales, Soil Biol. Biochem., 41, 1355–1379, doi:10.1016/j.soilbio.2009.02.031, 2009. Martin, M. P., Wattenbach, M., Smith, P., Meersmans, J., Jolivet, C., Boulonne, L., and Arrouays, D.: Spatial distribution of soil organic carbon stocks in France, Biogeosciences, 8, 1053–1065, doi:10.5194/bg-8-1053-2011, 2011. New, M., Jones, P. D., and Hulme, M.: ISLSCP II Climate Research Unit CRU05 Monthly Climate Data, doi:10.3334/ORNLDAAC/1015, 2011. Nishina, K., Ito, A., Beerling, D. J., Cadule, P., Ciais, P., Clark, D. B., Falloon, P., Friend, A. D., Kahana, R., Kato, E., Keribin, R., Lucht, W., Lomas, M., Rademacher, T. T., Pavlick, R., Schaphoff, S., Vuichard, N., Warszawaski, L., and Yokohata, T.: Quantifying uncertainties in soil carbon responses to changes in global mean temperature and precipitation, Earth Syst. Dynam., 5, 197–209, doi:10.5194/esd-5-197-2014, 2014. Nishina, K., Ito, A., Falloon, P., Friend, A. D., Beerling, D. J., Ciais, P., Clark, D. B., Kahana, R., Kato, E., Lucht, W., Lomas, M., Pavlick, R., Schaphoff, S., Warszawaski, L., and Yokohata, T.: Decomposing uncertainties in the future terrestrial carbon budget associated with emission scenarios, climate projections, and ecosystem simulations using the ISI-MIP results, Earth Syst. Dynam., 6, 435–445, doi:10.5194/esd-6-435-2015, 2015. Ostle, N. J., Smith, P., Fisher, R., Woodward, F. I., Fisher, J. B., Smith, J. U., Galbraith, D., Levy, P., Meir, P., McNamara, N. P., and Bardgett, R. D.: Integrating plant-soil interactions into global carbon cycle models, J. Ecol., 97, 851–863, 2009. 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, Soil Sci. Soc. Am. J., 51, 1173–1179, 1987. Prince, S. D. and Zheng, D. L.: ISLSCP II global primary production data initiative gridded NPP data, doi:10.3334/ORNLDAAC/1023, 2011. Ramankutty, N. and Foley, J. A.: ISLSCP II historical croplands cover, 1700–1992, doi:10.3334/ORNLDAAC/966, 2010. R Core team: R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, 2013. Scharlemann, J. P., Tanner, E. V., Hiederer, R., and Kapos, V.: Global soil carbon: understanding and managing the largest terrestrial carbon pool, Cabon Manag., 5, 81–91, doi:10.4155/cmt.13.77, 2014. Schimel, D. S., Braswell, B. H., Holland, E. a., McKeown, R., Ojima, D. S., Painter, T. H., Parton, W. J., and Townsend, A. R.: Climatic, edaphic, and biotic controls over storage and turnover of carbon in soils, Global Biogeochem. Cy., 8, 279–293, doi:10.1029/94GB00993, 1994. Scholes, E. and Brown de Colstoun, E.: ISLSCP II global gridded soil characteristics, doi:10.3334/ORNLDAAC/1004, 2011. Shao, P., Zeng, X., Sakaguchi, K., Monson, R. K., and Zeng, X.: Terrestrial carbon cycle: Climate relations in eight CMIP5 earth system models, J. Climate, 26, 8744–8764, doi:10.1175/JCLI-D12-00831.1, 2013. Sierra, C. A. and Müller, M.: A general mathematical framework for representing soil organic matter dynamics, Ecol. Monogr., 85, 505–524, doi:10.1890/15-0361.1, 2015. Sierra, C. A., Trumbore, S. E., Davidson, E. A., Vicca, S., and Janssens, I.: Sensitivity of decomposition rates of soil organic matter with respect to simultaneous changes in temperature and moisture, J. Adv. Model. Earth Syst., 7, 335–356, doi:10.1002/2014MS000358, 2015. Tarnocai, C., Canadell, J. G., Schuur, E. A. G., Kuhry, P., Mazhitova, G. and Zimov, S.: Soil organic carbon pools in the northern circumpolar permafrost region, Global Biogeochem. Cy., 23, GB2023, doi:10.1029/2008GB003327, 2009. Tian, H., Lu, C., Yang, J., Banger, K., Huntzinger, D. N., Schwalm, C. R., Michalak, A. M., Cook, R., Ciais, P., Hayes, D., Huang, M., Ito, A., Jain, A. K., Lei, H., Mao, J., Pan, S., Post, W. M., Peng, S., Poulter, B., Ren, W., Ricciuto, D., Schaefer, K., Shi, X., Tao, B., Wang, W., Wei, Y., Yang, Q., Zhang, B., and Zeng, N.: Global patterns and controls of soil organic carbon dynamics as simulated by multiple terrestrial biosphere models: Current status and future directions, Global Biogeochem. Cy., 29, 775– 792, doi:10.1002/2014GB005021, 2015. Todd-Brown, K. E. O., Randerson, J. T., Hopkins, F., Arora, V., Hajima, 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 century, Biogeosciences, 11, 2341–2356, doi:10.5194/bg-11-23412014, 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. Tuomi, M., Thum, T., Järvinen, H., Fronzek, S., Berg, B., Harmon, M., Trofymow, J. A., Sevanto, S., and Liski, J.: Leaf litter decomposition – Estimates of global variability based on Yasso07 model, Ecol. Model., 220, 3362–3371, doi:10.1016/j.ecolmodel.2009.05.016, 2009. ˇ Tupek, B., Ortiz, C. A., Hashimoto, S., Stendahl, J., Dahlgren, J., Karltun, E., and Lehtonen, A.: Underestimation of boreal soil carbon stocks by mathematical soil carbon models linked to soil nutrient status, Biogeosciences, 13, 4439–4459, doi:10.5194/bg13-4439-2016, 2016. Verdin, K. L.: ISLSCP II HYDRO1k Elevation-derived Products, doi:10.3334/ORNLDAAC/1007, 2011. Wang, Y. P., Chen, B. C., Wieder, W. R., Leite, M., Medlyn, B. E., Rasmussen, M., Smith, M. J., Agusto, F. B., Hoffman, F., and Luo, Y. Q.: Oscillatory behavior of two nonlinear microbial models of soil carbon decomposition, Biogeosciences, 11, 1817– 1831, doi:10.5194/bg-11-1817-2014, 2014. Wang, Y. P., Jiang, J., Chen-Charpentier, B., Agusto, F. B., Hastings, A., Hoffman, F., Rasmussen, M., Smith, M. J., Todd-Brown, K., Geosci. Model Dev., 10, 1321–1337, 2017 www.geosci-model-dev.net/10/1321/2017/ S. Hashimoto et al.: Data-mining analysis of the global distribution of soil carbon 1337 Wang, Y., Xu, X., and Luo, Y. Q.: Responses of two nonlinear microbial models to warming and increased carbon input, Biogeosciences, 13, 887–902, doi:10.5194/bg-13-887-2016, 2016. Wieder, W. R., Allison, S. D., Davidson, E. A., Georgiou, K., Hararuk, O., He, Y., Hopkins, F., Luo, Y., Smith, M. J., Sulman, B., Todd-Brown, K., Wang, Y.-P., Xia, J., and Xu, X.: Explicitly representing soil microbial processes in Earth system models, Global Biogeochem. Cy., 29, 1782–1800, doi:10.1002/2015GB005188, 2015. Wieder, W. R., Bonan, G. B., and Allison, S. D.: Global soil carbon projections are improved by modelling microbial processes, Nat. Clim. Change, 3, 909–912, doi:10.1038/nclimate1951, 2013. Wieder, W. R., Boehnert, J., and Bonan, G. B.: Evaluating soil biogeochemistry parameterizations in Earth system models with observations, Global Biogeochem. Cy., 28, 211–222, doi:10.1002/2013GB004665, 2014. Wutzler, T. and Reichstein, M.: Soils apart from equilibrium – consequences for soil carbon balance modelling, Biogeosciences, 4, 125–136, doi:10.5194/bg-4-125-2007, 2007. Zaehle, S.: Terrestrial nitrogen-carbon cycle interactions at the global scale, Philos. T. Roy. Soc. B, 368, 20130125, doi:10.1098/rstb.2013.0125, 2013. Zaehle, S., Jones, C. D., Houlton, B., Lamarque, J.-F., and Robertson, E.: Nitrogen availability reduces CMIP5 projections of twenty-first-century land carbon uptake, J. Climate, 28, 2494– 2511, doi:10.1175/JCLI-D-13-00776.1, 2015. www.geosci-model-dev.net/10/1321/2017/ Geosci. Model Dev., 10, 1321–1337, 2017