Global spatiotemporal distribution of soil respiration modeled using a global database

The flux of carbon dioxide from the soil to the atmosphere (soil respiration) is one of the major fluxes in the global carbon cycle. At present, the accumulated field observation data cover a wide range of geographical locations and climate conditions. However, there are still large uncertainties in the magnitude and spatiotemporal variation of global soil respiration. Using a global soil respiration data set, we developed a climate-driven model of soil respiration by modifying and updating Raich’s model, and the global spatiotemporal distribution of soil respiration was examined using this model. The model was applied at a spatial resolution of 0.5and a monthly time step. Soil respiration was divided into the heterotrophic and autotrophic components of respiration using an empirical model. The estimated mean annual global soil respiration was 91 Pg C yr (between 1965 and 2012; Monte Carlo 95 % confidence interval: 87– 95 Pg C yr) and increased at the rate of 0.09 Pg C yr. The contribution of soil respiration from boreal regions to the total increase in global soil respiration was on the same order of magnitude as that of tropical and temperate regions, despite a lower absolute magnitude of soil respiration in boreal regions. The estimated annual global heterotrophic respiration and global autotrophic respiration were 51 and 40 Pg C yr, respectively. The global soil respiration responded to the increase in air temperature at the rate of 3.3 Pg C yr C, andQ10 = 1.4. Our study scaled up observed soil respiration values from field measurements to estimate global soil respiration and provide a data-oriented estimate of global soil respiration. The estimates are based on a semi-empirical model parameterized with over one thousand data points. Our analysis indicates that the climate controls on soil respiration may translate into an increasing trend in global soil respiration and our analysis emphasizes the relevance of the soil carbon flux from soil to the atmosphere in response to climate change. Further approaches should additionally focus on climate controls in soil respiration in combination with changes in vegetation dynamics and soil carbon stocks, along with their effects on the long temporal dynamics of soil respiration. We expect that these spatiotemporal estimates will provide a benchmark for future studies and also help to constrain process-oriented models.


Introduction
The carbon balance of terrestrial ecosystems is the result of the balance between carbon uptake by plants and carbon loss by plant and soil respiration (Beer et al., 2010;Luyssaert et al., 2007;Malhi et al., 1999;Le Quéré et al., 2009, 2014;Trumbore, 2006).The value of the balance, i.e., whether terrestrial ecosystems act as sinks or sources of carbon, has been a subject of considerable interest for studies of climate change.Accurate evaluations of each sink/source component and their response to environmental factors are thus essential for understanding future changes in the terrestrial carbon balance.
Published by Copernicus Publications on behalf of the European Geosciences Union.
The carbon dioxide (CO 2 ) flux from the soil to the atmosphere (called soil respiration, R S ) is one of the major fluxes in the global carbon cycle (Schlesinger and Andrews, 2000).R S primarily consists of heterotrophic respiration (soil organic matter decomposition) and autotrophic respiration (root respiration) (Bond-Lamberty et al., 2004;Hanson et al., 2000).R S is the main contributor to the total ecosystem respiration (Malhi et al., 1999); hence, R S plays a role in determining the carbon balance of terrestrial ecosystems.R S is sensitive to environmental factors (e.g., temperature and precipitation) (Davidson et al., 1998;Hashimoto et al., 2011b;Raich and Schlesinger, 1992;Raich et al., 2002;Schlesinger and Andrews, 2000), and future climate change is expected to increase the rate of R S at the global scale (Bond-Lamberty and Thomson, 2010b;Hashimoto et al., 2011a;Raich et al., 2002).Even a small change in the global R S rate will have a strong impact on the terrestrial carbon cycle and may accelerate the increase in atmospheric CO 2 concentration (IPCC, 2001(IPCC, , 2007)).
Observations of R S have a long history; in particular, the amount of collected field data increased rapidly in the 1990s, and there are now thousands of records of observed data (Bond-Lamberty and Thomson, 2010a;Chen et al., 2014).Recently, a global data set of observed R S was established by collecting data from available studies published in peerreviewed scientific literature (Bond-Lamberty and Thomson, 2010a).The use of the accumulated data for field observations will improve the estimation of global R S .
The number of global estimates of R S is, however, quite limited compared to the estimates of other terrestrial carbon fluxes (e.g., gross and net primary production (GPP and NPP, respectively)), or other greenhouse gas fluxes (e.g., methane and nitrous oxide).For instance, based on a literature survey (Ito, 2011), there are at least 251 estimates of global NPP.For R S , to the best of our knowledge, there are only six global estimates, ranging from 68 to 98 Pg C yr −1 and, thus, characterized by a large uncertainty (Hashimoto, 2012).Another example that indicates the large uncertainty in estimating R S are the large variations in estimates of soil carbon stocks and heterotrophic respiration simulated by the stateof-the-art Earth system models of the Coupled Model Intercomparison Project Phase 5 (CMIP5; http://cmip-pcmdi.llnl.gov/cmip5/) (Exbrayat et al., 2014;Todd-Brown et al., 2013), which is a model intercomparison project that provides scientific knowledge to the Intergovernmental Panel on Climate Change (Taylor et al., 2012).These facts suggest that further efforts should attempt to constrain the estimate of global R S by employing a model-data integration approach and field measurements.
The purpose of this study was to provide a new global estimate of the spatiotemporal distribution of R S based on the available observational data sets.Using a global R S data set (Bond-Lamberty and Thomson, 2010a), we developed a semi-empirical climate-driven model for R S .The temperature and precipitation functions of Raich's model were re-fined, and the parameters of the model were newly determined using over one thousand data points.We explored the spatiotemporal distribution of R S and examined the temperature sensitivity of R S .We further divided the estimated global R S into the heterotrophic and autotrophic components of R S using an empirical relationship between R S and heterotrophic respiration, and examined the distribution of each type of respiration.

Models
We developed a climate-driven model by updating Raich's model (Raich and Potter, 1995;Raich et al., 2002).The original model, Eq. ( 1), simulates R S as a function of temperature and water (precipitation), and its sensitivities are parameterized using three constants (F , K, and Q).The model is applied at a monthly time step and requires the monthly mean air temperature (T , • C) and monthly precipitation (P , cm): where moR S is the mean monthly R S (g C m −2 d −1 ), and F , and K (cm mol −1 ) are the parameters.The advantage of this model is its simplicity.Although there are numerous factors that affect R S (Chen et al., 2014), it is often recognized that temperature and precipitation are the two best predictors to represent the spatiotemporal variability of R S (Bond-Lamberty and Thomson, 2010b).In this study, the temperature and water functions of the original model were modified as follows: where F (g C m −2 d −1 ) and K (cm mol −1 ) are the parameters; a ( • C −1 ) and b ( • C −2 ) are the parameters for the temperature function; and α is the parameter for the precipitation function; P m−1 (cm) is the precipitation of the previous month.
First, we introduced a more flexible temperature function that has been reported to behave better than the simple exponential temperature function (Tuomi et al., 2008).This function peaks at T = a (2b) −1 , and the function can take either a convex shape or a simple exponential-like shape depending on the parameters a and b.The simple exponential temperature function has been widely applied to the modeling of the temperature sensitivity of R S , but a limitation is often pointed out in that the Q 10 value (the factor by which the respiration rate increases for a temperature interval of 10 • C) of the exponential function does not change across temperatures, while analyses of observed temperature sensitivities of R S suggest that the Q 10 value decreases with an increase in temperature (Kirschbaum, 1995;Tuomi et al., 2008) (but see Mahecha et al., 2010).The Q 10 value of the new temperature function can change across temperature ranges.
Second, we adopted the weighted average of the precipitation of both the current month and the previous month, instead of only using the precipitation of the current month.One of the limitations of the precipitation function of the original model was the so-called "zero-precipitation-zerorespiration" problem (Reichstein et al., 2003).In the original precipitation function, R S becomes zero when precipitation of the present month is zero; however, although zero precipitation occurs at times, even in temperate regions, R S can be maintained.However, this assumption of R S is reasonable when we focus on a global-scale evaluation and distinguish between very dry regions, such as deserts, and other regions.Including a soil water submodel to simulate the soil water conditions would be one solution, but we used a weighted average of precipitation here to avoid model complexity.By modifying the temperature and precipitation functions, the model has an increased flexibility, and global parameters for the model were estimated.

Data set
The R S data used in this study were obtained from the SRDB database (Bond-Lamberty and Thomson, 2010a) (version 3).The database covers a wide range of geographical and climatic spaces (Fig. S1 in the Supplement), although the availability is limited for certain regions (i.e., with low temperature, and with high temperature and low precipitation).For the purpose of modeling, the data from non-experimentally manipulated, non-agricultural ecosystems that had been measured using an infrared gas analyzer or gas chromatograph were extracted.The data with quality check flags, except for Q 01 , Q 02 , and Q 03 , were excluded.We further extracted the data with the location information (latitude and longitude) to support their combination with the monthly climate data from the global climate data set.Annual R S in the SRDB was used for data-model synthesis.Some of the data points in the SRDB are based on multi-year observations; however, the data were not weighted in this study.Each data point has the information of the year in which the study was performed or the middle year if the observation was conducted in multiple years, and we assumed that the data were obtained in a year of observation (or in the middle year if multiple years) and linked to the climate data.For each data point, we ran the model using a monthly time step and calculated the annual R S .The air temperature and precipitation were derived from the CRU3.21 climate data (University of East Anglia Climatic Research Unit (CRU), 2013).The spatial resolution of the climate data is 0.5 • .Using the latitude and longitude information and the year of observation in the SRDB, we extracted the monthly climate data from the climate data set.The number of data points used for model parameterization was 1638.
We examined other models that included leaf area index (LAI) and GPP for the parameterization of F in Eq. (1) (Mahecha et al., 2010;Migliavacca et al., 2011;Reichstein et al., 2003) (Table S1 in the Supplement).The model with LAI and GPP was characterized by a higher R 2 value than the simple climate-driven model (Table S1), which supports the hypothesis that vegetation substantially influences the variation in R S (Migliavacca et al., 2011;Reichstein et al., 2003;Wang et al., 2010).However, the number of data points in the database with LAI and GPP were limited, and including LAI and GPP resulted in losses of over 70 and 90 % of the data points, respectively (Table S1) (e.g., Bond-Lamberty and Thomson, 2010b).For the purpose of providing global estimates based on the accumulated observed data, we placed a higher value on relatively larger data points that cover wider geographical and climatic spaces rather than building additional mechanistic models.Hence, the climate-driven model, described above, was adopted for the estimation of global R S in this study.Similar to previous studies, the impact of landuse change was not included in this study.

Parameterization
We used a Bayesian calibration scheme to determine the best parameter sets and their uncertainty (Bates and Campbell, 2001;Hashimoto et al., 2011b;Van Oijen et al., 2005;Ricciuto et al., 2008;Zobitz et al., 2008).We assumed a uniform distribution for the a priori distribution of every parameter (F , a, b, K, α) and assumed a Gaussian modeldata error (standard deviation: σ ).To generate the a posteriori distribution, we performed a Markov Chain Monte Carlo simulation (MCMC) based on the Metropolis-Hastings (M-H) algorithm; the log-likelihood was used in practice.The MCMC program was coded in C, and the statistical analyses of the output were conducted using R versions 3.0.2and 3.1.0(R Core team, 2013).We conducted 100 000 iterations of sampling and discarded the first 20 000 iterations as the burn-in period.The maximum a posteriori estimates (MAP) were designated as the best-fit parameters.Geweke's Z-score was calculated for convergence diagnostics (Geweke, 1992); a Geweke's Z-score range of ± 1.96 indicates convergence (significance level of 5 %).

Global application
The R S was evaluated at a spatial resolution of 0.5 • and a monthly time resolution.The air temperature and precipitation were derived from the CRU 3.21 climate data (University of East Anglia Climatic Research Unit (CRU), 2013).The global land-use data in SYNMAP (Jung et al., 2006) were converted to 0.5 • for use in this study.We calculated R S for the period from 1965 to 2012.A Monte Carlo simulation was applied to evaluate the uncertainty of the estimates; the model was run 1000 times using the parameter uncertainties derived from the Bayesian calibration.(Raich and Potter, 1995;Raich et al., 2002).The grey area is the 95 % confidence interval of the estimated functions.
Table 1.A priori and a posteriori probability distributions of the parameters.The MAP is the maximum a posteriori estimate.A uniform distribution was assumed for every a priori distribution.CI indicates the confidence interval.F , a, b, K, and α are the model parameters, and σ is the standard deviation of the model-data error.

Parameters
Prior

Partitioning the total R S into the heterotrophic (R H ) and autotrophic (R A ) respiration components
The estimated annual R S was divided into heterotrophic respiration (R H ) and autotrophic respiration (R A ) using an empirical relationship derived by a meta-analysis (Bond-Lamberty et al., 2004).From that meta-analysis, a global relationship between the heterotrophic and autotrophic components of R S was established from the analysis of published data.We adopted this relationship: ln(anR H ) = 1.22 + 0.73 ln(anR S ). (3) The annual R H (anR H ) was estimated by substituting the calculated annual R S (anR S ) into the relationship described above.The annual R A was then calculated by subtracting the annual R H from the annual R S .

Comparison with Earth system models
The estimated R H was compared with the estimates from Earth system models provided by CMIP5.We calculated global R H using 20 Earth system models of the CMIP5 (Table S2) and compared the results with our estimate.

Statistical analysis
We defined tropical, temperate, and boreal regions based on the annual temperature (T < 2 ) after a previous study (Bond-Lamberty and Thomson, 2010b).Statistical analyses were conducted using R versions 3.0.2and 3.1.0(R Core team, 2013).The Mann-Kendall trend test was applied to test for the significance of trends (R package, Kendall version 2.2).

Results
Table 1 lists the a priori and a posteriori distributions of the parameters, and the estimated best parameters with their uncertainties and statistics.The temperature function and precipitation function developed in this study are depicted in Fig. 1, and the two original functions are also plotted.Regarding the temperature function, the three lines overlap at low temperatures (e.g., below 10 • C), but the differences among the three lines increased with an increase in temperature.The temperature sensitivity of the newly estimated function was attenuated at high temperatures compared to the simple exponential functions applied in the original temperature functions, for which the temperature sensitivity steadily increased.Depending upon the parameterization, the newly introduced function can either peak at a certain temperature or behave as a simple exponential function.Our parameterization did not result in a peak in this temperature range, but the temperature sensitivity decreased as the temperature increased.The newly estimated precipitation function was similar to that of the previous study (Raich and Potter, 1995); note that the precipitation used in this study is the weighted average of the precipitation of the current month and the previous month.The best value for the weighting factor α was 0.98, but α was characterized by a large uncertainty (0.03-0.99, 95 % confidence interval).
The estimated mean annual global R S was 91 Pg C yr −1 (1965-2012; Monte Carlo 95 % confidence interval: 87-95 Pg C yr −1 ), and the spatial distribution of R S is depicted in Fig. 2. The estimated R S was high in tropical regions and low in boreal regions, following a temperature-oriented gradient from near the Equator to higher latitudes, but the estimated values were low in dry regions as well (Fig. 2a, b).Latitu-dinally, the regions between 30 • S and 30 • N contributed the most to global R S , but the contribution of the region between 30 and 60 • N was also large (Fig. 2c).The contributions of the tropical, temperate, and boreal regions to global R S were 64, 24, and 12 %, respectively.The monthly global R S was lowest in February (5.7 Pg C mo −1 ) and greatest in July (9.4 Pg C mo −1 ) (Figs. S2 and S3 in the Supplement).The mean annual grid-cell R S was characterized by a broad distribution, ranging from 0 to greater than 1500 g C m −2 yr −1 (Fig. S4).
The estimated R S followed an increasing trend over time, with fluctuations, and the rate of the estimated increase was 0.09 Pg C yr −2 (P < 0.0001) between 1965 and 2012 and 0.14 Pg C yr −2 (P = 0.0015) between 1990 and 2012 (Fig. 3, Table S3).The lowest value of R S (88 Pg C yr −1 ) occurred in both 1965 and 1970, and the highest value (95 Pg C yr −1 ) occurred in 2010.The higher and lower values were mainly coincident with El Niño-Southern Oscillation events.The trends were examined for the tropical, temperate, and boreal regions: the annual variation was largest in tropical regions and was lowest in boreal regions (Figs. 4 and S5).In tropical regions, large fluctuations in R S occurred during the 1970s.In all regions, R S followed an increasing trend with time.The rates of increase in R S for the tropical, temperate, and boreal regions were 0.048 (P < 0.0001), 0.025 (P < 0.0001), and 0.021 (P < 0.0001) Pg C yr −2 , respectively; hence, the highest rates of increase occurred in the tropical regions.The proportional increases in R S of the tropical, temperate, and boreal regions were 0.08, 0.11, and 0.19 %, respectively; thus, the proportional increase was greatest for the boreal regions.The difference between the earlier and later period of the simulation is shown by latitude in Fig. 2d.R S increased at nearly all latitudes.There were large increases in R S between 0 and 30 • N and between 30 and 70 • N.
The relationship between the annual mean global temperature and the global R S is characterized by a slope of 3.3 Pg C yr −1 • C −1 (P < 0.0001) (Fig. 5  of 1.4 (derived by fitting an exponential function).Figure 6 presents the distribution of the Q 10 values at the grid scale, which was calculated using the temperature function estimated in this study and the mean temperature of each grid.
The Q 10 values varied between 1 and 2 and were lower in the regions near the Equator and higher in the regions at high latitudes with colder climates.
The estimated global R H and R A were 51 and 40 Pg C yr −1 , respectively.The spatial distributions of R H and R A are depicted in Fig. 7.Both the R H and the R A were high in tropical regions and low in cold and/or dry regions.The R H and R A were nearly equivalent to each other, but in the regions of high R S , R A was greater than R H ; and in the regions with low R S , R H was greater than R A .The distribution of the R A : R S ratio indicates that, in tropical and temperate regions, the R A component contributes approximately 40-50 % of R S , while R A accounted for less than 30 % of R S in cold and/or dry regions (Fig. 8).
Figure 9 compares the R H estimated by our model to that estimated using other Earth system models.The value of R H estimated by the Earth system models varied from 40 Pg C yr −1 to greater than 77 Pg C yr −1 .The mean of the results from the Earth system models (54 Pg C yr −1 , 1965-2004) was similar to our estimate.The latitudinal distributions of R H differed among the Earth system models (Fig. 10).In particular, the differences among models were large between 30 • S and 10 • N and from 40 to 70 • N. The distribution of the R H estimated using this model was primarily in accordance with the mean of Earth system models; however, a large difference was noted between 10 • S and 10 • N.

Spatiotemporal distribution of R S
Overall, the estimated R S was high in tropical regions and low in cold and/or dry regions.The model parameters derived from the parameterization (Table 1 and Fig. 1) indicate that R S increases under conditions of high temperature and high precipitation.Our modeling suggests that the spatial distribution of R S at the global scale is controlled by both precipitation and temperature (Fig. 2a and Fig. S6).These patterns basically agree with those reported in previous studies (Bond-Lamberty and Thomson, 2010b;Chen et al., 2010;Hashimoto, 2012;Raich et al., 2002).However, the estimated total global R S of this study (91 Pg C yr −1 ) differs from the results of the previous studies.Previous estimates can be roughly divided into two categories, the highest estimate of 98 Pg C yr −1 (Bond-Lamberty and Thomson, 2010b) and the other estimates (76 Pg C yr −1 , on average, N = 5) (Hashimoto, 2012;Raich and Potter, 1995;Raich and Schlesinger, 1992;Raich et al., 2002;Schlesinger, 1977) (Table S4).Our estimate is based on the same data set as that analyzed for the highest estimate (Bond-Lamberty and Thomson, 2010b), but the new estimate of this study was 7 % lower than that estimate.We speculate that one of the reasons for this difference might be the differences in model structure.A non-linear model was used in this study, while linear models were employed in the previous study.In particular, we assumed that R S is sensibly reduced when the sum of precip-   ; please see Table S2).The orange line represents the result of this study.The blue line indicates the mean of the results of the 20 Earth system models.
itation of the current month and previous month is zero.R S was very low in dry regions (e.g., deserts in Africa and Mongolia) (Fig. 2).The numbers of observations are quite limited for very dry regions and deserts; for this reason, although we considered that it is reasonable to assume approximately zero-respiration in these regions, we should consider the potentially high uncertainty in these estimates.However, the new estimate was higher than other previous estimates (i.e., all of the estimates other than Bond-Lamberty and Thomson 2010b).In particular, the new estimate was higher than that of Raich et al. (2002) despite using nearly the same model structure.We attribute this difference to the differences in the data sets analyzed for parameterization (Table S1 and Fig. S7).The global R S followed an increasing trend, and the rate of the increase was comparable to that estimated by a previous Figure 10.Latitudinal distribution of heterotrophic respiration.The orange line represents the results of this study.The grey lines are the results from the 20 Earth system models (please see Table S2).
A smoothing spline was fit to each result because of the variation in the grid sizes of the Earth system models.The solid blue lines and broken blue lines indicate the mean and standard deviation, respectively, of the results of the 20 Earth system models.
study (Bond-Lamberty and Thomson, 2010b).Our model did not include a detailed carbon cycle for the evaluated ecosystems; hence, it is not possible to argue that this increasing trend indicates a net loss of carbon from the soil to the atmosphere (Gottschalk et al., 2012;Smith and Fang, 2010).However, our analysis provides additional data to support an increasing trend for global R S , even though a new model was applied for this study, and supports the assumption that the soil carbon flux from soil to the atmosphere is increasing in response to climate change.

Heterotrophic and autotrophic respirations
Although the number of reports of R H is limited, our estimate of R H is comparable to those of previous studies (IPCC, 2001;Potter and Klooster, 1997;Sitch et al., 2015;Tian et al., 2015).In addition, the mean value of R H estimated using the Earth system models is comparable to our estimate (Fig. 9).This agreement might imply that the carbon cycles in the Earth system models are, to an extent, well constrained by the carbon influx terms (GPP and NPP), and there are, in comparison to R H , numerous global estimates for GPP and NPP.However, when we look at the results from each Earth system model, the differences among estimates are distinct in terms of magnitude and spatial distribution.Because the air temperature simulated by the models in CMIP5 is well correlated with CRU surface air temperature (Todd-Brown et al., 2013), the variation in R H might be attributable to the differences in the description of the terrestrial carbon cycle in each model.R H is a major carbon flux in an ecosystem carbon cycle; therefore, the large variation in R H indicates that there are large uncertainties in the overall flows of carbon in ecosystems (e.g., photosynthesis and respiration) associated with the Earth system model.In addition, the Q 10 value for R H in each Earth system model in CMIP5 ranged from 1.4 to 2.2 (Todd-Brown et al., 2014); thus, the range of Q 10 is wide enough and must contribute to the large variation in R H .In fact, there are large variations among estimates of soil carbon stocks and soil carbon responses to climate change generated using Earth system models (Carvalhais et al., 2014;Nishina et al., 2014;Todd-Brown et al., 2013).
The mean terrestrial NPP reported in previous studies was 56.2 ± 14.3 Pg C yr −1 based on a thorough literature survey (Ito, 2011) (most data included were published after 1990), and our estimated R H between 1990 and 2012 was 51.5 Pg C yr −1 .The residual, the so-called net ecosystem production, is then 4.7 Pg C yr −1 .The global terrestrial carbon sink for 1990-2009 was estimated to be 2.4 Pg C yr −1 (Sitch et al., 2015); when global fire carbon emission (2.0 Pg C yr −1 ; 1997-2009) (van der Werf et al., 2010) is taken into account, despite these figures being based on different approaches, the figures show surprising consistency.
Although previously reported NPP trends vary and are still debated (Table S5) (Ahlström et al., 2012) -and care must be taken to ensure that different climate data were used among the studies -comparing the trends of R H with those of NPP may imply possible changes in net global ecosystem carbon uptake.Before 2000, both NPP and R H showed increasing trends (Table S5); however, the reported trends of NPP were larger than those of R H estimated in this study, suggesting a possible increase in global ecosystem carbon uptake.In the 2000s, the increasing trend of NPP is likely to continue; however, one study reported the possible decline of NPP, which may imply the possible diminishment of increasing global ecosystem carbon uptake (Table S5).However, in this study, R H was estimated using a simple empirical relationship with R S , and the interannual changes in R S are mostly climatedriven and do not include process-based changes in the carbon cycle.Therefore, the trends in R H obtained in this study may be underestimated and must be carefully evaluated.
The estimated global-scale contribution of R A (root respiration) to total R S was 44 %.At the grid scale, there was considerable variation in the R A : R S ratio, which is in agreement with the reports based on compilations of previous laboratory and field studies (Hanson et al., 2000).However, although there are observational reports of R A : R S ratio greater than 0.5, such high R A : R S ratios were not observed in our modeling study because of the relationship between R S and R H applied in this study.Another reason might be that the compilation studies included data observed under various vegetation/soil conditions and seasons, while our study provides a spatiotemporal average.For example, the R A : R S ratio will be high in densely planted vegetation.

Contributions of tropical, temperate, and boreal regions
Our study, similar to previous studies, revealed that tropical regions contribute the largest proportion of global R S (Bond-Lamberty and Thomson, 2010b;Hashimoto, 2012;Raich et al., 2002).This finding is not surprising because R S responds to temperature exponentially and also because there are large amounts of litter input to soil in tropical regions.However, strikingly, the contribution of R S from boreal regions to the rate of increase in R S at the global scale for the study period was on the same order of magnitude as that of the contributions from tropical and temperate regions, despite the lower contribution of R S from boreal regions to the total global R S in terms of absolute magnitude.This relatively large contribution is attributed to the temperature sensitivity of R S (quasi-linear response) and the magnitude of the temperature increase in boreal regions, which was greater than the increase for tropical regions.At present, tropical regions are the most influential regions in terms of global R S (e.g., 64 % of the global R S based on our results).As suggested in previous studies, the importance of boreal regions in the global carbon cycle is increasing and will continue to increase because a large amount of carbon is stored in soils in boreal regions (Batjes, 1996;Dixon et al., 1994;Eswaran et al., 1993;Post et al., 1982).

Temperature sensitivity
R S is strongly influenced by temperature, and an understanding of the response of global R S to the change in global temperature is critical to understanding and predicting the response of the terrestrial carbon cycle to climate change.In our study, global R S responded to the increase in global air temperature, over the study period, at the rate of 3.3 Pg C yr −1 • C −1 (Q 10 = 1.4,based on the air temperature, not the soil temperature), which is in accord with the results of previous studies (Bond-Lamberty and Thomson, 2010b; Raich et al., 2002;Zhou et al., 2009).There are several estimates of the global Q 10 for R H (organic matter decomposition) or ecosystem respiration (the sum of plant and soil respiration), and these values range, approximately, from 1 to 2 and are distributed around 1.5 (Ise and Moorcroft, 2006;Jones and Cox, 2001;Kaminski et al., 2002;Mahecha et al., 2010;Todd-Brown et al., 2013;Zhou et al., 2009) (Table S6 in the Supplement).At the field scale, the observed Q 10 values of R S are typically in the range of 2.0-3.0, are characterized by high variability, and decrease with an increase in temperature (Bond-Lamberty and Thomson, 2010a;Kirschbaum, 1995;Wang et al., 2010;Wei et al., 2010), although the calculated Q 10 value depends on temperature range and on the analyzed temperature (air/soil temperature, depth of soil temperature).In regards to ecosystem respiration, the observed temperature sensitivity at the ecosystem level is seasonally confounded, and an unconfounded Q 10 value of 1.4 has been suggested, even among multiple biomes and independent of the mean temperature (Mahecha et al., 2010).In other words, the Q 10 value, which has a value of 1.4, is observed to be approximately 2.0-3.0 (with a high degree of variation) at the field scale, while the globally estimated value (1.5) is approximately equal to the unconfounded value mentioned above.These apparent differences in temperature sensitivity are curious and are probably attributed to the confounding effects of other ecophysiological factors (e.g., photosynthesis), and differences among analyses conducted at multiple spatiotemporal scales (e.g., Kirschbaum, 2010;Subke and Bahn, 2010).These apparent differences in temperature sensitivity have not yet been fully interpreted.Some studies have addressed this issue; for example, a modeling study (Kirschbaum, 2010) reproduced, in part, such changes in temperature sensitivity across-scale that are introduced by seasonal temperature variations.When process-oriented ecosystem models are applied at the field or grid scale and then scaled up to the global scale, comparisons of the global-scale temperature sensitivity of such scaling efforts with the results of our study may be useful for examining whether R S has been properly scaled.

Conclusions and future work
In this study, we estimated the spatiotemporal variation of global R S using a global soil database, SRDB, and a semiempirical model.The study scaled up the observed fieldscale R S values to a global-scale R S to provide a dataoriented estimate of global R S .The estimated mean annual global R S was 91 Pg C yr −1 (1965-2012; Monte Carlo 95 % confidence interval: 87-95 Pg C yr −1 ), which differs from those of previous studies.Our model does not include detailed processes for ecosystem carbon cycles, imparting both limitations and advantages to this study.For example, plant photosynthesis, belowground carbon allocation, soil carbon stock changes, land-use changes, and nitrogen transformations can affect R S , and in particular, these processes play important roles in long-term simulations of terrestrial carbon cycles.Estimation of R S by satellite remote sensing (e.g., normalized difference vegetation index, NDVI), which includes the vegetation information, may be a promising solution (Huang et al., 2013).In regards to boreal regions, the impact of permafrost melting, which is an important process in northern regions, was not explicitly considered in this study, although SRDB includes some data measured in permafrost regions.However, simple semi-empirical models are good at assimilating accumulated observed field data and providing data-oriented estimations.The relationship between R H and R S is derived from data observations for forest ecosystems, which could affect our estimate of R H .The resolution of our analysis is coarse compared to the scale of the field observations.
Our study has demonstrated that the accumulated data for R S can be used to develop simple, data-oriented models, but in the future, data sets that include other related processes/properties (e.g., LAI, and GPP) will be necessary to generate relatively more sophisticated, simple models and to further constrain process-oriented models.Nevertheless, our approach, the use of a simple model for the analysis of accumulated data resources, provides data-oriented estimates and can be used to bridge a gap between processoriented modeling and observed data sets.We expect that our data-oriented, spatiotemporal estimates will serve as benchmarks and also help to constrain process-oriented models and Earth system models.The gridded outputs are available at http://cse.ffpri.affrc.go.jp/shojih/data/index.html.
The Supplement related to this article is available online at doi:10.5194/bg-12-4121-2015-supplement.

Figure 1 .
Figure 1.Shapes of the temperature function (a) and precipitation function (b).The red line represents the results of this study, and the green-dashed and blue-dashed lines indicate the functions estimated by previous studies(Raich and Potter, 1995;Raich et al., 2002).The grey area is the 95 % confidence interval of the estimated functions.

Figure 2 .Figure 3 .
Figure 2. Spatial distribution of the estimated annual soil respiration (a), the latitudinal patterns of soil respiration components (b, c), and the difference between the earlier (1965-1989) and later (1990-2012) periods of the simulation (d).
) and a Q 10 value www

Figure 4 .Figure 5 .
Figure 4. Interannual variations of soil respiration for boreal, temperate, and tropical regions.The orange lines represent the 5-year moving averages.

Figure 6 .
Figure 6.Spatial distribution of Q 10 values was estimated using the temperature function of Eq. (2) (f T = exp(aT − bT 2 )) and the mean temperature of each grid (T M ).The Q 10 value was calculated by f T (T M + 5)/f T (T M − 5).

Figure 7 .Figure 8 .
Figure 7. Spatial distribution of heterotrophic respiration (a) and autotrophic respiration (b).(c) and (d) depict the latitudinal distributions of heterotrophic and autotrophic respiration per square meter and per 0.5 • , respectively.

Figure 9 .
Figure 9.Comparison global heterotrophic respiration.The grey bars are the results from the 20 Earth system models of the CMIP5; please see TableS2).The orange line represents the result of this study.The blue line indicates the mean of the results of the 20 Earth system models.