UvA-DARE

. Heathlands are cultural landscapes which are managed through cyclical cutting, burning or grazing practices. Understanding the carbon (C) ﬂuxes from these ecosystems provides information on the optimal management cycle time to maximise C uptake and minimise C output. The interpretation of ﬁeld data into annual C loss values requires the use of soil respiration models. These generally include model variables related to the underlying drivers of soil respiration, such as soil temperature, soil moisture and plant activity. Very few studies have used selection procedures in which structurally different models are calibrated, then validated on separate observation datasets and the outcomes critically compared. We present thorough model selection procedures to determine soil heterotrophic (microbial) and autotrophic (root) respiration for a heathland chronosequence and show that soil respiration models are required to correct the effect of experimental design on soil temperature. Measures of photosynthesis, plant biomass, photosynthetically active radiation, root biomass, and microbial biomass did not sig-niﬁcantly improve model ﬁt when included with soil temperature. This contradicts many current studies in which these plant variables are used (but not often tested for parameter signiﬁcance). We critically discuss a number of alternative ecosystem variables associated with soil respiration processes in order to inform future experimental planning and model variable selection at other heathland ﬁeld sites. The best predictive model used a generalized linear multi-level model with soil temperature as the only variable. Total annual soil C loss from the young, middle and


Introduction
Soil respiration represents an important source of CO 2 in the biosphere as it is the second largest flux after gross primary productivity in the global carbon (C) cycle, contributing 20-40 % of atmospheric annual C input (Raich and Schlesinger, 1992;Schlesinger and Andrews, 2000).Soil respired CO 2 originates from a number of partitioned belowground sources and this total soil respiration (R S ) can be broadly categorised into autotrophic respiration (R A : the activity of roots and rhizosphere organisms) and heterotrophic respiration (R H : bacteria and fungi decomposition of organic matter and soil faunal activity in the organic and mineral horizons) (Hanson et al., 2000).There has been increasing research attention directed towards quantifying C losses from soil respiration, both at a local ecosystem scale and at a global scale, with the aim of quantifying C balances and predicting C flux changes for the future.
Changes to soil C fluxes have been linked to anthropogenically induced conditions, such as the IPCC predicted climate change (IPCC, 2007), where increased soil warming has resulted in increased C efflux rates (e.g.Davidson and Janssens, 2006;Rustad et al., 2001;Schindlbacher et al., 2012), and prolonged drought periods have resulted in reduced C efflux rates (e.g.Selsted et al., 2012;Sowerby et al., 2008;Suseela et al., 2012).Changes in C fluxes can also be associated with anthropogenic land management regimes, such as the selected land use (e.g.grazing; Peichl et al., 2012), any subsequent land use change (Perez-Quezada et al., 2012); soil disturbances (Novara et al., 2012) and cyclical vegetation management practices like heathland burning Published by Copernicus Publications on behalf of the European Geosciences Union.
These overall changes in CO 2 efflux are often associated with changes to the underlying drivers of R S activity in an ecosystem.These drivers include abiotic factors, such as temperature and soil moisture, and biotic factors, such as gross primary productivity (Bahn et al., 2010a;Davidson and Janssens, 2006;Trumbore, 2006).These factors can interact with each other or can independently affect soil respiration from either R H or R A sources (Davidson et al., 2006).The R H is proportionate to the decomposition of soil carbon by microbial communities, which use mostly the recently produced organic matter as an energy source (Ryan and Law, 2005;Trumbore, 2006).In contrast, CO 2 lost from autotrophic activity is tied to the assimilation of organic compounds supplied by plant metabolism with a part of this carbon rapidly released from the soil (Horwath et al., 1994;Metcalfe et al., 2011;Ryan and Law, 2005).Live root respiration is typically quantified either by using an isotopic approach, such as repeated pulse labeling, continuous labeling, natural abundance (following change of land use/species); by vegetation removal techniques, such as tree girdling; or by using one of the root exclusion methods, such as root removal, trenching and gap analysis (Chemidlin Prévost-Bouré et al., 2009;Díaz-Pinés et al., 2010;Gomez-Casanovas et al., 2012;Graham et al., 2012;Hanson et al., 2000;Jassal and Black, 2006).
Once field measurements have been collected, respiration data has generally been interpreted through statistical analysis to determine any treatment effects.Many studies then additionally processed their observations using the known exponential relationship between organic matter decomposition and temperature (Davidson and Janssens, 2006;Sierra et al., 2011) to determine Q 10 values and investigate the sensitivity of R S to temperature within their studied treatments (e.g.Sowerby et al., 2008;Suseela et al., 2012;Webster et al., 2009;Xiang and Freeman, 2009).However, a much fewer number of studies calculated a continuous CO 2 efflux time series for either the length of the study period or predicted for a projection into the future, to allow the annual C loss from R S (or R H and R A ) to be estimated.Where these continuous efflux series were modelled, the soil temperature relationship was consistently used in model equations while a smaller number of studies included additional measures of ecosystem processes in the model equations.These additional variables have commonly included soil water content or precipitation, as organic matter decomposition and plant activity are affected by moisture availability (Davidson et al., 2006;Raich and Schlesinger, 1992).Increasingly, direct measures of plant activity, such as plant metabolism or litter production, have also been included to link the aboveground processes with the belowground processes that occur within ecosystems (Bahn et al., 2010b;Metcalfe et al., 2011;Ryan and Law, 2005).The degree to which soil respiration models could be modelled with a priori parameter values or required calibration was often dependent on both the spatiotemporal scale at which the models were to be applied and the available environmental data (Keenan et al., 2012).
Where modelling was used to generate annual soil C loss estimates (rather than to generalize the results of an experiment or survey), most studies assessed their selected model using measures of fit for the calibration-data (e.g.Kutsch et al., 2010;Selsted et al., 2012), with many fewer studies evaluating models through a (cross-)validation procedure on separate observation datasets (Caquet et al., 2012;Webster et al., 2009).Furthermore, relatively few studies considered the evaluation of structurally different models and a complete variable selection procedure (Chen et al., 2011;Webster et al., 2009).Recently, several review studies have discussed progress in the modelling of soil respiration and proposed better model-data integration with more rigorous and critical procedures to test respiration models (Keenan et al., 2012;Vargas et al., 2011).Interestingly, soil respiration trial measurements have often been collected repeatedly in time (i.e.longitudinal) and clustered in space but this method has generally not been discussed within the context of soil respiration models.This type of data should ideally be analyzed by hierarchical (multi-level) model framework.In this framework, the measurements which are collected for the same observation unit are explicitly assumed to be dependent, which leads to a more realistic estimate of the effective degrees of freedom (and consequently more realistic confidence bounds) than when assuming independent observations.However, only a few soil respiration studies adopt a multi-level modelling approach (Bernhardt et al., 2006), whereas multi-level modelling is commonplace in many other areas of ecology and the environmental sciences (Qian et al., 2010).In this study, we aimed to follow these guidelines to implement good modelling practices and build predictive models for R S , R A and R H for a managed heathland site.The ultimate goal of this research was to evaluate soil respiration fluxes for the heathland at different vegetation development phases, which would allow for future calculation of a C balance.
Heathlands are cultural landscapes in which cyclical management practices, such as cutting, burning or grazing are undertaken (Webb, 1998).It is known that the structure of the dominant heathland plant (Calluna vulgaris) changes with increasing plant age, from a "net biomass gain" phase up until 15 yr of age, to a "net biomass loss" phase after this time (Gimingham, 1985).It was hypothesized that the younger vegetation community would have the highest plant activity, resulting in greater allocation of C to the roots and therefore a greater R A (and subsequently greater R S ) than on the older communities.Community age was not expected to influence R H as there was no significant difference in the quantity of microbial energy source (carbon stock) between the vegetation ages prior to treatment application (Kopittke et al., 2012).Therefore, based on the known relationships between microbial respiration, soil temperature and soil moisture it was hypothesized that soil temperature and soil moisture would be significant variables for the R H model.In addition, based on the contribution of plant metabolism to root processes, it was hypothesized that soil temperature, soil moisture and a measure of plant activity would contribute significantly to the R S models for all three community ages.

Study site
The investigation was undertaken at a dry heathland, located approximately 25 meters a.s.l. at Oldebroek, the Netherlands.The dominant vascular species at the site is Calluna vulgaris (L.) Hull which grows to a maximum height of 75 cm and provides approximately 95 % of the groundcover, with some Deschamspia flexuosa and Molinia caerulea.The dominant nonvascular species is Hypnum cupressiforme Hedw.with two ecological phenotypes, one growing under Calluna protection and the other adapted to more light between Calluna plants.
The trial was established within a 50 m × 50 m area, at the convergence of three Calluna communities of different ages.Each community age was considered to be a treatment.Replication of these treatments was not possible due to the inherent nature of the site.Therefore, a quasi-experimental design was used, in which groups were selected upon which the variables were tested but where randomization and replication processes were not possible (Campbell and Stanley, 1966).
The oldest heathland area (the old community) was approximately 28 yr of age at the conclusion of the investigation, while the vegetation on the south-eastern third of the research site was approximately 19 yr of age (the middle community).The southern portion of the site was last cut in the year 2000 as part of the creation of a fire break and was 12 yr old (the young community) at the conclusion of the study.
The site is relatively flat in the west and rises in the east and north-east onto a gentle slope with a south-western aspect.The soil is a nutrient-poor, well drained, acid sandy Haplic Podzol (van Meeteren et al., 2008).The soil has an organic horizon which ranged between 1.4 and 8 cm thick, with the mean thickness of 3.9 cm ± 0.04 (Kopittke et al., 2012).The carbon stock of the soil (organic layer and to 25 cm depth of mineral soil) was 8.01 ± 0.6 kg m −2 on the young community, 7.61 ± 0.5 kg m −2 on the middle community and 6.18 ± 0.4 kg m −2 on the old community and were not significantly different to each other (Kopittke et al., 2012).Further information about the site location, species composition and climate is provided in Table 1.

Experimental design
To measure soil respiration and calibrate soil respiration models, eight experimental plots (60 cm × 60 cm) were es-tablished within each heathland age in April 2011 (n = 24).Four of these plots were used to measure heterotrophic respiration on each community age (henceforth called 'trenched' plots; n = 12) and the other four were used to measure total respiration on each community age ("untrenched" plots; n = 12).In this study, the terminology "total soil respiration" and "heterotrophic soil respiration" refers to the observed field data from the untrenched plots and trenched plots, respectively.Due to the inherent nature of the site, randomization of the factor "community age" is not possible in our experiment.However, collinearity of weather data with the distribution of the three age classes is highly unlikely since the area is small compared to variations in weather-variables.Furthermore, soil data (including soil temperature and soil moisture) appear not to vary much between the age classes.The terminology "R S ", "R H " and "R A " refers to the modelled total soil respiration, modelled heterotrophic soil respiration and modelled autotrophic respiration, respectively.
The plots were placed in pairs (one trenched in combination with one untrenched plot) that were 1.5 m apart, but the exact location of the individual plot as well as the location of the pairs were randomly allocated within each vegetation age (Fig. 1).In May 2011, the aboveground biomass was harvested from the four trenched plots within each age group and a narrow trench was excavated to 50 cm depth around the 60 cm × 60 cm plot area.This depth extended below the main rooting zone, but was above the water table and did not encounter any impermeable layers, all of which may have affected CO 2 concentration productions at depths (Jassal and Black, 2006).A nylon mesh (Plastok Associated Ltd., Birkenhead, Wirral, UK) of 41 µm was placed in the trench to prevent the new roots growing into the plots during the experiment.The soil horizons were backfilled in the order of removal to keep soil disturbance to a minimum.Any subsequent vegetation regrowth was periodically removed but the remains left in the plot on the soil surface.The remaining four untrenched plots in each vegetation age were not disturbed and were used as a control treatment.
For the purposes of soil respiration model validation, an additional four plots ("trenched validation" plots) in each heathland age group were trenched using the described method (n = 12) and data collected for the purposes of validating the derived R H model.A further nine untrenched plots ("untrenched validation" plots) were established in the old vegetation and the collected data was used for validation of the derived R S model.

Site meteorological and treatment soil conditions
Site meteorological conditions were recorded on an hourly basis (Decagon Devices Inc.; DC, USA).Air temperature and relative humidity measurements were obtained from 20 cm above ground surface at a central location on the site.Rainfall was measured using a Vaisala tipping bucket rain gauge www.biogeosciences.net/10/3007/2013/Biogeosciences, 10, 3007-3038, 2013 a Water extraction of 1 : 5 for organic horizons and 1 : 1 for mineral horizons.
b Obtained following a rainfall event and reported as a percentage (g per g dry weight soil).
Treatment soil conditions were recorded on an hourly basis (Decagon Devices Inc.; DC, USA).Soil moisture (m 3 m −3 ) and soil temperature ( • C) measurements were obtained from 4-7 cm below ground surface in two trenched plots, two untrenched plots, and two trenched validation plots in each heathland age group (5TM Sensor, Decagon Devices Inc., DC, USA).The same measurements were obtained from the three untrenched validation plots on the old community.In total, 21 soil probes were installed, with six being in the young community, six in the middle community and nine in the old community.

Soil respiration measurements
Respiration collars of 10 cm diameter and 6 cm height were inserted approximately one centimeter into the soil surface in each plot, maintaining a minimum buffer zone of 10 cm from the plot boundary.In the untrenched plots, moss was removed from inside these collars, to ensure that only soil respiration was measured.Moss was not present on the trenched plots as it had been removed during trenching activities.Soil respiration measurements were obtained using a portable gas exchange and fluorescence system (LI-6400XT; LICOR Biosciences, Lincoln, NE USA) in combination with a soil CO 2 flux chamber (LI-6400-09; LICOR Biosciences) which fitted onto the collars.
Soil respiration measurements using this methodology commenced in May 2011, three days after trenching occurred, and continued until August 2012.A total of 29 measurement events occurred post-trenching on the three ages of vegetation.A common effect of a trenching methodology is a flush of CO 2 within the first weeks or months after trenching which originates from decomposing roots (Hanson et al., 2000).To minimise this effect of root decomposition, the first four months of CO 2 efflux measurements were excluded from the study and only observations after 21 September 2011 are included in the analyses.In addition, to determine if there had been significant root biomass loss from the trenched plots (ie decomposition) during the study period, the root biomass in the trenched and the untrenched plots was assessed one year after trenching activities.There were 19 soil respiration measurement events from September 2011 until August 2012.
Soil respiration measurements using the above methodology were also obtained from the trenched validation plots to validate the R H model and from the untrenched validation plots to validate the old vegetation R S model.

Photosynthesis measurements
The gross photosynthetic rate provided a measure of plant activity for the three heathland ages.This gross photosynthetic rate (P G ) was calculated as the net ecosystem exchange (NEE) rate of CO 2 flux minus the ecosystem respiration (ER) rate of CO 2 flux (µmol CO 2 m −2 s −1 ).This photosynthetic rate has a negative sign.A loess smoother curve was applied to the P G data to obtain daily estimates of plant activity for use in the soil respiration models.
Three permanent sampling locations were selected in each vegetation age.A metal base frame (60 cm × 60 cm) was permanently installed using small, narrow sandbags to provide a seal between the frame and the soil surface and fixed with metal pins.The CO 2 fluxes of the vegetation were measured with the same LI-6400 infrared gas analyzer as used for the soil respiration measurements (LICOR, Lincoln, NE, USA) but in this case attached to a 288 L ultra-violet light transparent Perspex chamber (60 cm × 60 cm × 80 cm).Full details of the NEE and ER method used at this site are provided in Appendix A.

Plant and microbial biomass
The biomass harvested from the trenched plots in April 2011 was separated into Calluna and moss layers.These components were oven dried at 70 • C and the dry weight recorded (n = 12).
Microbial biomass and root biomass were sampled in May 2012, approximately one year after trenching activities.Soil sampling was undertaken using a soil corer of 5cm diameter and intact soil samples were obtained from the organic horizon and 0-5 cm mineral soil.Three cores were obtained and were bulked by soil horizon from each trenched plot (n = 12) and untrenched plot (n = 12).The soils were kept refrigerated during preparation.All the soil was sieved and roots were separated, washed, oven dried at 70 • C and the root dry weight calculated for the organic and the mineral horizon.
In the organic horizon, each sample was divided into three subsamples of each 10 g.One part was analyzed for water content by drying at 70 • C and bulk density was then calculated.Samples were ground and the C concentrations were analyzed on a CNS analyzer (Vario EL Analyzer, Elementar).Another subsample was fumigated with the chloroform-fumigation method and extracted for 1 h in 50 mL 0.5 M K 2 SO 4 (Jonasson et al., 1996).The third soil fraction was extracted for 1 h without prior fumigation for initial content of carbon and nutrients.The extractions were frozen until shortly before analysis.Upon defrosting, analysis of total organic C (TOC) was undertaken on a vario TOC cube (Elementar).Microbial C was estimated as the difference between the concentration of TOC in the fumigated and un-fumigated extract.An extractability constant of K EC = 0.45 was used for microbial C (Jensen et al., 2003).Microbial C (mg) of the organic horizon is reported per gram of substrate C.

Data analysis
The data analysis workflow approach is described in the following sections and is summarised in Fig. 2. Initially, the observational data was analyzed to determine if there were statistically significant differences between community ages (an age effect) or between trenched and untrenched plots (a methodological effect).This indicated how the datasets should be grouped in the later modelling phase; for example, if there was no soil respiration difference between trenched plots on the three community ages and there was no hypothesized environmental reason as to why there should be a R H difference, then the three age datasets were grouped for the modelling phase.
Once the observation data had been statistically analyzed, a number of plausible model formats and explanatory variables were chosen for calibration and validation (Jørgensen and Bendoricchio, 2001).The explanatory variables were chosen around the major drivers of R S and R H , being abiotic factors, such as temperature and soil moisture, and biotic factors, such as gross primary productivity (Bahn et al., 2010a;Davidson and Janssens, 2006;Trumbore, 2006).A number of drivers were considered for inclusion as explanatory variables but the final decision was based on the observation data available, the outcome of the statistical analysis, the variables used in other studies and the outcome of a preliminary fitting of the models.Preliminary model fitting indicated that no model could account for the extreme values recorded on 21 March 2012, which were associated with an extreme meteorological episode (freeze followed by thaw).In addition, the misfit on this day dominated the overall performance criterion.These extreme values are most likely associated with the death of fine roots and microbial populations during a late winter (February), extreme freeze period (−20 • C), followed by the rapid recovery of microbial populations as daytime air temperatures reached > 15 • C (March) which all lead to shortterm fluxes of CO 2 from the soil (Matzner and Borken, 2008;Sulkava and Huhta, 2003).Although these CO 2 releases occur, there is strong evidence that these events have little effect on soil C losses at an annual time scale (Matzner and Borken, 2008), therefore it was decided to omit this specific extreme event in the modelling process.This allowed the model to be calibrated and validated more accurately on the observations in which non-extreme processes are believed to be dominant.
The models were calibrated and validated, using the procedures described in the following sections.Based on these results, a model was selected and soil respiration rates were predicted.These values were used to estimate annual R S , R H and R A C loss for each community.

Observational data analysis
The effect of community age on the single occasion measurements of biomass (plant leaves/shoots, plant roots and microbes) was investigated by a linear model ANOVA.If a treatment effect was identified, then a pairwise t tests (using the Bonferroni correction factor) was undertaken whereby an effect is considered as significant if it is associated p value is smaller than 0.05.The effect of vegetation age on the repeated measurements (soil respiration and on photosynthetic activity) was investigated using a linear multi-level model (Pinheiro and Bates, 2000).Where the response variable in the linear multi-level model was the CO 2 efflux measurement (a repeated measurement per location), the vegetation ages formed the fixed effects and the measurement locations formed the random effects.
Where mean results are referenced, the standard errors of the mean (SEM) are provided in both text and graphics.For all statistical analyses, the R statistical computing program was used (R Development Core Team, 2008).

Soil moisture model
A zero-dimensional finite difference soil moisture model (i.e. a "bucket model"), with a daily time resolution and model inputs of rainfall plus air temperature, was constructed and calibrated on the observed soil moisture data.When compared to observed data, the soil model gives an unbiased prediction and explains approximately 70 % of the variance (the details of this model are given in Appendix B).The soil moisture information in this study is used as a potential explanatory variable for respiration.A soil moisture model, rather than observed soil moisture, was used for two reasons.Firstly, a dynamic model is an appropriate method to integrate the soil moisture values per sensor to an average soil moisture value per treatment and this integration is necessary Biogeosciences, 10, 3007-3038, 2013 www.biogeosciences.net/10/3007/2013/because not all plots were equipped with a soil moisture sensor.Secondly, it overcomes problems of missing data, such as when a respiration model is used at other sites for predictive purposes: in these cases, the soil moisture data is usually not available, whereas daily rainfall and temperature are commonly present.The result of the soil moisture model was treated as if it were an observation for an individual observation plot (without assuming any additional variable uncertainty associated with modelling process).

Soil respiration model calibration and validation
A model comparison framework was used to assess the R S models and R H models (Burnham and Anderson, 2002).
A number of plausible models were calibrated and only the models with significant parameter values were retained.These models were ranked according to the root mean squared error for the calibration data (RMSE C ) and the models with low RMSE C were considered suitable for further validation and discussion.
Validation of the suitable models was done with soil respiration data obtained from the validation plots.The models were fitted and validated to data in accordance with Table 2, where validation was conducted on different observation data over the calibration time period (Validation Type I) and for a different time period (Validation Type II).A third validation method (cross-validation) was also used in the model selection procedure.The cross-validation results did not alter the outcomes of the model selection procedure; therefore, the full details on method and results are not discussed further here but are provided in Appendix C.
For each of the validation datasets, a root mean squared error (RMSE V ) was calculated.The RMSE is specified in Eq. ( 1).

RMSE
where Ri is the predicted respiration at time i, R i is the observed respiration at time i and n is total number the number of observations.The general equation is identical when applied to calibration or validation data, as well as for R S and R H .The group of plausible models was built-up as follows.First, an existing soil respiration model was selected from a study undertaken on a comparable Calluna vulgaris heathland located in Denmark (Selsted et al., 2012).This model (henceforth denoted as the Selsted model) is used in this study as a null model for both R S and R H .It is a nonlinear model with three explanatory variables (temperature, soil moisture and biomass) and four parameters that need to be calibrated (further details follow below).The model selection procedure calibrated and validated not only the full model with three explanatory variables, but also the more parsimonious variants with two variables (temperature and soil moisture or temperature and biomass) and with one variable (temperature).
Next, a linear multi-level model (LMM) with the same variables as the Selsted model was calibrated and validated.The multi-level structure is required to deal with the repeated measurements on individual locations.Furthermore, a generalized linear multi-level model (GLMM) with a Gammadistributed error and a log link function (again with the same variables) was calibrated and validated.Generalized linear models extend linear models that involve non-normal error distributions or heteroscedasticity and may also require a transformation to become linear.Linear functions of the predictor variables are obtained by transforming the right side of the equation by a so-called link function.In this case the shape of the relationship is exponential, so by taking its logarithm it becomes linear.The data are then fit in this transformed scale (using an iterative routine based on least squares), but the expected variance is calculated on the original scale of the predictor variables.The Gamma distribution describes that the error is right-skewed at low values of the predictor variable and becomes symmetric at higher values.In our case, the mean and variance of the model error are equal (McCullagh and Nelder, 1989).
In a next step, the soil moisture and biomass variables were transformed into quadratic variables and the LMM and GLMM models using these variables were also calibrated and validated (these models are denoted by LMM2 and GLMM2).These quadratic forms of the models were successfully applied in the study by Khomik et al. (2009).
Following the approach by Selsted et al. (2012), soil moisture as well as biomass was scaled to represent relative soil moisture and relative biomass.Equations ( 2) and (3), respectively, provide the details of these transformed variables.
where M is the relative soil moisture content (a fraction between approximately 0.1 and 1), θ is the volumetric soil moisture content (in this study output from a dynamic soil moisture model, Sect.2.7.2), θ fc is the soil moisture content at field capacity.An estimate for θ fc was available per treatment from the soil moisture model (Sect.2.7.2).
where B is the relative biomass (a fraction between approximately 0.3 and 1), "Biomass" is the aboveground Calluna biomass in g m −2 for a given observation plot and "Max Biomass" (a value of 2.2) for the plot with the greatest quantity of aboveground biomass.Moss was also harvested from the plots, however only the Calluna biomass was used in this calculation as the Calluna root systems were expected to contribute to R A but the moss layer lacks a rooting system and www.biogeosciences.net/10/3007/2013/Biogeosciences, 10, 3007-3038, 2013  2012), peak biomass was estimated using nondestructive techniques.In the current study, the biomass initially harvested from the trenched plots within each nested replicate was used as an estimate of aboveground biomass for the untrenched plots in the same nested replicate.However, as harvested biomass does not give a dynamic measure of plant activity throughout the year and the changes of seasons, a measure of photosynthetic activity (Sect.2.5) was included in the model testing process as an alternative variable for Calluna biomass.A value for relative photosynthetic activity was calculated as follows: where P G is the gross photosynthesis measured per plot in µmol CO 2 m −2 s −1 , and Maximum P G is the maximum CO 2 consumption rate measured during the study, as described in Sect.2.5.The absolute values for Maximum P G were 23.2, 12.2 and 8.1, respectively, for young, middle and old communities.
In the first modelling cycle, the soil temperature at 5 cm depth (T soil ) was used, as it is a common component of soil respiration models.However, in a second modelling cycle, air temperature (T air ) was also tested as a substitute for soil temperature, as it is often a more commonly recorded variable across ecosystems.The equations for the Selsted, LMM, LMM2, GLMM and GLMM2 models using the T , M, B and P variables (see Eqs. 2-4) as predictor variables are shown in Table 3.
In addition to the variables detailed above, a number of other variables were tested in an early explorative phase that occurred prior to the formal model identification process.These other variables included Photosynthetically Active Radiation (PAR) values used as a substitute for the P variable, the microbial biomass as a substitute for the B variable and the root biomass as a substitute for the B variable in both R S and R H models.In that explorative phase, it was found that the RMSE C and RMSE V values for the models involving these variables were higher or close to those variables shown in Table 3.Therefore, the results of these alternative variable combinations were not tested further.

Soil respiration model selection and generation of predictions
The final models for R S and R H were selected using the following rationale.Firstly, the calibrated models in which all coefficients were significant were identified and retained for further consideration.Secondly, only models in which parameters were feasible according to literature values and experience were retained.The reasonableness of these parameters were defined for basal respiration rate: (R 0 0 to 0.5 for the Selsted and GLMM models), a > 0 and c > 0 (GLMM models) or a < 0 and c < 0 (LMM models).For the R H models, a complete set of validation data for each vegetation age was available.Therefore, the subset of R H models with significant parameter values were further assessed by their RMSE V1 values, and those with the lowest values were considered most suitable.
In the R S models, the validation data and therefore, the RMSE V1 and RMSE V2 s, were only available for the old community.Consequently, the RMSE C provided a better measure of model performance across each age of vegetation.Hence, the R S models with significant parameter values and the lowest RMSE C were selected while the values for RMSE V1 and RMSE V2 were of secondary importance (these should lie in the lower to mid-range of all RMSE values).
Following the selection of the model, R S and R H were predicted for the length of the study period using a single hourly soil temperature dataset from the untrenched treatment.The mean annual C loss from R S and the 95 % confidence intervals of model predictions were calculated using a bootstrap procedure with 1000 replications (Davison and Hinkley, 1997).

Vegetation characteristics
Destructive vegetation sampling indicated that mean Calluna aboveground biomass was lowest on the young community and greatest on the middle community.The difference between young and middle communities was Biogeosciences, 10, 3007-3038, 2013 www.biogeosciences.net/10/3007/2013/ R 0 e kT e aM R 0 e kT e aM T B R 0 e kT e cB -T P R 0 e kT e cP -T MB R 0 e kT e aM e cB -T MP R 0 e kT e aM e cP -GLMM2 * * T M R 0 e kT e a(M−1) 2 R 0 e kT e a(M−1) 2 T B R 0 e kT e c(B−1) 2 -T P R 0 e kT e c(P −1) 2 -T MB R 0 e kT e a(M−1) 2 e c(B−1) 2 -T MP R 0 e kT e a(M−1) 2 e c(P −1) 2 just above the 0.05 significance level after the Bonferroni correction (p = 0.059; Fig. 3a), hence we do not interpret this as a difference.The biomass of the moss layer was almost double on the young community (0.43 ± 0.09 kg m −2 ) than the moss biomass on either the middle community (0.27 ± 0.04 kg m −2 ) or old community (0.26 ± 0.04 kg m −2 ; results not shown).Photosynthesis, as a measure of plant activity throughout the year, was greatest in the summer months, least in the winter months and was significantly different between communities (F = 25.1, p < 0.001; Fig. 3b).In winter months, there was no significant difference between mean photosyn-thesis on the young (−2. communities. However, in summer months there was significantly greater photosynthesis on the young community (−16.0 ± 1.4 µmol CO 2 m −2 s −1 ) than on either the middle community (−5.7 ± 1.5 µmol CO 2 m −2 s −1 ) or the old community (−5.2 ± 1.0 µmol CO 2 m −2 s −1 ).The old community was significantly different to the middle community in summer (F = 48, p = 0.049) but there were no other seasonal differences between the mean photosynthetic rates of the middle and old communities during the study period.

Soil respiration
The age of the vegetation had a significant effect on soil respiration (F = 5, p = 0.035) and in every season of the year, total soil respiration was significantly greater on the young community than on the old community (winter p = 0.034, spring p = 0.0144, summer p = 0.007, autumn p = 0.006).
The greatest mean total soil respiration was recorded in summer months on all three communities, ranging from a mean of 2.8 ± 0.2 µmol CO 2 m −2 s −1 on the young community to 2.1 ± 1.9 µmol CO 2 m −2 s −1 on the old community (Fig. 4a).
The lowest mean soil respiration values were recorded in winter, although soil respiration was still significantly greater than zero (t = 14.7, p < 0.001) in these colder conditions.The differences between the communities were greatest in spring with total soil respiration on the young community (1.9 ± 0.2 µmol CO 2 m −2 s −1 ) exceeding respiration on the middle community by a factor of 1.6 and exceeding the old community by a factor of 1.7.
There was no effect of community age in any season for heterotrophic soil respiration on the trenched plots (Fig. 4b).Therefore, the heterotrophic data was not split into age treatments for further analyses, but was treated as a single dataset.Mean heterotrophic soil respiration was least in winter months (0.4 ± 0.05 µmol CO 2 m −2 s −1 ) and greatest in summer months (1.7 ± 0.09 µmol CO 2 m −2 s −1 ).
A peak was observed in both total soil respiration and heterotrophic soil respiration on 21 March 2012.The elevated respiration results were observed on both trenched and untrenched plots and, although the CO 2 flux was variable between measurement locations, the largest fluxes were generally observed on the young community.The maximum respiration observed on this day for the trenched plots was 10.28 µmol CO 2 m −2 s −1 (young community) and for the un- trenched plots was 5.11 µmol CO 2 m −2 s −1 (also young community).

Treatment effect
Soil temperature at 5 cm below ground surface was significantly different between the trenched plots and the untrenched plots over the study period (Fig. 5a).The mean soil temperature in winter was significantly lower on the trenched plots (3.8 ± 0. (4.8 ± 0.10 • C) (t = 4.24, p = 0.016, 3 plots and using averages of soil temperature per plot for the complete winter).However, the reverse occurred in summer, where mean soil temperature was significantly greater on the trenched plots (16.9 ± 0.28 • C) than on the untrenched plot (15.5 ± 0. Microbial C was not significantly different between the trenched plots and the untrenched plots in the organic horizons for any of the young, middle or old vegetation ages (Fig. 6).On untrenched plots, the organic horizon microbial C was significantly greater in the young vegetation than the middle or old vegetation but there was no significant difference between the middle and the old vegetation.
Root biomass (summed for organic horizon and 0-5 cm mineral horizon) was not significantly different between the trenched and untrenched plots on the young, middle or the old vegetation (Fig. 6).Additionally, the root biomass in the untrenched plots was not significantly affected by the vegetation age.However, there was a significantly greater root biomass in the organic horizon than in the 0-5 cm mineral horizon for all vegetation ages (data not shown).

Calibration of the model for total soil respiration (R S )
All model predictions of soil respiration generally followed the seasonal soil temperature patterns, where the lowest respiration was recorded in winter (in February).However, not all models predicted the highest respiration equally, with some models predicting peak values in June, while others predicted peak values in August.
Stepwise application of variables into the different models using the untrenched datasets produced models with absolute RMSE C values that ranged from 0.30 to 2.32 (Fig. 7, left panel).When soil temperature (T soil ) was assigned as the T variable, the RMSEs were generally lower than when air temperature was used (T air ).The lowest RMSE C values were obtained using the Selsted, GLMM and GLMM2 models and therefore, the LMM and LMM2 models are not further discussed within this results section.A selection of models and the RMSE C values are provided in Table 4.
Within the GLMM and GLMM2 model formats the use of the explanatory variable T soil resulted in lower mean RMSE C values (0.31 to 0.49) than where T air (0.35 to 0.68) was included, with the exception of the T +M+P models.When all three variables T + M + P were used in the GLMM format, the model over-predicted soil efflux and resulted in very high errors (0.68 to 2.32); thus, these were excluded from further consideration.This did not occur with the GLMM2 format.When T (both for T soil or T air ) was the only variable used, the model parameters were significant for all three young, middle and old dataset.
The GLMM models in which all parameters were considered significant occurred on 18 occasions.However, the GLMM model was only significant for all three vegetation communities when the T variable (either T soil and T air ) was used alone (Fig. 7, left panel).The parameters in all of these significant models were considered reasonable.Table 5 lists the parameter values for the GLMM models T soil , T soil + M and T soil + P .It appeared that adding soil moisture to a model with only temperature only lowered R 0 while hardly influencing the parameter value associated with temperature (k), whereas adding photosynthesis had the reverse effect (it lowered the k-parameter, associated with temperature, and did not influence R 0 ).The stepwise addition of the Selsted equation resulted in RMSE C values that were very similar to the GLMM and GLMM2 models.However, there were fewer occasions for the Selsted models (10 occasions) than for the GLMM models in which all parameters were significant (Fig. 7

Calibration of the model for heterotrophic soil respiration (R H )
Stepwise application of variables into the different models using the trenched data produced models with absolute RMSE C values that ranged from 0.3 to 0.44 (Fig. 7, right panel).The RMSE C values were lower on the heterotrophic models in which T soil was used as the T variable, rather than T air .Similarly, the GLMM models and Selsted models resulted in RMSE C values that were lower than the LMM mod-els.Therefore, only Selsted, GLMM and GLMM2 T soil models are further discussed within this results section.The significant GLMM models were those which applied T soil variable singly and also T soil in combination with M. This was not the case for the Selsted model, where only T soil applied alone resulted in a model in which all parameters were significant.The GLMM T soil +M model had the lowest RMSE C , while the Selsted and GLMM T soil models had the second lowest RMSE C .All parameters were considered to be reasonable for these significant models.www.biogeosciences.net/10/3007/2013/Biogeosciences, 10, 3007-3038, 2013 Table 5 lists the parameter values for the R H GLMM models T soil and T soil + M. Similar to the models for R S , adding soil moisture to a model with only temperature lowers R 0 while not influencing the parameter value associated with temperature (k).

Model validation
The calibrated models were used on the validation data for period one (September 2011-August 2012) and period two (November 2010-August 2011, see Table 2).The resulting RMSE validation values (RMSE V1 and RMSE V2 ) were then compared to the RMSE calibration values (RMSE C ).A selection of these results, are shown in Table 4.The R S models which had the lowest RMSE C and RMSE V values used the GLMM format with T soil as a single variable, with T soil + P and with T air + P .Of these, the GLMM T soil model and the GLMM T air + P model were the only ones where all parameters were significant for all vegetation ages.The R S models which performed the worst in the validation phase also used the GLMM format and included the T variable (both T soil and T air ) in combination with M + P .A large part of the unexplained variance in the models with T + M + P appears to be due to location-effects (when the error of the multi-level models are evaluated with location as fixed effects, the misfit is in fact quite small).
The average ratio of RMSE V1 : RMSE C in the R S models was 1.5 and the ratio of RMSE V2 : RMSE C was 1.3.The ratio of validation error : calibration error measures the degree at which a model can generalize the results for a specific site (or experiment) to other locations or conditions.If the ratio is large, it indicates that the calibration data is unrepresentative or that the model for which the ratio is calculated is over-parameterized.In our experience, ratios smaller than 2 are quite acceptable and we therefore think that the calibration dataset is representative and that the models that were applied are not over-parameterized.The ranges for RMSE V1 (approx.0.5 to 2.8) and RMSE V2 (approx.0.45 to 3.5) were comparable, with the same four of the fifty-seven R S models (LMM and GLMM, using T , M and P , for both air and soil temperature) leading to very high values for RMSE V1 as well as for RMSE V2 .For the R S models, there was very a high correlation (> 0.99) between RMSE C , RMSE V1 and RMSE V2 .It should be noted that the validation was done only for the old vegetation.
The R H models produced relatively low RMSE V values for all combinations and formats (< 0.49).The R H models which produced the lowest RMSE V values were the GLMM format with T soil + M, the GLMM and Selsted format with T soil alone, and the LMM format with T soil + M. The ratio of RMSE V1 : RMSE C in the R H models was, on average, 1.15.There was a very strong correlation between RMSE C and RMSE V1 (Pearson correlation coef. of 0.997).

Model selection
Following the rationale described in the methodology to select the best predictive models, both the Selsted and GLMM models using only T soil or using T soil + M are selected as the best predictive models for R H .However, the GLMM models provide more realistic confidence bounds (by taking the property of repeated measurements in our data into account).Therefore, we preferred to use the GLMM model for predictive purposes.Furthermore, the RMSE, parameter values and predictions of the T soil and T soil + M models were similar.So the most parsimonious model (using only T soil as predictor) was selected to be used for further predictions.
For the R S models, the best predictive models were the Selsted and GLMM models, using only T soil or using T soil + P .As with the R H models, we have selected the GLMM models for prediction rather than the Selsted model.The differences between the GLMM T soil and GLMM T soil +P models (with respect to RMSE C , RMSE V1 and RMSE V2 ) were minor.So choosing the most parsimonious model for R S also leads to a model with T soil as the only predictor variable.The GLMM T soil models for R H and R S were used to predict soil respiration over the length of the study period (Table 5 and Fig. 8).

Autotrophic soil respiration
Autotrophic soil respiration was determined by subtracting the model predicted heterotrophic soil respiration results from the total soil respiration results in each vegetation community (R S − R H = R A ; Fig. 9).Soil R A was approximately zero on the middle and old communities in winter.The greatest R A was predicted to occur on the young community in the summer months, with a maximum in July when approximately 55 % of soil respiration was attributable to autotrophic sources.In this same time period, approximately 45 % and 37 % of soil respiration on the middle and old communities, respectively, was attributable to autotrophic sources.

Annual soil carbon loss estimates
Based on model predictions, annual R S C loss was significantly greater on the young community (650 g C m −2 yr −1 ) than either the middle (462 g C m −2 yr −1 ; p = 0.048) or the old (435 g C m −2 yr −1 ; p = 0.029) communities (Fig. 10).There was no significant difference between annual R S C loss on the middle and old communities (p = 0.39).The annual soil C losses from R A and R H were approximately equal in the young vegetation (50 % was R A ), but it was calculated that there was greater soil C loss from R H than from R A sources in both the middle and the old communities (30 % and 26 % was R A , respectively).The soil C loss was plotted against community age, using a "time for space" chronosequence approach to approximate changes in soil C loss over a 30 yr period.Year zero represents the bare soil which would be expected following a vegetation cutting cycle.In this case, all soil respiration would be expected to originate from R H , as no plant roots are respiring and the lack of vegetation would result in more variable soil temperatures, as observed in the bare Trenching plots.Therefore, soil C loss in year zero was predicted using the more variable trenched soil temperatures (350 g C m −2 yr −1 ).Soil temperatures were less variable under plant cover and so the untrenched temperatures were used in the model to predict annual R H C loss (322 g C m −2 yr −1 ) where plant cover was present.

Discussion
Carbon loss from soil respiration was greatest on the young community and root-associated respiration contributed approximately equally to the annual C sum as was contributed by microbial respiration.As the community age increased, the annual C loss from soil respiration decreased and this change was driven by the decreasing contribution of root respiration.
The following sections have been grouped around discussion of the soil respiration, of the trenching effects, the modelling process and finally a discussion of the annual soil C loss model predictions.

Soil respiration
Heterotrophic respiration rates were not statistically different between the three communities and this was consistent with the original hypothesis.In general, CO 2 effluxes from microbial decomposition are determined by the quantity and quality of available substrate, the soil temperature and other conditions that control decomposer activity (Kirschbaum, 2006).This was consistent with trial observations, as there was no difference between the quantity of available substrate in the different communities prior to trenching, that is, soil C stocks to 10 cm soil depth (Kopittke et al., 2012), and no soil temperature or soil moisture pattern differences between the trenched plots.However, the quality of the organic matter and recently deposited litter (prior to trenching) was not known.The proportion of lignin in the litter could be expected to increase with increasing community age, as woody stem growth increases with increasing plant age (Gimingham, 1985).Increasing the lignified material in organic mat-ter results in slower decomposition rates (Filley et al., 2008;Kalbitz et al., 2003).However, as no differences in respiration were observed, it is possible that the rapid decomposition of the labile organic matter masked any underlying differences (if indeed present) in the more recalcitrant pools.
The differences observed between total soil respiration on the community ages was not associated with heterotrophic respiration and therefore by elimination (R S − R H = R A ), was associated with autotrophic respiration.The greater total soil respiration on the young community indicated that the young Calluna plant roots were more actively respiring than on the middle or old communities.These higher rates corresponded to a higher P G and supported the hypothesis that the youngest plants, which were in a "net biomass gain" phase of growth (Gimingham, 1985), had the highest plant activity with greater allocation of carbon to the roots.
However, Calluna biomass was not the only contributor to P G .Mosses also contributed to P G , with almost double the moss biomass on the young community than on the middle or old communities.Although moss did not contribute directly to R A , as it lacks a root system, this mismatch in aboveground and belowground rates is likely to have introduced additional bias when including P G as a variable in the R S models.This study did not quantify the separate P G contributions of moss and Calluna.However, based on the preliminary data from a trial in May 2012, the young Calluna plants were approximately 2.5 times more photosynthetically active than the middle and old Calluna; therefore, P G would still provide a measure of the plant activity for each community.
The peak respiration values recorded in March 2012 corresponded to the first warm period in which air temperatures exceeded 15 • C, following from a severe frost (−20 • C) in February 2012.These extreme values were most likely associated with the death of fine roots and microbial populations, followed by the rapid recovery of microbial populations which lead to short-term fluxes of CO 2 from the soil (Matzner and Borken, 2008;Sulkava and Huhta, 2003).In addition, Calluna litter fall measurements on the old vegetation have shown peak fall rates occur approximately in January and old flowers are the dominant litter type (unpublished data from the adjacent long-term trial).This unlignified litter is likely to provide a rapidly decomposable energy source for microbial populations and may have contributed to the general CO 2 efflux peak that was observed in spring.
The observed total soil respiration rates were comparable to other Calluna heathland communities, such as in Brandbjerg, Denmark and a hydric Calluna heathlands in the Northern Pennines, England (Heinemeyer et al., 2011;Selsted et al., 2012).The mean summer total soil respiration rates in Brandbjerg ranged between 1.2 and 2.9 µmol CO 2 m −2 s −1 (2008 and 2006, respectively) and this was within the same range observed at Oldebroek in the summer of 2012.
Total soil respiration of other heathlands far exceeded the observations recorded at the Oldebroek study site.In the mesic heathland at Mols in Denmark, mean summer to-tal soil respiration rates were 16 µmol CO 2 m −2 s −1 in 2003 (Sowerby et al., 2008), which was approximately 5.8 times the mean summer respiration observed on the young community at Oldebroek in 2012.This large difference is most likely associated with the age of the vegetation and possibly differences in vegetation composition rather than soil differences.The soil type at Mols was similar, but the heathland experienced a heather beetle attack in 1999, which mainly resulted in Deschampsia regrowth and young Calluna plants (four years old).Similarly, total soil respiration on a hydric Calluna heathland at Clocaenog in Wales was also consistently greater in every season than the young community, even when the peak values of 5.6 µmol CO 2 m −2 s −1 (young community) and 7.6 µmol CO 2 m −2 s −1 (Clocaenog) were compared (Emmett et al., 2004).

Trenching effect
The soil temperature difference observed between trenched and untrenched plots is likely to be a function of the Calluna plants providing shade and the thick moss layer providing insulation at the soil surface.These two factors are hypothesized to have regulated soil temperature in the untrenched plots but not in the trenched plots where the aboveground vegetation had been removed.Since the temperature difference between trenched and untrenched plots was considerable (temperature on trenched plots minus temperature on untrenched plots: −1.0 • C, p = 0.016 in winter, and +1.4 • C, p = 0.021 in summer), total soil respiration and heterotrophic respiration cannot be directly compared, and R A cannot be obtained accurately by difference.Soil respiration models are capable of compensating such experiment related temperature artefacts to obtain meaningful partitioning into R H and R A .
Soil moisture patterns were also observed to differ between the trenched and untrenched plots, where the trenched plots were drier than the untrenched plots in non-rainfall periods.This is contrary to other studies in which trenching was observed to result in higher soil moisture than the control plots (Hanson et al., 2000).It is hypothesized that vegetation removal led to a loss of shade cover and this resulted in the organic layer and litter layer being exposed to greater evaporation rates.This hypothesis is supported by visual observations of a drier and cracked organic layer on the trenched plots.The respiration models being tested incorporated a soil moisture parameter so that any moisture effect could be assessed.

Model evaluation
All models followed generally the same pattern in the prediction of minimum effluxes in the winter, maximum effluxes in the summer and the highest autotrophic respiration for the young community (see Fig. 8, showing only the results for GLMM).However, the specific fit to the observations (as www.biogeosciences.net/10/3007/2013/Biogeosciences, 10, 3007-3038, 2013 summarised by RMSE varied between the different models (see Fig. 7).
The RMSE values for all models using T soil were consistently lower than those using T air .Additionally, the Selsted and GLMM models led to lower RMSE values and a lower spread in RMSE between the different vegetation ages than the other models.However, only the models using T soil alone were significant for all community ages.These results indicated that the complex parameterization of soil moisture and biomass effects in the Selsted model were not suitable for our site.This difference in model fitness may be due to site differences, such as the % grass cover or topsoil thickness, between the Brandbjerg heathland in Denmark (for which the Selsted model was developed) and Oldebroek.
Both for the R H and R S models, the RMSE values were very similar and highly correlated between the calibration and validation phases.Therefore, these models were considered stable and it can be assumed that the model predictive uncertainty was mainly due to parametric uncertainty.Also, the very high correlation between model prediction errors for calibration and validation indicates that the calibration and validation data contain data with a very similar information content.When the model misfit is analysed in greater detail (see Appendix D), several structural deviations of the residual are seen over time (the model residual is not uncorrelated but contains information which is not captured by the model).This misfit is not apparent with regard to temperature.In our view, the most plausible explanation for the structure in the model residuals is that one or more important covarying variables are lacking in the models that were parameterized so far.
From the variables assessed and available for inclusion in our model selection process, very similar fits of the observed data was provided by models using soil temperature (possibly in combination with soil moisture or plant activity).However, the model with soil temperature as only variable contained less parameters and was therefore preferred for predictive purposes.The application of only a temperature function to model soil respiration data has previously been questioned since, as already discussed, other factors such as soil moisture limitation of microbial processes and the C allocation via plant roots are all reported to influence soil respiration rates (Davidson et al., 2006;Rustad et al., 2000).However, our results indicated that soil moisture and plant activity (Calluna biomass, P G , microbial biomass and root biomass) were not significant variables for our site.To examine this further, it is first considered if it is possible that some of the measured variables would have been significant, if the data had been measured differently.Secondly, other variables are considered that have been used in similar soil respiration studies and may have improved model fit.
Soil moisture was not found to be a significant variable.Whilst the use of a soil moisture model may have introduced additional noise or bias, our interpretation of results is that it did not introduce additional artefacts.For the R H model, soil moisture has been shown to impact microbial respiration (and therefore R H ) only at extremely low water contents when desiccation stress becomes important for microbial substrate supply (Davidson et al., 2006).It is possible that in our study the soil did not reach these desiccation stress levels, thus resulting in a non-significant soil moisture parameter for the R H models.For the R S model, Calluna plants appear to be resilient to water stress and heathlands can withstand quite severe summer droughts, if annual rainfall is high enough to compensate for the drought (Loidi et al., 2010).Additionally, the Oldebroek heathland is established on a free-draining, sandy soil that has relatively low stored soil moisture in the mineral soil.The majority of the Calluna roots were identified within the organic layer of the soil and this is also where the largest proportion of the soil moisture is stored (see Table 1).However, continuous soil moisture measurements in the organic layer are very difficult due to instrumentation constraints (Schaap et al., 1997).Because of this, it is likely that a large proportion of the soil respiration response to reductions in soil moisture occurred in the organic horizon, and this was not able to be quantified with the current technology.Therefore, continuous soil moisture measurements in the organic horizon may have improved model fit.
Other variables from published soil respiration models that could be considered have included using relative PAR with soil temperature and soil moisture (Caquet et al., 2012).In our study, PAR was included in the initial model screening process as a single predictor variable and as a predictor variable together with temperature.However, neither of these models resulted in a better fit than soil temperature alone and therefore, PAR was not included in further model testing.
Alternatively, another plant variable which has been considered in other respiration models is the rate of litter decomposition (Kutsch and Kappen, 1997;Kutsch et al., 2010).However, it is unclear from these studies whether the addition of litter decomposition to the soil temperature and moisture model resulted in a better model fit, as the parameter significance was not reported.Soil temperature has been found to generally have a good relationship with organic matter decomposition rates (Davidson and Janssens, 2006) and therefore it is hypothesized that a litter decomposition variable would not explain significantly more variability than already explained by soil temperature.Other plant litter variables, such as litter fall rates, are also often included in dynamic models as they provide an important feedback into the carbon cycle and substrate available for decomposition (Keenan et al., 2012).Litterfall results were not available for the young and middle communities, although litter data was collected on the old community validation plots between March 2011 and February 2012 (unpublished results).The maximum litterfall rate occurred in January (8 g m −2 month −1 ) and the minimum in February (2 g m −2 month −1 ) with gradually declining rates recorded from March to November (7 to 5 g m −2 month −1 ).This pattern did not correspond to Root maintenance (as a function of root nitrogen concentrations) and root growth have also been included in other soil respiration models.In a study in Tennessee USA, a model with root variables was able to describe more of the biological dynamics than the other models tested although it was still not capable of capturing all the data variation across the different study treatments (Chen et al., 2011).Root dynamics provides a direct measure of root activity and, if it had been measured at the Oldebroek site, may have explained more variance that the photosynthetic rates.
A further layer of complexity to the discussion is that model results may be influenced by a suboptimal measurement integration volume or integration time, as well as the alignment in space and time of different measurements.Problems of this kind ("scale problems") are common in the natural sciences and are an important source of model error, thus are considered as the most important challenge in ecology (Blöschl and Sivapalan, 1995;Wiens, 1989).An example of a data alignment problem in our study was the collection of soil respiration measurements on different days than the photosynthesis measurements, which required intermediate data processing for photosynthesis (viz.Fig. 3).Also, soil temperature was measured at a depth (5 cm), whereas the soil respiration was an integral measurement over a soil column (e.g.Reichstein and Beer, 2008).There may also be a lag time present within the data, where plant growth on one day does not immediately correspond to root respiration (Gomez-Casanovas et al., 2012;Kuzyakov and Gavrichkova, 2010), which our non-continuous data would not have been able to detect.These trial design aspects may have resulted in small data mismatches and is possible that the model calibration and validation results would have improved if the resolution and alignment of the data had improved.
The model selection process resulted in a model that used T soil alone, which is arguably the simplest variable.However, this "simple" result does not negate the use of a detailed selection procedure, as the process also highlighted that the current variables measured were not adequate to model all the variation observed in the R S (and therefore R A ) data.This is an important outcome of this study, as many studies include processes that are theoretically associated with soil respiration but the model variables are not assessed for significance and may not explain any additional data variation.This practice leads to a publication bias (Dieleman and Janssens, 2011).The use and reporting of full data pre-processing and modelling workflows that apply sound scientific procedures, which also report the "negative" or "less interesting" results, helps to avoid such a publication bias.
It is not possible to measure all ecosystem processes on an experimental trial due to practical constraints and it is not always possible to know which measured process would improve the model fit, although pre-planning field experiments based on the models we wish to use may assist in this process.This finding supports the discussion presented by Subke and Bahn (2010) on the ability to use the immeasurable to predict the unknown.In our study, although one or more important covarying variables were lacking (and model fit would have improved had these been measured), it is worth considering that soil temperature was likely to also be related to seasonal plant activity and may simply be the overwhelming driver of soil respiration in this system.Therefore, in the absence of other variables, the T soil variable was sufficient to explain most of the seasonal variation of R S .Similar findings have been reported in other studies, where site differences in R S were largely determined by plant productivity but since both R S and P G fluxes increased with temperature, it was concluded that the soil temperature typically sufficed to explain R S in non-drought ecosystems (Bahn et al., 2010a;Janssens et al., 2001;Reichstein et al., 2003).
We think that the findings from this empirical study (on the basis of static models) can also be used to investigate or test dynamic soil respiration models (which are typically parameter-rich and often model more than only soil respiration in isolation).First of all, there are many dynamic soil respiration models which effectively contain a respiration equation similar to those used in this study (e.g.Keenan et al., 2012;Kutsch et al., 2010).In those models, the more complex equation could easily be replaced by a similarly appropriate, but simpler, equation, which leads to less parameters requiring calibration or more stable model behaviour.Otherwise, if the way in which respiration is modelled is incompatible with the static equations in this study (and results cannot directly be translated), discussion needs to be given to the model evaluation methods used for soil respiration research.The process by which different models or model components are evaluated on calibration and validation datasets (i.e. the method promoted in this study) deserves attention -not because it is new, but because it is currently uncommon in this soil respiration research area.A lack of critical model evaluation limits progress.
In nature, many interactions can occur and when our field trials don't test these interactions, it is not possible to incorporate them into long-term model predictions.Therefore, it is necessary to develop field trials which incorporate this increased complexity, as suggested by Dieleman et al. (2012).However, if early consideration is not given to the models that we later want to fit to the data (and the data required to rigorously test the models), then increasing the complexity of field experiments will not necessarily provide us with better predictions of these interactions.Therefore, attention should be given to the trial layout, variable selection, measurement intensity and model selection process prior to the start of a trial to determine if they will provide the appropriate data for www.biogeosciences.net/10/3007/2013/Biogeosciences, 10, 3007-3038, 2013 model predictions.Consideration also needs to be given to the cost associated with obtaining the appropriate measurements, in terms of collection frequency, method accuracy and overall outcomes of the project.In some cases, it may be that using a proxy such as soil temperature (or even air temperature for rough estimations) with the soil respiration observations is a suitable substitute in models in the absence of suitable and significant variables.

Annual soil C loss and links to global change
Our model interpolations identified an annual R S C loss that was at the lower end of the range identified on the Danish heathland ecosystem of 672-719 g C m −2 yr −1 (Selsted et al., 2012).To place this within a broader European context, the heathland soil respiration is within the same range as temperate forest ecosystems, which have been reported between 430 g C m −2 yr −1 (Belgium) and 859 g C m −2 yr −1 (Germany) (Bahn et al., 2010a;Khomik et al., 2009;Raich and Schlesinger, 1992).In contrast, the heathland is at the lower end of the scale for annual soil respiration in comparison to temperate grasslands, which ranged between 729 g C m −2 yr −1 (Germany) and 1988 g C m −2 yr −1 (Switzerland) (Bahn et al., 2010a).
The study also identified a change in soil respiration with an increasing age of heathlands.Soil R A provided the largest change over time, from a complete absence on bare soil to a maximum at the 12 years and then decreasing up to the maximum studied age of 28 years.A similar relationship between soil respiration and vegetation age has been previously found for forest stands, where the younger stands had significantly higher respiration rates than the more mature sites (Saiz et al., 2006;Wang et al., 2011).
Within the last 50 yr, the cutting, burning and grazing cycles on heathlands have not occurred as frequently or as regularly as during the intensive agricultural periods of past centuries (Webb, 1998).Management of heathlands is required to maintain these cultural landscapes and in the past this management occurred on a 3-4 yr cycle (Webb, 1998).Currently, this cycle length has extended or is non-existent (Diemont and Heil, 1984;Wessel et al., 2004).From the perspective of optimizing C uptake and minimizing C output, having an understanding the C dynamics of these ecosystems allows us to determine the optimum time to cut the vegetation, thus contributing to global C emission mitigation measures.

Conclusions
This study showed that soil respiration models are required not only for prediction but also to correctly interpret field measurements of total and heterotrophic soil respiration.That is, the effect of experimental design on soil temperature (and subsequently soil respiration) can be corrected for through the use of models to predict both R S and R H over the same temperature range, thus allowing R A to be calculated.
Soil temperature appeared to be the best predictor for R S as well as R H in the heathland ecosystem of our study, while soil moisture, photosynthesis, plant biomass, PAR, root biomass and microbial biomass did not significantly improve model fit when added to soil temperature.However, the temperature model still contained temporal autocorrelation in the residuals, suggesting that alternative variables which were not considered in this study (like litter decomposition rates or root growth) could be important predictors.For our future experimental work, this model based on soil temperature may act as a null model against which the performance of other models can be compared.

Details of the P G measurements
The gross photosynthetic rate provided a measure of photosynthetic activity for the three heathland ages.The gross photosynthetic rate (P G ) was calculated as the net ecosystem exchange (NEE) rate of CO 2 flux minus the ecosystem respiration (ER) rate of CO 2 flux (µmol CO 2 m −2 s −1 ).This photosynthetic rate has a negative sign.A loess smoother curve was applied to the photosynthesis data to obtain daily estimates of plant activity.
The CO 2 fluxes of the vegetation were measured a LI-6400 infrared gas analyzer (LI-COR, Lincoln, NE, USA) attached to a 288 L ultra-violet light transparent Perspex chamber (60 cm × 60 cm × 80 cm) using the method described by Larsen et al. (2007).The chamber was installed with a fan as well as a soil temperature probe (LI-6400-09 temperature probe) and a PAR sensor (LI-COR quantum sensor).
Three permanent sampling locations were selected in each vegetation age.A metal base frame (60 cm × 60 cm) was permanently installed using small, narrow sandbags to provide a seal between the frame and the soil surface and fixed with metal pins.Measurement of CO 2 fluxes commenced immediately prior to the Perspex chamber being placed on the frame so as to capture the point at which the chamber was sealed and NEE occurred entirely within the chamber.The LICOR measurement program ran for 180 seconds, however, the results obtained while the chamber was being fitted were later discarded so that only data obtained from the sealed chamber (approximately 150 s) were utilized for calculation of NEE rates.After the NEE measurements, the chamber was vented and measurements of the ER rate were obtained by covering the chamber with a fitted blackout-cloth, in which the outer layer was white and the inner lining was black, to minimise any heating effect within the darkened chamber.
In most cases, NEE decreased from the first to the third minute of measurement, indicating an effect of the chamber by the decreasing CO 2 concentration as photosynthesis progressed.Therefore a linear regression did not provide a good fit for all measurements.To overcome this problem, the HMR procedure was used (Pedersen et al., 2010).This procedure was developed for soil-atmosphere trace-gas flux estimation with static chambers and tests the fit of both loglinear and linear regression models to the NEE or ER data at each measurement.If linear regression provided the best fit, the flux value was determined by the slope of the regression line.If non-linear regression gave the best fit, the flux was determined by the slope at t = 0 s.The HMR procedure is implemented in an R-package (Pedersen, 2011) and this implementation was used in our study.

Details of the soil moisture model
The soil moisture model used in this study is a zerodimensional finite difference model using a daily time resolution of rainfall data and air temperature data as model inputs.It was constructed and calibrated on approximately one year of observed soil moisture, rainfall and temperature data for 12 individual soil moisture sensors.The model comprises the following equations: In the equations, t refers to a day.Equation (B1) calculates drainage (Drain t , in mm day −1 ) as a linear reservoir with soil moisture (Smoist t−1 , in mm) above a threshold (depth • fc) as the driving force.Drain t refers to the drainage of soil moisture from the soil layer under consideration (i.e. the top of the mineral soil down to depth mm); Smoist t−1 refers to the soil moisture in the soil layer under consideration, and depth, www.biogeosciences.net/10/3007/2013/Biogeosciences, 10, 3007-3038, 2013 fc (field capacity, as a fraction of the soil volume) and df (drainage fraction) are model parameters.The depth parameter is set to 100 mm, while the values for fc and df were identified by model calibration.Equation (B2) calculates the soil moisture available for evapotranspiration (AvSmoist t , in mm) and the parameter wp (as a fraction of the soil volume) represents the wilting point below which only a negligible rate of evapotranspiration occurred.The value for wp was found by model calibration.
Evapotranspiration (ET t ) is calculated in Eq. (B3).Evapotranspiration is a modelled linear reservoir with either the air temperature or the available soil moisture as the driving force, depending of which factor is limiting.The parameter tf is set to 1 mm ( • C) −1 , and the value for the parameter ef was identified by model calibration.
The effective rain, i.e. the rainfall which enters the soil layer under consideration (EfRain t , in mm), is calculated in Eq. (B4).EfRain t is proportional to a soil saturation factor which contains two parameters: soil porosity (poros) and a rainfall factor (rf).The porosity is calculated by taking the maximum observed soil moisture content over the measurement period, while the rainfall factor is calculated by model calibration.
In Eq. ( B5), an update of the soil moisture is calculated by a balance equation, whereby it is assumed that any rainfall which cannot be stored in the soil layer under consideration is lost as surface runoff.
The water balance model thus contains eight parameters, three of which have fixed values (depth = 100 mm, poros = max all t (Smoist t /depth), and tf = 1 mm ( • C) −1 , and five of which were found via calibration (df, ef, fc, rf and wp).Calibration was undertaken by minimizing the root mean squared error between observed and predicted soil moisture, using the optimization routine by Byrd et al. (1995), as implemented in the standard R function "optim".
The fit of the soil moisture model for the different treatments is shown in the diagnostic plots of Fig. B1.The plots illustrate that there is still quite some room for improvement in the soil moisture model.For each of the cases, the explained variance in the observed versus predicted plot is approximately 0.7.
To further check the adequacy of the modelled soil moisture, we compared the cross-correlograms for daily air temperature and modelled soil moisture versus air temperature and observed soil moisture.These appeared to agree well and, on the basis of this, we concluded that the modelled soil moisture did not lead a different correlation structure with respect to air temperature than observed soil moisture.
In subsequent modelling, we treated the result of the soil moisture model as if it were an observation for an individual observation plot.
Table C2.The measures of model fit for a selection of R S models using the Calibration data, Validation Type I data, Validation Type II data, Validation Type III data.See Table C1 for Type definitions and Table 3 in the main article for parameter explanations.Calibration Type I models which included non-significant parameters are not shown in the table.Grey italicized text indicates that not all variables in this Calibration Type II model were significant.

Fig. 3 .
Fig. 3. Measures of plant activity for the young, middle and old communities, showing (a) mean Calluna biomass (kg m −2 ) obtained in April 2011 during trenching activities (n = 12); and (b) C uptake by Photosynthesis (µmol CO 2 m −2 s −1 ) obtained between August 2011 and August 2012 (n = 9) with observations represented by symbols and the mean curves (loess curves) represented by lines.The SEM bars are shown for each mean value.The three means are not significantly different at the 0.05 level (adopting a Bonferroni correction).

Fig. 4 .
Fig. 4. Soil Respiration (µmol CO 2 m −2 s −1 ) on three ages of vegetation for the (a) total soil respiration as represented by the untrenched plots; and (b) heterotrophic soil respiration, as represented by the trenched plots from September 2011-August 2012 (n = 4 per age per sampling event).For plot (b), the young community SEM bar in March 2012 extends outside the graphical boundaries to 6.79 µmol CO 2 m −2 s −1 .

Fig. 5 .
Fig. 5. Environmental parameters for September 2011-August 2012, showing (a) hourly temperatures ( • C) of the air at 20 cm above ground surface and of the soil at 5 cm below ground surface; and (b) mean daily soil moisture (m 3 m −3 ) at 5 cm below ground surface for the trenched plots and the untrenched plots.Periods of frozen soil moisture are indicated by shading.
18 • C) (t = −4.11,p = 0.021, 3 plots and using averages of soil temperature per plot for the complete winter).Mean air temperature at 20 cm above ground surface was 3.0 • C in winter and 15.7 • C in summer for both trenched and untrenched plots.Soil moisture was significantly different between the trenched and untrenched plots with lower soil moisture values observed on the trenched plots than the untrenched plots in non-rainfall periods (Fig. 5b; 1.5 vol.% less soil moisture on trenched, plots, t = 31.78,p < 0.001 based on 3 plots).

Fig. 6 .
Fig. 6.A comparison of the trenched plots (n = 12) and untrenched plots (n = 12) for (a) mean microbial C biomass (mg C g C −1 ) in the organic horizon and (b) mean root biomass (g m −2 ) in the summed (organic + 0-5 cm mineral) horizons shown for the three ages of heathland vegetation.Different letters above the bars represent statistical significance and the SEM bars are shown for each mean value.

Fig. 7 .
Fig. 7. Comparison of RMSE C values for models of (a) total soil respiration data (untrenched plots) and (b) heterotrophic soil respiration data (trenched plots).The models tested are listed on the left side of the figure.The explanatory variables within each model are listed on the y axis and are abbreviated as T = temperature (soil or air • C as indicated), M = soil moisture, B = relative biomass, P = relative photosynthesis.The " * " indicates that all model parameters were significant for one of either "Y" (young), "M" (middle) or "O" (old) vegetation community models.The SEM bars on the total soil respiration means were calculated from the RMSE C s of the three community ages.SEM bars could not be calculated for the heterotrophic models.Two means are outside the plotted range: GLMM T air + M + P (1.326) and GLMM T soil + M + P (1.390).

Fig. 8 .
Fig. 8.Predicted and observed soil respiration on the young, middle and old community "untrenched" plots (total soil respiration: µmol CO 2 m −2 s −1 ) and the "trenched" plots (heterotrophic soil respiration: µmol CO 2 m −2 s −1 ) calculated with the GLMM model using the T soil explanatory variable.The observed values from 21 March 2012 are excluded from these plots.

Fig. 9 .
Fig. 9. Predictions of soil respiration for (a) young community, (b) middle community and (c) old community calculated using the GLMM model with the T soil explanatory variable.

Fig. 10 .
Fig. 10.Estimated total annual soil respiration (R S ), heterotrophic soil respiration (R H ) and autotrophic soil respiration (R A ) as predicted by the GLMM T soil model.Year 0 is represented by respiration from bare soil, year 12 by the young community, year 19 by the middle community and year 28 by the old community.Mean prediction values are provided with the bars representing the 95 % Confidence Intervals.
Fig. B1.Observed versus predicted soil moisture for the four different vegetation communities in this study.

Table 1 .
Description of the Oldebroek trial location.

Table 2 .
Description of the data used for model calibration and validation.

Table 3 .
The models to estimate R S and R H in µmol CO 2 m −2 s −1 .The explanatory variables are T (models using air temperature at 20 cm above ground surface and soil temperature at 5 cm below ground surface are evaluated) and M, B, P as defined in Eqs.(2) to (3).The model parameters are R 0 , k, a, b and c and the units vary per model.R 0 is always in µmol CO 2 m −2 s −1 .Parameter k is in • C −1 for Selsted, GLMM and GLMM2 models, and in µmol CO 2 m −2 s −1 • C −1 for LMM and LMM2.The parameters a, b and c are dimensionless for Selsted, GLMM and GLMM2 models and are in µmol CO 2 m −2 s −1 for the LMM and LMM2 models.

Table 4 .
The residual standard deviation with RMSE C and RMSE V (Validation Type I and II) values for a selection of GLMM models.Only models in which all variables were significant are shown.See Table2for definition of validation types, Table3for model definitions and Table5for parameter values of the models shown here.Results for cross-validation and with additional error metrics are provided in Appendix C.

Table 5 .
Optimal parameter values of R S and R H models for the young, middle and old communities (GLMM T soil , T soil + M and T soil + P Models), with 95 % confidence intervals for the parameters in brackets.See text and Table3for parameter explanations.For R H , model T soil + P does not exist.The values of the different parameters are given with different numbers of significant digits to reflect the uncertainty in the corresponding variable.The bold cells indicate parameter values that are not significant at the 0.05 significance level.Diagnostic plots for the models listed in this table are shown in Appendix D.

Biogeosciences, 10, 3007-3038, 2013 www.biogeosciences.net/10/3007/2013/ the
observed soil respiration rates and suggests R H is more closely associated with temperature than with litterfall patterns.However, if sufficient litterfall data had been available for inclusion as a variable with soil temperature, it may have improved R H model fit by explaining additional data variation.