Temperature response functions introduce high uncertainty in modelled carbon stocks in cold temperature regimes

Models of carbon cycling in terrestrial ecosystems contain formulations for the dependence of respiration on temperature, but the sensitivity of predicted carbon pools and fluxes to these formulations and their parameterization is not well understood. Thus, we performed an uncertainty analysis of soil organic matter decomposition with respect to its temperature dependency using the ecosystem model LPJGUESS. We used five temperature response functions (Exponential, Arrhenius, Lloyd-Taylor, Gaussian, Van’t Hoff). We determined the parameter confidence ranges of the formulations by nonlinear regression analysis based on eight experimental datasets from Northern Hemisphere ecosystems. We sampled over the confidence ranges of the parameters and ran simulations for each pair of temperature response function and calibration site. We analyzed both the long-term and the short-term heterotrophic soil carbon dynamics over a virtual elevation gradient in southern Switzerland. The temperature relationship of Lloyd-Taylor fitted the overall data set best as the other functions either resulted in poor fits (Exponential, Arrhenius) or were not applicable for all datasets (Gaussian, Van’t Hoff). There were two main sources of uncertainty for model simulations: (1) the lack of confidence in the parameter estimates of the temperature response, which increased with increasing temperature, and (2) the size of the simulated soil carbon pools, which increased with elevation, as slower turn-over times lead to higher carbon stocks and higher associated uncertainties. Our results therefore indicate that such projections are Correspondence to: A. Wolf (annett.wolf@env.ethz.ch) more uncertain for higher elevations and hence also higher latitudes, which are of key importance for the global terrestrial carbon budget.


Introduction
Anthropogenic CO 2 emissions from fossil fuel consumption, cement-manufacturing and deforestation are leading to an increase in atmospheric CO 2 concentrations, thus inducing considerable changes of the climate at global, regional and local scales (Solomon et al., 2007).Atmospheric CO 2 concentrations are also strongly affected by changes in the major global natural carbon reservoirs.For example, at present significantly more carbon is stored in the world's soils than in the atmosphere (Schlesinger, 1997).Climatic changes have a direct impact on global soil carbon stocks, but their quantification is subject to considerable debate and disagreement (Davidson and Janssens, 2006;Kirschbaum, 2006;Hakkenberg et al., 2008).If significant amounts of carbon currently stored as organic matter belowground are transferred to the atmosphere by a warminginduced acceleration of decomposition, a positive feedback to climate change may occur (Kirschbaum, 2000).Conversely, if increases of plant-derived carbon inputs to soils exceed increases in decomposition, the feedback would be negative.The estimate of long-term soil carbon stocks, and their changes in time, therefore is a major source of uncertainty in the projections of soil carbon dynamics (Denman et al., 2007).Despite much research, a consensus has not yet emerged on the climate sensitivity of soil carbon decomposition.
Published by Copernicus Publications on behalf of the European Geosciences Union.
H. Portner et al.: Uncertainty of temperature response functions Soil respiration is commonly divided into two components: root respiration with associated mycorrhizal respiration and soil organic matter (SOM) decomposition.We focus on SOM decomposition here, for which the alternative terms of heterotrophic or microbial soil respiration are often used.SOM has turnover times ranging from years to decades and even centuries.It is often conceptualised as several distinct pools with increasing residence times (Knorr et al., 2005;Kirschbaum, 2004;Eliasson et al., 2005) or as a continuous entity with gradually changing decay rates ( Ågren and Bosatta, 1987;Bosatta and Ågren, 1999).Decomposition of SOM is highly complex, as it is driven by a combination of factors such as temperature (Berg and Laskowski, 2005a), moisture conditions (Cisneros-Dozal et al., 2006) and its chemical quality (Berg and Laskowski, 2005b;Weedon et al., 2009;Cornwell et al., 2008).Many biogeochemical models have been developed and applied to study the response of the carbon cycle to past, current and future changes in climate.While the process of carbon uptake (photosynthesis) is represented in a fairly detailed manner in these models (e.g.BiomeBGC (Thornton et al., 2002), IBIS (Kucharik et al., 2000), LPJ-DGVM (Sitch et al., 2003), LPJ-GUESS (Smith et al., 2001), CLM (Oleson et al., 2004) or Triffid (Cox et al., 2000)), the equally important process of carbon release by soil respiration is represented in a comparatively simple manner (Cramer et al., 2001;Friedlingstein et al., 2006).Some models have been specifically developed to study soil carbon dynamics, but their representation of aboveground productivity and hence litter input is usually highly simplified (Parton et al., 1987;Jenkinson, 1990).Interestingly, there is no agreement on the choice of the form of the response function that is used to describe the sensitivity of soil carbon decomposition to temperature, although the temperature relationship of Lloyd and Taylor (1994) and Q 10 are often used.
In this study, we focus on the sensitivity of LPJ-GUESS (Smith et al., 2001) to a range of possible formulations for the temperature dependency of soil organic matter decomposition in order to evaluate their assets and drawbacks.Earlier studies have addressed parameter uncertainties of LPJ with respect to net primary productivity and net ecosystem exchange (Zaehle et al., 2005) and of LPJ-GUESS with respect to vegetation dynamics (Wramneby et al., 2008).We complement these studies by focusing on soil respiration, and we expand them by considering the relative importance of model formulation vs. parameter uncertainty.Specifically, we address the following questions: 1. Is it possible to identify one temperature response function which can be recommended across different sites as measured by its match with observations and low uncertainties in its parameterization?
2. How does the choice and parameterization of temperature response functions affect uncertainty of carbon flux estimates?
3. How does this translate into uncertainties in simulated long-term carbon storage, compared to the short-term fluxes?

Methods
We chose a comprehensive approach by considering not only the raw fits of candidate functions to calibration datasets, but also the number of parameters, the confidence ranges of parameter estimates and the sensitivity of model output variables.We placed a special focus on the identification of a suitable model formulation that fitted well to experimental data and also led to acceptable uncertainty in the output variables when employed in LPJ-GUESS.
We built upon the well-established LPJ-GUESS model (Smith et al., 2001) and chamber measurements of soil respiration from a range of Ameriflux and CarboEuropeIP sites (Hibbard et al., 2005(Hibbard et al., , 2006)).Instead of performing a complete model-intercomparison of the soil carbon dynamics in different vegetation models, we used only one ecosystem model, LPJ-GUESS.This avoids further uncertainties that would be unavoidable due to the different representations of other processes in different ecosystem models, such as vegetation carbon uptake, plant respiration and transpiration.These differences in model structure are a well-known factor that complicate the interpretation of model inter-comparisons (Cramer et al., 2001;Morales et al., 2005).
In biogeochemical models, the relationship between SOM decomposition and soil temperature is often described by one out of a set of related functions.We tested five candidate formulations: a simple Exponential function with a constant Q 10 , the Arrhenius, the Gaussian, the Van't Hoff and the Lloyd-Taylor function.The Exponential and Arrhenius functions are simplifications of the one proposed by Van't Hoff (1901).Lloyd and Taylor (1994) proposed a modified Arrhenius function and Tuomi et al. (2008) andO'Connell (1990) suggested a Gaussian formulation.The details of the five variants are described in Sect.2.1.2.
To evaluate the uncertainty in the soil carbon related variables simulated by LPJ-GUESS, we preferred to use a virtual climatic gradient (arranged along elevation a.s.l.) rather than the Ameriflux and CarboEuropeIP sites, which represent just scattered points in climate space.Therefore we based our study on the Ticino catchment in southern Switzerland.We evaluated the sensitivity of the model to the uncertainty in model parameters with respect to different process formulations and calibration datasets along the climate gradient.We compared simulated ecosystem properties with observed data from sites with comparable climate to establish the overall consistency of model behavior, but this is not intended to be a model validation.Particularly, there are no observed data of soil carbon dynamics for the catchment our elevation gradient is based on, but this is not a problem for the purpose of this study.

The LPJ-GUESS model
We used the dynamic ecosystem model LPJ-GUESS (Smith et al., 2001;Sitch et al., 2003).The model framework incorporates process-based representations of plant physiology, establishment, competition, mortality and ecosystem biogeochemistry.LPJ-GUESS has been successful in predicting vegetation distribution, net primary production and net ecosystem exchange in many different ecosystems (Smith et al., 2001;Morales et al., 2007).

LPJ-GUESS soil module
Soil carbon in LPJ-GUESS is divided into three distinct pools: litter, fast SOM and slow SOM.The temporal dynamics of the carbon stock (C i ) of each individual pool (i) are modeled on a daily basis and follow first-order kinetics with a decay rate k i (Eq.1).The decay rate itself depends on soil temperature and soil moisture, expressed as the product of the decay rate k i,T ref at a given reference temperature T ref , the temperature response function R T and the moisture response function R M (Eq.2).The decay rate k i,T ref is the reciprocal of turnover time τ i,T ref .
Litter from leaves, roots and tree stems is added to the litter pool at the end of each simulation year.As the model does not discriminate between different litter qualities, the simulated tree species composition plays a minor role, as vegetation influences the soil carbon pools through the amount of litter input only.Additionally, we ensured that litter input is identical for all simulations on the same elevation level.Each of the three soil carbon pools, i.e. litter, fast and slow SOM, has its own specific turnover time (τ i,T ref ) at the reference temperature T ref =10 • C and ample soil moisture: 2.85 y, 33 y and 1000 y for the litter, fast SOM and slow SOM pools, respectively (Meentemeyer, 1978;Foley, 1995a).The mineralized litter is divided into three parts: 70% are respired, whereas 0.45% are transferred to the slow and 29.55% to the fast SOM pool (Foley, 1995a).This fractionation is given by the two parameters f →a = 0.7 (fraction of mineralized litter entering atmosphere) and f →f = 0.985 (fraction of remaining mineralized litter entering the fast pool).All soil carbon pools then undergo decomposition independently, i.e. without feedbacks to the other pools.As there are no longterm experiments of soil carbon dynamics (decades to centuries), the soil dynamics are implemented in a simple way and mainly based on findings from short-term experiments.Nevertheless, process based models such as LPJ-GUESS that are specifically built to make long-term projections based on processes taking place on shorter time scales have been successfully used at a number of different sites (Smith et al., 2001;Morales et al., 2007;Hickler et al., 2004).
Throughout the paper, we refer to heterotrophic soil respiration when talking about soil carbon dynamics and soil carbon fluxes.

Temperature response functions implemented in LPJ-GUESS
Five potential response functions were implemented in the model (Table 1).The Exponential response (E) features a constant Q 10 value.It is motivated by Van't Hoff's rule, stating that the rate of a reaction increases two-to threefold for an increase in temperature by 10 • C ( Van't Hoff, 1901).The Arrhenius function (A) is based on the concept of the activation energy for chemical and biological reactions.However, realizing that the change of the rate of reaction is not constant across temperatures, Van't Hoff suggested a more complex formula (V).Importantly, the Exponential and Arrhenius formulations are direct derivatives of the Van't Hoff formulation, obtained by setting the parameters A = B = 0 and C = B = 0 (Table 1), respectively.The temperature response in the standard implementation of LPJ-GUESS is based on Lloyd and Taylor (1994) (L).It is a variant of the Arrhenius formulation suggested by Lloyd and Taylor (1994) because it often leads to better fits to empirical data by allowing for a decrease in activation energy with increasing temperature.It must meet the condition T > T 0 .The Gaussian function (G) in turn is based on Lloyd-Taylor by taking into account the first three terms of the Taylor series expansion of the exponent of the expression by Lloyd-Taylor (Tuomi et al., 2008).Note that the Exponential, Arrhenius and Lloyd-Taylor response curves are monotonically rising with temperature, whereas the Gaussian and the Van't Hoff curves have a maximum.As the decay constant k i,T ref is valid only at the reference temperature T ref , the response equations were expressed relative to this temperature (Table 1).We thus reparameterized the functions by combining Eqs.3-4, leading to the general scheme of Eq. 5, where f abs , f rel and R T ref refer to the absolute and the relative temperature response and to the reference respiration at a given reference temperature T ref , respectively.
In the default version of LPJ-GUESS, autotrophic (root and mycorrhiza) and heterotrophic soil respiration (SOM decomposition) are modelled using identical temperature responses.
As we focused on SOM decomposition here, only the heterotrophic soil respiration was varied using the five alternative formulations introduced above.The autotrophic respiration influences net primary production and hence plant growth and litter production.To not confound our results by using different litter inputs, we did not vary the temperature www.biogeosciences.net/7/3669/2010/Biogeosciences, 7, 3669-3684, 2010  b The candiate formulations are: Exponential (E), Arrhenius (A), Gaussian (G), Van't Hoff (V) and Lloyd-Taylor (L).response of the autotrophic respiration between the different simulations.

Fitting of the temperature response functions
We used the database compiled by Hibbard et al. (2006), which contains datasets of soil respiration from a range of experimental sites in Europe and North America.We used total soil respiration data to estimate the overall temperature response of soil respiration, as more data sets were available for total soil respiration.Heterotrophic and autotrophic soil respiration have been shown to have a similar contribution to total soil respiration (Bond-Lamberty et al., 2004).Eight sites were selected for calibration (Table 2) to reflect forest vegetation types that span a broad environmental gradient (evergreen-needleleaf, mixed deciduous-evergreen, deciduous-broadleaf); we only used datasets that provided more than 30 measurements of temperature and soil respiration.Measurements were made on a daily basis, distributed over the whole year for time periods ranging from 1995 to 2002, depending on the site.
In nonlinear regression, the usual parameter confidence intervals cannot be used for sensitivity studies because the parameters show non-linear behavior.Therefore, we first linearized all five standardized equations using the method of expected-value parameters (Ratkowsky, 1990).Models in expected-value parameterization are close to linear models in terms of the statistical properties of their parameter estimates, i.e. the confidence intervals of the parameters are comparable, and thus a follow-up uncertainty analysis will give unbiased results.
The functions were linearized by replacing the initial parameters by parameters reflecting the expected-value of the function output at a given position of the curve (Appendix Table A1).We linearized for all parameters but R T ref .
The confidence intervals are provided in the Appendix (Table A3).In order to make the temperature response formulations comparable across the different sites, they were normalized (R T norm ) such that the reference respiration R T ref at reference temperature T ref =10 • C is equal to 1 for each site and equation (Eq.6).
For the sensitivity study, we used all five response formulations at all eight sites and performed nonlinear fits for each dataset-function pair using nonlinear least-squares estimates in the statistics software package R (R Development Core Team, 2008).
To fit the Van't Hoff relationship, we had to introduce an additional data point in each data set at (−40 • C, 0 µmol C m −2 s −1 ) to find a solution of this equation that has no respiration at low temperature (< −40 • C).We determined the 99% confidence intervals for each parameter of each function and the correlation matrix of the parameters for each individual fit.The goodness of each fit was quantified by the Bayesian information criterion (BIC) introduced by Schwarz (1978).
We used the SIMLAB software from the European Joint Research Center (Saltelli et al., 2004) to generate the parameter sample sets.For each fit, we generated a latin hypercube sample (N = 20).We sampled uniformly over the confidence intervals of the parameters and included the parameter dependencies through the correlation matrix obtained in the fitting procedure based on the method of Iman and Conover (1982).We used the 99% confidence intervals of the parameters, and created a sample of parameter sets over their corresponding confidence range for each response function-site pair.
We analyzed the variations when (1) only uncertainty in the temperature response (R T ), was included (2) variable turnover times (R T + τ ) for the litter, fast and slow soil carbon pools τ l , τ f , τ s were incorporated, and (3) variable turnover times and fractionation (R T + τ + F ) of the mineralized litter were included.
As the confidence interval of the turnover times could not be estimated with the given experimental datasets, we followed Parton et al. (1987) and set the range of turnover times for the three carbon pools to 1-5 y for the litter pool, 20-40 y for the fast SOM and 200-1500 y for the slow SOM.The corresponding confidence intervals for the two fractionation parameters were set to 0.633-0.767for f →a and 0.980-0.989for f →f (Foley, 1995b).We thus assumed implicitly that the turnover times and the litter fractionation parameters depended neither on each other nor on the other parameters of the temperature response.

Interpolation of climate data
LPJ-GUESS is driven by daily weather input, including mean temperature, precipitation sum, percentage sunshine and atmospheric CO 2 concentration.The climate data were based on interpolated weather data of a large elevation transect in the Ticino catchment in the Southern Swiss Alps ranging from 300 to 2300 m a.s.l., sampled at 200 m intervals, resulting in a total of 11 individual sites from which we choose three representative sites (at 300 m, 1300 m and 2300 m) to report the estimates from.Annual mean temperatures varied widely, ranging from 11.5 to −1.0 • C along this virtual elevation gradient.
Climate data for the period of 1901-2006 were compiled from different sources.Daily mean temperatures and daily precipitation sums for the period of 1960-2006 were obtained from a spatially explicit climate data set of Switzerland with a spatial resolution of 1 ha.The data were derived using the DAYMET model (Thornton et al., 1997), which was developed specifically for complex terrain such as mountain ranges (data source: Land Use Dynamics, Swiss Federal Institute for Forest, Snow and Landscape Research, Switzerland).For each elevation level we calculated the mean daily temperature and precipitation of 100 adjacent grid points (using a 10×10 grid) at a south-facing slope.
Temperature and precipitation data for the period of 1935-1959 were based on the nearest automated meteorological station Locarno-Monti (distance 24 km), which served as a reference to derive the daily anomalies relative to the longterm climatology of this station.The daily anomalies of the Locarno-Monti station for the years 1935-1959 were applied to the climatology of the years 1960-1970 for each elevation.This prolonged the climate input for each elevation level back to the year 1935.Lastly, the climate for the period 1901-1934 was based on monthly data from the Climate Research Unit (CRU TS 1.2, Mitchell et al. (2003)).For this period, the daily climate anomalies were taken from 35 randomly chosen years out of the Locarno-Monti dataset.The CRU dataset was sampled along the virtual elevation gradient and the daily anomalies were applied to these samples.The dataset for percentage sunshine was based on the reference station Locarno-Monti (1960-2006) and the CRU TS 1.2 dataset for the period of 1901-1959.The same dataset was used for all elevation levels, assuming that mean daily cloud cover did not differ within the valley.

Simulation experiments
Simulations were run for all sites along the virtual elevation gradient for a total of 1106 years.In LPJ-GUESS, each modelled stand is represented by independent replicate patches (N = 30).It is assumed that the patches experience the same climate and have the same soil type.Stochastic processes in the vegetation dynamics, like establishment and mortality may result in different dynamics in different patches, the mean over the patches however approaches an average value.The first 1000 years were used for a model spin-up, whereas the subsequent 106 years corresponded to the calendar years 1901-2006.The spin-up period was based on a constant long-term climate, but considered interannual variations to estimate the equilibria for both soil carbon pools and vegetation composition (Sitch et al., 2003).During the spin-up period, the long-term equilibria of the litter, fast and slow SOM pools were estimated by solving the differential flux equations analytically assuming that the annual litter inputs www.biogeosciences.net/7/3669/2010/Biogeosciences, 7, 3669-3684, 2010 from the years 700 to 900 represent steady state litter inputs, which is legitimate as vegetation composition and productivity reached their equilibrium before simulation year 700 (results not shown).Therefore, in our simulations the length of the model spin-up period has no influence on the steadystate soil carbon stock or its associated variability.An uncertainty analysis was performed for each pair of response formulations and sites separately.We analysed the model output of the year 2006.As the key variable to assess uncertainty, we chose the sum of the three carbon pool sizes at the beginning of August as a proxy for mean annual pool size.The summed soil carbon pool fluxes were also evaluated as monthly sums.We used the month of August because soil respiration was generally highest at that time within the year.

Fit of the functions
We divided the response functions into three groups sharing similar characteristics: (1) Exponential and Arrhenius, (2) Gaussian and Van't Hoff and (3) Lloyd-Taylor.
The Exponential and Arrhenius equations overestimated soil respiration at temperatures below 10 • C in all datasets (Fig. 1).Lloyd-Taylor generally performed better not showing an overestimation at lower temperatures.At five sites (Figs.1a-1e), the Gaussian and Van't Hoff equations yielded a maximum in the temperature range of 15-25 • C, but they provided the best estimates below 10 • C. Because the maximum was located at rather low temperatures, they tended to underestimate respiration at high temperatures, where only few measurements were available.
Most parameters were significant when predicting the soil respiration from temperature, with the exception of the first parameter of the Van't Hoff equation (Appendix Table A2).
The only parameter estimate directly comparable between the different temperature responses was the reference respiration, which ranged from 1.06-1.15µmolCm −2 s −1 at the site BEP to 3.49-3.63µmolCm −2 s −1 at the site THA, respectively (cf.Appendix Table A3; site acronyms are provided in Table 2).
The ranking of the performance of the temperature response functions depended on the criterion used: When the sum of squared residuals was used (Table 3), Van't Hoff performed best (7/8), Gaussian dominated the second rank (5/8) and Lloyd-Taylor dominated the third rank (5/8), but it showed the best fit at the site MEO.When the data for all sites were combined, thus comprising a larger variability of environmental conditions than any site-specific dataset, Lloyd-Taylor showed the best overall fit.The Exponential and Arrhenius formulations generally showed an inferior fit compared to any of the other three equations.2.
Based on the Bayesian information criterion, i.e. when also considering the number of parameters employed in a given formulation, the performance of the Van't Hoff equation was lower because it features the largest number of parameters (Table 4).It was ranked the second best model at four sites.The Gaussian model was best at five sites, the Lloyd-Taylor model at two sites and the Arrhenius model at one site.When assessed using the sum of squared residuals, a BIC: Bayesian information criterion (Schwarz, 1978).Best (lowest) values for each site shown in bold numbers.All is the compound dataset consisting of all eight individual datasets.
Lloyd-Taylor showed the best performance when all the data were analyzed together.It was best at two sites, second best at another two sites and third best at the remaining four sites (Table 4).
The uncertainty of the response function increased with increasing temperature for all sites (Fig. 2, results only shown for site HOW).As expected, uncertainties increased with the number of parameters used: the Exponential and Arrhenius formulations had the lowest ranges of variability (Figs.2a-2b), Gaussian and Van't Hoff the highest (Figs.2c-2d), and Lloyd-Taylor was characterized by intermediate uncertainty ranges (Fig. 2e).

General model behavior along the elevation transect
Vegetation: The model somewhat overestimated biomass along the elevation gradient (13-39 kgCm −2 , results not shown) compared to measurements of forests in southern  Switzerland (12-17 kgCm −2 , Swiss national forest inventory, Speich et al., 2010) and Northern Italy (4.2-15.9kgCm −2 above ground, Rodeghiero and Cescatti, 2005).The higher estimates of the model, however, can partly be explained by the intensive land use in this region in the past (Tinner et al., 1998), as many forests are young and still regrowing after abandonment of pastures and orchards (Baldock et al., 1996).This is not reflected in the model LPJ-GUESS, which simulates potential natural vegetation (Smith et al., 2001) and does therefore not consider management or land use history.The shifts from deciduous trees to needle leaved trees that is simulated to occur at around 1100-1300 m fit with expectations of natural vegetation in Europe (Ellenberg et al., 2009).
Soil respiration: The simulated yearly sums of total soil respiration for the years 2000-2006 varied between 0.3 and 1.0 kgCm −2 y −1 (results not shown) and are in accordance with measurements in a deciduous forest at 800 m in Northern Switzerland (0.49 kgCm −2 y −1 , Ruehr et al., 2010).With an average contribution of heterotrophic to total soil respiration of 50% (Hanson et al., 2000), the simulated range compared well with estimations of total soil respiration made over a narrower elevation gradient (220-1740 m) in Northern Italy (0.5-1.2 kgCm −2 y −1 , Rodeghiero and Cescatti, 2005).

Short-term soil carbon flux
Below, the results for the Exponential and Arrhenius response formulations are combined and referred to as E and A because of the high degree of similarity in their behavior.The results for the Lloyd-Taylor formulation are reported separately (L), but the expressions of Gaussian and Van't Hoff are combined and referred to as G and V due to their similar behavior.
The total soil carbon fluxes to the atmosphere are presented for the case where only the response functions and their parameters were varied.The differences to the case with varying turn-over times and to the case with varying litter fractionation parameters were negligible (results not shown).Unless stated otherwise, units of monthly carbon fluxes in August are given in kgCm −2 month −1 .
Elevation 300 m: Soil carbon fluxes ranged between 0.06 and 0.11 (Fig. 3a), whereby the range was somewhat smaller for the E and A functions.The uncertainty ranges of G and V and Lloyd-Taylor were 1.4 and 1.5 times larger relative to the range of E and A.
Elevation 1300 m: On 1300 m elevation the median values were rather similar ranging from 0.087 to 0.161 (Fig. 3b), although the range of uncertainty was larger for the Gaussian and the Lloyd-Taylor formulations.
Elevation 2300 m: While carbon fluxes increased from 300 to 1300 m, they decreased again (for E and A) up to 2300 m, and three distinct subgroups were identified: E and A with a range of 0.076-0.105,G and V with a range of 0.082-0.159,and Lloyd-Taylor with a range of 0.078-0.145(Fig. 3c).This resulted in uncertainty ranges for G and V and Lloyd-Taylor that were 2.7 and 2.3 times the range of E and A.
Changes with elevation: The medians of monthly respiration for each individual temperature relation followed a bellshaped curve over all 11 simulated sites of the elevation gradient (results not shown), starting with low values at 300 m (Fig. 3a), inflecting at around 1300 m (Fig. 3b) and then decreasing again up to 2300 m (Fig. 3c).Although the medians always were in the range of 0.1 ± 0.02 kgCm −2 month −1 , the ranges of uncertainty increased steadily with elevation (Fig. 3a-3c), particularly for the temperature responses G and V and Lloyd-Taylor, leading to ranges at 2300 m that were 1.5 and 1.7 times larger than the range at 300 m.

Long-term soil carbon stock
Looking at the carbon stock estimates of 2006, the response functions could be divided into the same groups as found in the regression analysis, both according to their median and the magnitude of their uncertainty range (Fig. 4a-4i).If not stated otherwise, the units of carbon pools are kgCm −2 .(d-f) additionally varying turnover times (R T + τ ) and (g-i) litter fractionation (R T + τ + F ). Pairs of response formulations and sites have been grouped according to the temperature response used.The box plots show the median, the lower and upper quartiles and span over the 95% confidence interval.Models are separated by the dashed lines into three distinct groups with similar means and uncertainty ranges.Abbrevations as in Fig. 1.At 300 m and 1300 m (a,b,d,e,g,h) the same ordinate scale is used, whereas at 2300 m (c,f,i) a different scale is used.(d-f) additionally varying turnover times (R T + τ ) and (g-i) litter fractionation (R T + τ + F ). Pairs of response formulations and sites have been grouped according to the temperature response used.The box plots show the median, the lower and upper quartiles and span over the 95% confidence interval.Models are separated by the dashed lines into three distinct groups with similar means and uncertainty ranges.Abbrevations as in Fig. 1.At 300 m and 1300 m (a,b,d,e,g,h) the same ordinate scale is used, whereas at 2300 m (c,f,i) a different scale is used.(d-f) additionally varying turnover times (R T + τ ) and (g-i) litter fractionation (R T + τ + F ). Pairs of response formulations and sites have been grouped according to the temperature response used.The box plots show the median, the lower and upper quartiles and span over the 95% confidence interval.Models are separated by the dashed lines into three distinct groups with similar means and uncertainty ranges.Abbrevations as in Fig. 1.At 300 m and 1300 m (a,b,d,e,g,h) the same ordinate scale is used, whereas at 2300 m (c,f,i) a different scale is used.(d-f) additionally varying turnover times (R T + τ ) and (g-i) litter fractionation (R T + τ + F ). Pairs of response formulations and sites have been grouped according to the temperature response used.The box plots show the median, the lower and upper quartiles and span over the 95% confidence interval.Models are separated by the dashed lines into three distinct groups with similar means and uncertainty ranges.Abbrevations as in Fig. 1.At 300 m and 1300 m (a, b, d, e, g, h) the same ordinate scale is used, whereas at 2300 m (c, f, i) a different scale is used.
Elevation 300 m: If only the temperature response and their parameters were varied, soil carbon stock estimates for E and A ranged from 9.2-13, for Gaussian and Vant't Hoff from 6-15.7 and for Lloyd-Taylor from 8-14.1.The ranges of uncertainty of G and V and Lloyd-Taylor were a factor 2.5 and 1.6 higher than those of the E and A formulations (Fig. 4a).When turnover times were varied as well (Fig. 4d), uncertainty ranges generally increased.The differences between the groups decreased, however, as the medians were more similar.In addition, the range of uncertainty differed less between the groups G and V vs. Lloyd-Taylor, amounting to 1.4 and 1.2 times the uncertainty range of the E and A formulations, respectively (Fig. 4d).The group E and A showed a strong increase in the spread of soil carbon stock, when aslo the uncertainty in the turnover times of the carbon pools was considered.When turnover times and litter fractionation parameters were varied, the variation compared to the E and A formulations did increase slightly to 1.9 times for G and V and to 1.4 times for Lloyd-Taylor (Fig. 4g).

H. Portner et al.: Uncertainty of temperature response functions
Elevation 1300 m: When the response functions and their parameters were varied, the E and A resulted in soil carbon stocks in the range of 14.8-20.2,whereas G and V as well as Lloyd-Taylor showed a larger range of carbon stock estimates: 14.1-23.7 and 15-21.5, respectively (Fig. 4b).The uncertainty ranges of G and V and Lloyd-Taylor amounted to 1.8 and 1.2 times the range of E and A (Fig. 4b).When the variability in turnover times was additionally considered, the uncertainty ranges were much larger (2.0, 1.4 and 2.0 times) for E and A, Gaussian and Van't Hoff and Lloyd-Taylor, respectively (Fig. 4e).If the litter fractionation parameters were also varied, the uncertainty was 2.0, 1.4 and 2.0 times higher for E and A, Gaussian and Vant't Hoff and Lloyd-Taylor compared to the case with only varying temperature responses, which means that it increased only slightly compared to the former case (Fig. 4h).
Elevation 2300 m: When only varying the temperature response formulations and their parameters, soil carbon stocks were generally largest at the highest elevation and showed a much larger range compared to lower elevation sites.Projections ranged from 17.7-38, from 21.4-80.4and from 18.5-64.6for E and A, G and V and Lloyd-Taylor, respectively (Fig. 4c).For the case with varying turnover times we found ranges of 13.6-37.7, 15.8-75.8 and 15.1-59.7,respectively (Fig. 4f).
For the case with varying litter fractionation the ranges amounted to 16.7-37.7,22.2-74.8and 18.2-68.5(Fig. 4i).In contrast to the other two elevations, the range of carbon stock predictions was only slightly affected by the variation in turnover times and litter fractionation parameters.
Changes with elevation: The range of uncertainty increased with increasing elevation for all three subgroups (Fig. 4), whereby the largest uncertainties were found at the 2300 m elevation site for all model formulations (Figs.4c, f,  and i).

Discussion
The reliability of model outputs heavily depends on the uncertainty associated with the choice of functional dependencies in the model and the data sets used to derive parameter values.Often, regression analysis is employed based on experimental data.The uncertainty inherent in the parameter estimates will propagate through the model and lead to a corresponding variation in model output.This has been shown by Jones et al. (2003) who analyzed the temperature sensitivity (Q 10 ) of soil respiration in a fully coupled global circulation model.

Fit of the functions
The temperature response functions could be assigned into three groups: Exponential and Arrhenius, Gaussian and Van't Hoff and Lloyd-Taylor.Using Exponential or Ar-rhenius responses led to an overestimation of respiration at low (<10 • C) temperatures, which resulted in an insufficient fit overall, thus corroborating the results of earlier research (Lloyd and Taylor, 1994).The Exponential response, which is based on a constant Q 10 value is not adequate as the Q 10 value has been shown to decrease with increasing temperature (Kirschbaum, 1995).This has been confirmed by Schindlbacher et al. (2010), who measured respiration in incubated soil samples and showed that the Lloyd-Taylor and Gaussian formulations had an increasing temperature sensitivity along elevation gradients in Spain and Austria.Nevertheless, in our study the Exponential formulation was included in the analysis because the usage of Q 10 values is common (Qi et al., 2002).
For the other three response functions, the rankings differed depending on the criterion employed.As expected, the Van't Hoff equation ranked best when considering the summed square residuals, as it has the largest number of parameters.However, when we used the Bayesian information criterion, which evaluates the model fit relative to the number of parameters, the Gaussian and Lloyd-Taylor formulations performed better.The good performance of the Gaussian function is in line with results from agricultural and forest soils in Finland and Sitka spruce plantations in Scotland (Tuomi et al., 2008).The Lloyd-Taylor formulation has been reported to give good results for a variety of soil types (Lloyd and Taylor, 1994) and it is widely used in soil and ecosystem models (Adair et al., 2008;Kucharik et al., 2000;Thornton et al., 2002).
Although the Gaussian and Lloyd-Taylor equations feature the same number of parameters, the Gaussian formulation outperformed the Lloyd-Taylor function in this study, which is in line with findings by Tuomi et al. (2008).However, when all individual sites were combined, Lloyd-Taylor outperformed both Gaussian and Van't Hoff with respect to a ranking based on both the summed-squared-residuals and the Bayesian information criterion.This is in agreement with Fang and Moncrieff (2001), who found Arrhenius, Lloyd-Taylor and Gaussian function all perform well on data of incubation measurements, whereby the Lloyd-Taylor formulation performed slightly better than the Arrhenius function.
The decrease of respiration rates at high temperatures, found for the Gaussian and Van't Hoff formulations was mainly an artifact of model parameterisation.Although, a decline in respiration rates would be expected at higher temperatures due to microbial protein denaturation, the declines found in our datasets started at too low temperatures (Larcher, 2001).Especially at the sites in the colder temperature regime (BEP, HAR, HES, HOW), Gaussian and Van't Hoff inflected below 20 • C and their parameterization should therefore not be used if the temperature regime frequently includes temperatures higher than the inflection point temperature (Friedlingstein et al., 2006).For temperatures above 35-40 • C, the decrease in respiration rates at higher temperatures is not an artifact, but a process that needs to be considered.
In such a case the Exponential, Arrhenius or Lloyd-Taylor functions would have to be complemented by an additional curve and parameters describing this decline in respiration rates.The data sets used here, however, do not have enough measurements at higher temperatures to provide reliable estimates of the inflection point of the Van't Hoff and the Gaussian curves or for an additional declining curve for the Exponential, Arrhenius or Lloyd-Taylor formulations.
The larger the number of parameters in a given expression, the larger the variability of the overall parameter space.Although each additional parameter improved the curve fit significantly, it also contributed to the total uncertainty inherent in a given response function.
Due to the heteroscedasticity in the data set, the ordinary least square regression we used can cause the uncertainty of the function parameters to be underestimated, but the estimations are neither biased nor inconsistent (Greene, 2002).Due to the risk of the underestimation of the uncertainty, our analysis should be seen as a conservative estimate of the actual spread of the simulation results.The uncertainties of the different parameters are directly comparable however, as we have first linearized the equations prior to the regression, avoiding the bias due to the otherwise non-linear behavior of the parameter estimates (Ratkowsky, 1990).
As the Gaussian and Van't Hoff functions showed better fits at low temperatures one would be tempted to favor them over the other temperature responses.However, as the functions were optimized using a dataset that comprises temperate test sites only, they would need to be verified over a larger temperature range first.Hence, when applying these expressions with the parameters estimated here, in the context of global vegetation modelling efforts, they are likely to have an unsatisfactory performance at warmer future conditions and in tropical and subtropical regions.In our test region, even the site with the highest annual mean temperature (at 300 m on our virtual elevation gradient), soil temperatures of 20 • C were exceeded on average on only 10% of the days per year.For sites at higher elevations and hence lower temperatures, soil temperatures never reached the values where the response curves had the highest uncertainty.Hence, the high degree of uncertainty at higher temperatures has only small or even negligible consequences for the variation in model output in regions where soil temperature normally does not exceed values of 20 • C, for instance in forests at high elevations or high latitudes.
We have to bear in mind, however, that measured data at each individual site may be influenced by additional factors apart from the temperature response function, such as soil moisture conditions (Rodrigo et al., 1997;Cisneros-Dozal et al., 2006), litter chemistry (Berg and Laskowski, 2005b) and soil quality (Conant et al., 2008).Still, the regression analysis based on the compound data set shows that the default response equation of Lloyd-Taylor in LPJ-GUESS can still be used for further work.These findings are in agreement with those by Adair et al. (2008), who found that the proposed relationship of Lloyd-Taylor performed best with a three-pool model on the Long-term Intersite Decomposition Experiment Team (LIDET) data set.The good performance of Lloyd-Taylor, when short-term carbon fluxes are considered, was also shown for range land sites (Del Grosso et al., 2005) and at different flux tower sites (Richardson et al., 2006).Our findings, however, are in contrast to those by Tuomi et al. (2008), who found the Gaussian formula to fit best incubation measurements from different sources.A Gaussian relationship may be generally preferable, but requires reliable calibration across a broad range of temperatures as previously stated by Bauer et al. (2008).This can be illustrated by looking at the variability of the flexing point at higher temperatures that determines the uncertainty of the long-term results.

Short-term soil carbon flux
The short-term soil carbon fluxes in the month of August 2006 as an example output showed a diverse picture along the virtual elevation gradient for both the response function and the size of the soil carbon stock.The medians of the projections under all temperature relationships showed a bellshaped behavior from low to high altitudes, the highest values being found at 1300 m.
The modelled respiration is positively correlated with both soil carbon pool size and soil temperature.As the soil carbon pool size increases with elevation, modelled respiration initially increases, too.The temperature however decreases with elevation.This counteracts the increasing trend of the modelled respiration due to pool size and finally reverses it at higher elevations.The soil carbon fluxes therefore changed from being more limited by the carbon pool size at low elevations to being more limited by the rate of decomposition (i.e.temperature) at high elevations.This is analogous to Atkin and Tjoelker (2003), who found that the temperature dependence of plant respiration is limited by the turnover rate (enzyme activity) at low temperatures and by substrate availability (pool size) at high temperatures.
We could show in our study that beside the size of the soil carbon pool, the choice of a particular soil temperature function had a significant impact on the estimations of the short-term soil carbon turnover, whereas uncertainties in turnover times and litter fractionation had little direct impact on the short-term carbon flux.This is in accordance with Bauer et al. (2008), who have applied different temperature reduction functions to the modified SOILCO2-ROTHC model and found deviations on the simulated estimates of the cumulative CO 2 fluxes.Also Tang and Zhuang (2009) applied a global sensitivity analysis for the Terrestrial Ecosystem Model (TEM) and found the temperature dependency of soil respiration to be one of the important key factors for the uncertainty in estimates of net ecosystem productivity.

Long-term soil carbon stock
We found that at low elevations and hence high temperatures, carbon pools turned over relatively quickly and therefore large carbon stocks did not accumulate.Carbon pools at higher elevations tend to be higher, due to the slower turnover rates.These findings are in agreement with experimental measurements by Rodeghiero and Cescatti (2005); Zinke and Stangenberger (2000), but not with Perruchoud et al. (2000) who found little evidence for a significant influence of climate on soil carbon stocks in Swiss forests.
The uncertainty bounds of total soil carbon stocks generally increased with elevation, i.e. they decreased with increasing mean temperature for all response equations and sites.At first sight, this may appear counter-intuitive as the uncertainty of the response formulation itself was found to increase with temperature.This apparent paradox is caused by the fact that the high spreading of the response function at high temperatures does not result in a high uncertainty of long-term carbon stocks because the carbon is readily decomposed and no large soil carbon pools are formed.It is important to take into account that the accumulation of uncertainty was larger as the average decomposition rate became slower.This was illustrated by the result that the influence of the alterations in turnover times and litter fractionation parameters diminished with increasing elevation.An additional change in an already very low decomposition rate did have only minor effects on the estimations of carbon storage.But the turnover times and litter fractionation may well play a role when considering global estimates as has been shown by Yurova et al. (2010), who did a limited sensitivity analysis of the LPJ soil carbon dynamics module coupled to their climate-C cycle model INMCM (Institute of Numerical Mathematics Climate Model).Of the three parameters studies (litter pool decomposition rate and fractionation parameters), the proportion of decomposed litter allocated to the slow soil carbon pool had the biggest impact on their estimates of soil carbon stocks.
At higher temperatures and thus at lower elevations, uncertainty in long-term soil carbon stocks resulted from the uncertainties in the temperature response itself.Due to high turnover rates, only little carbon accumulated and therefore variation in carbon stock estimations was comparatively low.This may nevertheless be important when comparing ecosystems within the tropics and subtropics, as Holland et al. (2000) showed.They have derived lower and upper confidence limits of the exponential temperature response from measurements and found substantial differences in modelled carbon fluxes and pools, analogous to our results.
The elevation range used in our simulations resulted in a temperature range that covers the temperature range of the calibration data sets, but extents it to lower temperatures (i.e. higher elevations).However, the extension over the lower end of temperatures is not problematic for the inter-pretation of the results, as the uncertainty in the temperature response curve was very low for low temperatures.

Summary
We have found two main sources of uncertainty for model simulations of both short-term and long-term soil carbon dynamics.On one hand there is the variation in the parameter estimates of the temperature responses and on the other hand there is the uncertainty in carbon pools that turn over slowly.The answers to our initial questions are: 1.The equation of Lloyd-Taylor fitted the compound data set of observed soil respiration best and did not add disproportionate to the spreading in parameter estimates.As this equation is already in use in the model LPJ-GUESS, we can confirm its applicability for the tested temperature range.The alternative formulations where not as favorable, because they either resulted in poor fits (Exponential, Arrhenius) or were not applicable when extrapolating to temperatures beyond the calibration datasets (Gaussian, Van't Hoff).
2. The uncertainty of the estimates of the short-term soil carbon dynamics was mainly due to the variation in the parameter estimates of the underlying temperature response functions.This uncertainty decreased with increasing elevation.
3. The uncertainty of the simulated estimates of the longterm soil carbon stocks, however, increased with increasing elevation.The soil carbon at low elevations was readily degraded due to faster turn-over times, whereas at higher elevations, the slower turn-over times lead to higher carbon stocks and as a consequence higher associated uncertainties.However, changing elevation here not only reflects changing of height above sea level but also incorporates more abstract qualities such as differing temperature, species distribution, litter input and soil carbon stock.This increased uncertainty in the size of carbon pools with slow turn-over rates has implications for the uncertainty in the projection of the change of soil carbon stocks driven by climate change.
The increased variation for higher elevations and, when taken as an analog for the higher latitudes in respect to low-temperature regimes, contributes to a high uncertainty when estimating the global carbon budget.

)a
Functions expressed relative to reference temperature T ref = 10 • C with reference respiration R T ref normalized to 1 at mean reference respiration R T ref .

Fig. 1 :Fig. 1 :
Fig. 1: Best non-linear fit for the soil respiration as a function of soil temperature for all sites are shown (E: Ex Arrhenius, G: Gaussian, V: Van't Hoff, L: Lloyd-Taylor).The abbreviations of the sites are explained in Tab. 2.

Fig. 1 .
Fig. 1.Best non-linear fit for the soil respiration as a function of soil temperature for all sites are shown (E: Exponential, A: Arrhenius, G: Gaussian, V: Van't Hoff, L: Lloyd-Taylor).The abbreviations of the sites are explained in Table2.

Fig. 2 :Fig. 2 :Fig. 2 :
Fig. 2: Uncertainty bound for each candidate temperature response formulation spanned out by the sampled functio range sets for the site HOW.The abbreviation of the site is explained in Tab. 2.

Fig. 2 .
Fig. 2. Uncertainty bound for each candidate temperature response function spanned out by the sampled function parameter range sets for the site HOW.The abbreviation of the site is explained in Table 2.

Fig. 3 :
Fig. 3: Uncertainty in short-term soil carbon flux with climate of August 2006 as an example output for the case with only varying temperature response functions (R T ) on (a) 300 m , (b) 1300 m and (c) 2300 m of elevation.Pairs of response functions and sites have been grouped according to the response formulation used.The box plots show the median, the lower and upper quartiles and span over the 95% confidence interval.Models are separated by the dashed lines into groups with similar means and uncertainty ranges.Abbrevations as in Fig. 1.

Fig. 3 .
Fig. 3. Uncertainty in short-term soil carbon flux with climate of August 2006 as an example output for the case with only varying temperature response functions (R T ) on (a) 300 m, (b) 1300 m and (c) 2300 m of elevation.Pairs of response functions and sites have been grouped according to the response formulation used.The box plots show the median, the lower and upper quartiles and span over the 95% confidence interval.Models are separated by the dashed lines into groups with similar means and uncertainty ranges.Abbrevations as in Fig. 1.
Fig.4: Uncertainty in long-term soil carbon stocks with climate of August 2006 as an example output for the cases with (a-c) only varying temperature response functions (R T ), (d-f) additionally varying turnover times (R T + τ ) and (g-i) litter fractionation (R T + τ + F ). Pairs of response formulations and sites have been grouped according to the temperature response used.The box plots show the median, the lower and upper quartiles and span over the 95% confidence interval.Models are separated by the dashed lines into three distinct groups with similar means and uncertainty ranges.Abbrevations as in Fig.1.At 300 m and 1300 m (a,b,d,e,g,h) the same ordinate scale is used, whereas at 2300 m (c,f,i) a different scale is used.
Fig. 4: Uncertainty in long-term soil carbon stocks with climate of August 2006 as an example output for the cases with (a-c) only varying temperature response functions (R T ),(d-f) additionally varying turnover times (R T + τ ) and (g-i) litter fractionation (R T + τ + F ). Pairs of response formulations and sites have been grouped according to the temperature response used.The box plots show the median, the lower and upper quartiles and span over the 95% confidence interval.Models are separated by the dashed lines into three distinct groups with similar means and uncertainty ranges.Abbrevations as in Fig.1.At 300 m and 1300 m (a,b,d,e,g,h) the same ordinate scale is used, whereas at 2300 m (c,f,i) a different scale is used.
Fig. 4: Uncertainty in long-term soil carbon stocks with climate of August 2006 as an example output for the cases with (a-c) only varying temperature response functions (R T ),(d-f) additionally varying turnover times (R T + τ ) and (g-i) litter fractionation (R T + τ + F ). Pairs of response formulations and sites have been grouped according to the temperature response used.The box plots show the median, the lower and upper quartiles and span over the 95% confidence interval.Models are separated by the dashed lines into three distinct groups with similar means and uncertainty ranges.Abbrevations as in Fig.1.At 300 m and 1300 m (a,b,d,e,g,h) the same ordinate scale is used, whereas at 2300 m (c,f,i) a different scale is used.

Fig. 4 .
Fig. 4. Uncertainty in long-term soil carbon stocks with climate of August 2006 as an example output for the cases with (a-c) only varying temperature response functions (R T ),(d-f) additionally varying turnover times (R T + τ ) and (g-i) litter fractionation (R T + τ + F ). Pairs of response formulations and sites have been grouped according to the temperature response used.The box plots show the median, the lower and upper quartiles and span over the 95% confidence interval.Models are separated by the dashed lines into three distinct groups with similar means and uncertainty ranges.Abbrevations as in Fig.1.At 300 m and 1300 m (a, b, d, e, g, h) the same ordinate scale is used, whereas at 2300 m (c, f, i) a different scale is used.
given in P-Values for all the parameters of nonlinear model fits for each pair of temperature response formulation (as given in Appendix TableA1) and calibration site.a Significance codes for P-values: 0 *** 0.001 ** 0.01 * 0.05 + 0.1 ns 1.All is the compound dataset consisting of all eight individual datasets.

Table 1 .
Temperature response functions

Table 2 .
Site characteristics a MAT: Mean annual temperature in • C. b N: Number of data points.

Table 3 .
Summed squared residuals of nonlinear model fits SSR: Summed Squared Residuals.Best (lowest) values for each site shown in bold numbers.All is the compound dataset consisting of all eight individual datasets. a

Table 4 .
Ranking of nonlinear model fits