Parameter interactions and sensitivity analysis for modelling carbon heat and water fluxes in a natural peatland, using CoupModel v5

In contrast to previous peatland carbon dioxide (CO2)model sensitivity analyses, which usually focussed on only one or a few processes, this study investigates interactions between various biotic and abiotic processes and their parameters by comparing CoupModel v5 results with multiple observation variables. Many interactions were found not only within but also between various process categories simulating plant growth, decomposition, radiation interception, soil temperature, aerodynamic resistance, transpiration, soil hydrology and snow. Each measurement variable was sensitive to up to 10 (out of 54) parameters, from up to 7 different process categories. The constrained parameter ranges varied, depending on the variable and performance index chosen as criteria, and on other calibrated parameters (equifinalities). Therefore, transferring parameter ranges between models needs to be done with caution, especially if such ranges were achieved by only considering a few processes. The identified interactions and constrained parameters will be of great interest to use for comparisons with model results and data from similar ecosystems. All of the available measurement variables (net ecosystem exchange, leaf area index, sensible and latent heat fluxes, net radiation, soil temperatures, water table depth and snow depth) improved the model constraint. If hydraulic properties or water content were measured, further parameters could be constrained, resolving several equifinalities and reducing model uncertainty. The presented results highlight the importance of considering biotic and abiotic processes together and can help modellers and experimentalists to design and calibrate models as well as to direct experimental set-ups in peatland ecosystems towards modelling needs.


Introduction
Understanding and quantification of interactions between different processes and between different parameters is required for reducing uncertainty in prognostic modelling in carbon (C) cycle research.Undisturbed peatlands act as carbon sinks and have accumulated at least 550 Gt of C, which is equivalent to twice the C stock in the forest biomass of the world (Gorham, 1991;Parish, 2008).A more recent estimate for exclusively northern peatlands amounts to 436 Gt of C (Loisel et al., 2014).Management or climate change can cause this carbon to be released as CO 2 emissions as has been shown from measurements (Maljanen et al., 2010;Drösler et al., 2013;Petrescu et al., 2015).Process oriented models are necessary to transfer the knowledge gained from measurements to different locations, management or future climate scenarios.Furthermore, such models can help to understand the processes underlying the observations.Although only a few of the parameters used in process models are known as site-independent, unambiguous constants from laboratory experiments.All others need to be either assumed, or gained from calibration procedures (e.g.Kennedy and O'Hagan, 2001;Wang and Chen, 2012), but not all parameters have a strong impact on model output and performance (i.e.fit between modelled and measured variables, whereas in this manuscript, variable always refers to a time series that is ei-Published by Copernicus Publications on behalf of the European Geosciences Union. ther the output of the model or the measurement to which the model output is compared).Monte Carlo-based sensitivity analyses are used to identify key parameters for both the performance and the impact on various major model outputs (e.g.Verbeeck et al., 2006;Van Oijen et al., 2011;Santaren et al., 2014).
Many studies investigated single processes and their parameters, whereas only a few consider different biotic and abiotic processes and multiple calibration variables; several modelling studies have explored peatland hydrology (e.g.Dimitrov et al., 2010;Dettmann et al., 2014) and heat fluxes in peatlands (e.g.Granberg et al., 1999;Keller et al., 2004), whereas others concentrate on carbon fluxes and pools (e.g.Frolking et al., 2002;Verbeeck et al., 2006;Wu et al., 2013) where the focus is sometimes only on heterotrophic respiration (e.g.Abdalla et al., 2014).However, many processes are involved in the C cycle of peatlands; net ecosystem exchange (NEE) is the balance of photosynthesis and autotrophic respiration from plants as well as heterotrophic respiration from microbes.All NEE component fluxes are strongly inter-connected in several ways with the amount of plant biomass, temperature, radiation, nutrients and moisture availability (e.g.Clymo, 1984;Lindroth et al., 2007).Photosynthesis, soil temperature (Ts) and moisture depend, among others, on incoming radiation, transpiration and plant coverage.Heterotrophic respiration further depends on quality and quantity of plant litter (e.g.Yeloff and Mauquoy, 2006).In addition, phenological events such as the timing of snowmelt are important for soil temperature dynamics, biologic activity and peatland CO 2 fluxes (Aurela, 2004;Peichl et al., 2015).Different biotic and abiotic processes are realized in some modelling studies on peatlands, though, only the sensitivity to carbon fluxes or pools was tested (e.g.Yurova et al., 2007;St-Hilaire et al., 2010;Quillet et al., 2013;Webster et al., 2013;Wu and Blodau, 2013;Kim et al., 2014).Also, models are continuously extended or coupled with other models (e.g.Wang et al., 2005;Prentice et al., 2007;Giltrap et al., 2010;Hidy et al., 2012;Jansson, 2012;Tang et al., 2015), developing into more and more holistic models, accounting for plant and soil carbon processes, water and energy flows and biochemistry.However, often only parameters of the new module are tested (e.g.Belassen et al., 2010;Wania et al., 2010;Zhu et al., 2014;Tang et al., 2015), ignoring possible interactions between processes.
Another limitation of previous peatland modelling studies is the use of local sensitivity analyses, changing only one parameter or one input driver at a time (e.g.Hilbert et al., 2000;Yu et al., 2001;Zhang et al., 2002;Wania et al., 2009;Frolking et al., 2010;Tang et al., 2010;St-Hilaire et al., 2010).This approach does not account for possible interactions and non-linearity in equations (e.g.Saltelli et al., 2008;Quillet et al., 2013), but peatland processes are often non-linear and interact in many ways (Belyea, 2009).Therefore, we performed a global sensitivity analysis, calibrating parameters simultaneously and accounting for interactions.This allows for inter-correlation between the different parameters, which complicates the parameter constraint to an unambiguous solution; several combinations of different parameter values can lead to a similar good fit of model output to measured variables, which is defined as equifinality (Beven and Freer, 2001).The model sensitivity to such parameters might be hidden if equifinalities are not considered.Constraining a model based on multiple observation variables can help to resolve or reduce equifinalities (Carvalhais et al., 2010).The profit of using multiple constraints for model calibration and the importance of interactions between parameters and across different processes has been shown by sensitivity analyses on, e.g., forest ecosystems (Carvalhais et al., 2010;Santaren et al., 2014;Tian et al., 2014).Unlike previous modelling studies on peatlands, we therefore investigate the sensitivity to parameters from several different modules simultaneously, in their effect on not only on NEE but also on latent heat flux (LE), sensible heat (H ), net radiation (Rn), leaf area index (LAI), Ts, water table (WT) and snow, and identify parameter interactions.
However, criteria based on multiple variables imply a subjective weighting of variables and performance indices.Fitting the model to a certain variable might improve or worsen the performance in another variable (Carvalhais et al., 2010) and might therefore have implications for the parameter range judged as valid (e.g.Schulz and Beven, 2003).In this study, the effects of selecting a certain criteria on the resulting parameter range will be investigated.We avoided using a Bayesian approach, which was tested by Van Oijen et al. (2011), with several models including the CoupModel, using a data set of more than one variable.The single probability of the model as the summation of many different variables requires a detailed understanding of an error model that is typically not available in field measurements that represent many different errors for each set of variables.
The detailed ecosystem model CoupModel (Jansson and Karlberg, 2010) was used in this study for the following reasons; it is a well-established and widely used model (Jansson, 2012).Its structure is flexible and allows for simulation of different abiotic and biotic processes based on wellestablished physical equations, which can be selected by the user.The CoupModel includes all main components expected to have an impact on the carbon cycle: (i) a detailed module for simulation of heat and water fluxes in the soil and at the interface to the atmosphere; (ii) plant growth from photosynthesis, limited by water availability and temperature; (iii) plant respiration and litter fall; and (iv) a module for soil organic carbon (SOC) decomposition.A user-defined time step allows one to use the full information contained in measurements with high temporal resolution (i.e.hourly) on site scale.

Objectives
The aim was to identify and explore the connections within and between biotic and abiotic processes and parameters that are relevant for modelling NEE in a natural open peatland.Therefore, 54 parameters of the CoupModel v5 from different plant, decomposition, energy and water flux processes were calibrated to several different output variables, and several different sets of criteria for selecting acceptable runs were tested.The specific objectives were 1. to identify which processes impact which measured variable, by testing the sensitivity of model performance to the parameters of the different processes; 2. to evaluate the dependence of model performance and resulting parameter ranges on the performance index, the measured variable and the time period of the variable that are chosen as criteria; 3. to identify and describe equifinalities between parameters from different processes simulating carbon, energy and water fluxes; 4. to test the potential of all available observation data for model constrain and identify missing measurement variables by identifying sensitive or interacting parameters that cannot be constrained by the available data.
The answers to these questions will be crucial for model development and future calibrations of carbon models on peatlands: they will represent the most valuable information for selecting processes that need to be taken into account, for selecting parameters and their value ranges and considering parameter connections, as well as selecting sites and observed variables.They further help experimentalists to decide on the measurement of variables to make their site suitable for modelling.
2 Materials and methods

Site description
Degerö Stormyr (64.182016 • N, 19.55663 • E) is an oligotrophic, minerogenic mire located on a highland (270 m a.s.l) in the county of Västerbotten, Sweden.A detailed description of the site and the measurements can be found in Peichl et al. (2014) and references therein."The mire catchment is predominantly drained by the small creek Vargstugbäcken in the north-west.The depth of the peat is generally between 3 and 4 m, but depths up to 8 m have been measured.[...] The micro-topography is dominated by mainly carpets and lawns, with only sparse occurrences of hummocks" (Peichl et al., 2014).The plant community of the mire is dominated by cotton grass (Eriophorum vaginatum L.), tufted bulrush (Trichophorum cespitosum L. Hartm.) and Sphagnum mosses (Nilsson et al., 2008;Laine et al., 2012).Total above-ground biomass (moss capitula and vascular plants) is 141 ± 45 g m −2 (Laine et al., 2012).Seasonal maximum leaf area index of vascular plants was estimated at 0.8 m 2 m −2 in 2012 (Peichl et al., 2015).
The 30-year (1961-1990) mean annual precipitation and air temperature are 523 mm and +1.2 • C, respectively, while the mean air temperatures in July and January are +14.7 and −12.4 • C, respectively (Alexandersson et al., 1991).The snow cover normally reaches a depth of up to 0.6 m and lasts for approximately 6 months (Peichl et al., 2014).The peatland was continuously a sink for atmospheric CO 2 during 12 years of eddy covariance (EC) measurements, with a 12-year average (± standard deviation) NEE of −58 ± 21 g C m −2 yr −1 (Peichl et al., 2014).

Data used in this study
Hourly values of global radiation, air temperature, relative humidity, precipitation and wind speed were used as meteorological input data (Table 1).They were measured at the same tower on which the EC sensors were mounted.For gap filling (due to instrument failure) of the input data, as well as for the pre-evaluation period 1991-2000, daily data from the nearby (13 km away) standard climate station at the Svartberget field station were obtained.In the case of air temperature and relative humidity, seasonal regression relationships were applied to account for temperature and humidity differences between the site and the standard climate station.
An overview of the data used for calibration can be found in Table 2, a more detailed description is provided by Peichl et al. (2014) and references therein.Measured carbon concentrations per soil layer were used for estimation of pool sizes as described in Sect.2.3.5.The model was calibrated based on measured NEE, LE, H , WT, Rn, soil temperatures in −2 cm (Ts 1 ) and −42 cm (Ts 2 ) depth, snow depth and LAI of vascular plants (Table 2).NEE, LE and H were measured using the eddy covariance technique, and details for data processing were previously described in Peichl et al. (2014).In this study, only the measured values of NEE, LE and H were used for calibration (i.e.gap-filled values were omitted).Negative NEE values indicate net CO 2 uptake by the ecosystem from the atmosphere while positive NEE values indicate emission from the ecosystem to the atmosphere.All calibration data were averaged to hourly values, except snow depth and LAI values, which had a daily and biweekly to monthly resolution.

Model description and application to Degerö Stormyr
The CoupModel v5 from 12 December 2014 was used for simulations.The current version can be downloaded from the CoupModel home page (CoupModel, 2015).A detailed description can be found in Jansson and Karlberg (2010)  Every 4 cm between 0 and −32 cm, and every 12 cm between −60 and −338 cm 0 to −338 cm a Measurement resolution was the same or higher, except where mentioned differently.b The method description of meteorological input data applies to the climate station at the site.For gap-filling and for the pre-evaluation period, the data were obtained from the nearby standard climate station (Svartberget field station).
CoupModel allows the user to select between different submodels, different equations and different complexities of the used equations.The following sections describe the configuration as applied in this study.The model represents the ecosystem with a description of C and N fluxes in the soil and in the plants.It includes the main abiotic fluxes, such as soil heat and water fluxes that represent the major drivers for regulation of the biological components of the ecosystem.For application to Degerö Stormyr, the vegetation canopy was defined as two layers: vascular plants and mosses.The soil profile was divided into 16 layers with an increasing layer depth from 4 cm for the upper 9 layers to 60 cm in the lowest layer, resulting in a total depth of 3.4 m.The model internal time step was half-hourly for abiotic processes and hourly for nitrogen and carbon-related processes.The simulations were started 10 years prior to the evaluation period, so the system could adapt to the site conditions and become more independent of initial values.
The most important equations and the corresponding calibrated parameters can be found in Tables S1 and S2 in the Supplement.The major model assumptions related to the model application of the peatland are described below.Detailed assumptions with respect to fixed parameter values can be found in Table S3 in the Supplement.

Radiation interception, evapotranspiration and snow
An interception model for radiation and precipitation, a snow model and a surface pool model was used to provide boundary conditions at the soil surface.Cloud fraction was calculated from global radiation input and latitude.Incoming radiation was partitioned between one part, which was absorbed by the plant canopy and another part, which reached the soil according Beer's law (see Impens and Lemeur, 1969).Radiation absorbed by the canopy was partitioned between the two plant layers (Fig. 1), depending on their height and surface cover, whereas it was assumed that leaves are uniformly distributed within the total height of the canopy.Interception and plant evaporation depended on the simulated leaf area index of the vegetation as well as the degree of area coverage.
Transpiration depended additionally on the simulated plant water uptake.Soil evaporation was derived from an iterative solution of the soil surface energy balance of the soil surface, using an empirical parameter for estimating the vapour pressure and temperature at the soil surface.Vapour pressure deficit was calculated from the relative humidity input.Snow fall was simulated from precipitation and air temperature, while snowmelt was estimated from global radiation, air temperature and simulated soil heat flux.Fluxes were calculated by the EcoFlux software (In Situ Flux AB, Ockelbo, Sweden) according to the EUROFLUX methodology (Aubinet et al., 1999;Sagerfors et al., 2008;Nilsson et al., 2008) 1

Soil temperatures and heat fluxes
Surface temperature was simulated based on an energy balance approach, where the radiation reaching the soil equals the sum of sensible and latent heat flux to the air and heat flux to the soil.The same approach was used for the snow surface temperature.Heat flow between adjacent soil layers were calculated based on thermal conductivity functions accounting for the content of ice and water.The heat flow equation is based on a coupled equation also accounting for freezing and thawing in the soil (Jansson and Halldin, 1979).Convection heat flows were not accounted for.The lower boundary temperature was calculated based on a sine variation including parameters for the annual mean temperature and amplitude at the site.

Soil hydrology
Soil water flows and water contents were calculated for each of the 16 soil layers.Soil water depended on infiltration to the soil, soil evaporation, water uptake by plants and groundwater flow.Soil moisture represented as liquid water content was calculated based on the water storage and temperature in the corresponding soil layer.Water flows between adjacent soil layers were calculated based on Richards' equation (Richards, 1931), considering hydraulic conductivity, water potential gradient and vapour diffusion.Saturation conductivity was assigned depending on mean-measured dry bulk density values of the corresponding layers (see Päivänen, 1973).With respect to hydrologic characteristics, the soil profile was divided in two horizons representing the acrotelm and the catotelm (see Ivanov et al., 1981), whereas the boundary between these horizons was positioned at −30 cm as suggested for Degerö Stormyr, based on visual differences in the soil profile and water table depth measurements (Granberg et al., 1999).The soil water characteristics were described by the Brooks and Corey equation (Brooks and Corey, 1964) and unsaturated conductivity by the Mualem function (Mualem, 1976).When the current simulated groundwater table is above the assumed drainage level, outflow of saturated layers above that level was simulated, based on a linear model.
Surface runoff was controlled by a surface pool of water that covers various fractions of the soil surface.During periods of a fully saturated soil profile the flow of water in the upper soil compartment could be directed upwards, towards the surface pool.Surface runoff was calculated as a function of the amount of water in the surface pool.

Vegetation
Two plant layers were simulated, representing vascular plants and mosses.They differed in their parameters for size, shape, carbon allocation, litter fall and temperature response for assimilation and respiration.A detailed description of the carbon pools of the two plant types and the partitioning of assimilates to the pools can be found in the Supplement.Vascular plants additionally had a pool for mobile reserves, which was filled during litter fall.They were assumed to have a maximal height of 50 cm compared to 2 cm for mosses.
Plant development was temperature-sum and day-length dependent.Senescence and litter fall for vascular plants depended on growth stage, temperature sum and day length.In the case of mosses, litter fall was proportional to assimilation.Litter from above-ground carbon pools went through a surface litter pool and then to the upper soil litter pool, and litter from below ground to the corresponding soil layer.In the case of mosses, litter fall occurred only in below-ground parts.For both plant types, assimilation was simulated using the light-use efficiency approach (see Monteith, 1972), where total plant growth is proportional to the net of global radiation absorbed by the canopy but limited by unfavourable temperature and limited soil water.The response to soil water was defined from the ratio of actual to potential transpiration.Potential transpiration depended on vapour pressure, temperature, wind speed and aerodynamic resistance of the plant.Actual transpiration was assumed to equal water uptake from soil layers, depending on the relative amount of roots, the specific response to soil water potential and the soil temperature of each layer.Both plant layers were assumed to be well adapted to wet conditions (see Keddy, 1992;Steed et al., 2002) and therefore experiencing water stress only due to dry conditions, which was supported by pre-study modelling results.Plant respiration was assumed to be proportional to assimilation (growth respiration) and to the amount of biomass (maintenance respiration), whereas maintenance respiration also depended on temperature through a simple Q 10 (constant factor changing the rate with a 10 • C increase in temperature) approach.

SOC decomposition
The organic substrate was represented by three C and N pools for each of the 16 soil layers: one representing more stable, partly decomposed material (SOM s ), one representing fresh or little decomposed moss litter (SOM m ) and one representing fresh or little decomposed litter from vascular plants (SOM v ).Initial conditions were selected to fulfil the measured total carbon per layer and partitioned into the pools in a way that they were approximately in equilibrium for a certain parameter combination that produces a reasonable fit to NEE (prior calibration).Decomposition followed first-order kinetics with pool-specific rates that were reduced under unfavourable soil temperature and moisture conditions.Temperature dependence was described by the Ratkowsky function, which was originally developed for bacteria (Ratkowsky et al., 1982) but has also been applied to fungal growth by Bazin and Prosser (1988).Soil moisture response was zero at moisture contents below the wilting point, rising to 100 % between two threshold moisture contents and falling to a certain level under saturated conditions.
Decomposition products from the SOM m and SOM v pools were partitioned into CO 2 that was released to the atmosphere and C that is partly moved to the SOM s pools and partly returned to the SOM m and SOM v pools.Decomposition products from the SOM s pools were partly released as CO 2 and partly returned to the SOM s pools.Under saturated conditions, carbon could leave the pools as methane (CH 4 ), which was later oxidized to CO 2 or transported to the atmosphere via plants or through ebullition.Nitrogen-and methane-related processes were considered by a model including the most important pathways and fluxes, but no emphasize on the calibration of these processes were made in this study.
Peat depth growth during the simulation period was considered by the following; the initial organic concentration was preserved for each layer except for the lowest in the profile.Instead, the difference in the total amount of C in all pools in one layer between the beginning and end of each year was moved to or from the layer below, to simulate growth or decrease the peat depth.Thereby, carbon was taken from the different pools according to the relative abundance of each pool in the source layer and inserted to the corresponding pool in the target layer to allow for dynamic changes in litter quality.The lowest layer (−2.8 to −3.4 m below the surface) represented the entire depth change of the whole profile, but was excluded from a constant concentration to avoid adjustments of the number of layers.

Calibration procedure
A Monte Carlo calibration including acceptance criteria was performed to identify process and parameter interactions.The resulting parameterizations were analysed for correlations between different parameters, between parameters and model performance and between performances in different variables.A total of 50 000 runs were performed to calibrate 54 parameters from different processes.Parameter values were randomly assigned from a uniform distribution within assumed prior ranges (i.e.all values had the same probability of being used).The parameters were selected as candidates to demonstrate the role of various regulating processes, which we group into eight different process categories that describe (1) plant growth, (2) decomposition, (3) radiation interception, (4) soil temperature, (5) aerodynamic resistance, (6) transpiration, (7) soil hydrology and (8) snow.Prior ranges for calibrated parameters were selected according to literature values or experiences from previous model runs, in most cases a certain range around the default value (Table S1 in the Supplement).Many parameters were still considered with fixed single values (Table S3 in the Supplement).Model outputs were compared with measured field data including many variables in high temporal resolution, spanning up to 12 years of observations (Table 2).Several combined criteria were defined to select runs (behavioural models) with an acceptable performance (see Sect. 2.4.2) in different variables.Resulting parameter value ranges of the accepted runs were then compared with the prior ranges and between the different criteria selections to examine the effect of the criteria selection.Correlations between parameter values and model performance in the different measurement variables were analysed, as well as between accepted values of different parameters.Parameters were ranked on their effect on model performance, their correlation with other parameters and their constraint ability from the available data.

Splitting of calibration variables into sub-periods
Additional to the calibration data for the whole period, we introduced further sub-variables for certain sub-periods and times of the day.NEE was separated into night-time values (22:30-02:30 CET), representing ecosystem respiration, and daytime values (09:30-15:30 CET), representing the sum of the respiration component and the assimilation component.Additionally, springtime values were considered separately for NEE and snow depth, and spring-and wintertime values for Rn, Ts, H and LE.This is justified as low values with few dynamics during winter, and the critical transition of plant emerge and snowmelt in spring might not be properly accounted for if only the whole period was considered.WT was calibrated and analysed in the whole profile and additionally in lower soil layers (one sub-variable for WT depths > −0.15 m and one for > −0.2 m).This was encouraging because as WT in the upper soil layers showed high fluctuations in the modelled and also partly the measured WT, and while our interest was to achieve a good overall water table with good representations of dry summer periods.

Performance indices
The selection of runs and evaluation of model performance were based on three indices: coefficient of determination (R 2 ) assess how well the dynamics in the measurementderived values are represented by the model; mean error (ME) is the difference between the average of the simulated compared to the average in the measured (i.e. it shows the error in the magnitude); and the Nash-Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970) accounts for both deviation of dynamics and magnitude (i.e. it ranges from −∞ to 1, whereas 1 means the best fit of modelled to measured data).Values < 0 indicate that the mean-measured value is a better predictor than the simulated value (Moriasi et al., 2007).As the NSE may be understood as a combination of the coefficient of determination (R 2 ) and ME, it was only evaluated if the R 2 and ME alone did not narrow the parameter range.
The decoupling of turbulent transport and biological activity during night-time may introduce spike-type fluxes in NEE, LE and H if accumulated concentrations during calm night-time conditions are released during the onset of turbulence in the morning hours.To attenuate the effect of these spikes, the simulated and measured values were transformed to cumulated total amounts, starting from the beginning of the observation period.An additional R 2 value was calculated for the cumulated values (AR 2 ).

Criteria for posterior selection
Criteria were applied in two steps.In the first step, a basic set of 1285 behavioural models was selected.Out of these, several sets of 50 runs each were selected in the second step in two different ways: one for sensitivity analyses and parameter ranges, which was based on single criteria, and the other for identification of equifinalities, based on multiple criteria.

Basic selection
The basic selection was applied, as the lowest summer water levels and a reasonable representation of the plant was assumed to be crucial for most of the processes of interest.Criteria were on performance in WT and vascular plant LAI (Table 3).The criteria on water level below 0.2 m was chosen, as a correct representation of summer drought conditions was of higher interest in this study than a correct water level during, e.g., frozen conditions in winter, causing the water table to drop down to 0.15 m.The criteria on LAI ME of ±0.2 m 2 m −2 was a relatively wide range, as the mean of measured values was 0.4 m 2 m −2 ; i.e. a underestimation of LAI by −0.2 m 2 m −2 would result in a maximum LAI of 0.2-0.4,which was close to the minimum for being able to re-establish new biomass after a low productive year.A wide www.geosci-model-dev.net/9/4313/2016/Geosci.Model Dev., 9, 4313-4338, 2016 range of daytime NEE ME was additionally applied to exclude outliers due to numerical problems, which reached an ME in NEE up to 8 × 10 27 g CO 2 -C day −1 m −2 in the prior.

Single criteria to identify parameter range
For sensitivity analyses and to test if, and how, parameter ranges depend on the selected criteria, the best 50 behavioural models for each performance index of each variable were selected out of the basic selection.Thereby, best means highest in the case of the R 2 and NSE, but closest to zero in the case of ME.We defined posterior parameter range as the interval between the 5th and the 95th percentile of the distribution of parameter values of the runs selected.Posterior parameter ranges were compared with the ranges resulting from the basic selection.If the upper or lower limit of a posterior parameter range of the final selections differed by ≥ 10 % from the upper or lower limit of the posterior range of the basic selection, the parameter was assumed to be sensitive to the selected criteria and further analysed.
The same was done for each of the best 200 behavioural models, but as the results were similar, they were only plotted with respect to parameter ranges.Furthermore, all parameters were plotted against all performance indices of each variable and checked visually for discrepancies with the resulting ranges (results are not shown).

Multiple criteria to identify parameter correlations
For identification of equifinalities, a set of multiple criteria for each variable (Table 3) was applied to select sets of 50 behavioural models each.Again, these selections were based on the basic selection.Parameter ensembles of these accepted behavioural models were then analysed to identify covariance between parameters.A pair of parameters was considered to interact if their values correlated with a R 2 of at least 0.1 in the basic selection, respectively, 0.2 in the final selection.If a pair showed correlations in several criteria sets, the highest R 2 value was reported in the results.

Evaluation and measures
To rank the parameters in their concern, several measures were used to quantify parameter sensitivities and constrainabilities, as well as equifinalities.The sensitivity (S) of a parameter to each performance index of each variable was quantified by the sum of the differences between the posterior range and prior range (range reduction).If a parameter was sensitive to more than one period of each variable, the highest value for each variable was chosen for further analysis.To identify trade-offs and supporting effects between different criteria, correlations of the performances between different variables and indices were plotted and visually analysed.Due to limited computer capacity, this was based on a random set of 3200 runs.Furthermore, the parameter value ranges resulting from the different criteria were compared with each other and determined how well they were overlapping, i.e. how unambiguously they could be constraint.Overlap (O) for each parameter was defined as the difference between the minimum of the upper limits of the posterior ranges of the different criteria minus the maximum of the lower limits of posterior ranges and, therefore, become negative if ranges were not overlapping.Furthermore, it was compared how well overlapping ranges differed between performance indices within the same variable and between different variables.The overlapping range of each parameter was normalized by dividing it by the average of the posterior ranges of this parameter; therefore, a value of 1 would be reached if all posterior ranges of that parameter would be identical for all performance indices and variables.Equifinalities were quantified by the R 2 value of a simple linear regression through the values of the interacting parameter pair in the accepted runs.Parameter concern (P ) was defined based on three components: the sensitivity of the parameter, how unambiguously it could be constraint and the sum of correlation coefficients of equifinalities with other parameters: . (1) Thereby, sensitivity was the sum of the range reduction for the R 2 and for ME, respectively NSE in case where no sensitivity was detected for the R 2 and ME but for the NSE.
The sensitivity was multiplied by the factor 1 minus the normalized overlapping range, so that the sensitivity of parameters, which could be unambiguously constrained are down weighted; therefore, with high uncertainty due to different results for different performances, indices or variables are up weighted.Equifinalities were considered by the sum of the R 2 values for each correlation of that parameter with another parameter, displayed in exponential form and weighted so that strong correlations were emphasized and the contribution of equifinalities were in a comparable scale to the sensitivity measures.

Results
Processes as well as parameters were strongly interacting, which was reflected in sensitivities of each variable to several different process categories, correlations between the performance in different variables, and in equifinalities between parameters of different process categories.About half of the parameters were sensitive to model performance in one or more variables, but only a few had a distinct range (Sect.3.1).Instead they affected several processes, causing trade-offs not only in model performance between the different measurement variables and between the different performance indices, but also several supporting effects could be identified (Sect.3.2).A lot of equifinalities were identified between parameters.Parameters were correlated with up to seven other parameters, often from different process categories.Therefore, a good performance often requires certain combinations of parameter values, rather than specific parameter values (Sect.3.3).Each of the available measurement variables (NEE, LAI, sensible and latent heat fluxes, net radiation, soil temperatures, water table depth and snow depth) constrained parameters from several different process categories, without any variable being redundant (Sect.3.4).Nevertheless, large uncertainty remained in especially the unsaturated water distribution (ψ a ) in the soil (Fig. 2), which affected all considered processes and hindered further parameter constrain.This might be solved by additional measurements of, e.g., soil hydraulic properties.Other important parameters that could not be constrained were the following: define aerodynamic resistance, radiation interception (in particular moss albedo), timing of snowmelt, and in the case of NEE mostly the leaf-litter fall rate of vascular plants during the growing season (Fig. 2).A detailed description of the key parameters for each process and the detected interactions can be found in Sect.3.5.Results for model fits to the different variables can be found in Fig. S1 in the Supplement.

Parameter sensitivity
Model performance was sensitive to parameters across the different process categories: out of 27 sensitive parameters 21 affected model performance in more than one variable.
For 15 of the sensitive parameters, resulting value ranges differed strongly (less than 50 % overlapping range), depending on both, the variable and the performance index (Fig. S2 in the Supplement).Performance in Ts and WT was determined by 12 key parameters belonging to seven and six different process categories, respectively (Fig. 3).In contrast, snow depth and LAI depended mainly on parameters from their own process categories.Radiation and LAI refer to the simplest processes with respect to number of connected parameters (Fig. 3).However, radiation, together with snow depth, was the variable with the strongest average disagreement in parameter value ranges between the different selection criteria (Fig. 4).Four parameters were sensitive to at least half of the considered variables (Fig. 2); the parameter defining the water retention curve and unsaturated soil hydraulic conductivity (ψ a ) affected model performance in variables of all eight considered variables.The moss transpiration coefficient (g max,moss ), vascular plant respiration coefficient (k gresp,vasc ) and litter fall rate (l Lc1 ) were important parameters for not only LAI and NEE, but also H , LE and WT.Furthermore, g max,moss and k gresp,vasc were also important for Ts.The sensitivities of the single parameters are described in more detail in Sect.3.5.The full table of the correlation coefficients between parameters and performance can be found in the Supplement (Table S4).

Confounding and supporting effects of interacting processes
The performances of several variables were connected in supporting and co-founding ways (Figs. 5 and 6).Trade-offs existed not only between the performances of different variables but also within a variable, depending on chosen performance index or sub-period.This was also reflected in the large differences in resulting accepted ranges.accepted value ranges overlapped with 35 % between different performance indices and between different sub-periods of the same variable, and additionally with 6 % if the differences between different variables were considered (Fig. 4).In the case of 11 parameters, the accepted ranges did not overlap at all (Fig. S2 in the Supplement).
Strong connections existed especially between ME of LE and WT, but also the ME of LAI had an impact on the per- formance in many other variables; the magnitude of vascular plant LAI was strongly correlated with magnitude of LE, WT, H and NEE, especially if daytime and night-time values were considered (Fig. 5).Thereby the lowest ME in daytime and night-time NEE, as well as the ME and dynamics of H , went along with a slight underestimation, and for LE and WT with a slight overestimation of vascular plant LAI.The best performance for WT dynamics was reached if the magnitude of vascular plant LAI was correct (Fig. 6).A noticeable existence of the vascular plants (LAI ME > −0.4) increased the fit in the NEE R 2 to at least 0.2, but this was not a necessary precondition for good NEE performance (Fig. 6).The highest performance in dynamics of WT, H and Ts in the upper layer coincided with a good fit in NEE magnitude (Fig. 6).This relationship was even stronger if these variables were compared to ME in NEE night-time and NEE daytime.A correct representation of WT dynamics and depth coincided with high performance in H dynamics and a correct or slightly underestimated H (Figs. 5 and 6).A small ME in H correlated with high performance in WT dynamics.Performances in soil temperatures of different layers were strongly correlated with each other in both dynamics and magnitude.
Underestimation of LE was connected to an overestimation of H , but also to better dynamics in H (Fig. 5).The ME in net radiation was positively correlated with the ME in H .A good fit between modelled and observed snow depth did not correlate with the performance in any other variable.The only exception was a negative correlation between the dynamics in snow depth and H if performance during springtime exclusively was considered (Fig. S3 in the Supplement).
Especially for snow, Rn and in the case of some parameters also for Ts, accepted ranges were contradictory depending on whether the R 2 or ME was chosen.In the case of moss albedo (a pve,moss ) and aerodynamic resistance dependency on LAI (r alai ), the ranges also strongly depended on the season during which the variable was considered.For two aerodynamic resistances and one soil parameter (z 0M,snow , c H 0,canopy , s k ) ranges differed between the R 2 of actual values and the R 2 of accumulated values.In addition to the uncertainty from unambiguous parameter ranges, further uncertainty results from equifinalities between parameters.

Equifinalities
Parameters were strongly inter-correlated, often with several parameters, and often across different process categories.Equifinalities can hinder the identification of sensitivities, which was especially true for the basic selection; despite reducing the number of runs by 97.5 %, posterior and prior ranges hardly differed (Table S5 in the Supplement).Instead certain-value triples for photosynthetic efficiency (ε L,vasc ) with the respiration coefficient (k gresp,vasc ) and with the storage fraction for plant regrowth in spring (m retain ) were crucial for the survival of the vascular plant layer.Certain-value pairs for the moss transpiration coefficient (g max,moss ) with the shape parameter of soil water retention (ψ a ) were crucial for a reasonable water table depth.
Equifinalities existed not only between parameters from the same process categories, but even more often between parameters from different process categories (Fig. 7).Parameters defining radiation interception, soil temperature, aerodynamic resistance, transpiration, and soil hydrology exclusively correlated with parameters from different process categories.Parameters defining radiation interception were mostly correlated with parameters defining aerodynamic resistance.Only in the case of plant and SOC decomposition parameters, equifinalities existed mainly between parameters of the same process category.Except for ρ smin , all sensitive parameters and other further parameters were detected to correlate with up to five other parameters in the final selections, ψ a correlated with even seven others (Fig. 2).Two parameters had very strong correlations (R 2 ≥ 0.3) with two other parameters each, which belong to different process categories (ψ a with c H 0,canopy and g max,moss and a pve,moss with z 0M,snow and r alai ) (Table S6 in the Supplement).

Usefulness of measurement variables
All available measured variables (NEE, LAI, LE, H , Rn, Ts, WT and snow depth) were helpful in constraining parameter ranges (Fig. 2).None of the supporting effects was strong enough, to allow one variable to be fully replaceable by another.Even for the strongest correlation between soil temperatures of the different layers, the remaining uncertainty in one temperature when knowing the other would be of the magnitude of 0.5 • C, which corresponds to more than 25 % of the total uncertainty resulting from the tested parameter ranges (Fig. 5).In the case of 15 variables, the usage of several variables revealed that constrained ranges were not robust.In total, 12 parameters could be unambiguously constrained to a more narrow range, as their resulting ranges had at least 50 % overlap, or affected only one variable (Fig. S2 in the Supplement).Each variable constrained parameters from several different process categories (Fig. 3).The highest number of correlations was detected for the performance in WT and Ts, which constrained 12 parameters from different process categories.Also the available data for LE, H and NEE constrained many parameters.Nevertheless, large uncertainty remained due to equifinalities and differences in accepted ranges; the largest uncertainty was caused by a parameter defining the shape of the water retention curve (air entry, ψ a ).As this was the only calibrated parameter of the water retention curve, it determined the unsaturated hydraulic conductivity of the soil.ψ a was sensitive to all considered variables and had many strong interactions with other parameters, while it was not possible to constrain it to an unambiguous value range (Fig. S2 in the Supplement).Therefore, it would be of great value to be able to deduce such parameters from additional measurements.This also applies to the following parameters, which could not be constrained unambiguously; leaf-litter fall rate of vascular plants during the growing season (l Lc1 ) was the second most sensitive parameter, affecting the performance in NEE, H , LE and WT.Moss albedo (a pve,moss ), aerodynamic resistance dependency on LAI (r alai ) and transpiration coefficients (g max,vasc , g max,moss , g maxwin ) had a similar concern, due to their equifinalities to other parameters.Plant respiration (k gresp,vasc ) had strong sensitivity but could be constrained unambiguously by the available data.

Detailed description of sensitivities and interactions per process
Detected sensitivities, connections between performances, and equifinalities all showed strong interactions between the different processes and parameters of different process categories.Connections existed between all variables and process categories, but most strongly inter-linked were LE with WT and Rn with H and Ts (Fig. 2).H , LE and WT were also linked to each other and to NEE.The impact of the plant is further reflected in the correlations between performances in LAI with performances in many other variables (Fig. 5).The implications on the performance for each considered variable will be described in the following sections.

Water level depth and soil moisture conditions
Performance in water level depth was determined by 12 key parameters (Table S4 in the Supplement).It was most strongly connected to the shape of the soil water retention curve (ψ a ) as well as to the transpiration coefficients for mosses and winter transpiration (g max,moss , g maxwin ).
The transpiration coefficient from vascular plants played a smaller role due to the high sensitivities of parameters defining the growth and therefore the magnitude of the vascular plant (i.e.k gresp,vasc , m retain , l Lc1 ).Equifinalities existed between several of these parameters.ψ a had a strong effect on the performance of all variables and several strong equifinalities, in particular with parameters defining aerodynamic resistance and transpiration.On the other hand, ψ a could not be constrained to an unambiguous range and was therefore the parameter causing the largest overall uncertainty (Fig. 2).Performance in WT was further sensitive to parameters defining aerodynamic resistance, i.e. r alai and c H 0,canopy .Both parameters had equifinalities with ψ a and moss albedo (a pve,moss ) as well as with timing of snowmelt (m T ) and thermal conductivity of snow (s k ).In addition, the distance between drainage (d p ) showed some sensitivity.

Transpiration and evaporation
The nine most important parameters for WT performance were also key parameters for LE (ψ a , g max,vasc , g max,moss , g maxwin , k gresp,vasc , m retain , l Lc1 , r alai , c H 0,canopy ).This explains the strong correlation between the performance in WT and LE ME (Fig. 5) and shows the connections with plants, WT and H . Another parameter, sensitive to LE was the roughness length of snow (z 0M,snow ), belonging to the aerodynamic resistance process category and correlating with moss albedo, hinting to the connections between LE-and R-associated processes.Dynamics in WT and LE, as well as magnitude of H , was improved if the transpiration coefficient was on its lower range in the case of mosses and on www.geosci-model-dev.net/9/4313/2016/Geosci.Model Dev., 9, 4313-4338, 2016 its upper range in the case of vascular plants (Fig. S2 in the Supplement).Despite the lower values for mosses, transpiration prior to criteria selection was dominated by mosses, due to their higher LAI and coverage (Fig. S4 in the Supplement).Crucial for LE performance was also a parameter defining the aerodynamic resistance of the canopy under stabile conditions (c H 0,canopy ): a very small value improved the R 2 of LE and spring LE, but downgraded the R 2 of accumulated LE and of winter radiation.Spring LE was overestimated in most of the runs (see Fig. S1 in the Supplement).
The strongest sensitivity on spring LE was by the coefficient for winter transpiration (g maxwin ): the higher the g maxwin the better the R 2 and ME.Together with (z 0M,snow ) this was also the most important parameter for winter LE.

NEE and LAI
Seven of the nine parameters, which were common for LE and WT, were also among the most effective parameters for NEE (ψ a , g max,moss , g max,vasc , k gresp,vasc , m retain , l Lc1 , r alai ) and belong to four different process categories: plant, transpiration, soil hydrology and aerodynamic resistance (Table S4 in the Supplement).However, the most sensitive parameter for NEE was the rate coefficient for heterotrophic respiration (k l1 ), which was especially important for night-time NEE.Further sensitive parameters for night-time NEE were the growth respiration coefficient for mosses (k gresp,moss ) and the temperature dependency coefficient for heterotrophic respiration (t min .The rates of photosynthesis and its temperature dependence (ε L,vasc , ε Lmoss , p mn,vasc ) were key parameters for LAI, NEE magnitude or temporal NEE dynamics, respectively.Many strong interactions existed between plant parameters, which were especially visible in the basic selection (see Sect. 3.3).The rate of leaf-litter fall during the growing season l Lc1 was one of the parameters with the highest concern, due to its sensitivity on many different processes, its equifinalities and as it could not be constrained to an unambiguous solution (Fig. 2).Resulting ranges for l Lc1 differed not only between the different performance indices within NEE and within LAI, but also between NEE and LAI (Fig. S2 in the Supplement).

Sensible heat fluxes, soil temperatures and net radiation
Many inter-connections existed between H , Ts and Rn but all three were also linked with LE, WT, snow and NEE.A snow parameter, determining the timing of snowmelt (m T ) was the most crucial parameter for heat fluxes, not only in springtime but also for the whole year period.Furthermore, m T was important for Ts in springtime (see Sect. 3.5.5).The shape of the soil water retention curve (ψ a ) was the second most sensitive parameter for both variables.The aerodynamic resistance dependency factor on LAI (r alai ) was the most sensitive parameter for Ts, and affected also LE, WT and night-time NEE, while it strongly correlated with moss albedo (a pve,moss ), the third most sensitive parameter for H and most sensitive parameter for Rn.The accepted ranges for r alai contradicted within the soil temperature variables, depending on the chosen performance index and considered season; high values were important for the Ts ME and R 2 during winter, but low ones improved the Ts R 2 during spring and during the whole period.Therefore, r alai was the parameter causing the largest overall uncertainty after ψ a .This was followed by a pve,moss , which had low values for accepted ranges in the case of H , Rn and Ts during the whole period but high values in the case of winter H and Rn.It further showed strong equifinalities with the roughness length of snow (z 0M,snow ), which was the second most sensitive parameter for Rn, but also affected H and LE.The coefficient for thermal conductivity of snow (s k ) affected Rn and Ts but not H .The thermal conductance coefficient of soil organic material (h 2 ), the lower boundary mean temperature (T amean ), the snowmelt dependency to radiation coefficient (m Rmin ) and the density of new and old snow (ρ smin , S dw ) affected only soil temperatures, the latter two also snow depth.Parameters defining moss and winter transpiration (g max,moss , g maxwin ) and the growth respiration coefficient of vascular plants with its effect on vascular plant biomass and LAI (k gresp,vasc ) were sensitive to Ts, g max,moss and k gresp,vasc as well as to H .The most important parameter for LE, c H 0,canopy was another key parameter for Rn and H .

Snow
The temperature coefficient in the snowmelt function (m T ) was the most important parameter for ME in snow and determined the timing of snowmelt.However, resulting parameter ranges did not overlap between the different performance indices within the snow-depth variable and between different other variables.A longer lasting snow cover (low m T < 3) was crucial for spring H and reduced mean error in snow depth, but lowered R 2 values in spring Ts and snow depth.m T interacted with another snow parameter (T RainL ) as well as with parameters from the temperature and transpiration process category (T amean , g maxwin ).The density coefficients for old (S dw ) and new snow (ρ smin ) had medium effect on snow-depth performance and also affected spring and winter soil temperatures in all layers, but the latter could be unambiguously constrained by the available data.

Discussion
Unlike many previous sensitivity studies for carbon modelling that often focus on only one or a few calibration variables and parameters of the associated process category, we considered many different abiotic and biotic measurements (NEE, LAI, Rn, Ts, H , LE, WT and snow depth) to investigate the interactions between various process categories Geosci.Model Dev., 9, 4313-4338, 2016 www.geosci-model-dev.net/9/4313/2016/(SOC decomposition, plant growth-related processes, radiation interception, soil temperature, aerodynamic resistance, transpiration, soil hydrology and snow) in a peatland ecosystem.Similar to results from a forest modelling study using the DRAINMOD-Forest model (Tian et al., 2014) and a N 2 O study using the CoupModel on a drained peatland forest (He et al., 2016), we found that processes were sensitive to parameters from several different process categories.Together with the discovered supporting effects between model performances in different variables, this confirms the connections and dependencies between different processes as implemented in the model (see "Model description and equations", Sect.2.3, Table 2 in the Supplement and Janson and Karlberg, 2010).The many equifinalities within and between different process categories reveal the dependency of constrained parameter ranges as well as parameter sensitivities to model structure, calibration set-up and parameters with fixed values; a deviation in one of these factors leads to different optimal value ranges, whereas a non-sensitive parameter might become sensitive if an interacting parameter is set constant.This implies a limited transferability of parameter values between models in general and even between studies using the same model in a different configuration.Resulting parameter ranges were moreover affected by the applied criteria for selecting runs.Yet, it is quite common practice to adopt at least some parameter values from other modelling studies (e.g.Frolking et al., 2002, Yurova et al., 2007;St-Hilaire et al., 2010;Wania et al., 2010;Gong et al., 2013;Kim et al., 2014;Kurnianto et al., 2014;Zhu et al., 2014), which includes the usage of model default values that were estimated under a different model configuration.
The strong interactions across different process categories also emphasize the importance of measurements of ancillary data additionally to the variable of interest and model input data (meteorological and SOC data).Measurements of NEE, LAI, LE, H , Rn, Ts, WT and snow were all found to be valuable for constraining parameters from several different process categories and can therefore reduce uncertainty in model predictions.Further constraint of the parameters in this study would be possible if especially additional water content or soil hydraulic properties were measured.
Beside parameter uncertainty, also uncertainty in model structure and in measured input and calibration data contribute to model uncertainty (Thorsen et al., 2001;Beven and Freer, 2001).This was tested for other peatland models (e.g.model structure: Tang et al., 2015;input drivers: Wania et al., 2009;St-Hilaire et al., 2010;Grant et al., 2011;Kim et al., 2014), but goes beyond the scope of this study.Here, only one model and one site were investigated.A previous study using CoupModel investigated the differences of parameter ranges between several different peatland sites (Metzger et al., 2015).

Parameter sensitivity
The sensitivity of variables to parameters from many different processes revealed the importance of process interactions.Especially abiotic processes were strongly inter-linked, but also biotic variables showed sensitivities to parameters from up to seven different process categories, suggesting that parameter sensitivities and model performance of a certain process depend on which other process categories are considered in a model and in a calibration.This is an important finding, as many studies investigate the sensitivity of often only a few parameters from mainly the same process category as the output variable (e.g.Yu et al., 2001;Frolking et al., 2002;Belassen et al., 2010;Wania et al., 2010;Morris et al., 2012;Wu and Blodau, 2013;Zhao et al., 2014;Zhu et al., 2014), which might lead to sensitivities and resulting ranges that are not robust.The identified interactions can help modellers to develop or select an appropriate model including the parameters, processes and process categories that need to be considered together, depending on the variable of interest.
Parameter sensitivity analyses can also help to simplify future calibrations (Saltelli et al., 2000), by focussing on the most striking parameters and narrowing the ranges for parameter that could be successfully constrained.Although while the existence of interactions between the processes and their parameters is supposed to be less dependent on site conditions and model structure, the exact shape of the connections, constraint parameter ranges, as well as the relevance of the specific processes and the specific interactions might strongly depend on these factors.Still, one or more of the following parameters that we identified as most influential, correspond to key parameters in other studies using other models and partly different ecosystems; the respiration rate coefficients, radiation use efficiency, transpiration coefficients or the soil water retention capacity were among the most sensitive parameters for NEE, its components or yield, respectively, in, e.g., the PCARS (Frolking et al., 2002) and the GUESS-ROMUL (Yurova et al., 2007) model on peatland, the SiB v2.5 model on a forest area including some wetlands (Prihodko et al., 2008), the LPJ-GUESS model on forest and herbaceous vegetation (Pappas et al., 2013), the EPIC model on cropland (Wang et al., 2005), the BIOME-BGC model for different tree species (Tatarinov and Cienciala, 2006), or the ACASA (Staudt et al., 2010), the 3-PG (Esprey et al., 2004;Xenakis et al., 2008), the FORUG (Verbeeck et al., 2006) or the DRAINMOD-FOREST (Tian et al., 2014) model on forest.These sensitivities seem to be therefore quite independent of model structure, included processes and parameters used for calibration and apply to different types of ecosystems.The resulting value ranges of these parameters should be compared between ecosystems and models to find out to what extent they can be related to site conditions and therefore used for predictions and upscaling.They might be connected to the environmental scenario (Hidy et al., 2012;Ben Thouhami et al., 2013;Sulman et al., 2013) and the chosen prior distributions of the parameters (e.g.Tatarinov and Cienciala, 2006).Furthermore, our results have shown that the parameter ranges depend on model structure, on the selection of parameters for calibration and on the selected acceptance criteria.Thereby, not only the selected variable but also the selected sub-period was relevant, as has been shown by other studies as well (e.g.Prihodko et al., 2008;Van Huisteden et al., 2009;Safta et al., 2015).

Confounding and supporting effects of interacting processes
Criteria selection is a subjective choice of the modeller if multiple output variables are available.The identified supporting effects and trade-offs between the performances in different variables allow modellers to assess the implications of a certain criteria on model performance and parameter ranges and to choose criteria according to the processes of interest; however, some of them might be ecosystem or model specific.Trade-offs existed not only between different variables but also within the same variable, depending on whether ME, R 2 of actual or R 2 of accumulated values were chosen and which season was considered.This implies that the problems of a subjective criteria selection also exist if only one time series variable is considered.Even if a standardized multi-criteria optimization algorithm like Bayesian calibration or a more sophisticated performance index combining several performance measures is used, the choices and the corresponding weightings are moved to the developer of the algorithm or index, but still remain subjective.More than half of the sensitive parameters in this study could not be constrained to an unambiguous range.Constraining such a parameter by only one variable and one index would result in a range that is not robust.Using several measurement variables and several indices can therefore help to test the robustness of calibrated parameters.A parameter that is robust might better represent a physical constant, whereas controversial resulting ranges might hint at a not well represented system; there is no value for this parameter that leads simultaneously to the best performance for dynamics and magnitude in all variables and during all periods.Instead of a physical constant this parameter might correspond to a dynamic process.Beside model inadequacy, mismatching ranges could be caused in some cases by an inappropriate performance index (see discussion in Sect.4.5.4) or measurements that do not truly represent the modelled variable.For example, with the EC technique, NEE is not directly measured as the CO 2 exchange between biosphere and atmosphere at a certain point, but rather results from calculations of the turbulent exchange of vertical fluxes measured several metres above the ground.Moreover, fluxes may originate from a footprint area that changes diurnally and seasonally and thus may include different soil conditions and vegetation.
Usually, LE is assumed to be closely connected to NEE due to the coupling of transpiration and carbon assimilation in vascular plants (e.g.Schulze, 2006), but has also been shown to correlate for mosses (e.g.Robroek et al., 2009).Our study reveals much stronger relations between parameters defining H and NEE, than between LE and NEE.Tradeoffs between performance in LE and NEE were also found by Staudt et al. (2010) and Prihodko et al. (2008) in a forest and a forest complex including wetlands.However, only the effect of parameters, not the effect of model input (i.e.meteorological input data), on these processes were tested in both studies, as well as in ours.Such a confounding effect might also result from a parameter value compensating for a process not implemented in the model.For example, parameter values that lead to an overestimation of NEE in spring result in higher transpiration and therefore better LE, whereas the reason for the underestimated LE during mid-April to mid-June (Fig. S1 in the Supplement) might in fact be caused by evaporation from open water bodies that form on the peatland during spring and early summer, a process not implemented in the applied version of CoupModel.
The detected supporting effects indicate that some measurement variables can partly compensate absence or low resolution of a connected variable, even though they were not strong enough to make one variable fully redundant.For example, LAI measurements could reduce uncertainty in model predictions of the magnitudes of NEE, LE, H and WT on locations where these variables are not available.Tight relationships between plant and LAI, soil hydrology, C fluxes and soil temperatures have been found by other model sensitivity studies as well (e.g.Ben Thouhami et al., 2013;Quillet et al., 2013;Tian et al., 2014;Sándor et al., 2016) and strong correlations between LAI and NEE (Lund et al., 2010) and between NEE and water availability (Reichstein et al., 2007) have also been found by data syntheses of eddy covariance sites.These relationships can be explained by the many dependencies between LAI and, e.g., photosynthesis, transpiration, heat insulation and water uptake (Schulze, 2006), of which several are also implemented in the model (see model description and equations, Sect.2.3, Table S2 in the Supplement and Jansson and Karlberg, 2010).Other examples for detected supporting effects indicate that if H fluxes are available, the model is constrainable to produce improved WT dynamics, even if WT measurements were missing.High temporal resolution of soil temperature measurements in one layer are sufficient to model good temperatures if just the magnitude of soil temperature in an upper and a lower layer is known, e.g., due to short-time or low-resolution measurements.The knowledge on supporting effects helps modellers in their site selection and in uncertainty estimation of model predictions depending on available ancillary data.It further can help experimentalists to decide which variables should and which need to be measured if the site should be usable for model constraint.

Equifinalities
The fit of model output to measured data in complex models is often not driven by a particular parameter but instead by interactions among parameters (e.g.Beven and Freer, 2001), which was also the case for several parameters in our study, hindering the constraint of parameters to a more narrow range.Furthermore, other carbon modelling studies found that parameter values and sensitivities depend on the values of other parameters (e.g.Tatarinov and Cienciala, 2006;Verbeeck et al., 2006;Quillet et al., 2013).This implies that especially if only a few parameters and processes are calibrated (as in e.g.Yu et al., 2001;Wania et al., 2010, Zhu et al., 2014;Kim et al., 2014, Tang et al., 2015), resulting constrained ranges might not be comparable and transferable between models differing in their constant parameter values.Many equifinalities were identified, not only between parameters from the same process category but also across different process categories.This means that the problem of limited transferability also applies if parameters from only one process category are calibrated (as e.g. in Wang et al., 2005;Belassen et al., 2010;Wania et al., 2010;Sándor et al., 2016) or if models differ in the structures and implementations of their modules.The knowledge on equifinalities is needed for a better parameter constraint in future calibrations as it allows for calibration of the connected parameters dependent on each other.Another way to respond to identified equifinalities is to calibrate only one of the connected parameters.However, the resulting range will then not be transferable to other models using different values for connected, constant parameters.
Some equifinalities included several parameters, making their visualization impossible and simple regression an insufficient tool for fully detecting and describing them (see Saltelli et al., 2008).These equifinalities need to be further investigated in additional calibrations that incorporate those parameter interactions and constrained ranges, which were unambiguous, to achieve a higher number of acceptable runs.This is needed, because the numbers of accepted runs in the final selections (50) did not allow for a much more detailed analysis in such a complex model, as was apparent in comparison with the basic selection; a R 2 threshold value of 0.1 was sufficient to identify equifinalities in the basic selection of 1286 accepted runs, but with just 50 accepted runs in the final selections; this threshold value could easily be exceeded by a random distribution, even when a higher threshold value of 0.15 was used.A threshold of 0.15 was, on the other hand, already too high to detect, for example, the strong relationships between the plant parameters that were only clearly visible in the basic selection.Nevertheless, the six equifinalities with a R 2 of higher than 0.30 are unambiguous in this application of the CoupModel and those with lower values are still very useful to design future calibrations to further investigate and describe these equifinalities.

Usefulness of measurement variables
Models can be improved and their uncertainty reduced by calibrating their parameters to measurement data (e.g.Friend et al., 2007;Wang et al., 2009;Williams et al., 2009).We tested the potential of several measurement variables (NEE, LAI, LE, H , Rn, Ts, WT and snow depth) and found that all contributed to a better parameter constraint.
Thereby none of the variables could be fully replaced by another.Due to the strong interactions and as parameters of each process category were constrained by several different variables, ancillary variables are valuable even if only one certain process is of interest.In the case of snow, our results suggest that data on snow cover might be sufficient if snow depth is not available.
In a forest site simulation with the ORCHIDEE model, H and Rn were found to be redundant for constraining energy balance parameters if NEE and LE were available (Santaren et al., 2007).In contrast, some energy balance-related parameters in our study were constrained exclusively by Rn and H , or additionally by LE, but with different resulting ranges.This reveals the usefulness of Rn and H measurements for model constraints, and shows that variables which might have been identified as redundant in one study could be of high importance on another ecosystem or for another model calibrating a different parameter selection.
Several influential parameters could not be unambiguously constrained or showed equifinalities and need additional measurements to be further investigated.This includes soil water content or soil water retention properties, as well as canopy albedo and leaf-litter fall during the growing season.Except for water retention properties these variables are needed as time series throughout the year.A more detailed discussion of the benefit of such measurements can be found in the following sections.

Detailed discussion of sensitivities and interactions per process
The parameters that were identified as most influential or that showed the strongest equifinalities were related to soil hydrology and water content, to a stable representation of the plant, to radiation, temperature and heat fluxes or to snow.As only one parameter per equation was calibrated, a high sensitivity to this parameter means a high sensitivity to the corresponding process.Some of such process sensitivities might also be interesting for other models and similar ecosystems.The introduced index to measure parameter concern includes subjective choices such as weighting factors, the choice of considered calibration variables and their subperiods as well as the chosen performance indices.However, several tested variations in especially the weighting did not noticeably change the results; ψ a was always the most important parameter, followed by the group of parameters with www.geosci-model-dev.net/9/4313/2016/Geosci.Model Dev., 9, 4313-4338, 2016 medium importance, which differed slightly in their ranking among each other.

Unsaturated water distribution and soil moisture conditions
Our results suggest that model uncertainty could be greatly reduced if data for either soil hydraulic properties, water content or plant transpiration characteristics were available; despite available data of detailed WT and LE in our study, large uncertainty remained in simulated water content due to the combined uncertainty in estimates of soil hydraulic properties (ψ a ) and plant water uptake (g max,vasc , g max,moss , g maxwin ).Their sensitivity to many variables and the high number of equifinalities hindered the constraint of other parameters and therefore the uncertainty reduction in all involved processes.For example this might explain why the water response functions for neither plant assimilation nor soil respiration could be constrained.The shape parameter of the water retention curve (ψ a ) was among the top two most sensitive parameters for NEE, WT, LE, H , Ts, and the third and fifth most sensitive parameter in the case of Rn and snow.That confirms the importance of the implemented interactions of soil moisture with water and heat fluxes, soil temperature, assimilation and respiration processes, as reported from empirical studies (Kim and Verma, 1996;Bridgham et al., 1999;Tezara et al., 1999;Kellner, 2001;Flangan and Johnson, 2005;Lafleur et al., 2005;Schulze, 2006;Belyea, 2009).Furthermore, the transpiration coefficients (g max,vasc , g max,moss , g maxwin ) were among the top 10 most important and influential parameters.In the case of vascular plants, they correspond to the stomatal conductance parameter in other models, which was shown to be crucial for modelling NEE, biomass, LE or H in other studies (Esprey et al., 2004, for forest stand volume;Tatarinov and Cienciala, 2006, for NEE and carbon pools;Staudt et al., 2010, for NEE, LE and H ;Hidy et al., 2012, for carbon fluxes and LE;Bonan et al., 2011, andTian et al., 2014, for LE andH ).The control of stomatal conductance on transpiration and photosynthesis has also been emphasized by several empiric studies (e.g.Jarvis and Morison, 1981;Quick et al., 1992;Tezara et al., 1999;Yordanov et al., 2000).The strong sensitivity of ψ a , g max,vasc , g max,moss , g maxwin for many processes is especially remarkable as parameters and parameter combinations could only vary to such an extent that the water level fit the measurements as restricted by the basic selection.
The importance of the water table on NEE fluxes has been widely mentioned (e.g.Silvola et al., 1996;Yurova et al., 2007;Kurbatova et al., 2009;Dušek et al., 2012) but our results point out that the knowledge on WT alone is not sufficient for model calibration and reliable predictions.In addition, also measurements of soil hydraulic properties are crucial for model calibration.The usefulness of water retention properties for modelling carbon dynamics was also found by other sensitivity analyses on peatlands as well as on min-eral soils (e.g.Wang et al., 2005;Pappas et al., 2013;Quillet et al., 2013).Nevertheless, many of the available peatland sites in current databases (e.g.European Fluxnet Database Cluster, http://gaia.agraria.unitus.it)still do not contain information on water retention properties or water content.We strongly recommend experimentalists to include water retention measurements in their experimental set-up.Thereby, the horizontal and vertical variability in peat hydraulic properties needs to be accounted for (Baird et al., 2012;Waddington et al., 2015).Such measurements might also help to resolve the strong equifinalities of ψ a with transpiration coefficients and a parameter in the calculation of aerodynamic resistance of the plant canopy, defining the minimum exchange under stabile conditions (c H 0,canopy ).

C balance of vascular plants
A stable vascular plant that establishes a reasonable amount of biomass every year throughout the simulation period could only be achieved by certain-value combinations for the photosynthetic efficiency (ε L,vasc ), the respiration coefficient (k gresp,vasc ) and the storage fraction for plant regrowth in spring (m retain ).Despite their high impact in the basic selection, neither equifinalities nor sensitivities of these parameters reached high measures in final selections, probably because several parameters were interacting simultaneously.This indicates the need for either calibrating these parameters dependent on each other or setting at least one of them to a constant value, as the available data were not sufficient to resolve these equifinalities.Many studies on other ecosystems have found NEE or biomass to be strongly sensitive to a parameter corresponding to photosynthetic efficiency (ε L,vasc ) (Esprey et al., 2004;Verbeeck et al., 2006;Prihodko et al., 2008;Staudt et al., 2010;Bonan et al., 2011;Pappas et al., 2013;Tian et al., 2014;Xenakis et al., 2008), but were performed without a simultaneous calibration of parameters related to plant respiration and storage for regrowth.Pappas et al. (2013) discussed a possible overestimation of model sensitivity to photosynthetic efficiency due to processes that are not implemented like the active simulation of plant growth including growth limitations.A strong negative correlation between two of the parameters (plant respiration and photosynthetic efficiency) was also found in a sensitivity analysis using the LPJ model (Zaehle et al., 2005).Despite their effect on model performance, ε L,vasc , k gresp,vasc and m retain had a low rank in parameter concern, as ranges for these parameters could be narrowed unambiguously due to well overlapping ranges between the different variables.Nevertheless, these parameters would be of high importance for predictions if none of the constraining variables are available.
Compared to a previous application of the CoupModel on five different open peatlands including different management intensities (Metzger et al., 2015), vascular plants had to have a much more effective C household to produce the measured leaf area given a limited amount of assimilates.This can be realized by low respiration and litter fall losses and a large storage pool for regrowth in spring.Even if respiration losses from vascular plants were 1/10 of the ones used at the sites in Metzger et al. (2015), the model tended to either underestimate vascular plant LAI or overestimate CO 2 uptake (Fig. 2).A possible explanation for the differences in parameter value combination of vascular plants might lie in the vegetation communities.Although Metzger et al. (2015) included several different types of treeless peatland vegetation communities, none of these sites had a similar vegetation community typical for nutrient poor habitats, consisting of mainly mosses and Eriophorum vaginatum, as at Degerö Stormyr.Eriophorum vaginatum is known to be much more effective in maintaining C compared to other sedges and having a highly efficient remobilization from senescing leaves (Shaver and Laundre, 1997;Jonasson and Chapin, 1985).Uncertainties in measurements and the distribution of modelled respiration over the hours of the day might accelerate or diminish this effect.Explanations by differences in model structure can be excluded, as the same effect was observed when using exactly the same structure (unpublished data).To identify the difference between the sites, which causes the deviations in the combined parameter value ranges, the model needed to be applied to further open peatland sites differing in vegetation community, nutrient status and plant productivity.This might allow for finding trends in parameter ranges, which is a necessary precondition for estimation and reducing model uncertainty in predictions on other peatland sites.
Another plant parameter, which was important for a stable vascular plant layer, and was ranked as one of the overall most important parameters was the rate coefficient for the leaf-litter fall during the growing season (l Lc1 ).Probably due to the high number of correlations with other parameters, these correlations did not exceed the threshold value.l Lc1 is directly connected to the filling of the storage pool, but also for maintaining C in the leaves.The strong sensitivity of LAI to l Lc1 affects transpiration and thereby water uptake, which explains the strong sensitivity to WT depths below −0.2 m and the equifinalities with a transpiration parameter and a parameter describing the response of heterotrophic respiration to water.In Metzger et al. (2015), a value of l Lc1 = 0.01 day −1 could be used site independent.This contradicts the much lower ranges of l Lc1 in our study, necessary for acceptable performance in several variables, in particular the R 2 of LAI, WT depths below −0.2 m and ME of springtime NEE.However, species in nutrient-poor habitats are associated with longer-lived leaves than those of nutrientrich habitats (Ryser, 1996) and fast growing species (Reich et al., 1992), whereas Eriophorum vaginatum in particular is known for long-lived leaves and therefore have a very low litter-fall rate (Jonassson and Chapin, 1985).Less complex models such as the GUESS-ROMUL model, which was also applied to this site, use annual accumulated NEE as estimate for litter fall (Yurova et al., 2007) which is therefore directly dependent on site productivity.Only one site in Met-zger et al. (2015) had lower annual NEE compared to Degerö Stormyr, but this is probably a result of the shorter vegetation period at that site, whereas a site with similar annual NEE was formerly drained, so that the soil respiration contribution to NEE is much larger, compensating for the larger productivity.A high sensitivity of litter fall rate to plant biomass and soil carbon pools was also found by Xenakis et al. (2008) using the 3-PG model on forest.
Further investigations including model applications to additional sites are needed to resolve the differences in resulting ranges and equifinalities with other parameters.Thereby, measurements of leaf-litter fall throughout the year would be of high value.

Sensible heat fluxes, soil temperatures and net radiation
The large number of strong connections between H , Ts and Rn and the equifinalities between their determining parameters indicate the importance to consider, model and calibrate the related processes together.However, the constraint of two of the most important parameters (aerodynamic resistance dependency on LAI, r alai and moss albedo, a pve,moss ) failed not due to different ranges between variables but due to the differences depending on which performance index and season was considered.This emphasizes the importance of the subjective criteria choice, even if only one variable is considered.
Accepted values for r alai were exceptionally high (200 s m −1 for the Ts R 2 and 550 to 800 s m −1 for Ts 1 ME, whereas a r alai of 200 multiplied with the moss LAI of 1.8 leads to an aerodynamic resistance of 360 s m −1 ).Mosses might form a well-insulating layer, but still the values are much higher than the aerodynamic resistance estimates for this site (approximately 50 s m −1 , Peichl et al., 2013) or of a bog in southern Sweden (60 s m −1 , Kellner, 2001).Price (1991) reported very high resistance when moss surface moisture is low, e.g., during dry periods, but these values were still lower than ours.A possible explanation might be an interaction with a non-calibrated, fixed parameter.A highaerodynamic resistance causes better temperature insulation leading to higher summer soil temperatures with lower diurnal oscillations.Furthermore, it leads to strongly reduced soil evaporation and therefore reduced LE, even though this is partly compensated for by slightly higher transpiration from mainly mosses, which profit from the higher water contents in the upper soil layers.This explains the sensitivities to WT and LE, which also supported a higher r alai value.The main cause for the much lower optimum range for dynamics in Ts compared to magnitude in Ts is probably an overestimation of the diurnal amplitude.A lower moss LAI can reduce this overestimation, but the corresponding parameter was not calibrated to avoid further equifinalities: r alai showed already strong interactions with a pve,moss and z 0M,snow .The correlations of the conductivity of organic material (h 2 ) with plant, LE and WT parameters might be explained by the dependency of thermal conductivity from peat wetness (Kellner, 2001).
Seasonal differences in moss albedo (a pve,moss ) could be expected as their radiation reflection properties vary with moss water content (Graham et al., 2006).However, higher values would be expected in summer, when the moss surface is dry and lighter, but our calibration resulted in higher values during spring and winter.These values were much higher (> 22 %) compared to literature values (11-16.5 %, Berglund and Mace, 1972;16.4 %, Zhao et al., 1997;11 %, Kellner, 2001) and therefore rather compensate for values of interacting parameters (in particular z 0M,snow and r alai ) or not implemented processes.Especially the effect on winter H and Rn might result from the strong interaction with z 0M,snow , as the mosses in winter are covered with a thick snow cover, so that their albedo should not show any sensitivity in winter.Furthermore, H in spring tended to be overestimated, which would be compensated by a high albedo during this time, but might be caused in the real world by open water over frozen soil, which was not realized in the model.Interestingly, albedo of vascular plants did not show any sensitivities, neither during vegetative stage (a pve,vasc ) nor after start of senescence (a pgrain ) when a higher value would have been expected due to leaf yellowing.Direct measurements of plant albedo were not available in this study.A time series observation of those would be very helpful for clarification, as this parameter is known to vary substantially within and between peatlands (Belyea, 2009).

Snow
The model performance in simulating snow depth was not connected to performance in any other variable, except to performance in H if exclusively springtime values were considered.This was surprising, as the uncertainty for timing of snowmelt ranged for about 2 weeks but determined the start of temperature rise, water table dropping and biotic activity.A possible explanation might be the poor ability of the snow-depth R 2 and ME to assert a good fit in duration of snow cover.This is supported by the fact that the most important parameter for timing of snowmelt (m T ) strongly affected performance in dynamics of H , NEE and Ts during springtime.Parameters defining timing of snow depth might be better constrained if future calibrations include an additional variable with a stronger conclusiveness to the timing of snowmelt, e.g., by a Boolean time series indicating if snow cover is present or not.It needs to be tested if this could also help to solve the disagreements in value ranges between the performance indices in the case of the density coefficient of old snow (S dw ), which in combination with m T caused the low average overlap within snow-depth sensitive parameters.
According to Jansson and Karlberg (2010), a high value for m T (4-6 kg • C −1 m −2 day −1 ) could be expected for open fields.A possible explanation for the low accepted values (< 3 kg • C −1 m −2 day −1 ) of m T in the case of criteria on H in contrast to the high values if criteria were on Ts could be that high values compensate for overestimated springtime H (see Fig. S1 in the Supplement).However, the overestimation of spring H might be connected to different reflection properties of mosses during springtime or to missing consideration of radiation reflection and evaporation from open water, which might be formed during snowmelt on still frozen soils.The latter is further supported by the underestimation of LE during April and May (Fig. S1 in the Supplement), which cannot be connected to underestimated plant transpiration, as the model even tended to overestimate CO 2 uptake during this period.

Conclusions
CO 2 models are commonly calibrated on NEE as only measurement variable.Here, we investigated the interactions between different abiotic and biotic processes and their parameters, as well as the implications and usefulness of data on not only NEE but also LAI, sensible and latent heat fluxes, radiation, water table depth, soil temperatures and snow depth for model calibration on a boreal peatland.Processes and model performance in the different observation variables were strongly interlinked across process categories.This means parameter ranges that result from calibration depend on model structure, included processes, other parameter values and calibration set-up, and might therefore not be transferable between studies.It further implies that a study aiming to understand and interpret parameter values needs to calibrate processes and parameters of many different process categories, using a wide range and multiple criteria on various observation variables.
The key parameters identified will help to simplify future model calibrations by selecting only the most influential parameters for the variable of interest and using a narrower range for the constrained parameters.This means a simpler calibration and faster computation and, in turn, allows for the inclusion of a more detailed investigation of a process of certain interest.Furthermore, it helps model developers to include the most sensitive processes for simulating a certain variable.
Parameter interactions were found to be more important than parameter value ranges, revealing the need for accounting for equifinalities, also across different biotic and abiotic processes: either by calibrating correlated parameters dependent on each other or by calibrating only one of the correlated parameters.The latter will lead to a narrower constrained range, but this range might not be transferable to other sites and other models.
The gained knowledge on trade-offs will be useful to avoid modelling studies with too many purposes and helps model users assessing the implications of their criteria choice.The validity of calibrated models is always restricted and robust-ness of obtained parameter ranges should be questioned.The identified supporting effects between some variables indicated that some measurement variables can partly compensate absence or low resolution of the connected variable.This information tells experimentalists which measurement variables are helpful and which are obligatory if a certain process should be understood from the underlying regulating principles.It further helps modellers to decide if a site has enough available data for model calibration and to estimate uncertainties in model predictions depending on available ancillary data.
All observed calibration variables (NEE, LAI, sensible and latent heat fluxes, net radiation, soil temperatures, water table depth and snow depth) helped for model constraint and interpretation.Ancillary variables are in particular important for evaluating the robustness of calibrated parameter ranges.They should therefore be measured on sites used for calibration of complex process oriented models.Additional measurements of, in particular, soil hydraulic properties or water content would largely reduce uncertainty and help for a better parameter constraint.

Code and data availability
The model and extensive documentation can be downloaded from the CoupModel home page http://www.coupmodel.com/ (CoupModel, 2015).The source code can be requested for non-commercial purposes from Per-Erik Janson (pej@kth.se).The simulation files including the model and calibration set-up, the used parameterization and corresponding input and validation files can be requested from Christine Metzger (cmetzger@kth.se).They cannot be made publicly available, as they include climate and site data that require authorization from the data owners.
The flux data and ancillary data are available from the European Flux Database Cluster (http://www.europe-fluxdata.eu/ (Nilson et al., 2014) The Supplement related to this article is available online at doi:10.5194/gmd-9-4313-2016-supplement.

Figure 1 .
Figure 1.Energy flux partitioning and related soil water flows in the CoupModel as applied to a peatland using two plant canopies and root systems.Rn: Incoming radiation; LE: latent heat fluxes (sum of actual transpiration, interception evaporation and soil evaporation); H : sensible heat fluxes.

Figure 2 .
Figure2.Parameter concern is shown on the y axis as sum of equifinalities (hatched) and sensitivities that could not be constrained unambiguously (solid).The x axis shows the parameters that belong to the process category of the background colour.

Figure 3 .
Figure 3. Connections between processes and parameters of different process categories.The y axis shows the count of parameters from the different process categories (colours) that are sensitive to model performance in the various variables (x axis).

Figure 4 .
Figure 4. Average overlap of accepted ranges per parameter within each process and between processes, i.e. how unambiguously the parameters could be constrained.Negative values indicate the distance between accepted ranges when ranges did not overlap at all.

Figure 5 .
Figure 5. Correlations between performance indices in the prior distribution (3200 random runs): R 2 vs. R 2 (upper panel); mean error (ME) vs. ME (lower panel).Each of the dots represents a parameter set.Grey lines indicate the axes through zero.

Figure 6 .
Figure 6.Correlations between performance indices in the prior distribution (3200 random runs): R 2 (columns) vs. mean error (ME) (rows).Each of the dots represents a parameter set.Grey lines indicate the axes through zero.

Figure 7 .
Figure7.Process category belonging to parameters that correlated with parameters of a certain process category.
; site name: Degerö; site code: SE-Deg), with open data access for the years 2001-2006, and restricted data access (the principal investigator of the site has to authorize the data request) for the years 2007-2015.

Table 1 .
Measurement data used as model input.

Table 2 .
Measurement data used for model calibration.

Table 3 .
Different criteria sets for the selections of accepted runs.