Interactive comment on “ Sensitivity analysis of a wetland methane emission model based on temperate and Arctic wetland sites

Comments This manuscript evaluated the sensitivity of a well-known CH4 emission model by applying the GLUE (Generalized Likelihood Uncertainty Analysis) methodology. The response of CH4 flux to model parameters (microbial reaction rate, vegetation and soil properties) is well addressed. This type of study is necessary for estimating regional and global wetland CH4 emission and is of scientific significance. I suggest publishing this paper after minor revision. (1) Determination of CH4 flux was made by chamber/gas chromatograph system before May 2004 (Horstermeer and Ruwiel sites). Thereafter CH4 analysis was performed in the field using an Innova 1312 photoacoustic gas analyzer. Please make sure the CH4 fluxes determined using these two systems are comparable. Otherwise, it would obscure the sensitivity analysis. (2)


Introduction
Together with water vapour and carbon dioxide (CO 2 ), methane (CH 4 ) is an important greenhouse gas, because of its strong global warming potential of 23•CO 2 on a 100-year time scale.The atmospheric mixing ratio of CH 4 has increased with 151±25%, since pre-industrial times.About 60% of the global CH 4 emission is of antropogenous origin.From the natural sources (wetlands, termites, oceans, methane seeps and hydrates), the wetland environments are the major natural source of atmospheric methane (IPCC, 2001).Moreover, the atmospheric methane concentration appears to be strongly linked to climate change during the last 800 000 years (Loulergue et al., 2008).
Understanding of feedbacks between climate and wetland CH 4 emission, in particular in boreal/arctic regions, is a problem for predicting future climate change (Denman et al., 2007).Wetland CH 4 emission is also influenced by land management (e.g., Van Huissteden et al., 2006;Hendriks et al., 2007).With the need to reduce greenhouse gas emissions, the relation between wetland CH 4 emission and wetland management may become an important question in the future.Predictive models may contribute to a better understanding of feedbacks between climate and CH 4 emission, or the effects of wetland management on CH 4 emission (e.g.Petrescu et al., 2009b;Hendriks et al., 2007;Petrescu et al., 2009a).However, to reduce modelling uncertainty extensive Published by Copernicus Publications on behalf of the European Geosciences Union.sensitivity testing and uncertainty analysis is required, in particular when models are scaled up from a local to regional or global scale.As yet, existing CH 4 emission models have not been subjected to rigorous uncertainty analysis going beyond simple model-data comparisons.Here, we present an uncertainty analysis of a wetland CH 4 emission model, based on the GLUE (Generalized Likelihood Uncertainty Estimation) methodology (Lamb et al., 1998;Beven, 2001 and references therein).
Methane emission from wetland soils is essentially the net result of a balance between CH 4 production by methanogenic bacteria in anaerobic soil zones, and CH 4 oxidation by methanotrophic bacteria in aerated soil zones and in plants.Several process models of wetland soil methane emission have been designed (Walter, 2000;Segers and Leffelaar, 2001;Granberg et al., 2001;Segers et al., 2001;Wania, 2007).These papers and references therein give an overview of the processes involved.CH 4 is generated by methanogenic bacteria in anaerobic parts of the soil, when other electron acceptors for organic matter oxidation are exhausted or unavailable (nitrate, sulfate, Fe and Mn oxides).The substrate for methanogenesis is mainly derived from labile organic compounds, produced by the roots of the wetland vegetation.In wetlands very rapid transfer (1-2 days) of photosynthesis products to CH 4 has been observed (King and Reeburgh, 2002).The two major reaction pathways for methanogenesis are CO 2 reduction and acetate splitting (e.g., Bréas et al., 2001).transport of CH 4 from soil to atmosphere occurs along three pathways: diffusion in soil pores, bubbles rising to the surface (ebullition), and transport through plant roots and stems.In particular during diffusive transport in the soil and plant-mediated transport, CH 4 is subject to oxidation by methanotrophic bacteria (e.g.Whalen et al., 1996;Raghoebarsing et al., 2005;Van Huissteden et al., 2008).
Several models have been developed to model methane fluxes either at a plot scale (Walter, 2000;Granberg et al., 2001;Segers and Leffelaar, 2001;Segers et al., 2001) or on a larger scale, ranging from regional to global scale (Petrescu et al., 2009a, b), usually coupled to climate models (e.g., Cao et al., 1996;Kaplan, 2002;Gedney et al., 2004;Van Huissteden, 2004;Wania, 2007).Process-based modelling of methane fluxes from wetland environments is difficult because of the complicated interactions between soil biochemistry, vegetation and soil chemical and physical processes; most of these processes require parameters that are difficult to measure and generally not available (e.g., Walter, 2000).In fact, many 'process' based models therefore contain rather course bulk parameterizations of key processes, and the more detailed the process formulation in a model, the higher the parameter requirements of the model.This may result in overparameterized models, containing parameters that do not contribute significantly to a better fit of the model to field data.Careful parameter sensitivity analysis is therefore necessary to assist model improvement.Usually CH 4 emission or soil respiration models are tested only by varying a few key model parameters and input data and determining the resulting variation in model output, without further analysis of the model uncertainty.Van Huissteden et al. (2006) tested the PEATLAND-VU model on sensitivity to climate and water table input data and a limited number of model parameters.Granberg et al. (2001) consider in the data-model comparison also the standard error of the data, which can be large compared to the measurements in case of CH 4 emission.Wania (2007) tests the senitivity of her model by regressing output on parameter values for a range of model parameters.Berritella and Van Huissteden (2009) tested a large scale CH 4 flux model with varying complexity.However since they modelled paleo-wetland CH 4 fluxes a rigorous data-model comparison was impossible.
Plot-scale models have the advantage that they can be validated against site CH 4 flux measurements under a variety of conditions (e.g., Walter, 2000;Granberg et al., 2001;Petrescu et al., 2008) and can make use of detailed on-site measurements of key parameters of soil physical and chemical conditions.Larger scale modelling of CH 4 fluxes always requires aggregated and simplified information on vegetation and soil and are more difficult to validate.However, to properly understand interactions of wetland CH 4 emission with climate or wetland management, large scale modelling of these emissions and coupling to climate or hydrological models is highly important (Petrescu et al., 2009b).For that purpose it is necessary to know at which level of detail processes need to be modelled to represent the interactions between climate or management correctly.CH 4 fluxes are known to be spatially highly variable on a small scale (e.g., Van Huissteden et al., 2005;Hendriks et al., 2009a).Water table position is the most important variable, but also difference in vegetation and soil properties have been shown to be influential (Hendriks et al., 2009a).
An approach that can give information on the required model complexity for large scale modelling, is to test the parameter sensitivity of the more detailed, plot-scale models (Beven, 2001).If a model parameter has a strong influence on the modelled fluxes on the plot scale, it is likely that it also has a large influence in an upscaled version of this model.Depending on model structure, this may hold also for other models that use the same or similar parameters.In that case large-scale modellers should focus on obtaining correct values of this parameter, or at least obtaining a good proxy estimate, for example from remote sensing of vegetation cover.Conversely, model parameters that do not contribute significantly to model-data fit on the plot scale will neither contribute to large-scale modelling of fluxes.To distinguish influence from specific local conditions, models need to be tested with several data sets, from spatially and environmentally different locations.In particular for largescale modelling also simpler, reduced complexity modelling approaches should be considered (e.g., Berritella and Van Huissteden, 2009).For model testing the following questions should be asked: 1. What is the uncertainty in the model results, given the uncertainty in input parameters?
2. What is the sensitivity of the model results to variations in parameter values, in particular those that vary spatially?
3. What is the interaction with other model parametersdoes variation in one parameter affect the sensitivity of another parameter?
4. Can optimum parameter sets be found?Are there unique parameter sets or multiple sets that produce realistic model simulations?
A widely used plot-scale process model is that of Walter (2000); modified versions have been incorporated in PEATLAND-VU (Van Huissteden et al., 2006;Petrescu et al., 2008;Petrescu et al., 2009a), and WetlandDNDC (Zhang et al., 2002).The model of Walter (2000) includes methane generation by bacterial consumption of labile soil organic matter, bacterial methane oxidation, and transport of methane to the atmosphere by ebullition, diffusion and fluxes through plants.This model contains several site-specific parameters that are difficult to quantify properly.We tested the version of this model included in PEATLAND-VU (Van Huissteden et al., 2006), which is a slightly modified version of the Walter (2000) model (for a description of the model and modifications of the original model of Walter (2000) see below).We use the GLUE (Generalized Likelyhood Uncertainty Estimation) methodology (Lamb et al., 1998;Beven, 2001 and references therein), with validation data from three different sites, including a natural and a managed temperate wetland and a permafrost tundra wetland.
GLUE is an approach that includes a combined evaluation of model result uncertainty and parameter sensitivity.It has been applied extensively to hydrological models (Beven, 2001) and overcomes several problems that usually arise with model calibration and sensitivity analysis in complex environmental models.In simple cases, for model calibration the outcome of a model for a given parameter set is compared to observation data on the modelled system using an objective function.This objective function indicates with a goodnessof-fit measure to what extent the model results agree with the observed values.The parameter set that yields model results with the best agreement between model and observations (lowest value of objective function) is chosen as the optimal parameter set.This is a straightforward approach when a clear optimal value of the objective function exists and the number of model parameters is small and their value range well constrained.However, in complicated systems like hydrological and soil systems the number of relevant parameters that have to be considered may be prohibitive.Also there may be considerable interaction between the parameters.In such cases widely different parameter sets may yield similar model results (equifinality).This situation also has been observed for the PEATLAND-VU model (Van Huissteden et al., 2006).For instance, PEATLAND-VU can generate CH 4 flux time series with a good model-data fit using different combinations of microbial CH 4 production rates and plant oxidation rates.
GLUE makes no assumptions about the nature of the optimal parameter set of the model.The approach specifically recognizes the occurrence of non-unique solutions of model optimization.GLUE is based on a large number of model simulations with randomly generated parameter sets.Each parameter can vary within a specified range; multiple parameters are changed at each model run.For each run also an objective function value is generated.Although there will be one most optimal value of the objective function among the simulations, there may be many that are nearly as good and may represent also valid parameter sets.By studying the distribution of the objective function values for all model simulations that are well-behaved, not only optimal parameter sets can be found but also conclusions can be drawn on the parameter sensitivity, parameter interaction and predictive uncertainty of the model (Freer and Beven, 1996;Lamb et al., 1998;Beven, 2009).

The model
PEATLAND-VU is a process-based, plot scale model of CO 2 and CH 4 emission from peat soils at various climate scenarios.The model has been used by Van den Bos et al. (2003) and Petrescu et al. (2009b) for regional scale simulation of CO 2 and CH 4 fluxes in the Netherlands.Van Huissteden (2004) and Berritella and Van Huissteden (2009) employed the model for simulation of paleo-CH 4 fluxes from wetlands in Europe during the last glacial.Petrescu et al. (2009a) have used the model for global scale simulation of present-day boreal and arctic wetlands, by coupling the model to a global hydrological model.PEATLAND-VU consists of four submodels: a soil physics sub-model to calculate temperature, water saturation and ice content of the soil layers, a CO 2 submodel, a CH 4 sub-model and a soil organic matter (SOM) production sub-model.For a complete model description we refer to Van Huissteden et al. (2006).Here only recent modifications are discussed.The CH 4 sub-model is based on Walter (2000).The model of Walter (2000) (Van Huissteden et al., 2006).
According to Walter (2000), the production factor for methane from labile organic compounds in the soil (termed R 0 in Walter's model description) should be regarded as a tuning parameter to adapt the model to different sites and climatic conditions.In PEATLAND-VU R 0 has been made dependent also on soil pH, using an empirical linear relation derived by Dunfield et al. (1993) (Van Huissteden et al., 2006).Additionally, the model appears sensitive to parameters influencing the soil-atmosphere CH 4 transport through plants (Van Huissteden et al., 2006;Petrescu et al., 2008), in particular to the fraction of CH 4 that is assumed to be oxidized during plant transport P ox ).Also parameters related to primary production and distribution of labile organic compounds in the SOM production submodel potentially influence modelled fluxes: the Net Primary Production (NPP), and the fraction of NPP transferred to labile organic compounds.This fraction in turn is determined in PEATLAND-VU by the fraction of below-ground organic production f roots and the fraction of f roots that is transferred to rhizodeposition (dead root material and exudates, f dep ).
With respect to the the model description by Van Huissteden et al. (2006) and Walter (2000), modifications have been made to the model.Field observations suggest that after a dry period in wetland soils, a time lag occurs between a rise of the water table at the onset of rain and the increase of CH 4 fluxes (Hendriks et al., 2007;Hendriks et al., 2009a).This time lag is caused by the decrease of redox potential in pore water due to progressive oxidation of labile organic compounds (e.g.Segers et al., 2001).Within the PEATLAND model lowering of redox potential is not explicitly modelled by modelling the successive redox processes since it would require addition of extra soil chemical parameters.However, it can be mimicked by assuming an exponential increase of CH 4 production to its maximum rate, depending on the time lag l sat (days) after the onset of completely water saturated conditions in a soil layer that has previously been unsaturated, and the availability of labile organic carbon C labile (µmol C kg −1 ).This is modelled as follows: where k delay (range 0.01-0.05) is a constant defining the restoration rate of maximum CH 4 production.Second, our field observations at the Horstermeer site (see below) suggest that in dense, partly oxidized fen peats CH 4 production also may occur above the water table at low stands in summer, presumably due to the presence of anaerobic microsites in the soil.Wagner and Pfeiffer (1997) have found viable methanogenic bacteria above the water table in similar marsh soils.In PEATLAND-VU CH 4 production above the water table is modelled as a fraction (p anaer ) of the the production below the water table.This fraction depends linearly on the pore water saturation fraction.The slope of this relation, f anaer , is the model parameter that determines p anerobe : 3 Study sites and field methods

Sites
Horstermeer (52 • 14 30 N, 5 • 5 E) is located SE of Amsterdam, in a drained lake.The water level in the ditches is at approximately 3.5 m below sea level, and up to 2 m below that of surrounding areas.The area is subject to strong seepage, in particular in the drainage ditches.The soil consists of 2 m of clayey gyttja (organic lake sediment), erosively overlying eutrophic fen peat on Pleistocene sand.Until 1997 the area was a grazed pasture, thereafter the water level has been raised to 0.2-0.4m below the surface, to create a wetland nature reserve.The present vegetation, a degraded pasture, is not harvested or managed otherwise.Dominant species in the wetter parts are Holcus lanatus, Equisetum palustre, Glyceria maxima and Typha latifolia; dryer patches are dominated by Urtica dioica and Phalaris arundinacea.At the site ten chamber flux measurement stations have been installed, of which two are located on ditches, the others on the land surface.Data have been collected from May 2003 until August 2008 with monthly to weekly intervals.The average annual air temperature is 9.8 • C and an average precipitation of 793 mm yr −1 .The site was extensively described by Hendriks et al. (2007).Ruwiel (52 • 10 30 N, 4 • 56 30 E) is a small nature reserve (Armenland Ruwiel) with a high water table.Climatic conditions are the same as those of Horstermeer.It is a species-rich, mesotrophic hay pasture, dominated by sedges (Carex sp.) and Eriophorum angustifolium, and has never been manured or fertilized.It is mown only once a year.The water table is kept artificially 0.3-0.5 m higher than that of the surrounding agricultural land.Within the reserve, the water level varies between 0 to 30 cm below the surface, outside the reserve it varies between 20 and 60 cm below the surface.The soil is a clayey fen peat.Four measurement stations have been installed in the winter of 2003-2004 in the reserve.CH 4 flux chamber measurements were taken with a bi-weekly to monthly interval from 22 January, 2004 till 20 December, 2005.
Kytalyk is a tundra wetland site, located in Northeastern Siberia, in the Indigirka lowlands near Chokurdakh (70 • 48 N, 147 • 26 E, elevation 48 m).The climate is arctic, with an annual average temperature measured at the Chokurdakh airport weather station of −14.3 • C, the warmest month being July, the coldest January.The research site consists of two different morphological units: a river floodplain, and the bottom of a former thaw lake, both underlain by continuous permafrost with a network of ice wedge polygons.The area is characterized by silty soils with a peaty topsoil.The CH 4 flux measurements have been made on both the thaw lake bottom and the river floodplain.The sites at the river floodplain are situated in Carex/Eriophorum or Arctica fulva vegetation and show very high fluxes.The vegetation on the thaw lake bottom is more varied, with hummocks and pools dominated by Sphagnum, Carex/Eriophorum meadows and vegetation dominated by Betula nana or Eriophorum hummocks on higher parts.Compared to the river floodplain, the fluxes are modest, being lowest in the Sphagnum vegetations, despite high water table.Air temperature, precipitation and snow data are based on local site measurements in summer, supplemented with data from the Chokurdakh airport weather station.CO 2 flux measurements using chambers and eddy covariance started in 2003, CH 4 flux measurements using chambers in the summer of 2004.From 2004 till 2006, CH 4 flux was measured only once a year in short (4-6 days) field campaigns, from 2007 onwards measurement campaigns included the months of July and part of August, with a higher measurement frequency.The measurement stations sample the entire range of wetland vegetation types in the area.The site and its CO 2 and CH 4 flux measurement methodology have been described extensively by Van Huissteden et al. (2005) and Van der Molen et al. (2007).A first attempt at modelling the CH 4 fluxes was undertaken by Petrescu et al. (2008).

Field methods and error sources
The flux measurements were carried out using closed chambers (non-transparant PVC, of different sizes; in Kytalyk a smaller sized chamber was used).The measurement procedure has been described in detail by Hendriks et al. (2007), andVan Huissteden et al. (2005).For each flux measurement, at least five gas concentration measurements were taken at regular time intervals per chamber per flux measurement.Before May 2004 (Horstermeer and Ruwiel sites) CH 4 concentrations were determined from syringe samples taken from the chambers and analysed on a gas chromatograph.Thereafter CH 4 analysis was performed in the field using an Innova 1312 photo-acoustic gas analyser, fitted with a CO 2 (sodalime) and H 2 O (silica gel) filter to prevent interference of high concentrations of these gases with the CH 4 analysis.A cross-check was made with flux measurements at the Ruwiel site to check the agreement between the two methods, no significant differences were detected.Thereafter, the Innova1312 has been frequently calibrated according to the recommendations of the manufacturer (Van Huissteden et al., 2005;Hendriks et al., 2007).
For all sites analysis of soil organic matter content and dry bulk density was available as input for soil profile informa-tion for the model.For the Horstermeer and Ruwiel sites also soil pH was available, and pF curve estimates.The pF curve estimates for Kytalyk have been based on average pF curve data from peat profiles (Petrescu et al., 2008).
For data-model comparison, also error sources in the data should be considered.In the case of the chamber flux measurements, these consist of: 1.The statistical error in the flux measurements which is inherent to the method of flux calculation.This consists of calculating the gradient of CH 4 concentration vs. time using regression.This gradient is subject to statistical error, which is specified as a standard deviation on the flux.
2. The flux calculation method.Here, we assumed that the time-CH 4 concentration relation is linear, which is a common approach and is valid when the measurement period is kept as short as possible.However, the relation may not be linear, for instance as a result of a decreasing soil-chamber concentration gradient during the measurement.In that case a linear approximation causes underestimation of the fluxes (Kutzbach et al., 2007).
3. Other technical errors of the flux measurements (concentration analysis errors, chamber leakage, other disturbances of the measurement) may result in faulty measurements.In particular on extremely wet sites with soft soils, excessive CH 4 flux by ebullition is easily induced by site access.This results in overestimation of fluxes.With the Innova 1312 such events are detected by high starting concentrations of the measurement, otherwise these errors can be detected by plotting the time-CH 4 concentration relation for every measurement and checking for irregularities.However, it cannot be excluded that faulty measurements remain unnoticed.
4. Spatial and temporal variability of the CH 4 fluxes.Although observed fluxes are generally related to water table, soil temperature and vegetation, the variability of fluxes within measurement points with similar soil type, vegetation and water table position is usually high (e.g., Van Huissteden et al., 2005;Hendriks et al., 2009a).This small-scale spatial variation is probably related to unquantified differences in vegetation characteristics and soil.Also small-scale (daily and shorter) temporal variation in CH 4 fluxes occurs.This variation has been observed with CH 4 flux measurements using eddy covariance, but is unnoticed with daily chamber flux measurements.This temporal variation may be caused by air pressure variations and variations in near-surface turbulence (Hendriks et al., 2008;Wille et al., 2008).Also the diurnal variation of CH 4 fluxes may be strong, as is the case at the Horstermeer site (Hendriks et

Procedure
The GLUE method is based on Monte Carlo simulations of the model with randomly chosen parameter values (Beven, 2009).Monte Carlo simulations are large sets of model simulations, for each single simulation within this set the value of one or more parameters of the model is chosen randomly within a pre-defined value range.The parameter values have been sampled from a uniform distribution, assuming no prior knowledge of the correct parameter value.The results of each model run are compared with the data from the study site being considered.The performance of the model run is summarized by an objective function value, derived from the differences between data and model.Different types of objective functions can be chosen, depending on the desired features of the data to which the model should fit best (see below).We used 5000 model runs for each site/data set separately.
After completion of all model runs, the distribution of the objective function values over the value range of the parameters is used to analyse the sensitivity of the model.In particular the difference of this distribution for "behavioural" (model runs that that fit well to the data) and "non-behavioural" (poorly fitting model runs) indicates to which parameters the model outcome is sensitive, and the range of parameter values that contribute to a good modeldata fit (Hornberger and Spear, 1981;Young, 1983).The selection of behavioural models is based on the objective function value, different criteria can be used.In this study we selected the 2% of all runs with the highest objective function value as "behavioural", which allows to study the parameter sensitivity for all sites irrespective of the maximum objective function value.To plot the results of the behavioural model runs the 1% best runs have been selected.
A large difference between the cumulative distribution of the behavioural runs and that of all model runs indicates a strong sensitivity for the parameter in question.The Kolmogorov D statistic for testing of differences in distribution functions is a measure of the parameter sensitivity.The D value should be seen as a qualitative measure of the difference between the distributions, since for large numbers of simulations the statistical test for D is not robust (Beven, 2001).

Objective functions
We tested three different objective functions for use in subsequent analysis.In the case of CH 4 fluxes, data-model comparison fluxes can be performed in different ways.Comparing the model results with a single measurement station appears obvious, but because of the high spatial variability of fluxes mentioned above, it can be argued that the model should reproduce the average flux of a group of measurement points with similar soil, hydrology and vegetation (grouped sites hereafter), rather than the measurements of single stations.We tested both approaches.Also it may be desirable to account for errors of the measurements in the data-model comparison.In case of grouped sites, the within-group variance can be taken as a statistical error on the flux.If the model results are compared with individual sites, the statistical error of the flux measurements is taken.This results in the following choice of objective functions: 1.The Nash-Sutcliffe efficiency (NS hereafter, Nash and Sutcliffe, 1970) is often used for model-data comparison (Beven, 2001).It is defined as where σ 2 e is the error variance, in which ŷt is the predicted value at time t, and y t the observed value, and σ 2 o the variance of the observations.E has the value of 1 for a perfect fit, and values close to, or below 0 when the error variance is of the same magnitude or larger than the variance of the observations.In that case the model performs not better, or worse than a flux estimate simply based on the average of the data.

Regression Comparison (RC). Since the methane flux
can also easily be modelled by regression of local flux data on water table and soil surface temperature (e.g., Van Huissteden et al., 2005), we also tested a variant of the NS objective function, that compares the model results with an estimate from a regression model.In this Regression Comparison (RC) function, σ 2 o is replaced by the variance of the residuals of a multilinear regression with water table position and soil temperature as independent variables.A value close to, or below 0 indicates that the model performs not better, or worse than the regression equation.
3. Summed Z score accounting for data error.The NS and RC objective functions do not account for the statistical error in the data outlined above.Therefore an objective function that accounts for this error has also been tested, based on the summed z scores of the deviation between modelled and measured flux: Here z t is the absolute standardized model-observation deviation at time t, and σ t the standard error of the flux measurement.To combine these in a single measure the z scores are summed and divided by the number of observations.To convert this to an objective function value which increases with better model-data fit and to scale between 0 and 1, the exponential of the result is computed, with S as a shape parameter.This scales the objective function value between 1 and 0. Depending on S, Z rises rapidly with low values of T t=1 z t , allowing good discrimination of best fitting model runs.This objective function can be easily adapted for datamodel comparison using grouped measurement points, by comparing with the average flux of the group.In this case of the Z score objective function, denoted below as Z group , the standard deviation σ t is the group standard deviation.This has the effect that the data error of individual measurements is not included in the objective function, but the within-group spatial variation instead.
The NS and RC objective functions also can be used for grouped measurement points as well as single points.Since NS and RC are based on comparison of variances, significance of the objective function values can be evaluated using an F (variance ratio) test, with the degrees of freedom determined by the number of observations on which the objective function value is based.

Model parameters and analysis procedure
The model parameters that potentially influence the CH 4 flux and hence are tested, can be grouped into 1) microbial reaction rate parameters, 2) vegetation parameters and 3) soil physical and chemical parameters (horizon thicknesses and properties).We used for all sites a soil profile definition consisting of two horizons, in the case of the Horstermeer and Ruwiel sites this was a aggregation of a detailed soil profile with more horizons (Table 1).The ranges of the parameters are based on Walter (2000) and Van Huissteden et al. (2006), or for the soil parameters, on parameter ranges measured at the sites.Parameter ranges are derived from literature and discussed by Walter (2000).For a first approximation, parameters in group 1 and 2 and group 3 have been tested separately to select the parameters that have a significant influence on the model output.Next, a combination of sensitive parameters of all groups (soil and non-soil parameters) has been tested.
The following test procedure has been applied: 1.For the Ruwiel site, analysis of CH 4 production/oxidation and vegetation parameters, to compare the objective functions described above, and to compare the effects of using data from individual measurement stations versus grouped stations.Based on these tests, an objective function is selected for subsequent analysis.The Ruwiel site has been chosen because previous modelling experiments showed a good model-data fit for this site (Van Huissteden et al., 2006).
2. All sites: CH 4 production/oxidation and vegetation parameters, to test parameter sensitivity and its consistency among different wetland sites, and effects of using data time series of different length.
3. All sites, combining soil, CH 4 production/oxidation and vegetation parameters, to study the effects of soil parameters.

Objective function selection
The tests for the Ruwiel site shows the effects of selection of the objective function.From the site, three measurement points with similar vegetation (species-rich grasses and sedges) and water table (frequently at or above soil surface) have been selected for data-model comparison, in grouped and single station mode.The number of observations for each measurement station is 26.Monte Carlo simulations have been made for all CH 4 production/oxidation and vegetation parameters in Table 1. Figure 1 shows the objective function values of NS plotted against parameter value for each parameter, Fig. 2 the plots of the parameter distributions of the behavioural model runs.Figure 1 indicates that there is strong equifinality.Behavioural model runs that exceed the F test p = 0.1 probability limit (NS > 0.3937) are realized with quite different sets of parameters; for all parameters the entire parameter range is covered.In Fig. 2, the deviations of the parameter distributions of the behavioural model runs from the original parameter distribution, measured with the D statistic, indicates the parameter sensitivity.Here we take the behavioural model runs as the best 2% runs (NS > 0.6408) in stead of the F test criterion, for comparison with the other sites.The sensitivity is highest for maximum root depth Z root , with a value for the D statistic of 0.30.Other sensitive parameters are the plant transport factor V transp , Q10 for CH 4 production and the maximum primary production of the vegetation with values for D of respectively 0.29, 0.26 and 0.24.Model runs with relatively high Q10, shallow root depth, low primary production and high plant transport factor tend to produce a better fit to the data.; blue: non-behavioural runs below the probability limit.For parameter explanation, see Table 1.

Fig. 1.
Objective function values (Nash-Sutcliffe) for 5000 runs of PEATLAND with randomly chozen parameters using a uniform distribution over the parameter range.The model results have been compared with data from the Ruwiel measurement site 1, 2 and 3. Red: model runs of which the objective function value exceeds the 0.1 probability limit (F-test, see text); blue: non-behavioural runs below the probability limit.For parameter explanation, see Table 1.  ; blue: non-behavioural runs below the probability limit.For parameter explanation, see Table 1.  1.
Similar plots for the other objective functions are not shown, but Fig. 3 and Table 2 summarize the results for all objective functions.The maximum value for the N S efficiency in Table 2, 0.75, is quite high and indicates that the model explains the data significantly better than an estimate based on the mean of the data.The F test shows that the variance of the model-data residuals is significantly smaller (p < 0.1) than the variance of the data.Also the maximum value for RC, 0.30, is positive, indicating that the model performs better than a regression on water table and soil temperature.However, the value does not indicate a significant difference between the variance of the regression residuals and the variance of the model-data differences (F test, p = 0.18 > 0.1).Although revealing on the performance of the model with respect to a simple regression model, it will not be discussed further here since it behaves similar as the NS function.
The Z objective function is a stricter requirement for model-data fit, since it does not test on variances over the entire data range but requires a good fit for each individual data point, weighed against the known data error.If the model is compared with the average fluxes of the three sites, a maximum value for Z group of 0.28 results (shape parameter S = 1).However, when compared with individual sites, the results become much worse, resulting in near-zero Z values.The model clearly cannot follow the individual data points of an individual site.The large difference between Z and Z group also results from the fact that the within-group variance of a site group is much larger than the statistical measurement error of the individual measurements.2) for all objective functions for the Ruwiel site.For parameter explanation, see Table 1.Above: D statistic for grouped Ruwiel measurement sites 1, 2 and 3; below; the same, evaluated for the measurement points individually.NS: Nash-Sutcliffe efficiency; R: Regression comparison; Z: Z statistic (see text); Z grouped: Z group statistic for grouped sites (see text).
The tests for the three objective functions generally indicate high sensitivity (highest values for the D statistic) for the Q10 and T ref parameters for CH 4 production and the parameters related to vegetation biomass and plant transport of CH 4 , P max , P ox and V transp (Fig. 3).However, there are also conspicuous differences for the objective functions.The Z group function indicates a much higher sensitivity for P ox , Q10 and T ref than NS and RC and a lower sensitivity for V transp .The parameters P ox and V transp may affect plant CH 4 transport rate antagonistically, in the absence of interaction with other parameters, a high V transp may be compensated by a high P ox .Apparently, different objective functions select one or the other parameter as the most sensitive one.Q10 and T ref may affect differences in CH 4 flux due to soil temperature and may result in a better fit of variations in CH 4 flux to seasonal and shorter term differences in temperature; the Z group function weighs these differences more strongly than NS and RC.Also the tests for individual measurement points show differences in the parameter sensitivities.Here, also the f anaer parameter proves to be sensitive.Apparently, part of the spatial variation between individual measurement points is explained by the occurrence of CH 4 production in anaerobic microsites in the soil at lower water tables in the quite dense, clayey peat of this site.
For individual measurement points, the model cannot capture flux differences that are related to small-scale spatial and temporal variation (Table 2 and Fig. 3).This is most clearly shown by the low values of the Z function for individual sites, but also the maximum NS value is lower for the individual sites than for the grouped sites.Figure 4 shows the 1% best model runs for the NS and Z group functions, compared with the data.For grouped sites, the NS function results in a very slight positive bias with respect to the data, the Z group function shows a slightly lower bias.However, multiday temporal fluctuations in fluxes are captured much better when the NS function is used.For both objective functions, some model runs show high flux peaks.For the highest measured fluxes in the second summer season, these peaks may be realistic, but cannot be checked against the data because the data density is too low to reject unrealistic peaks that fall between to measurement dates.
A possible cause of the observed high flux peaks may be ebullition events that are induced by air pressure variations (Whalen, 2005), in particular for the high water table sites.The representation of ebullition in the model is very simple and depends only on the pore water CH 4 concentration (Walter, 2000).However, for the Horstermeer site a statistically significant relation between air pressure and CH 4 flux measured by eddy covariance is absent (Hendriks et al., 2009b).
Concluding, the NS objective function performs best as it results in model runs that follow better the yearly and withinyear variations in the fluxes, and performs well also for individual measurement points.The test on the other research sites have been restricted to the NS function.

Comparison between study sites
A comparison of the sensitivity for the different sites (Ruwiel, Horstermeer, Kytalyk) shows the parameter sensitivity for sites that differ in geography and wetland type.Only grouped measurement points have been considered, using the NS objective function.For Ruwiel, the same sites have been selected as above (number of observations n = 26).The Kytalyk site has been split into two contrasting measurement point groups: the river floodplain with sedge and grass vegetation (n = 30, Kytalyk Floodplain hereafter), and the oligotrophic terrace with submerged Sphagnum vegetation (n = 27), mainly located in ice wedge polygon centres (Kytalyk Terrace hereafter).These are highly contrasting sites, both with high water table but different vegetation and strongly different CH 4 fluxes.Likewise, two site groups have been tested at Horstermeer: a site group consisting of sites with varying water table (high in winter, up to 35 cm below surface in dry spells in summer) and vegetation that is well capable of CH 4 transport (Holcus lanatus grass, Equisetum palustre, Glyceria maxima, Juncus effusus, point nrs.3-5, n = 64), and two extremely marshy sites along ditches dominated by Typha angustifolia and Glyceria maxima, where the water table is at a constant level throughout the year (nrs.7-8, n = 24, Horstermeer Wet hereafter).For the varying water table sites at Horstermeer, a longer time series is available, the tests have been run for both a shorter time series (n = 24) for compatibility with the other sites, and the longer time series to study the behaviour of the model over longer runs (Horstermeer 1 and Horstermeer 2 hereafter for respectively the short and long time series).
The results (Fig. 5, Table 3) show clear differences in the ability to model the CH 4 fluxes for the sites.The model performs best for Ruwiel; for Horstermeer 1 and for Kytalyk Terrace the model also performs significantly better than an estimate based on average measured fluxes.For Horstermeer 2, and for Kytalyk Floodplain, also positive objective function values are produced, but these do not exceed the significance limit.For Horstermeer Wet only negative NS values were calculated, meaning that he model does worse than an average of the data.
The relatively poor fit of the model to the longer time series at Horstermeer is caused mainly by high flux peaks observed in the third and fourth year.The model simulates flux peaks, but not exactly at the same dates as the observations (Fig. 6).The same holds for the Horstermeer Wet and Kytalyk Floodplain (not shown).Remarkably, the model performs less well for the eutrophic high water table sites.However, more tests on more sites would be necessary to confirm whether this is a consistent feature of the model.For Horstermeer Wet, the model completely fails to simulate the high measured fluxes.Measurement error cannot be excluded here, since these sites are extremely sensitive to disturbance during flux chamber measurements, despite precautionary measures such as boardwalk construction and careful  5.For parameter explanation, see Table 1.  3. For parameter explanation, see Table 1.
analysis of the measurement data.However, if suspected data points are deleted, the result only slightly improves to a still negative maximum NS value of −0.10.
Differences in the parameter sensitivity (D statistic) between the site tests are found for most parameters (Fig. 5), although in general the same parameters that were insensitive for Ruwiel are also insensitive for the other sites.Exceptions are R 0 and P ox which show a high sensitivity for all sites except Ruwiel.The sensitivity of Q10, Z root , T ref , P max , f anaer , R 0,peat , V transp also vary among the sites.Within the Horstermeer site, the sensitivity of T ref , S, Q10 ox and Z root varies depending on the length of the data time series against which the model is tested.With the more difficult model fit for the longer times series these parameters also contribute to a better model fit, while they contribute insignificantly for the shorter time series.With exception of Z root , all these parameters influence the temporal variation of CH 4 emission throughout the year.The sites also differ markedly in sensitivity to T ref , the reference temperature for the Q10 relation of CH 4 production.Highest sensitivity is found for the Kytalyk Terrace and Horstermeer 2 data; for the first, low values of T ref give the highest objective function values, for the second high T ref values.Also for the Kytalyk Floodplain sites low T ref values result in a better model fit, but the effect is strongest for the Terrace sites where the active layer is thinner and soil temperatures generally lower.This agrees well with the expected differences in microbial communities between arctic and temperate wetland sites.Microbial populations in arctic soils tend to show metabolic activity also at low temperatures (Rivkina et al., 2007).
Also within Kytalyk differences in parameter sensitivity arises between the Floodplain and Terrace sites, in particular for the plant oxidation P ox transport V transp and P max primary production parameters.The floodplain and terrace points markedly differ in biomass and probably also net primary production of the vegetation, with highest biomass occurring on the floodplain.    5.For parameter explanation, see Table 1.parameter value for the P ox and V transp parameters are shown in Fig. 7.For the terrace points, the model performs best with a higher plant oxidation rate and a lower transport rate factor, whilst for the floodplain this pattern is reversed.The sensitivity pattern confirms the inferred processes that are responsible for the spatial difference in fluxes between river floodplain and terrace.On the terrace, the transport rate of CH 4 through the dominating Sphagnum vegetation is low, while oxidation is high due to the presence of symbiotic methanotrophic bacteria in the plants (Raghoebarsing et al., 2005; confirmed by measurements on Sphagnum samples, N. Kip, personal communication, 2008).On the floodplain, transport rate through Carex and Eriphorum species is high (Van Huissteden et al., 2005), with low oxidation rate.The results do not confirm that the net primary production on the floodplain is high on the floodplain.For the Kytalyk Floodplain sites, low P max values result in a higher objective function value.

Soil parameters
Next to the parameters above, also soil parameters may influence model results strongly and could explain the spatial variability of CH 4 fluxes.The soil properties tested for each horizon are the water-filled porosity θ sat , organic matter percentage O, thickness of the upper soil horizon H 1 , soil pH, and for Kytalyk, the shape parameter of the relation between temperature and frozen water content k freeze has been added.These parameters are combined in the test with the most sensitive vegetation and microbial population parameters: R 0 , Q10, P ox , Z roots , P max , f anaer , k delay , R 0,peat and V transp .All sites have been tested, for the Ruwiel site also the individual measurement stations have been tested (Fig. 8).
With exception of the thickness of the upper soil horizon, all soil parameters appear to be less sensitive than the vegetation parameters.Only soil pH has a somewhat higher D statistic value, in particular for Kytalyk Floodplain.This site shows the strongest sensitivities to soil parameters.A remarkable feature is the different sensitivities of pH for the Horstermeer long and short time series; apparently more model parameters need to be adjusted for making the model fit to the longer time series.The sensitivity of the k freeze parameter is inconsistent, it is more sensitive for Kytalyk Floodplain than for Kytalyk Terrace.This may be a spurious effect introduced by the generally low model fit for Kytalyk Floodplain.Concluding, the model is not very sensitive to uncertainty in soil characteristics.The thickness of the upper soil horizon is the most critical soil parameter.
For the Ruwiel individual measurement stations, the soil parameters are hardly sensitive.Apparently variability in soil properties does not contribute here to the observed smallscale spatial variability in CH 4 fluxes of individual measurement points within the site.Comparing Figs. 5 and 8, there are no large differences in sensitivity of the vegetation and microbial parameters if tests are done with and without soil parameters.The same parameters that showed high D statistic values without soil parameters also show high values with soil parameters included.The only exception is the CH 4 production factor from peat, R 0,peat , which becomes less important for the Horstermeer site when soil parameters are added.For other vegetation and microbial parameters D tends to be higher when soil parameters are included, in particular for the Horstermeer 2 and Kytalyk Floodplain data.We infer that model fit problems arising from soil parameters can be compensated by adjustments of the vegetation and microbial parameters, in particular for sites where the model fit in general is rather poor.
Adding the soil parameters does not improve the model fit (Table 4).In particular for Ruwiel and Kytalyk Floodplain, the maximum objective function value with the soil parameters included is lower than without the soil parameters.This might be caused by the deletion of some of the other param-    4, including soil parameters.For parameter explanation, see Table 1 and text

Recently added parameters
The parameters k delay and f anaer were newly added to the model, respectively to simulate time delay in restoration of anaerobic conditions in the soil at rapid water table rise, and CH 4 formation above the water table.f anaer is sensitive only for the Horstermeer 2 dataset, k delay for both Ruwiel and Horstermeer 2. For the other sites, f anaer influences the model fit more strongly when the soil parameters are added.
We conclude that both parameters may be useful depending on site conditions, in case of a poor model-data fit these parameters may improve the model fit to some extent.However, these parameters never appear to have a strong overall influence on model fit.
A parameter of the SOM production submodel, that is not included in the original model by Walter (2000) is the correction factor on stronger exudate production in spring S (Van Huissteden et al., 2006).In all tests, this factor attains only low values for the D statistic, so it does not influence the modelled CH 4 emission significantly.Soil pH also has been added to the original Walter (2000) model as a factor influencing CH 4 production in the model (Van Huissteden et al., 2006).It proves an effective parameter in some cases (Kytalyk, Horstermeer).

Conclusions
With regard to the overall model performance, we conclude that the PEATLAND-VU model is capable of simulating CH 4 fluxes in temperate and arctic wetlands, under different type of site conditions.However, not in all cases the model improves prediction of emissions, compared to a sim-ple emission factor approach based on averages of measurement data.In three of the six data sets the model results were significantly better than an estimate of the fluxes based on averaged data.In two data sets, the model still performed better but the difference was not large enough to classify it as significant.For one data set, the model did not perform well, but in that case data error cannot be excluded.For one data set (Ruwiel) the model also has been compared with a multilinear regression model derived from regression of the flux measurement data on soil temperature and water table.Although the objective function values indicate better performance for the PEATLAND-VU model, it does not significantly outcompete the regression model.However, a regression model is less relevant for upscaling purposes since larger scale spatial upscaling based on regression results of individual sites depends strongly on local data availability compared with a process model.
Details of spatial and temporal variation are poorly reproduced by the model.Analysis of the model-data deviations shows that the model is not capable of simulating shortterm temporal variation that may occur on a daily time scale.However, the model simulates longer term temporal variation (seasonal and weekly-monthly) correctly.Since longer term variations and the average yearly cycle is more important for temporal upscaling than timing of the peaks, this does not have to be a problem.We experimented with different types of objective functions in applying the GLUE method.One type of objective function accounts for short-distance spatial variability of the fluxes by comparing the model results with averages of groups of points with homogeneous vegetation/soil characteristics, another objective function for the flux measurement standard error.The first approach provided the best results since it averages out some of the small scale spatial variability inherent to CH 4 fluxes.
The parameters to which the model is most sensitive are vegetation parameters and temperature sensitivity of the methanogenic microbial population.In particular parameters related to transport of CH 4 through plants (transport rate, oxidation during transport and root depth) determine the model sensitivity.This is not surprising since in wetland sites plant transport is usually the dominant soil-atmosphere transport mechanism for CH 4 .By contrast, the model is not very sensitive to soil parameters, which are an important source of spatial variability in input data.For all sites, adding soil parameters to the GLUE analysis resulted at best in very small improvements of the model results.Uncertainty with respect to water table and soil temperature input has not been tested in this study.
The parameter sensitivity and the parameter values resulting from the GLUE optimalisation agree well with a priori knowledge on the parameters.For the arctic site, a lower reference temperature for CH 4 production temperature sensitivity resulted, compared with the temperate climate sites.This agrees well with observed temperature sensitivity of microbial populations in arctic soils.Also, at the site where oxidation of CH 4 by symbiotic methanotrophs was observed in Sphagnum vegetation, the optimalisation correctly resulted in higher values for the plant oxidation parameter and lower values for plant transport rate.We conclude that GLUE analysis may enhance insight in the local relevance of processes included in the model.
There is considerable parameter interaction within the model.In particular parameters for vegetation and microbial population strongly interact.The reason is that some of the parameters act antagonistically in the model.For instance a higher plant transport oxidation rate suppresses the effects of high plant transport rate.Tuning both parameters in an opposite way results in a good model fit for a large range of parameter values.For the purpose of model tuning it may be useful to summarize these two parameters in one, in particular since these have been defined in the original model as bulk parameters without physical background.
The GLUE analysis suggests future improvements for wetland CH 4 emission modelling.The vegetation parameters contribute strongly to model uncertainty.These are spatially highly variable, and several of the relevant vegetation parameters above are difficult to quantify.Therefore it is highly important to invest in improvement of vegetation data, in particular data on methane emission functionality of wetland species and wetland vegetation units, and wetland vegetation mapping.Fortunately the results for the Kytalyk site suggest that some of these parameters for vegetation can be derived from general vegetation characteristics (e.g.dominance of arenchymous tissue in wetland plants, oxidation of CH 4 in Sphagnum) and can be constrained by model fitting.

Fig. 1 .
Fig. 1.Objective function values (Nash-Sutcliffe) for 5000 runs of PEATLAND with randomly chozen parameters using a uniform distribution over the parameter range.The model results have been compared with data from the Ruwiel measurement site 1, 2 and 3. Red: model runs of which the objective function value exceeds the 0.1 probability limit (F-test, see text); blue: non-behavioural runs below the probability limit.For parameter

Fig. 1 .
Fig. 1.Objective function values (Nash-Sutcliffe) for 5000 runs of PEATLAND with randomly chozen parameters using a uniform distribution over the parameter range.The model results have been compared with data from the Ruwiel measurement site 1, 2 and 3. Red: model runs of which the objective function value exceeds the 0.1 probability limit (F-test, see text); blue: non-behavioural runs below the probability limit.For parameter

Fig. 2 .Fig. 2 .
Fig. 2. Cumulative distributions of the parameters in the model runs of Figure 1.Green: distribution of behavioural runs; red: distribution of all runs.Behavioural runs are the best 2% of the Monte Carlo simulations.D is the Kolmogorov-Smirnoff D statistic for comparison of distribution functions.For parameter explanation, see Table1.

Fig. 3 .Fig. 3 .
Fig. 3. Bar graphs showing the D statistic for the tested model parameters (see also Figure 2) for all objective functions for the Ruwiel site.For parameter explanation, see Table 1.Above: D statistic for grouped Ruwiel measurement sites 1, 2 and 3; below; the same, evaluated for the measurement points individually.NS: Nash-Sutcliffe efficiency; R: Regression comparison; Z: Z statistic (see text); Z grouped: Zgroup statistic for grouped sites (see text).

Fig. 5 .
Fig. 5. Bar graphs showing the D statistic (see also Figure 2) for all sites indicated in Table5.For parameter

Fig. 5 .
Fig. 5. Bar graphs showing the D statistic (see also Figure 2) for all sites indicated in Table5.For parameter

Fig. 7 .Fig. 7 .
Fig. 7. Objective function value plot for the plant transport parameters Pox (oxidation of CH4 during transport) and Vtransp (Plan transport factor) for Kytalyk floodplain and terrace.

Fig. 8 .
Fig. 8. Above: Bar graphs showing the D statistic (see also Figure 2) for all sites indicated inTable 5, including

Fig. 8 .
Fig. 8. Above: bar graphs showing the D statistic (see also Fig. 2) for all sites indicated inTable 4, including soil parameters.For parameter explanation, see Table 1 and text.Below: the same for Ruwiel site, individual measurement stations.
Fig. 8. Above: bar graphs showing the D statistic (see also Fig. 2) for all sites indicated inTable 4, including soil parameters.For parameter explanation, see Table 1 and text.Below: the same for Ruwiel site, individual measurement stations.
Walter (2000)the day fluxes at Horstermeer are higher than during the night, the daytime flux chamber measurements result in an overestimation of the fluxes.These sub-daily variations and processes are not included in the PEATLAND-VU/Walter (2000)model.

Table 1 .
Model parameters of PEATLAND-VU which have been included in the GLUE analysis of the model, with value ranges of the parameters.

Table 2 .
Summary of objective function values for Monte Carlo runs for the Ruwiel site.The probabilities for the NS and RC objective functions are based on an F variance ratio test (see text).The cutoff value for behavioural runs is the lowest value of the best 100 (2%) runs, and has been used for calculating the D statistic in Fig.2.
The plots of objective function value vs.

Table 3 .
Summary of objective function values for Monte Carlo runs for all sites.The probabilities for the NS objective function is based on an F variance ratio test (see text).The cutoff value for behavioural runs is the lowest value of the best 100 runs, and has been used for calculating the D statistic in Fig.2.

Table 5 ,
including soil parameters.For parameter explanation, see Table1 and text.Below: the same for Ruwiel site, individual

Table 4 .
Summary of objective function values for Monte Carlo runs for all sites.The probabilities for the NS objective function is based on an F variance ratio test (see text).The cutoff value for behavioural runs is the lowest value of the best 100 runs, and has been used for calculating the D statistic in Fig.2."X": no behavioural model runs.