Influences of observation errors in eddy flux data on inverse model parameter estimation

Eddy covariance data are increasingly used to estimate parameters of ecosystem models and for proper maximum likelihood parameter estimates the error structure in the observed data has to be fully characterized. In this study we propose a method to characterize the random error of the eddy covariance ﬂux data, and analyse error 5 distribution, standard deviation, cross-and autocorrelation of CO 2 and H 2 O ﬂux errors at four di ﬀ erent European eddy covariance ﬂux sites. Moreover, we examine how the treatment of those errors and additional systematic errors inﬂuence statistical estimates of parameters and their associated uncertainties with three models of increasing complexity – a hyperbolic light response curve, a light response curve coupled to wa-10 ter ﬂuxes and the SVAT scheme BETHY. In agreement with previous studies we ﬁnd that the error standard deviation scales with the ﬂux magnitude. The previously found strongly leptokurtic error distribution is revealed to be largely due to a superposition of almost Gaussian distributions with standard deviations varying by ﬂux magnitude. The crosscorrelations of CO 2 and H 2 O ﬂuxes were in all cases negligible ( R 2 below 0.2), 15 while the autocorrelation is usually below 0.6 at a lag of 0.5 hours and decays rapidly at larger time lags. This implies that in these cases the weighted least squares criterion yields maximum likelihood estimates. To study the inﬂuence of the observation errors on model parameter estimates we used synthetic datasets, based on observations of two di ﬀ erent sites. We ﬁrst ﬁtted the respective models to observations and then added 20 the random error estimates described above and the systematic error, respectively, to the model output. This strategy enables us to compare the estimated parameters with true parameters. We show that the correct implementation of the random error standard deviation scaling with ﬂux magnitude signiﬁcantly reduces the parameter


Introduction
The availability of carbon dioxide and water vapour flux measurements between ecosystems and the atmosphere around the world offers various opportunities to improve our knowledge about processes connected with the global carbon cycle (Friend et al., 2007;Baldocchi et al., 2001).The interplay of models and data gives us insights into the performance of models, our level of understanding the system, but also into the quality of data and the information content therein about the processes represented in the model.Classically, parameters were often derived from experiments at leaf or plant scale or from expert judgement.If nonlinear relationships are involved the parameters are scale-dependent and cannot be easily transferred to but also not observed on larger scales.An alternative option to obtain parameter estimates is the inversion of models against data.In this case a cost function describing the misfit between model output and observations is minimized by varying the parameters.The inversion of models against Eddy-Covariance (EC) data leads to parameter estimates at ecosystem scale, our scale of interest, thus EC data are increasingly used for model inversions.EC data contain information about the actual ecosystem flux, but a measured quantity is always the sum of the "true" value and errors, these errors need to be addressed in an adequate way.The measurement errors can be distinguished into random errors, fully systematic errors and selective systematic errors (Moncrieff et al., 1996).Fully systematic errors appear constantly and arise for instance from inaccurate calibration or consistently missing high or low frequency components of the cospectrum, while selective systematic errors appear only during special temporal periods, for instance at night under unfavorable micrometeorological conditions.The random error of EC data arises from the measurement instruments, the stochastic nature of turbulence and varying footprint (the area that influences the measurement, it depends primarily on atmospheric stability and surface roughness).Quantification of the random error is a prerequisite for statistical comparisons between models and data and model-data synthesis as it expresses our confidence in the data (Raupach et al., 2005).The characteristics of the errors play an important role for the parameter estimation, the error distribution, error cross-and autocorrelations or inhomogeneous variance can bias the parameter retrieval if not accounted for (Tarantola, 1987).The study of Trudinger et al. (2007) showed that how data errors and uncertainties are treated in the optimization criterion will have a significant impact on the retrieved parameters.Studies using EC data in inverse modelling often assume constant error variance (Reichstein et al., 2003;Owen et al., 2007;Wang et al., 2007), use the standard deviation of the model residuals (Sacks et al., 2006;Braswell et al., 2005) or an adhoc fraction of the observations (Knorr and Kattge, 2005).During the last few years approaches for the quantification of random errors of EC data came up, they used paired observations, first spatially separated measurements (Hollinger et al., 2004), but as there are only few appropriately distanced towers available, Hollinger and Richardson (2005) developed a methodology using daily differenced measurements with equivalent environmental conditions that allowed to characterize the univariate distribution for several sites (Hollinger and Richardson, 2005;Richardson et al., 2006).However, the auto-and crosscorrelation of the errors have so far not been systematically quantified and are assumed to be zero.Moreover, the systematic errors are still under investigation and challenging the scientific community (Wilson et al., 2002;Friend et al., 2007).Hence the aim of this study is and their uncertainties, -and to explore how selective systematic errors influence parameter estimates of models describing carbon and water exchange.
We carry out the parameter estimation experiments with synthetic data based on eddy covariance data from two European sites and with three models of different complexity, a hyperbolic light response curve, a light response curve coupled to water fluxes and the SVAT scheme BETHY, a process-based model that calculates the CO 2 , H 2 O and energy exchanges of soil, vegetation and atmosphere for the terrestrial land surface (Knorr and Heimann, 2001).

Analysis strategy
The first part of the study deals with the characterization of the random error.We estimate the random error for four different sites, Hainich, Loobos, Puechabon and Hyyti äl ä, using the gapfilling algorithm of Reichstein et al. (2005).We focus on the statistical properties important for parameter estimation, e.g.variance, distribution, autocorrelation, crosscorrelations.The kurtosis is a measure of peakedness and can be used as an indicator for the type of distribution, the excess kurtosis used here is zero for Gaussian distributions and three for double exponential distributions.To reveal the influence of errors to parameter estimates we designed 20 synthetic data sets with random errors and 20 synthetic data sets with systematic errors for each model that are based on EC data from two sites.We optimized model parameters for three models to match ten periods consisting of 14-day EC data measured at Hainich and Loobos in 2005 from May to September to get a range of reasonable parameter estimates.The estimated parameters were used to create a reference model output.Then we added a random error and systematic error respectively.The random errors were estimated from the real data in the same way as for the first part of the study, the selective systematic error is a fixed percentage of the averaged observed night time flux subtracted from the modelled night time flux.Afterwards the parameters were reestimated using different ways to account for data uncertainty and error distribution.This strategy offers the advantage that the properties of the error are known and the model error is zero, but the dataset is still realistic.Knowing the true properties of the reference data we could compare estimated parameters with true parameters and model output to a reference model output to reveal the influence of the errors.

Data
We used half hourly EC and meteorological data from the CarboeuropeIP database.
In the statistical analysis of the random error we inlcuded data from four sites: Hainich in Germany, an unmanaged deciduous broad-leaf beech forest, Loobos in the Netherlands, a planted maritime coniferous forest, Hyyti äl ä (Finland), an evergreen needleleaf forest and Puechabon (France), an evergreen broadleaf forest.For the parameter retrieval experiments we chose two sites, Hainich and Loobos.The data sets were processed using the standardized methodology described in Papale et al. (2006); Reichstein et al. (2005).CO 2 fluxes are corrected for storage, low turbulence conditions are filtered using the u * criteria and spikes (outliers) are detected.Subsequently gap filling and fluxpartitioning is applied.For the parameter estimation only filtered and corrected high quality measurements are used.

Observation errors
In this study we assume that the measurement value consists of the actual value and an additive systematic and random error where δ is a systematic error and ǫ is a random error.The commonly used ordinary least squares (OLS) optimization assumes the random error standard deviation, Introduction

Conclusions References
Tables Figures Full Screen / Esc Printer-friendly Version Interactive Discussion e.g. the data uncertainty to be constant (homoscedasticity).A constant standard deviation, can usually be provided by the manufacturer of a measurement device or it can be determined with simple tests.For flux data the standard deviation of the random error is not constant in this case tests need to be performed for varying conditions, quantifying the changes of the standard deviation.One option is to perform measurements close to each other, temporally or spatially, provided that the conditions are the same or very similar, then the actual value is equal and the variation is caused by the random error.For the flux data meteorological conditions, the state of the vegetation and if spatially seperated the footprint and topography have to be the comparable.To get an estimate of the random error we used the gapfilling algorithm of Reichstein et al. (2005).This tool computes the expected value of the flux using data measured under the same meteorological conditions in a time window of ±7 days.The small time window is necessary to ensure a similar condition of the ecosystem.The residual of the gap filling algorithm can be used as a random error estimate (Moffat et al., 2007), it is comparable to the paired observations approach used in Hollinger and Richardson (2005), as shown in Richardson et al. (2007).For the parameter estimation an error standard deviation has to be assigned to each observation.For the parameter estimation experiments we compared the different estimates for the standard deviation of the random error: 1. constant weights, 2. the standard deviations of the observations with similar meteorological conditions within a time window of ±7 days is used directly from the gapfilling algorithm (std),this is equal to the standard deviation of the residuals between observations with similar meteorological conditions and expected value, 3. the standard deviations of the residuals of the gapfilling algorithm (res) were obtained grouping the data according to the flux magnitude in 30 groups with an equal number of data points, for each group the standard deviation was computed (see Fig. 1).Afterwards the standard deviation was related to the flux mag-Introduction

Conclusions References
Tables Figures Printer-friendly Version Interactive Discussion nitude using two linear regression lines to allow for a minimum for net ecosystem exchange of carbon (NEE) and one linear regression line for the latent heat (LE).
For the third method the modeled flux (here the flux derived from the gapfilling algorithm) has to be used to derive the dependency of error standard deviation on flux magnitude, because the relationship between the residuals and measured flux is biased (Draper and Smith, 1981).Furthermore an observation that is accidentally lower is given a higher weight than an overestimated value, which will lead to an underestimation of flux magnitude by the model (Evans, 2003).

Parameter estimation
The procedure of parameter estimation can be described as varying the parameters until the best fit between model and data is found.The fit or misfit between model and data is quantified via the costfunction: x d represents the data vector, x m the model output vector, C d the error covariance matrix.The best parameter set is found at the minimum of the cost function.T denotes that the vector is transposed.For uncorrelated errors the function simplifies, as all off diagonal elements of the matrix C d are zero, to: than data with low uncertainty (low error variance).For constant error variance the Eq. ( 3) simplifies to the OLS method summing up only the squared distances.Given a double exponential distribution as proposed by Hollinger and Richardson (2005) and Richardson et al. (2006), parameter estimation should be based on the sum of absolute deviations rather than on squares.To find the minimum of the costfunction we used the Levenberg-Marquardt algorithm implemented in the data analysis package "PV-WAVE 8.5 advantage" to (Visual Numerics, 2005) for the simple models.For the complex BETHY model a Bayesian approach was used to determine the a posteriori probability density function (PDF) of parameters including prior information and the Metropolis Markov Chain Monte Carlo (MCMC) technique was used to sample the PDF of parameters, which was then characterised by mean and 95% confidence intervalls (Knorr and Kattge, 2005).The Optimisation Intercomparison of Trudinger et al. ( 2007) compared different algorithms, including the two used here, for the optimisation of a simple coupled model.The optmisation algorithms were found to be comparable with respect to the parameter retrival.

Evaluation of the parameter estimation performance
The reestimation of the parameters was evaluated through the deviation from the original parameter value, the parameter uncertainty and the root mean squared error between model output computed with the reestimated parameters and the reference model output without noise.The uncertainty of the parameters determined with the Levenberg-Marquardt algorithm, was derived by bootstrapping (n=500), which is only based on the empirical sample not on assumptions about probability theory of the normal distribution (Wilks, 1995).As a measure of uncertainty for the paramters we used the 95% confidence intervall (=1.96• standard error) of the mean of the parameter distribution.When using the Metropolis algorithm the parameter uncertainties can be directly calculated from the sampling of the MCMC approach.The uncertainty reduction when using the Metropolis algorithm was computed as 1− α is an approximation of the canopy light utilization efficiency, β is GPP (Gross primary production) at light saturation and γ is the ecosystem respiration.Instead of R g photosynthetic active radiation (PAR) or photosynthetic photon flux density (ppfd) is often used, they are closely related to R g , but not measured at all EC sites.Using R g instead of PAR or ppfd changes only the value of α, as PPFD is approximately twice the R g .

Water use efficiency model
To increase the complexity of the model we coupled NEE with the latent energy (LE) using the HLRC and connecting it to LE via the water use efficiency (WUE), which is the ratio of gross primary production and latent heat.The WUE times water vapour deficit (WUE VPD) is considered constant (Beer et al., 2007).Using VPD as additional driver NEE and LE can be connected as follows: This model is inverted against NEE and LE.To make sure, that LE and NEE errors contribute to a similar extend to the cost function when using constant weights, we scaled the sum of the added synthetic errors to the sum of synthetic errors weighted with std for NEE and LE respectively: c denotes the constant weight, the same weight is used for the inversion of the BETHY model.The model underestimates LE, but as we use the model output as reference, model performance is of minor importance.We used this model to show, that the results derived with the simple model hold for models of various complexities and to mediate between the very simple HLRC and the quite complex SVAT scheme of BETHY.

BETHY
BETHY is a process-based model of the coupled photosynthesis and energy balance system to simulate the exchange of CO 2 , water and energy between soil, plant canopy and atmosphere (Knorr and Heimann, 2001).It computes absorption of PAR in three layers, while the canopy air space is treated as a single, well mixed air mass with a single temperature.Evapotranspiration and sensible heat fluxes are calculated from the Penman-Monteith equation (Monteith, 1965).Carbon uptake is computed with the model by Farquhar et al. (1980) for C3.The stomata and canopy model of Knorr (2000) simulates canopy conductance in response to PAR, VPD and soil water availability.In the version of BETHY applied here, autotrophic respiration is calculted as a temperature modulated fraction of photosynthetic capacity while heterotrophic respiration is based on a basal respiration modulated by soil water availability and air temperature.
The inversion set up was the same as in Knorr and Kattge (2005), inverting all 21 parameters simultaneously.The prior uncertainties of the parameters were set to 20% of the prior parameter value.For the synthetic datasets the prior parameter were the parameters used to generate the data.Introduction The standard deviation of the error has been derived from the residuals of the gapfilling model, e.g. standard deviation of the residuals depending on the flux magnitude (res), and using the standard deviation of the gapfilling algorithm directly (std).Figure 1 shows the relationship between flux magnitude and error standard deviation for NEE and LE.The standard deviation is not homogeneous, e.g. the errors are heteroscedastic and increase with increasing flux magnitude.Thus the residuals have to be weighted with the reciprocal of the standard deviation of the random errors as already suggested by previous studies (Richardson et al., 2006).The magnitude of the error variance is similar for the two methods of deriving the error variance described in the previous section, see Fig. 1, for res the observations needed to be grouped to derive the standard deviation.The res standard deviation for NEE ranges from 1 to 5, for LE from 5 to 40.
With the std method the ranges are wider because the data were not grouped.For NEE the standard deviation lies between 0.5 and 9.5, for LE between 2.8 and 85.For eddy covariance data it is known, that the error variance increases with increasing flux magnitude, Richardson and Hollinger (2005) showed that the error standard deviation of NEE not only scales with flux magnitude, but that wind speed also has a fundamental effect on the uncertainty.Thus not the whole variability of the standard deviation can be reproduced when only the flux magnitude is used.Another source for the higher scatter of the std results is the uncertainty in the estimation of the standard deviations derived directly (std).Introduction

Conclusions References
Tables Figures Printer-friendly Version Interactive Discussion

Distribution
Previous studies (Richardson and Hollinger, 2005) showed, that the error distribution of NEE is rather double exponential (Laplace) than normal and this is also found for the data used here (see Fig. 2e, f).For LE the distribution is even more peaked than the double exponential distribution.The normal distribution is characterized by the mean and the standard deviation.As error standard deviation increases with increasing flux magnitude the distribution of all error estimates is a superposition of normal distributions with varying standard deviation.If we group the data according to the flux magnitude, we find Gaussian distributions for high flux magnitudes (see Fig. 2a,  b), adding more data to the distribution plot, we find a rather double exponential distribution (see Fig. 2c-f).Another possibility to show the Gaussian distribution is shown in Fig. 2g, h , we normalized the errors with the standard deviation derived with the gapfilling algorithm (std) this transforms all error distributions to a standard deviation of unity.For NEE the normalized errors are slightly closer to a normal distribution than for LE.Thus the double exponential distribution is largely due to a superposition of Gaussian distributions and the least squares criteria can be used for eddy covariance data show here.Figure 3 shows the distribution of errors from the four different sites.The normalization of NEE resulted in a rather Gaussian distribution, for LE the distribution is inbetween Gaussian and Laplace distribution and is slightly skewed.This indicates, that the distribution of the error varies from site to site or that the error estimation does not perform well for all sites.One indicator for the peakedness of the distribution is the excess kurtosis (=kurtosis-3), it is zero for a normal distribution and 3 for a double exponential distribution, a high kurtosis indicates a strong peak.Figure 4 shows the kurtosis for ten two week periods for the four sites.For NEE the normalization of the errors decreases the kurtosis and changes the distribution to a less peaked shape.
The kurtosis is in general below the kurtosis for double exponential distributions, but some outliers indicate a much stronger peak (excess kurtosis=5.5).For LE the kurtosis decreases also for HAI, LOO and PUE.For HYY the kurtosis increases after normalization.The kurtosis shows a high sensibility to outliers, but also to the rule of the detection of outliers, excluding outliers will always decrease the kurtosis.As the outliers are important for the characterisation of the distribution and the data is already prefiltered (spike detection according to Papale et al., 2006) we did not exclude them.
For errors of fluxes with high magnitude, normal distribution is still found (excess kurtosis 0.3) and seems to be valid across sites.Overall, we conclude that random error characteristics should be considered on a site-by-site basis.When comparing different studies and different sites regarding their error distribution, a careful documentation of the influence and the treatment of outliers is strongly recommended.

Autocorrelation
Figure 5 shows the autocorrelation function of the random errors for the four eddy sites.The behaviour of the function is similar for all sites, the autocorrelation decays fast, after 10 hours there is no considerable change in the correlation.Figure 6 shows boxplots for the autocorrelation for a lag of 30 min, it is usually below 0.7, with one exception for Puechabon (0.82).Hyyti äl ä shows the highest autocorrelation for LE and NEE, Loobos the lowest for NEE and Hainich the lowest for LE.Although the gapfilling algorithm provides a reasonable estimate for the random error, the autocorrelation could partly be an artefact of the algorithm, if the deviation from the statistical expection value was not caused by a random error the following and previous value would deviate in a similar way and the actual autocorrelation of the random error would be lower.To make sure that error autocorrelation does not influence the parameter estimation one could prefer to use only every second or third value for the parameter estimation.Introduction

Conclusions References
Tables Figures Printer-friendly Version

Interactive Discussion
Crosscorrelation R 2 values for the crosscorrelation between NEE and LE errors of the four sites and ten data periods for each site are summarized in Table 4.In our study the correlation between NEE and LE errors is close to zero, thus the correlation between NEE and LE errors is of minor importance and does not need to be considered in the error covariance matrix.The highest R 2 was 0.24 for one period for Puechabon, for the same period the outlier of the LE autocorrelation for a lag of 30 min was found (see Fig. 6).
The measurements of NEE and LE are both based on the vertical wind velocity and errors introduced via the wind velocity measurement, such as errors due to turbulence sampling must show up as a correlation between NEE and LE errors.This indicates that the variation in the measured fluxes under similar meteorological conditions seems (i.e. the flux errors) to be rather caused by changes in concentrations of water and CO 2 than by the measurement of the vertical wind velocity.As auto-and crosscorrelation are low, the generalized least squares method (Eq.2) can be simplified to the weighted least squares method (Eq. 3) by setting off-diagonal elements in the error covariance matrix to zero.

Ordinary least squares vs weighted least squares
The parameters were estimated for three models of different complexities, the synthetic data is based on data from two different sites (Loobos and Hainich).We are comparing constant weights with two ways of estimating the standard deviation of the observation errors, which is then used to weight the data for the parameter estimation to account for the non constant error standard deviation.The standard deviation of the errors was estimated as the standard deviation of the observations measured under similar metereorological conditions (std) and as the standard deviation of the gapfilling algorithm residuals related to the modelled flux magnitude (res), see Fig. 1.

Conclusions References
Tables Figures Printer-friendly Version

Interactive Discussion
The results of the parameter retrieval experiments (Fig. 7) show, that the random error introduces no systematic error to the parameter estimates and the true parameters are usually within the parameter uncertainty (95% confidence interval) derived from bootstrapping.The mean of the parameter ratios is not significantly different from unity (α=0.05).The mean uncertainty of the parameters using a non constant estimate for the error standard deviation as weight is between 10 and 24% lower for the HLRC than using constant weights (see Table 2).Due to the stochastic nature of the procedure, these results are true for the mean results but there exist data periods for both sites in which the results are opposite.The std weights decrease the mean uncertainty more than res and therefore describes the error standard deviation better.The root mean squared error (rmse) between reference model output without noise and the model output with the reestimated parameters can be decreased using std as weights for the HLRC, for res it increases, again indicating, that a description of the data uncertainty only based on flux magnitude is not sufficient.
For the water use efficiency model the results of the model parameterization are similar, estimates of parameter uncertainty decrease between 5% and 60% and the RMSE between reference model output and model output of the reestimated parameters is decreased when using std, while res increases the value (see Table 3).For simplicity we focus on the comparison between constant weights and std, since std gave the best results.For the inversion of the BETHY model the distance between retrieved and true parameters can be decreased using std compared to the constant weights (see Table 4).The influence of using varying data uncertainty compared to constant data uncertainty with the MCMC algorithm is less pronounced as the absolute value of the data uncertainty is more important than the relative changes.Nevertheless the reduction of uncertainty for parameters is higher when using std and the rmse between reference and model output is decreased.Another advantage of weighting the data showed up during the initial fit to real data for the creation of the reference model output, the parameters estimated with std weights resulted in reasonable parameters, whereas using constant weights for some periods negative values for α were estimated (not shown).The random error complicates the parameter retrieval, it can increase the number of local minima or the minimization can become an ill-posed problem.Using weights representing the data uncertainty seems to improve the behaviour of the cost function and improves the extraction of information inherent to the data.This shows that the standard deviation provided by the gapfilling algorithm is a good measure for the eddy covariance data uncertainty, it improves the parameter retrieval and therefore model performance after optimization, at least for the sites used here.For skewed error distributions we would expect the parameter estimates to be biased.

Least squares vs absolute deviations
As the use of absolute deviations in the cost function was suggested previously by Richardson et al. (2006) we compare least squares and absolute deviations, to illustrate the effect to parameter estimation.Comparing the parameter ratio again shows no significant difference between the methods.For our sites, the parameter uncertainty increases using absolute deviations compared to the ordinary least squares method (see Table 2).The rmse increases compared to the OLS using constant weights and for the

Systematic error
Figure 8 shows the results of the parameter retrieval based on data with a selective systematic nighttime error of 10, 20 and 40%.The α and γ parameters of the HLRC show a systematic bias, estimated parameters underestimate the underlying true parameter.The bias is stronger for higher data error.For β the parameter bias seems to be not systematic, the retrieved parameter is for some periods lower, for some periods higher than the original parameter.β is GPP at light saturation, as NEE at light saturation does not change but only the night time flux, representing the respiration, is lower β should also be lower to sum up to the same NEE values during daytime.The effect on γ seems to be too low to show up in the comparably high values of β.For the water use efficiency model all parameters are biased, all estimates are lower than the true values.Through the interconnection of GPP and LE and the use of water and CO 2 fluxes to constrain the parameters the distance to the true value decreases for all parameters.This indicates the potential of using multiple constraints for inverse model parameter estimation.The parameter uncertainties increase the higher the error but the real value of the parameter is not within the uncertainty range of the estimated parameter.This means, that the real uncertainty of the parameter is underestimated, projection of the parameter uncertainty to model output will result in uncertainties for the fluxes that are too low.To get the real uncertainty for parameters and fluxes further knowledge about the systematic errors is needed and methods need to be developed to incorporate them into the estimation of uncertainty, if the systematic errors cannot be removed.

Conclusions
Previous work to quantify the random error structure of eddy covariance data (Hollinger et al., 2004;Richardson and Hollinger, 2005;Richardson et al., 2006)  Interactive Discussion magnitude of the error (i.e. its standard deviation) to the flux magnitude, and evaluating whether or not flux errors are Gaussian.Here we have built on these efforts by considering the auto-and crosscorrelation, introducing a new method to quantify the standard deviation of the random errors.We show the effect of the varying standard deviation to the distribution and investigate how random and systematic errors affect parameter estimates.The analysis of the error distribution shows that the apparently double exponential distribution can be almost entirely due to the superposition of Gaussian distributions with inhomogeneous variance.Whether this is the case for a special site can be affirmed by testing the normality of the normalized error distribution.If it cannot be affirmed one should consider using robust methods.The autocorrelation is low, but one might consider to analyse the autocorrelation function and use only every second or third flux for parameter estimation if there is enough data available.As a reason for the low but significant autocorrelation of errors we can not exclude artefacts of the gap filling tool.The crosscorrelation between LE and NEE is low and can be neglected.The assumption for ordinary least squares that is not met is the constant error standard deviation, thus the ordinary least squares method needs to be extended to weighted least squares, using the reciprocal of the standard deviation as weight in the costfunction.We propose a measure for data uncertainty, e.g. the standard deviation of the values used to compute the expected value, that can be used to weight the data in the costfunction.Weighting the data decreases the parameter uncertainty and the parameter retrieval is improved.We showed that this result holds true for a wide range of model complexities.We show that the impact of systematic errors varies by parameter, but the bias is systematic, therefore the interpretation of parameters derived from data with systematic errors might be misleading.The parameter uncertainty slightly increases when a systematic error is added, but the true parameter is not within the uncertainty range of the estimate.Not considered here but of similar importance is the model error, which was set to zero by using the model output as basis for the synthetic data.For the least squares optimization the model output random error is additive to the data ran-  -16.period 15.5. 31.5. 15.6. 30.6. 15.7. 31.7. 15.8. 31.8. 15.9. 30.9 to fully analyze the random error of EC water and carbon fluxes regarding the properties important for parameter estimation, i.e. beside the univariate distribution, also autocorrelation and multivariate correlations of CO 2 and H 2 O fluxes, -to elucidate the effect of the error model choice on model parameter estimates σ d is the standard deviation of the random errors, N the number of data points.In this study we use synthetic data based on a model output, therefore the model error is zero.To consider the uncertainty of flux measurements is necessary if the errors are hetereoscedastic, e.g.error variance (=squared standard deviation) increases with increasing flux magnitude, or if different data sources are used.From another point of view this means that data with high uncertainty (high error variance) get a lower weight weighted least squares.Since by normalizing the errors with the standard deviation we get a Gaussian distribution for our selected sites the absolute deviation minimization cannot improve the parameter retrieval.If the errors show a double exponential distribution as a result of the superposition of different Gaussian distributions, then least squares optimization should be applied.If the error distribution is more peaked due to outliers or a different data filtering, robust methods like the minimization of absolute deviations or robust regression tecniques, which exclude outliers, may be advantageous.Testing whether the normalized error distribution is Gaussian could support the choice of the costfunction.
has focused on describing the moments of the distribution of the error, particularly relating the expected

Table 1 .
depending on the point of view part of the data random error can also be seen as model errors, e.g.footprint heterogeneity.Even more severe problems can arise from model structural problems, which are comparable to systematic data errors, the model cannot reproduce the dynamics in the data because processes are missing or represented in an insufficient way.Using inverse model parameter estimation we can only accept the best fit although the model is not meant to reproduce the patterns emerging from these insufficiencies.Hence we conclude that potential systematic errors in flux data or models need to be addressed more thoroughly in data assimilation approaches since otherwise uncertainties will be vastly underestimated.Introduction Crosscorrelation between NEE and LE errors for ten two week periods between March and September 2005.

Table 3 .
Mean of retrieved normalized parameters and the mean uncertainty for the ten two week periods for the WUE-model using last squares minimization.