Carbon dioxide balance of subarctic tundra from plot to regional scales

We report here the carbon dioxide (CO 2) budget of a 98.6 km2 subarctic tundra area in northeast European Russia based on measurements at two different scales and two independent upscaling approaches. Plot-scale measurements (chambers on terrestrial surfaces, gas gradient method and bubble collectors on lakes) were carried out from July 2007 to October 2008. The landscape-scale eddy covariance (EC) measurements covered the snow-free period of 2008. The annual net ecosystem exchange (NEE) of different land cover types ranged from−251 to 84 g C m−2. Leaf area index (LAI) was an excellent predictor of the spatial variability in gross photosynthesis (GP), NEE and ecosystem respiration (ER). The plot-scale CO 2 fluxes were first scaled up to the EC source area and then to the whole study area using two data sets: a land cover classification and a LAI map, both based on field data and a 2.4 m pixel-sized QuickBird satellite image. The good agreement of the CO 2 balances for the EC footprint based on the different measuring techniques ( −105 to−81 g C m−2 vs.−79 g C m−2; growing season 2008) justified the integration of the plot-scale measurements over the larger area. The regional CO 2 balance based on area-integrated plot-scale measurements was −41 or −79 g C m−2 yr−1 according to the two upscaling methods, the land cover classification and the LAI map, respectively. Due to the heterogeneity of tundra, the effect of climate change on CO 2 uptake will vary strongly according to the land cover type and, moreover, likely changes in their relative coverage in the future will have great impact on the regional CO2 balance.


Introduction
The strong warming predicted for the Arctic by global climate models (Kattsov et al., 2005) has underlined the need to understand how carbon (C) fluxes in tundra will respond to climate change. Carbon storage in plant biomass is low in the Arctic (Bazilevich, 1993;Hugelius et al., 2011), but northern permafrost soils contain as much as 50 % of global belowground soil C pool (Tarnocai et al., 2009). There, higher temperatures and permafrost thaw will likely enhance the mineralization of soil organic matter in the future, thereby increasing release of C as carbon dioxide (CO 2 ) to the atmosphere (Dorrepaal et al., 2009;Schuur et al., 2009). At the same time, there will be changes in vegetation composition and productivity (Walker et al., 2006;Forbes et al., 2010). The increase in above-ground C stocks may partly compensate for the respiratory below-ground C losses, at least in the short term (Qian et al., 2010).
In order to predict the changes in the tundra C balance in the future, we need to accurately estimate the present-day C balance and understand its dependence on environmental factors. This is a real challenge, taking into account the high temporal (e.g., Kwon et al., 2006;Groendahl et al., 2007) and spatial (e.g., Heikkinen et al., 2004;Fox et al., 2008) variability of CO 2 fluxes in tundra, as well as the logistical difficulties when working in remote areas. Studies on tundra CO 2 exchange have been previously conducted using either the micrometeorological eddy covariance (EC) method (e.g.,  Kutzbach et al., 2007b;Lafleur and Humphreys, 2008;Humphreys and Lafleur, 2011;Parmentier et al., 2011;Lund et al., 2012) or chamber techniques (e.g., Christensen et al., 2000;Heikkinen et al., 2002Heikkinen et al., , 2004Shaver et al., 1998Shaver et al., , 2007. Although the EC and chamber methods are complementary of each other, there are only few studies where these two techniques have been applied in parallel in order to study CO 2 exchange in tundra (Soegaard et al., 2003;Zamolodchikov et al., 2003;Fox et al., 2008).
With high-frequency micrometeorological EC measurements it is possible to get continuous data set over a source area of 0.1-1 km 2 and catch the short-term variations in gas fluxes often missed by (manual) chambers (e.g., Pihlatie et al., 2010). On the other hand, chamber measurements give information on spatial variability in fluxes of different functional ecosystem types, often located at short distances in tundra landscape. This kind of data provides good grounds for extrapolating, or "upscaling" fluxes to a regional scale. However, upscaling based on chambers can be risky if the measuring frequency is low or the distribution of various surfaces is not well documented. Chamber measurements can also cause disturbance to the studied system and, thus, affect the gas flux rates observed (e.g., Kutzbach et al., 2007a). In the upscaling process, comparison of EC results and plotscale measurements in the EC source area serves as an intermediate step towards extrapolation to larger scale, making the regional estimates scientifically sound.
Due to the fragmented nature of tundra environment, it is critical to select a relevant spatial scale for upscaling efforts. Landsat satellite data with 30 m resolution, used in most studies on regional CO 2 balance of tundra (Soegaard et al., 2003;Heikkinen et al., 2004), are not detailed enough to detect the small scale differences typical for tundra (Laidler and Treitz, 2003). The first remote sensing studies on Arctic vege-tation using higher resolution than Landsat satellite images have been conducted only recently (Fuchs et al., 2009). In the present study, we use a land cover classification based on QuickBird imagery with 2.4 m resolution, appropriate for detecting the relevant spatial variability of the landscape investigated.
In this study, we aimed at quantifying CO 2 fluxes in a tundra landscape and investigating the controlling factors, in order to provide a sound basis for upscaling the fluxes to the regional level. The study was conducted in discontinuous permafrost region in Komi Republic, northern European Russia. We measured the CO 2 exchange over a whole snow-free period at two different scales: the chamber technique in terrestrial surfaces and the gas gradient method and bubble traps in lakes (plot scale), and the EC method (landscape scale). The results of these methods were compared in the EC source area using the footprint analysis. Plot-scale measurements, conducted over a prolonged period (two summers and one winter), give insights to the inter-and intra-annual variability in CO 2 exchange. The final goal was to estimate the CO 2 balance for a region of 98.6 km 2 based on the fluxes of different land cover types that were scaled up using two alternative methods. The upscaling methods were based on two high-resolution satellite derived data sets, a land cover classification and a leaf area index (LAI) map.

Study site
The study site, near the settlement Seida, is located in southern tundra with discontinuous permafrost in northeast European Russia (67 • 03 N, 62 • 56 E). The mean air temperature in the region is −5.6 • C and the mean annual precipitation Most of the landscape in the region is hilly upland with tundra heath vegetation ("upland tundra"; coverage 58 % based on a land cover classification described in Sect. 2.5.1) dominated by shrub-lichen-moss communities (Table 2). Also, peatlands with bog vegetation are typical for the region ("dry peatlands"; coverage 24 %). These peatlands include peat plateaus with up to several meters thick peat deposits, raised up by permafrost. The peat plateaus are spotted by unvegetated patterned ground features ("bare peat") and small thermokarst lakes ("lakes"). Narrow fens and willow stands ("wetlands"; 14 %) are located on low-lying parts of the landscape that act as water conduits in the terrain. The abovementioned land cover types, all studied at plot scale (Table 2), represent 97 % of the landscape. The residual landscape area (3 %) is covered by deciduous and coniferous forest stands, rivers, human-impacted tundra and sand. The terrestrial microsites and lakes studied at plot scale were located within the EC footprint area northeast from the EC tower (Fig. 1a). More details on the soils and vegetation at the Seida study site are provided by Hugelius et al. (2011) andMarushchak et al. (2011).

Auxiliary data
Auxiliary data were collected on air temperature, soil heat flux, net radiation, photosynthetically active radiation (PAR), precipitation, wind speed and wind direction by a weather station located at the study site. Parameters showing large heterogeneity across the landscape, such as soil moisture (θ v ), soil temperature, groundwater level and active layer depth, were monitored on various land cover types. The leaf area index of vascular plants (LAI) was measured in June-September 2008 at 284 field measurement points (Fig. 1b). More details on auxiliary data collection are available in the Supplement.

Chamber measurements at terrestrial microsites
Exchange of CO 2 was measured at the terrestrial microsites with a closed chamber technique (Heikkinen et al., 2002) in July-October 2007 and May-October 2008 (Table 2). Measurements were carried out mainly between 08:00 a.m. and 09:00 p.m. using a transparent polycarbonate chamber (30 × 60 × 60 cm) with an infrared gas analyzer (Li-840, LiCor). The chamber had a cooling unit with ice water circulation that kept the temperature in the chamber headspace close to the ambient. Data on CO 2 concentration, air temperature inside and outside the chamber (107 Thermistor Probe C/W, Campbell Scientific, UK) and PAR (SKP215, Skye Instruments, UK) were collected to a data logger (CR850, Campbell Scientific) at 2-s intervals during 2-min measurements. Net ecosystem CO 2 exchange (NEE) measurements were done in ambient and reduced light conditions (∼ 50 % PAR), created by shadowing the chamber with a net. For measuring ecosystem respiration (ER), the chamber was darkened with an aluminum lid. Each microsite had 3 replicate soil collars (depth 15 cm) that served as permanent measurement plots. At the willow site additional collars (height 60 cm) were used during the CO 2 flux measurements.
Carbon dioxide fluxes were calculated from the increase in CO 2 concentration in the chamber headspace using an exponential non-linear model by Kutzbach et al. (2007a) in MAT-LAB (R2008a) program. Residual standard deviation of the regression > 1.5 ppm was used as a filtering criterion, based on which 2 % of the flux data were rejected. Gross photosynthesis (GP) was calculated as a difference between consecutive NEE and ER measurements. We follow a sign convention where C loss from the ecosystem is defined as positive and C uptake as negative values.

Interpolating gross photosynthesis and ecosystem respiration over time
The two components of net CO 2 flux, GP and ER, were interpolated over the measuring period based on their dependence on environmental parameters using non-linear regression in SPSS 14.0 statistical software (Table 3). Response functions for ER were formulated individually for each chamber plot, while parameterization of the more complex GP functions was done at the microsite level (n = 3). Data were split by years, and for GP growing season 2008 was further divided into early and late season. The GP models included a Michaelis-Menten type equation for the light response, a Gaussian or linear temperature term and a linear LAI term. The LAI data from 2008 were used also for 2007, which was justified by the good model fit in 2007. For surfaces without any vascular plants ("bare peat") a soil moisture term was used instead of the LAI term. Ecosystem respiration was described with an Arrhenius-type temperature dependence (Lloyd and Taylor, 1994). Anomalously high respiration peaks at soil freezing and thawing (2 % of all data) were removed before the modeling (Supplement Fig. S1).

Carbon dioxide fluxes during snow period
During the snow period CO 2 fluxes were measured with a snow-gradient method as described by Marushchak et al. (2011). Fluxes were determined in January-June 2008 a total 2-5 times per plot depending the timing of the snow melt. The samples were stored in glass vials a maximum of three months before analysis on a gas chromatograph (HP 5890 series II, Hewlett-Packard, USA) with a thermal   (Nykänen et al., 1995). A leakage test with a gas standard (2500 ppm) showed that the reduction of CO 2 concentration in the sample vials over two months was ≤ 3 % (data not shown).

Carbon dioxide fluxes from lakes
Emission of CO 2 by diffusion and ebullition pathways was studied in three thermokarst lakes from July to August 2007 (11 samplings) and from June to October 2008 (19 samplings). The determination of CO 2 concentrations in the surface water and flux calculation using the thin boundary layer Q = Maximum GP; k = PAR level at which GP reaches half of Q; T 2 cm = Soil temperature at 2 cm; T mean = Average of air temperature and T 2 cm ; T opt = Temperature optimum of GP; T tol = Temperature tolerance of GP; a = Correction term for LAI; M = The volumetric soil moisture; R 10 = ER at 10 • C; E 0 = Activation energy of ER.
model followed Repo et al. (2007). Surface water samples were collected during daytime (08:00 a.m.-19:00 p.m.). Linearly interpolated daily CO 2 concentrations and hourly wind speed (measured at 2 m, normalized to 10 m using a logarithmic wind profile) were used to calculate hourly flux rates. Ebullitive CO 2 flux was monitored with permanently installed submerged funnel gas collectors (Repo et al., 2007). Each lake had 6-7 replicate gas collectors (Ø 0.35 m). Gas samples were stored and analyzed as described above.

Eddy covariance setup
Landscape-scale NEE was measured in May-October 2008 (139-280 days) at a frequency of 10 Hz using the micrometeorological EC method (Aubinet et al., 2000;Baldocchi, 2003). The system was set up with an R3 ultrasonic anemometer (Gill Instruments, UK) and an LI-7500 openpath CO 2 /H 2 O IRGA (LI-COR) mounted at 3.95 m height. The power was supplied by a fuel generator placed 40 m eastsoutheast (110 • ) of the mast.

Data processing
The raw data were processed using the Alteddy software (version 3.5, University of Wageningen, the Netherlands, http://www.climatexchange.nl/projects/alteddy/), which is based on EUROFLUX methodology (Aubinet et al., 2000). The means and variances of turbulent CO 2 fluxes (= NEE) were calculated in half-hour time steps.
The open-path gas analyzer bears the disadvantage of self heating of the instrument surface which may result in overestimation of the CO 2 uptake especially under cold conditions (e.g., Goulden et al., 2006;Ono et al., 2008). To account for this effect we applied the correction proposed by Burba et al. (2008). The quality of the data was assessed using quality flags according to Foken et al. (2004), and only data with quality flags between 1 and 6 were included in the final data set. Furthermore, data were rejected when measured during rain and fog events, when the wind was from the direction of the generator (90 to 130 • ) and when friction velocity (u*) was < 0.1 m s −1 . Out of the 6835 measured half-hourly NEE fluxes, 50 % passed these criteria. Subsequent gaps in the time series were filled using an online tool (http://www.bgc-jena.mpg.de/bgc-mdi/html/ eddyproc/) based on the algorithms by Falge et al. (2001) and Reichstein et al. (2005). A nighttime gap from 25 August to 6 October 2008 (238-280 days) was filled using a regression function based on soil temperature (Lloyd and Taylor, 1994;see below). During this period, 70 % of the nighttime values had to be discarded due to rain events, dew formation or low turbulence conditions, which made the use of the online tool unfeasible.

Nighttime fluxes and flux partitioning
Measured NEE was partitioned into its components, GP and ER, by modeling ER based on its temperature response. For this purpose, we used the nighttime NEE (PAR < 50 µmol m −2 s −1 ) that is equal to ER. Only flux data with a quality flag 1-3 were used. The regression was 442 M. E. Marushchak et al.: CO 2 balance of subarctic tundra performed on nightly average fluxes and soil temperatures by fitting data into an Arrhenius type function (Lloyd and Taylor, 1994): where R 10 is the ER at 10 • C and T s is the soil temperature at 5 cm in • C. GP was calculated as a difference between NEE and ER. The temperature response of ER is shown in Supplement Fig. S2.

Footprint analysis
The flux measured by the EC technique originates from a large number of ground level points located on various land cover types in the EC source area, the so-called footprint. The footprint area varies according to meteorological conditions. In order to upscale the results of the plot-scale measurements to the landscape scale we needed to know the contribution of each grid cell of the QuickBird image to the EC flux. The contribution of each 2.4 m × 2.4 m grid cell within a 2.5 km radius around the EC mast was calculated by spatial integration of a source weight function for hourly observation intervals. The footprint analysis was performed in MATLAB R2009a using the footprint model originally developed in 2-D by Gash (1986) and Schuepp et al. (1990) and expanded for use in 3-D by Soegaard et al. (2003). By superimposing these grid source weights upon the land cover classification and LAI map (see Sect. 2.5.), the contribution of each land cover type to the recorded flux and the mean LAI of the footprint could be estimated for upscaling purposes.
2.6 Upscaling of the fluxes using remote sensing data 2.6.1 Land cover classification The land cover classification, used as one method for scaling up the CO 2 balance, was based on a QuickBird satellite image covering 98.6 km 2 around the flux measurement site, acquired on 6 July 2007 (QuickBird © 2007, DigitalGlobe; Distributed by Eurimage/Pöyry) (Fig. 1a). Four channels were used in the classification procedure (blue, green, red and infrared (NIR); pixel size 2.4 m). Classifications were produced using a multiple level segmentation in the Definiens Professional 5.0 software. In the segmentation process, the neighboring pixels of the image are grouped together to form uniform regions using information on the spectral properties of the pixels. Vegetation descriptions made at 150 transect points and additional field notes and photographs were used as ground-truthing data. The classification was tested using vegetation descriptions from 130 randomly selected field points. The classification is described in more detail by Hugelius et al. (2011). From the total of 13 land cover types distinguished from the satellite image, 8 were studied with plot-scale techniques ( Fig. 1a; Table 2). Five land cover types corresponded directly to microsites studied at plot scale, while some incorporated more than one microsite type. The relative contribution of different microsites within these land cover types was estimated from photographs taken at transect points around the EC mast. The chamber fluxes of different land cover types were weighed by their relative area contributions in order to obtain CO 2 balance over the study area. For rivers, we used the CO 2 emission of 33 g C m −2 during the open water period (100 days), measured in the same region by Heikkinen et al. (2004). A zero CO 2 balance was assumed for forest stands, sand and impacted tundra.

Leaf area index mapping
An alternative upscaling data set, a LAI map, was developed based on the same QuickBird satellite image as was used for the land cover classification. The LAI map was based on a regression model, where mean LAI values for growing period measured in the field (Fig. 1b) were predicted by NDVI (normalized difference vegetation index: (NIR − red)/(NIR + red)) and individual channel reflectance values. Mean values from a 5 m radius around the LAI measurement points were calculated from the QuickBird data in order to eliminate spatial inaccuracies caused by GPS device. Different transformations were tried for both LAI and satellite image variables, and the best model fit was obtained when the LAI values were square root transformed. NDVI was the best single explanatory variable to explain LAI (30 % of variation), and by adding channels 1 (blue) and 2 (green) the explanatory power was increased to 36 % (n = 284): sqrt (LAI) = 4.437 · NDVI + 0.027 · Channel 1 − 0.014 (2) ·Channel 2 − 2.784.
The residuals of the model were normally distributed. The rather large unexplained variation can be explained by the fact that the LAI measurements were made at single points and vegetation cover around them was not uniform enough to perfectly match the satellite data. When the model was applied to the whole QuickBird image, predicted negative LAI values were reclassified as 0. Mean LAI of the terrestrial surfaces was calculated by subtracting the water bodies from the LAI of the whole region. Dependence between cumulative CO 2 fluxes and LAI, observed for terrestrial microsites (see Sect. 3.3 and Fig. 6), was used to upscale the CO 2 balance to landscape (EC footprint) and regional scales.

Uncertainty estimates for the CO 2 fluxes
Uncertainty of the CO 2 balance from area-integrated plotscale measurements was estimated by weighing the standard deviations flux measured by chambers (see Bubier et al., 1999). For EC fluxes, the gap-filling procedure was considered as the main uncertainty factor due to the high percent of gap-filled values (51 %). Artificial gaps were added to the data set and subsequently modeled by the usual gap-filling procedure described above. These artificial gap-filled fluxes and the original fluxes showed a good agreement with a R 2 of 0.85 and RSME of 1.3 µmol m −2 s −1 . The average daily difference between the artificial gap-filled and original measured fluxes was considered as the uncertainty of the EC flux.

Meteorology
Air temperature and precipitation in 2007 and 2008 measured at the study site and the long-term averages from the nearby Vorkuta station are shown in Table 1. Air temperatures were above the long-term means in December and January. July was anomalously hot during both years, especially in 2007. The precipitation during the snow-free period 2007-2008 was mostly equal to the long-term mean. The length of the thermic growing season, defined as a period when the daily mean air temperature is permanently above + 5 • C, was 80 days in 2007 and 79 days in 2008. The permanent snow cover lasted from mid-October to mid-May-early June.

LAI map
The distribution of LAI values predicted by the regression model based on spectral satellite data is realistic when assessed in relation to the land cover classification (Fig. 1a,b). Overall mean LAI for the region was 0.98 (max 5.36). The proportion of zero LAI values (1.8 %) matched very well to the coverage of water and non-vegetated soils in the land cover classification (1.6 %). There was a strong correlation between the LAI values measured at the chamber plots and those derived from the LAI map for the corresponding land cover types as a whole (p < 0.05; Fig. 2). The only land cover types deviating from 1 : 1 line are fen and Betula nana tundra heath. Fig. 3 shows the raw chamber data that were used for building the regression functions to obtain hourly CO 2 fluxes at the terrestrial microsites (Table 3). The regression functions were able to explain 86 and 74 % of the observed overall variability in the measured GP and ER fluxes. Details of the GP and ER models are presented in the Supplement (Tables S1  and S2). The spatial variability in soil conditions and vascular plant cover (Table 4) resulted in large differences in CO 2 fluxes across the landscape (Figs. 4-5). The growing season ER of terrestrial microsites varied from 64 to 226 g C m −2 , GP from −31 to −518 g C m −2 and NEE from −292 to 48 g C m −2 . The annual CO 2 balance ranged from −251 to 84 g C m −2 . A general trend of increasing C uptake with increasing wetness was observed for the three main land cover types (upland tundra < dry peatlands < wetlands), and wetlands also had the highest ER (Fig. 4).

Spatial variability in CO 2 exchange
The mean LAI was an excellent predictor of the spatial variability in GP and NEE, and even in ER. The cumulative CO 2 fluxes during the growing season were linearly dependent on LAI across all the terrestrial microsites studied (Fig. 6). Moreover, the predictive power of LAI was good also if the period under review was the whole EC measuring period (May-October), or even the full year (data not shown). These linear regression functions between LAI and CO 2 flux, observed at plot scale, were used together with the LAI map to upscale the CO 2 balance to landscape and regional scales (see Sect. 2.5.2.). The microsites with LAI < 0.3 (dry shrub tundra heath, dry lichen tundra heath, bare peat, Eriophorum dominated fen) had the lowest GP and were net sources of C to the atmosphere on an annual basis. All the other microsites showed annual net uptake of C.
The thermokarst lakes studied were supersaturated with CO 2 and, thus, atmospheric sources of this gas during the whole open water period (Fig. 7). The annual emissions of CO 2 from the three lakes ranged from 13 to 82 g C m −2 , averaging 43 ± 35 g C m −2 . Diffusion was the main pathway of CO 2 release, the importance of ebullition being negligible.

Seasonal variations
Plot-scale fluxes showed an initiation of the summer uptake period on day 170-178, while the first day with daily net C uptake observed by EC was day 174 (Fig. 4). Both nighttime respiration and daytime net uptake measured by EC reached the maximum values during the second half of July, being 4 µmol m −2 s −1 and −8 µmol m −2 s −1 , respectively. In autumn, the EC measurements indicated a shift to a daily source by day 240, while area integrated plot-scale measurements indicated this only on day 258 (Fig. 3). The seasonal amplitude of NEE was higher for EC fluxes than for the area integrated plot-scale measurements, i.e., EC showed higher net uptake rates during mid-summer and higher net CO 2 release during spring and autumn. Based on the plot-scale measurements, the non-growing season ER comprised 25-45 % of the annual ER. The measured values of wintertime respiration were 0-0.24 µmol m −2 s −1 in January and 0-0.22 µmol m −2 s −1 in March. The area-weighed ER for the EC footprint based on modeled ER of different land cover types was within the range of the measured respiration rates: 0.13 and 0.11 µmol m −2 s −1 for January and March, respectively. We observed a clear peak in ER at topsoil (2 cm) temperatures ∼ 0 • C in upland tundra and dry peatlands, but not in wetlands (Supplement Fig. S1). The peak was more pronounced in the spring than in the fall.

Interannual variations
Interannual comparison is limited to July and August, the only months with intensive chamber measurements during both study years (Fig. 8). Upland tundra (n = 12) had 25 % higher ER in July-August 2007 than in 2008 (Wilcoxon signed-ranks test; p < 0.05). 75 % of the cumulative difference in ER occurred in July, which had higher air and topsoil temperatures in 2007 than in 2008. As a result of the higher respiration, the net C sink on upland tundra was lower by 50 % in 2007, while the GP was similar during both years. No interannual differences were observed in the peak-season GP, NEE or ER on wetlands and dry peatlands. The CO 2 Fig. 4. Seasonal dynamics of CO 2 fluxes of the three main landscape types (upland tundra, dry peatlands and wetlands) based on chamber measurements and integrated landscape-scale CO 2 flux from eddy covariance. Data are daily mean values. emissions from thermokarst lakes were of similar magnitude in growing seasons 2007 and 2008 (20 ± 18 g C m −2 and 29 ± 22 g C m −2 , respectively).

Comparison of CO 2 balances by the two measuring techniques in landscape scale
The EC and area integrated plot-scale measurements resulted in similar CO 2 sink strength for the EC footprint during growing season 2008 (−81 ± 37 g C m −2 vs. −105 ± 27 to −79 ± 35 g C m −2 , respectively; Table 5). Outside the growing season EC showed a higher source than the plot-scale measurements (upscaled with land cover classification), the mean difference being 0.44 g C m −2 d −1 in the early season and 0.82 g C m −2 d −1 in the late season. Still, the CO 2 balances based on the two methods were within the error range of each other for the whole EC measuring period of 141 days. Upscaling with the LAI map resulted in a somewhat stronger sink than upscaling with the land cover classification, but the difference was not significant. The partitioning approach used for the EC fluxes (Fig. 4) allows the comparison of the two CO 2 flux components, ER and GP, between the two methods. Although the EC and area-integrated chambers resulted in similar cumulative NEE for the growing season, the plot-scale measurements showed lower fluxes to both directions. Cumulative ER and GP for the growing season from EC were 226 and −305 g C m −2 (−ER/GP = 0.74) and from plot-scale measurements 111 to 130 and −233 to −188 g C m −2 (−ER/GP = 0.56-0.59).

Seasonal and annual CO 2 balance for the study region
After verifying the plot-scale fluxes by comparing them with the EC results in the landscape scale, we scaled them up to the whole QuickBird area of 98.6 km 2 . The regional CO 2 balance during the growing season 2008 was −127 ± 30 to −94 ± 37 g C m −2 based on the two independent upscaling approaches. The annual CO 2 balance (NEE) for the study region was −79 ± 22 to −41 ± 57 g C m −2 , consisting of ER of 197 ± 35 to 212 ± 50 g C m −2 and GP of −294 ± 58 to −238 ± 39 g C m −2 . 33-39 % of the annual ER and 8-10 % of the annual GP occurred outside the growing season. The study region as a whole was a 16-21 % stronger CO 2 sink than the EC footprint area, which corresponds to a difference in LAI between the two scales (0.98 vs. 0.83 for the whole study region and EC footprint, respectively). The EC foot- print had, e.g., more tundra bog and lakes and less willow than the whole region (Supplement Fig. S3).

Plot-scale CO 2 fluxes
Despite the heterogeneity of the studied landscape, a clear trend was found to explain the spatial variability in CO 2 exchange. Vascular LAI explained ∼ 90 % of the variability in the growing season GP and NEE across the chamber microsites, and even for ER its explanatory power was 67 %. Given that the response functions for interpolating the CO 2 fluxes were determined independently for each microsite, this was a true dependence, not a modeling artifact. Also, in earlier studies on tundra and peatland ecosystems LAI has explained well the spatial variability in GP (Lund et al., 2010) or both GP and ER Mc-Fadden et al., 2003;Humphreys et al., 2006). Positive correlation between ER and LAI may be explained by dominance of autotrophic respiration over microbial decomposition of soil organic matter (heterotrophic respiration), or a tight link between these two respiration components through, e.g., rhizomicrobial respiration (Hutsch et al., 2002). Since non-destructive LAI measurements with a plant canopy analyzer are fast and easy and can be related to remote sensing data, the strong relationship between LAI and CO 2 fluxes is very promising from the upscaling point of view. Still, multiyear studies are needed to verify how strongly the shape and slope of the CO 2 flux-LAI dependence are affected by interannual variability in weather conditions.  We did not observe significant differences in GP between the study years at any of the land cover types. However, the peak summer ER was higher and net CO 2 sink smaller on upland tundra in 2007, which was the warmer of the two years. There was no similar increase in ER at the peatland sites, most probably due to wetter soil conditions that were limiting soil respiration over temperature. In the same region Heikkinen et al. (2004) measured higher ER from non-waterlogged soils during the warmer and drier of the two study years, and Zamolodchikov et al. (2000) observed a switch in net C flux from sink to source in shrub tundra communities when canopy temperature rose above +14 • C. Also, multiannual EC flux studies in other tundra regions have indicated good adaptation of certain vegetation communities to the present temperatures, seen as reduced net C sink during warm summers (Parmentier et al., 2011;Lund et al., 2012).

Differences in landscape-scale CO 2 fluxes between the two measuring techniques
Previous tundra CO 2 studies combining EC and chamber measurements are strikingly few (Table 6). This study demonstrates the suitability of this kind of nested study design for flux studies in heterogenous environments like tundra. The two methods resulted in very similar CO 2 balances for the growing season. However, the area-integrated plotscale measurements showed lower amplitude of the seasonal NEE cycle, and partitioning of the EC fluxes into ER and GP components revealed that the chamber-based CO 2 fluxes were lower to both directions. Further, chambers showed a slightly longer net uptake period than EC. Due to these differences, the plot-scale measurements resulted in a significant CO 2 sink for the EC measuring period from May to October, while the CO 2 balance from EC was also negative but did not differ significantly from zero.
A key reason for the differences between the methods is that the chamber fluxes used for the areal integration were modeling results whereas NEE from EC were mostly obtained by direct measurements. A feature of all models is that they are dominated by the data in the middle of the data range and tend to filter out the extremes. This has implications also for determination of the timing and length of the summer net uptake period. Here, more frequent measurements of LAI and fluxes during spring and autumn would have been helpful.
There are some problems related with the closed chamber technique that may cause underestimation of the gas flux, such as reduction of the concentration gradient between the soil and the atmosphere (Kutzbach et al., 2007a) or underestimation of the effective chamber volume when the air-filled soil pore space is not accounted for (Rayment, 2000). We took care of all necessary precautions to avoid chamber bias: the incubation time was kept short, and fluxes were calculated using non-linear regression. In turn, a known problem of the open path EC gas analyzer is overestimation of C uptake during the cold season due to heating of the sensor head (e.g., Goulden et al., 2006;Ono et al., 2008). Also, in this study a spring uptake up to 0.5 g C m −2 d −1 was observed while the landscape was still covered by snow, and the correction suggested by Burba et al. (2008) was applied on the data. Since the sensible heat flux in the optical path was not directly measured, the correction was based on the measured air temperature (Burba et al., 2008;Method 4). This may have caused overcorrection of the fluxes (Wohlfahrt et al., 2008;Bowling et al., 2010).
Biased selection of chamber plots (see Fox et al., 2008) or a failure to estimate correctly the coverage of different land cover types would also affect the flux estimates, particularly when the spatial variability in CO 2 fluxes is so high. In our study, we found a good match in LAI values between the chamber plots and the corresponding land cover types as a whole, showing that the chamber plots were well representative for the study region. Further, we can be sure that our QuickBird-based land cover classification represented better the small-scale variability in the patchy tundra vegetation than the previously used Landsat image classifications Heikkinen et al., 2004). While the size of one Landsat image pixel is 900 m 2 , the mean patch size of homogenous vegetation in this study was 816 m 2 , the average fen patch being only 260 m 2 .

Regional CO 2 balance
The regional CO 2 balance of the tundra area of 98.6 km 2 was obtained by upscaling the plot-scale measurements using two independent approaches. From October 2007 to October 2008, the studied tundra acted as a CO 2 sink of −79 to −41 g C m −2 yr −1 . The gross CO 2 fluxes (GP = −294 to −238 g C m −2 yr −1 and ER = 197 to 212 g C m −2 yr −1 ) were relatively large compared with the earlier studies summarized in Table 6. This could be due to the milder climatic conditions prevailing in the southern East-European tundra compared with high arctic and more continental sites (see Zamolodchikov and Karelin, 2001). In addition to large-scale variability driven by climate, large site-specific differences exist even within climatic zones. Lowlands with wet soils typically have higher C accumulation rates than upland tundra in the same region (Table 6; e.g., Vourlitis et al., 2000;Kwon et al., 2006), which is also seen as intrasite variability in this study.
The dependence of CO 2 on vascular LAI, discussed above, would provide a more generic explanation for the variability in CO 2 budgets in Table 6: The relatively high gross CO 2 fluxes at the Seida site could be related to high above-ground biomass compared with the reference studies. Although this hypothesis cannot be verified here since many of the reference studies do not provide data on LAI or biomass, there is published evidence on the potential of LAI to explain variability in CO 2 balance even across climatic zones (Lund et al., 2010).
The C budgets of northern ecosystems also show very large year-to-year fluctuation due to variability in local weather conditions (Table 6; Aurela et al., 2004;Groendahl et al., 2007;Lafleur and Humphreys, 2008;Lund et al., 2010Lund et al., , 2012. The two study years were warmer than the long-term mean, and particularly high temperatures were measured in July during the intensive growth period (3-5 • C higher than the long-term mean). This probably enhanced both flux components, GP and ER.
Importance of a specific land cover type for the regional CO 2 balance is a result of two factors: total coverage in the landscape and magnitude of the CO 2 fluxes per unit area. For example, the willow stands that cover only 8.7 % of the study region contributed to the growing season net CO 2 uptake by 27 %, being as important for the regional balance as the dominant land cover types, tundra bog and shrub tundra heath. The willows with high C sink capacity may become even more important with future warming, as suggested by the trend of increased willow growth observed in the region during the last decades (Forbes et al., 2010). Studies from northern Alaska show that also deciduous shrubs, dominant in upland tundra, will benefit from warmer temperatures (Tape et al., 2006). The positive dependence observed between LAI and NEE suggests that better growth of vascular plants would mean increased net C sink to the studied tundra ecosystem, at least in the short term, before the increased respiration due to soil C mobilization takes over as has been predicted (e.g., Qian et al., 2010).
An example of the changes expected in the landscape composition in the future is that the coverage of fen type peatlands will likely increase at the expense of peat plateaus with permafrost thawing (e.g., Johansson et al., 2006). This would increase the C uptake at regional scale because of the higher 450 M. E. Marushchak et al.: CO 2 balance of subarctic tundra CO 2 sink character of fens. In general, landscape reorganization will have a strong influence on regional CO 2 balance in the study region due to the large differences in CO 2 balance between different tundra surfaces.