Towards operational remote sensing of forest carbon balance across Northern Europe

Abstract. Monthly averages of ecosystem respiration (ER), gross primary production (GPP) and net ecosystem exchange (NEE) over Scandinavian forest sites were estimated using regression models driven by air temperature (AT), absorbed photosynthetically active radiation (APAR) and vegetation indices. The models were constructed and evaluated using satellite data from Terra/MODIS and measured data collected at seven flux tower sites in northern Europe. Data used for model construction was excluded from the evaluation. Relationships between ground measured variables and the independent variables were investigated. It was found that the enhanced vegetation index (EVI) at 250 m resolution was highly noisy for the coniferous sites, and hence, 1 km EVI was used for the analysis. Linear relationships between EVI and the biophysical variables were found: correlation coefficients between EVI and GPP, NEE, and AT ranged from 0.90 to 0.79 for the deciduous data, and from 0.85 to 0.67 for the coniferous data. Due to saturation, there were no linear relationships between normalized difference vegetation index (NDVI) and the ground measured parameters found at any site. APAR correlated better with the parameters in question than the vegetation indices. Modeled GPP and ER were in good agreement with measured values, with more than 90% of the variation in measured GPP and ER being explained by the coniferous models. The site-specific respiration rate at 10°C (R10) was needed for describing the ER variation between sites. Even though monthly NEE was modeled with less accuracy than GPP, 61% and 75% (dec. and con., respectively) of the variation in the measured time series was explained by the model. These results are important for moving towards operational remote sensing of forest carbon balance across Northern Europe.


Introduction
The global carbon balance is the result of fluxes of carbon into and out of ocean and land. In this balance, terrestrial carbon sinks result in high quantities of carbon being drawn from the atmosphere into land. The compensating effect of the carbon sinks in relation to the anthropogenic fossil fuel emissions -now being the main source of 5 atmospheric CO 2 (Keeling et al., 1996;Schulze et al., 2000) -makes the knowledge of sink distribution of utmost importance for understanding the biosphere's interaction with climate and its impact on future carbon levels. There is scientific evidence of a missing carbon sink in the global carbon cycle, probably located in the Northern Hemisphere (Tans et al., 1990;Keeling et al., 1996). With boreal forests covering large parts of 10 the Northern Hemisphere, knowledge of the carbon sink/strength of these regions is if great importance.
It is possible to measure the net exchange of CO 2 between atmosphere and biosphere for long periods of time using eddy covariance methods. These measurements have proven to be of great importance for studies concerning carbon budgets and their 15 seasonal patterns (e.g. Lindroth et al., 1998;Falge et al., 2002). However, measurements can be made only at a limited number of locations due to the high costs of implementing the technique, for example, less than ten eddy covariance towers are currently operating in the forested areas of Sweden which is not enough to represent the different forest age classes, species classes, etc., needed for making a precise 20 national carbon balance estimate (Lagergren et al., 2006).
Two main approaches can be identified for obtaining larger scale estimates of carbon balance: i) employment of physiological process-based models, and ii) direct estimates on the basis of remote sensing. The former approach simulates ecosystem processes using detailed data sets of biophysical and meteorological conditions as input (e.g. The method has been proven attractive to implement on the basis of remote sensing since it is possible to obtain these parameters from satellite (e.g. Sellers et al., 1994;Nichol et al., 2000;Lobell et al., 2002).
While the gross primary productivity (GPP) gives the total amount carbon fixated by photosynthesis, NPP -which is obtained by subtraction of the autotrophic respiration 10 from GPP -gives the input of carbon to the ecosystem. Depending on the definition of the LUE factor, both of these productivity measures can be successfully modeled on the basis of light-use efficiency and remote sensing. For example, Xiao et al. (2004a,b) modeled the LUE factor as a function of temperature, water and leaf phenology in a light-use efficiency-based GPP model. Running et al. (1999Running et al. ( , 2000 used a process-15 based ecosystem model for determining LUE factors for daily GPP calculations. A number of NPP models based on the light-use efficiency concept have been published (e.g. Ruimy et al., 1994;Lagergren et al., 2005). However, since the computation of the LUE-factor usually requires meteorological data, an operational LUE-based model may easily be limited by the spatial resolution and accuracy of the meteorological input 20 data sets (Sims et al., 2006).
The net ecosystem exchange (NEE) is obtained by subtracting ER from GPP. Since NEE gives the net amount of carbon uptake or release, this is a more appropriate measure than NPP and GPP of carbon sink strength (Schulze et al., 2002). However, using the LUE approach, or remote sensing techniques in general, for NEE modeling is 25 problematic due to the difficulties involved in obtaining information on the heterotrophic respiration from space (Valentini et al., 2000). The mechanisms behind decomposition of soil organic matter are not yet fully understood as decomposition depends not EGU only temperature but on a range of different environmental constraints (Davidson and Janssens, 2006). Accordingly, estimating soil decomposition is hard in general, and especially so using only data from satellite. Given the importance of ER in controlling the forest carbon balance, especially boreal forests at high latitudes (Valentini et al., 2000), there are strong arguments for including information on ER when studying car-5 bon balance from space, despite the difficulties involved.
Most studies aiming at estimating NEE using remote sensing employ physiological process-based models which are parameterized for a specific location using flux data together with meteorological data; the model output is then scaled up using satellite data. Veroustraete et al. (1996) estimated NEE for a deciduous forest in Belgium us-10 ing a physiological process model driven by the fractional absorption of PAR (FAPAR) from satellite and meteorological data. Chiesi et al. (2005) estimated monthly NEE for a forested site in Italy by parameterizing and calibrating FOREST-BGC (Running and Coughlan, 1998), including FAPAR derived from Landsat normalized difference vegetation index (NDVI). With correlation coefficients above 0.9, NEE was obtained with 15 high accuracy when comparing to ground-measured NEE. Hunt et al. (2004) used a LUE-based approach instead of a process model for obtaining NEE for two rangeland sites in the U.S. The model was driven by meteorological data together with satellitederived FAPAR. The authors also investigated the relationship between APAR and NEE for the rangeland sites and a coniferous forest site; a linear relationship was found for 20 the rangeland sites but not for the coniferous site.
NEE can be obtained with high accuracy for a parameterized site using a physiological process-based model, however, if the model requires a large set of input parameters (of, for example, meteorological data) which may be hard to obtain for large areas, the spatial representation of the study may be limited. 25 The RACES (Regional Arctic CO 2 Exchange Simulator) equations allow for calculation of ER and GPP using only three variables: NDVI and PAR for calculation of GPP and temperature and NDVI for calculation of respiration (Vourlitis et al., 2003). This limited input demand make the equations attractive from a remote sensing perspective EGU and, hence, allow for larger spatial representation. The model has been used for NEE studies in arctic ecosystems (Oechel et al., 2000;Vourlitis et al., 2000Vourlitis et al., , 2003. Turner et al. (2004) used remotely sensed data on land cover, stand age and harvesting in combination with BIOME-BGC (Law et al., 2001), coupled with a regional climate dataset, for monitoring of the carbon sequestration for a region in the state of 5 Oregon, USA.
A different approach was presented by Churkina et al. (2005) who used flux measurements from 28 sites in North America, Europe and Brazil for investigating the relationship between NEE and the carbon uptake period. The authors found a linear relationship between the two parameters and with the possibility of obtaining the car-10 bon uptake period from satellite vegetation indices, the approach has the potential of extrapolating NEE over large areas without relying on heavy input datasets. Rahman et al. (2005) investigated the relationship between MODIS enhanced vegetation index (EVI) and GPP; and MODIS surface temperature and respiration using data from ten flux tower sites across the U.S. The authors identified the scientific need 15 for an operational per-pixel production model and suggested that NEE can be obtained solely from satellite parameters. Sims et al. (2006) also identified the need for a fully satellite driven productivity model not relying on the traditional LUE approach, and investigated the use of MODIS EVI for GPP estimations using data from nine flux tower sites across the U.S. The authors 20 showed that a GPP model driven solely by EVI performs as good as or even better than the MODIS GPP product. The authors state that they are exploring the use of satellite derived air temperature for a robust carbon balance model entirely based on satellite data. The MODIS GPP product, which gives eight day averages of GPP globally at 1 km resolution, was in turn validated by Turner et al. (2003)  EGU in the possibilities of operational implementation (Running et al., 2000), and that it can output productivity with high accuracy if parameterized correct (Turner et al., 2005). However, a LUE-model usually relies on meteorological data sets for obtaining the LUE-factor. With some meteorological variables being only available at about 1 • ×1 • , a meteorological data set may be too coarse to take the heterogeneity in LUE into 5 account, and accordingly, may result in erroneous estimates of carbon uptake (Turner et al., 2002;Sims et al., 2006). Alternatively, parameterization of process-based model can give the full carbon balance (NEE) but this will require a detailed input data set that is rarely available on a regional scale, making it hard to implement operationally. For these reasons, the aim of this paper is to investigate the possibilities of imple-10 menting a model for obtaining NEE, GPP and respiration over large areas. The aim is to no rely on large input data sets but to drive the model using only a few parameters that can be obtained from satellite.  et al., 2003). A PAR-to-shortwave conversion factor of 0.43 was used (Olofsson et al., 2007).

5
AT is defined as the averaged daily temperature when the global radiation exceeds 1 W m −2 .
The reason for including Hainich is that data from only one deciduous site (Sorø) is present in the material, and by including Hainich, which is rather similar to the beech forest in Southern Scandinavia, the amount of deciduous data is doubled. 10 2.1.1 GPP and ER from measured NEE GPP and ER were derived from the NEE measurements. Respiration at Norunda was derived from a two month relationship of average night NEE and night temperature; while respiration at Skyttorp was derived from a half month relationship between air temperature and night NEE. At the other sites, respiration was obtained through a 15 short-term temperature response of night-time fluxes based on NEE. This processing had been performed beforehand since all data were obtained through the CARBOEU-ROPE project. A review of the different methods for separating ER and GPP from NEE measurements is given by Reichstein et al. (2005). Correct ground measurements are crucial when comparing against values derived from satellite measurements for 20 evaluation of model accuracy and robustness. Although no error analysis of the flux measurements was performed, averaging the data over long periods (from half-hourly to monthly) reduced the random sampling errors to relatively small values (Baldocchi, 2003

NDVI
The NDVI is the normalized quotient of the near infrared and red surface reflectance (ρ NIR and ρ R , respectively): The NDVI data set for Scandinavia is based on seasonally adjusted NDVI at 250 m 5 resolution from Terra/MODIS. Creation and validation of the data set is described and discussed in Olofsson and Eklundh (2007). MODIS NDVI for Hyytiälä and Hainich were obtained from the MODIS ASCII Subset project at 1 km resolution. The data were processed in the same manner as the data set in Olofsson and Eklundh (2007), which includes season adjustments by nonlinear least square fitting of local double logistic 10 model functions to the time-series, using the computer program TIMESAT. TIMESAT fits the function to the upper envelope of the VI data, thus, effectively reducing negatively biased noise due to remaining atmospheric influence. The least-squares procedure also, to some extent, eliminates single outliers such as the high values often occurring during early and late parts of the season, when Solar zenith angles are very 15 high (Jönsson andEklundh, 2002, 2004). This processing is done in order to filter out the noise in the VI time series, clearly visible in Figs. 1 and 2. Since Hyytiälä and Hainich were not included in the data set created by Olofsson and Eklundh (2007), 1 km NDVI time series at these two sites were obtained from The Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC).

EVI
The EVI was developed to optimize the vegetation signal but maintaining sensitivity in high biomass areas, while minimizing background and atmospheric influences (Huete et al., 2002). The EVI is defined as

EGU
where G is a gain factor; C 1 and C 2 are coefficients for correction of atmospheric influences utilizing the red and blue bands; and L is a background adjustment factor that compensates for the higher near infrared reflectance due to the red extinction through the canopy (Liu and Huete, 1995;Huete et al., 1997Huete et al., , 2002. EVI at 250 m and 1 km were acquired from Terra/MODIS and processed in the same 5 manner as the NDVI data (Olofsson and Eklundh, 2007), which includes the TIMESAT preocessing.
The EVI data cover Scandinavia and Finland (tiles H18V02, H18V03 and H19V02). Only 250 m data was used in the analysis for Sorø since the MODIS 1 km pixel represents other land covers than deciduous forest at this site (Olofsson and Eklundh, 10 2007).

VI issues
Time series of EVI, NDVI and FAPAR for a specific flux site were constructed using values from the single pixel centered on the flux site. Figure 1 shows raw and TIMESATadjusted EVI at both 250 m and 1 km resolution for the Knottåsen site (coniferous). It 15 is apparent that the 250 m EVI is unable to trace the seasonal development of the vegetation, mainly because of unrealistically high values in between the growing seasons. TIMESAT eliminates the most severe outliers but is unable to extract the seasonal dynamics. The 1 km data on the other hand displays the expected seasonal trend clearly. Although not plotted, there is a large difference between the 250 m and 1 km data also 20 at the other coniferous sites. The situation is strikingly different for the Hainich site (deciduous) as seen in Fig. 2. The NDVI is not exhibiting the same difference between resolutions, and as stated, there is no discrepancy in EVI for the deciduous pixels. Why there is such a difference in EVI, but not for NDVI, between the two resolutions for the coniferous sites is not fully clear. It is hypothesized that three different factors inter- EGU communication). A study is currently being planned aiming at investigating the reasons for the observed differences and evaluating the performance of the MODIS NDVI and EVI at 1 km, 500 and 250 m across Scandinavia.

FAPAR
An FAPAR data set based on MODIS 250 m NDVI across Scandinavia was created by 5 Olofsson and Eklundh (2007). The data set was evaluated at Sorø, Asa, Norunda, and Skyttorp, and produced an RMSE of the means at these sites between 2.5 to 6.6% at a daily basis. The data were monthly averaged for use in this study which is likely to have decreased the errors further. Since the data set was created for Scandinavia, FAPAR data for Hyytiälä and Hainich was lacking. To obtain these data, the NDVI-10 FAPAR relationship observed in Olofsson and Eklundh (2007) was applied to the 1 km NDVI at these sites. However, the NDVI-FAPAR relationship established is valid for the 250 m NDVI while the NDVI for these sites has a spatial resolution of 1 km. The 1 km and 250 m NDVI data, just like the EVI, do not appear to behave in the same manner, and a comparison between the NDVI time series at Norunda and Hyytiälä reveals that 15 the winter values for all years are lower at the latter site. Bearing this in mind, applying the NDVI-FAPAR relationship observed in Olofsson and Eklundh (2007) to the 1 km NDVI at Hyytiälä will generate a somewhat biased FAPAR. Exactly how large this bias is and how large the difference between 1 and 0.25 km NDVI is, will be investigated in a future study. 20

Models
Obtaining estimates over large areas requires NEE to be derived from parameters obtained at high areal and temporal resolution as the involved processes' responses to the environment are non-linear. Since NEE (F n ) is the difference between GPP (P g ) and ER (R e ), EGU both terms need to be modeled. The vegetation indices (VIs) are obtained from satellite; APAR is estimated with FA-PAR from satellite (see below) and measured PAR. AT is measured -these parameters are as input for predicting the biophysical parameters: NEE, GPP, and respiration. To solve this task, linear regression models were constructed for monthly values of GPP 5 and ER, for the two different forest types.
For GPP, a simple linear model with EVI and APAR as independent variables gave the best fit to measured values: where E (t) is EVI at month t, I p,a (t) the amount of absorbed PAR at month t, and 10 {ǫ(t)} is independent Gaussian white noise. APAR is calculated as the product of measured incident PAR and the FAPAR obtained from Terra/MODIS NDVI in Olofsson and Eklundh (2007) Since soil respiration is a highly complex biological process involving many variables which are not possible to derive using the limited number of satellite parameters we 15 have at our disposal (Davidson and Janssens, 2006), an external variable is most likely needed. It has been shown that the Lloyd-Taylor equation for calculation of soil respiration can be successfully used with the original parameterization for obtaining ecosystem respiration for many of the sites in this study (Anders Lindroth, unpublished results). The expression is given by where R 10 is the site-specific respiration rate at T =10 • ; and a 1−3 are model coefficients (Lloyd and Taylor, 1994). The original coefficients were derived from a large data set on soil respiration representing ecosystems all over the world, however, no boreal forests were included and a new set of coefficients based on the data in the present study was EGU Different sets of model coefficients were calculated for deciduous and coniferous stands, respectively, and in order not to use the same data for model parameterization and evaluation, about half of the data was used for estimation and the other half for validation. Three years from Hainich and three years from Sorø were used for estimating the deciduous model and the rest of the data (also 6 years) was used for evaluating the 5 performance of the model. All the data from Norunda and Hyytiälä except the last year at each site was used for estimating the model for the coniferous stands (8 years of data) while the rest of the data was used for evaluation (6 years of data, all stands represented). Since the data used for evaluation was completely excluded from the data used for estimating the model, the model performance gives the accuracy by which 10 these biophysical parameters can be obtained using this technique.
Given that the validation set differs from the calibration set, both in time and space, the validation tests for model robustness both in location and time. The model accuracy is reported by the root mean square error (RMSE) and the coefficient of determination r 2 (the squared correlation coefficient of the predicted and measured values, y andŷ), 15 respectively) which gives the amount of response variation explained by the model.

Operational remote sensing of APAR and AT
The models were estimated using measured values on AT and PAR, with APAR being the product of measured PAR and FAPAR estimated from satellite. In order for 20 the methodology to be used on an operational basis, these parameters have to be available from satellite at the same resolution as the VIs and FAPAR. Van Laake and Sanchez-Azofeifa (2004, 2005 presented a method which gives incident PAR by implementation of a simple radiative transfer model where the amount of PAR is a function of the solar constant, solar zenith angle and the atmospheric transmittance, which in turn 25 is calculated using daily atmospheric data from the MODIS sensor onboard the NASA ing incident PAR using reflectance data from MODIS. Surface and top-of-atmosphere reflectance were used for deriving both diffuse and direct instantaneous PAR. Daily values were obtained by regression analysis. With the use of look-up tables and by assuming known aerosol properties the method does not rely on atmospheric data which has proven to be hard retrieve from satellite. The results published so far are promising, and relative errors as low as 2% for daily values have been reported (Liang et al., 2006). Accordingly, incident PAR can be obtained at high accuracy from e.g. the MODIS sensor. Temperature can also be obtained from MODIS: the land surface temperature product gives temperature at 1 km resolution globally every day by the use of 15 seven thermal infrared bands (Wan and Li, 1997). The product have been validated against ground measured temperature and high accuracies have been reported, however, since thermal infrared signals from the earth surface are hard to register by a satellite sensor when clouds are present, temperature is only retrieved when the sky is clear (Wan et al., 2004). Furthermore, the MODIS sensor is carried by both the 20 Terra and Aqua platforms, resulting in two daily observations around noon instead of daily means. Assuming that AT, as it is defined in this study, can be obtained from satellite is therefore not a fully valid assumption. As an alternative to satellite observations, regional climate data sets is a solution. For example, the SWECLIM data set gives modeled air temperature with a three hour temporal resolution across the Nordic 25 countries (Rummukainen et al., 2000). EGU 3.2 Correlation between satellite and biophysical parameters Figure 3 shows the relationships between 1 km EVI and GPP (a) and daily maximum NEE (b) averaged over one month and for all coniferous and deciduous data, except for Sorø where 250 m data was used since the 1 km pixel is clearly influenced by thez surrounding crop lands. The linear correlation is stronger for the deciduous 5 sites, r(E, P g )=0.90; r(E, F n,max )=0.91, than for the coniferous sites (0.83 and 0.85, respectively). One of the reasons for this is that the deciduous data are collected from only two sites (12 years of data), which are quite similar in terms of species and structure, while the coniferous data are collected at five different sites, including both pure and mixed stands of Scots pine and Norway spruce, for a total of 15 years. Another 10 reason is that the EVI has a larger dynamic range for the deciduous sites than for the coniferous, with EVI values ranging between about 0.15 and 0.8 as compared to 0−0.5 for the coniferous (Fig. 4, also evident in Figs. 1 and 2), which in combination with higher uptake levels results in a stronger correlation. It was found that the productivity of the coniferous ecosystems is one month ahead the EVI (Fig. 5), or in other words, the EVI reaches the corresponding GPP level one month "too late", and the reported relationships were established using the EVI data from the one month following of the biophysical parameters. (This is clearly illustrated in Fig. 5 which shows EVI and GPP for the Knottåsen site.) The reason for this is not fully clarified. One issue that may help explain the phenomenon is the use of TIMESAT; the original EVI sometimes exhibits 20 a peak in the beginning of summer which tends to be eliminated by the adjustment. This is obvious in Fig. 1  EGU main advantages of the EVI over the NDVI is the higher sensitivity over regions with high biomass, making the the former less prone to saturation, a fact which is confirmed by the results of this study. As stated, the NDVI values are higher than the EVI values for almost all data points, with an average NDVI of 0.65 and 0.75 for the coniferous and deciduous sites, respectively, as compared to 0.27 and 0.41 for the EVI. The same 5 has been found in ecosystems considerably different from the ones in this study (Huete et al., 2002;Fensholt et al., 2006). However, NDVI correlates fairly well with the maximum uptake (r=0.86 and 0.73 for deciduous and coniferous data, respectively), but due to the exponential behavior induced by the NDVI saturation, the confidence and prediction intervals are wider compared to the EVI relationships, especially for the de-10 ciduous data. The correlation coefficients are given in Table 3, which also gives the correlation between the VIs and monthly values of NEE and respiration; again, EVI is superior to the NDVI and correlates better with all parameters for both forest types, and again, the saturation of the NDVI at higher levels of uptake (or release for respiration) is responsible for the lack of correlation. The saturation generates a characteristic 15 "wall" of NDVI values around 0.9, parallel to the y axis (Figs. 6a, b and 8), which is not present for the EVI (Figs. 3a, b and 7). There is low correlation between monthly NEE and the VIs, considerably lower than for GPP and ER, mainly because of a slight time discrepancy in the seasonal development between productivity and ER. This occurs because the altter depends on both air and soil temperature, and soil temperature lags 20 air temperature, which causes a decrease in the net carbon uptake sometime around late summer. With NDVI being less able to capture the vegetation peak, EVI is better correlated with these parameters. An example of this is given in Fig. 9  EGU and 0.89 at Norunda (Fig. 10a). The relationship is weaker for the deciduous stands (Fig. 10b), r[T, log(R)]=0.85 (0.91 for Hainich and 0.87 for Sorø). As seen by the slope of the linear fit, ER also seems to be responding to temperature quite differently, with the response being much sharper at Sorø indicating that an increase in temperature results in sharper rise in ER at Sorø than in Hainich.

5
3.3 Models of ER, GPP and NEE Different sets of models were constructed for the two forest types ( Table 4). The models were constructed using data from certain sites and evaluated in others. The strong correlations reported here (Table 3) indicate that the biophysical parameters, GPP and max. NEE above all, can be obtained using a simple linear regression model with only one independent variable. For example, when estimating P g (t) using E (t+1) as independent variable, the regression model explains 79% of the variation for the coniferous data. Using I p,a (t) instead gives an r 2 value of 90%. The corresponding figures using E (t) and I p,a (t) for the deciduous data are 76% and 78%, respectively. These results support the findings of Sims et al. (2006) who stated that EVI alone can 15 be used successfully for estimating GPP. Using APAR instead of EVI will give even higher accuracies, assuming that incident PAR can be obtained with high accuracy. Coniferous max. NEE was estimated using E (t+1) as independent variable with an r 2 value of 80%. Using APAR instead gave a lower accuracy which is a bit surprising since the full APAR time series correlated better with max. NEE than the EVI. Also de-20 ciduous max. NEE was more accurately modeled using E (t) as independent variable, with 80% of the variation in measured max. NEE explained. Accordingly, the maximum carbon uptake, calculated as the mean of the five observation around the highest daily recorded value, can be modeled on a monthly basis using EVI alone with an r 2 value of 80% for both forest types. 25 As expected, monthly coniferous NEE can not be modeled using neither APAR, EVI nor AT alone, nor combined in a multiple regression model. These models explained only about 50% or less of the measured variation in NEE. EGU study was more regular and displayed a higher degree of periodicity, resulting in higher correlations with e.g. APAR, thereby making it possible to estimate it directly from APAR alone. Even though such a model explains about 75% of the variation, it is likely that the inclusion of other stands would lower that figure, making it necessary to include e.g. temperature. However, more data is needed before such conclusions can be drawn.

5
To obtain NEE, individual multiple regression models for GPP and ER were created. It was found that using EVI and APAR as independent variables for GPP gave the highest accuracy, with APAR being a more important variable explaining more of the variation than EVI. No accuracy was gained by using E (t+1) instead of E (t), probably because APAR is the more important variable. As suggested by the observed correla-10 tions in Table 3, using EVI instead of NDVI gave a higher correlation, probably because of the larger dynamic range of the EVI. Monthly GPP was well described by the regression models, with a large proportion of the total variability in measured GPP being explained by the model: r 2 =86% (dec.) and r 2 =95% (con.), respectively (Figs. 11 and 12; and Tables 4 and 6). The lower accuracy for the deciduous stands is to a large 15 extent caused by an underestimation of the productivity peaks during the summers in Sorø. Since the coniferous GPP model is constructed using data from only two sites and evaluated at five sites representing four years, the r 2 value of 95% indicates that the model is robust both in time and space. With only two stands, the same conclusion cannot be drawn for the deciduous model. However, a high accuracy is achieved and 20 even though more validation is needed, the results indicate that deciduous GPP can be accurately modeled on an operational basis using only EVI and APAR as independent variables. It is evident from Fig. 11  EGU ER was modeled using Eq. (4), and just as with GPP, a high accuracy was achieved, especially for the coniferous data (r 2 =73% and r 2 =93%; Figs. 11 and 12; and Table 4 and 6). As evident in Fig. 12, modeled coniferous ER traces measured ER very well. The accuracy of the deciduous ER model is not able to capture the much higher respiration rates at Sorø compared to Hainich, especially during 2005. There is also a 5 weaker relationship between AT and ER for the deciduous data (Table 3), with the two stands exhibiting quite different responses in ER to changes in AT (Fig. 10b). Since AT is strongly correlated to logarithmic ER, an attempt was made to model ER as an exponential function of AT without relying on R 10 . Even though such a model will be able to explain much of the variation in ER at a specific site, the model -which was con-10 structed in the same manner as the GPP model -was unable to capture the respiration magnitudes for the different sites, confirming the findings of Janssens et al. (2001) who found AT to control changes in soil respiration within sites but not between sites. For example, the model was able to simulate ER at Hyytiälä 2005 and Asa which exhibits quite similar AT response (Fig. 10a), but failed to capture the high respiration rates 15 during summer at Knottåsen. As a consequence, an external factor which defines the site-specific magnitude in ER has to be introduced. Equation (4) uses the year-and site-specific ER rate at 10 • , R 10 , which is derived from the observed logarithmic ER and AT relationship, and hence, makes the model non-operational as there is no established approach to obtain this parameter from space. The observed R 10 are given 20 in Table 5. An attempt was made to relate these values to the VI's but no significant relationships were found. NEE is the difference between productivity and respiration, and modeled NEE was here defined as the difference between modeled GPP and ER. Much of the carbon fixed by photosynthesis in Northern European forests is lost through respiratory processes 25 (Valentini et al., 2000), which -from a computational aspect -means that the remaining NEE becomes small and therefore is more likely to generate a higher relative error in the modeled time series, and hence, weaker correlation. This is especially true for the boreal sites which have an annual NEE close to zero. At Norunda, for example, EGU an annual carbon loss has been observed, probably due to past soil drainage which enhance the respiratory processes at this site (Lindroth et al., 1998). Therefore, it is not surprising that the fit of modeled NEE to measured data was weaker compared to the GPP and ER models for both forest types, r 2 =61% (dec.) and r 2 =75% (con.) (Figs. 11 and 12;and Tables 4 and 6). In general, it is obvious that a small discrepancy 5 in the modeled GPP and ER can result in a large relative error in the final NEE. For example, the GPP model misses the growing season production peak in Hyytiälä by about 1 g m −2 d −1 which results in an average relative error of about 6%. Although rather small, this error, together with a very accurate estimation of ER, results in a 45% average relative error in final NEE for Hyytiälä. As a consequence, the accuracy of the 10 coniferous GPP and ER models, with coefficients of determination above 90%, results in 75% of the total variation in measured NEE is being explained by the models. These are still good estimates, proving the potential for operational remote sensing of boreal carbon balance. The modeled deciduous NEE is less accurate, with the carbon uptake at Sorø 2004 being completely missed due to the inaccurate 2004 GPP estimation.

15
In summary, GPP was modeled with higher accuracy for the coniferous sites than for the deciduous, in spite of the higher correlations observed between the satellite parameters and all ground measured biophysical parameters (Table 3). This could be attributed to the fact that only two deciduous site were used for parameterization and evaluation, with none of the sites being very variable from year to year. For the 20 coniferous data, two sites were used for parameterizing, sites that have a higher interannual variation, especially Norunda where, as mentioned, the higher variability may be attributed to soil disturbance caused by past drainage. Accordingly, the availability of measured biophysical data -GPP, ER, and in particular, NEE -is of utmost importance: the more data that is available, both in time and space, the more reliable 25 estimates can be obtained. For this study, it would be desirable to use more data in general, especially from other deciduous sites in Scandinavia, such as oak and birch. The latter species is abundant in forests in Northern Scandinavia.
Finally, annual modeled NEE, obtained by the averaging of the monthly estimates EGU over the years, was compared to averaged measured NEE. The result is seen in Fig. 13. The points for the coniferous data are clearly dispersed along the one-toone line (r 2 =79%), while the deciduous points deviate from the line. One explanation for the very low accuracy of the deciduous annual estimates compared to the monthly, is that the interannual NEE behavior is quite similar from year to to year with small 5 changes even between Hainich and Sorø. Again, the lack of deciduous data makes it hard to draw further conclusions on the annual NEE estimates. The coniferous sites exhibit more variation which in combination with a more accurate monthly modeling results in a high annual accuracy.

10
The results of this study illustrate the potential that remote sensing can be used for assessing the carbon balance of forested areas in Northern Europe. Linear relationships between parameters that can be made available from satellite (VIs, APAR, and AT) and measured biophysical parameters (NEE, max. NEE, GPP and ER) were found for both coniferous and deciduous data. APAR was highly correlated 15 with the biophysical parameters. There were no linear relationships between NDVI and the biophysical parameters found at any site, mainly because of a distinct saturation in the NDVI at higher uptake levels.