Exchange of CO2 in Arctic tundra: impacts of meteorological variations and biological disturbance

An improvement in our process-based understanding of carbon (C) exchange in the Arctic and its climate sensitivity is critically needed for understanding the response of tundra ecosystems to a changing climate. In this context, we analysed the net ecosystem exchange (NEE) of CO2 in West Greenland tundra (64 N) across eight snow-free periods in 8 consecutive years, and characterized the key processes of net ecosystem exchange and its two main modulating components: gross primary production (GPP) and ecosystem respiration (Reco). Overall, the ecosystem acted as a consistent sink of CO2, accumulating −30 gCm−2 on average (range of −17 to −41 gC m−2) during the years 2008–2015, except 2011 (source of 41 gCm−2), which was associated with a major pest outbreak. The results do not reveal a marked meteorological effect on the net CO2 uptake despite the high interannual variability in the timing of snowmelt and the start and duration of the growing season. The ranges in annual GPP (−182 to −316 gCm−2) and Reco (144 to 279 gCm−2) were > 5 fold larger than the range in NEE. Gross fluxes were also more variable (coefficients of variation are 3.6 and 4.1 % respectively) than for NEE (0.7 %). GPP and Reco were sensitive to insolation and temperature, and there was a tendency towards larger GPP and Reco during warmer and wetter years. The relative lack of sensitivity of NEE to meteorology was a result of the correlated response of GPP and Reco. During the snow-free season of the anomalous year of 2011, a biological disturbance related to a larvae outbreak reduced GPP more strongly than Reco. With continued warming temperatures and longer growing seasons, tundra systems will increase rates of C cycling. However, shifts in sink strength will likely be triggered by factors such as biological disturbances, events that will challenge our forecasting of C states.


Introduction
Quantifying the climate sensitivity of carbon (C) stocks of the terrestrial biosphere is a major challenge for Earth system science (Williams et al., 2005). In the Arctic, organic soil C storage has the potential for very large C releases following thaw (Koven et al., 2011) that could create a positive feedback on climate change and accelerate the rate of global warming. Recent reviews have estimated the Arctic terrestrial C pool to be 1400-1850 Pg C, more than twice the size of the atmospheric C pool (Hugelius et al., 2014;McGuire et al., 2009;Tarnocai et al., 2009) and approximately 50 % of the global soil organic C pool (AMAP, 2011;McGuire et al., 2009). Further, Arctic ecosystems have experienced an intensified warming tendency, reaching almost twice the global average (ACIA, 2005;AMAP, 2011;Callaghan et al., 2012c;Serreze and Barry, 2011). The projected Arctic warming is also expected to be more pronounced in coming years (AMAP, 2011;Callaghan et al., 2012a;Christensen et al., 4468 E. López-Blanco et al.: Exchange of CO 2 in Arctic tundra perature, precipitation and growing season length will likely increase in the Arctic (ACIA, 2005;Christensen et al., 2007Christensen et al., , 2004IPCC, 2007). Given this situation, an improvement in our process-based understanding of CO 2 exchanges in the Arctic and their climate sensitivity is critical (McGuire et al., 2009).
Measuring the interannual C exchange variability in the Arctic tundra is challenging due to extreme conditions and the patchy nature of the landscape linked to microtopography. Different eco-types are linked to different C exchange rates (Bubier et al., 2003). Synthesis studies have found a significant spatial variability in NEE Mbufong et al., 2014) between different tundra sites Lund et al., 2010) and also large temporal variability within sites (Aurela et al., 2004Christensen et al., 2012;Grøndahl et al., 2008;Lafleur et al., 2012). Minor variations in the key process of photosynthesis (gross primary production, GPP) and ecosystem respiration (R eco ) may promote important changes in the sign and magnitude of the C balance (Arndal et al., 2009;Elberling et al., 2008;IPCC, 2007;Lund et al., 2010;Tagesson et al., 2012;Williams et al., 2000). With continued warming temperature and longer growing seasons, tundra systems will likely have enhanced GPP and R eco rates, but long-term data with which to investigate and quantify these responses are rare. Further, the effects on net CO 2 sequestration are not known, and may be altered by long-term processes such as vegetation shifts and short-term disturbances like insect pest outbreaks, complicating the prognostic forecast of upcoming C states (Callaghan et al., 2012b;McGuire et al., 2012). Consequently, there is a need to understand how the C cycle behaves over timescales from days to years and the links to environmental drivers. There is a lack of reference sites in the Arctic from which full measurement-based data are available, documenting carbon fluxes at the terrestrial catchment scales. Here we investigate the functional responses of C exchange to environmental characteristics across eight snowfree periods in 8 consecutive years in West Greenland.
In recent decades, eddy covariance has become a fundamental method for carbon flux measurements at the landscape scale Lund et al., 2012;Reichstein et al., 2005). Eddy covariance measurements of landatmosphere fluxes or net ecosystem exchange (NEE), of CO 2 can be gap-filled and subsequently separated into the modulating components of GPP and R eco using flux partitioning algorithms (Reichstein et al., 2005). These techniques are critical for providing a better understanding of the C uptake vs. C release behaviour (Lund et al., 2010), but they also allow for an examination of the environmental effects on ecological processes (Hanis et al., 2015). However, large gaps in the measured fluxes may introduce significant uncertainties in the C budget estimations. Moreover, GPP and R eco estimates can be calculated in different ways. Some algorithms fit an instantaneous temperature-respiration curve to night-time data to calculate R eco and estimate GPP Reichstein et al., 2005); others calculate R eco from a light-response curve (Gilmanov et al., 2003;Lindroth et al., 2007;Lund et al., 2012;Mbufong et al., 2014;Runkle et al., 2013). Unfortunately, different interpretations of the flux gap filling and partitioning lead to different estimates of NEE, GPP and R eco as well as undefined uncertainties.
The main objectives of this paper are (1) to explore the uncertainties in NEE gap filling and partitioning obtained from different approaches, (2) to determine how C uptake and C storage respond to the meteorological variability, and (3) to identify how the environmental forcing affects not only the interannual variability, but also the hourly, daily, weekly and monthly variability of NEE, GPP and R eco . The intention of this paper is to elaborate on the information gathered in an existing catchment area under an extensive crossdisciplinary ecological monitoring programme in low Arctic West Greenland, established under the auspices of the Greenland Ecosystem Monitoring (GEM) (http://www.g-e-m.dk). Using a long-term (8-year) data set to explore uncertainties in NEE gap-filling and partitioning methods and to characterize the interannual variability of C exchange in relation to driving factors can provide our understanding of landatmosphere CO 2 exchange in Arctic regions with a novel input. Our overarching hypothesis was that both GPP and R eco would respond positively to warmer and longer growing seasons. However NEE response to warming would be more complex and variable (positive or negative) depending on subtle balances between plant and microbial climate sensitivity.

Site description
Field measurements were conducted in the low Arctic Kobbefjord drainage basin in south-western Greenland (64 • 07 N; 51 • 21 W) (Fig. 1a). The study area is located ∼ 20 km SE of Nuuk, the Greenlandic capital. Kobbefjord has been subject to extensive environmental research activities (the Nuuk Ecological Research Operations) since 2007 (http://www.nuuk-basic.dk). The lowland site is located 500 m from the south-eastern shore of the bottom of Kangerluarsunnguaq Fjord (Kobbefjord), and 500 m from the western shore of the 0.7 km 2 lake called "Badesø" (Fig. 1b). Three glaciated mountains, all above 1000 m a.s.l., surround the site. The landscape consists of a fen area surrounded by heath, copse and bedrock. The current fen vegetation is dominated by Scirpus cespitosus, whereas the surroundings are dominated by heath species such as Empetrum nigrum, Vaccinium uliginosum, Salix glauca and copse species such as S. glauca and Eriophorum angustifolium (Bay et al., 2008). Kobbefjord belongs to the "arctic shrub tundra" (bioclimate zone E) according to The Circumpolar Arctic Vegetation Map (CAVM Team, 2003;Walker et al., 2005). This map is based on the summer warmth index (SWI), which is the sum of the monthly mean temperature above 0 • C from May to September and the southernmost bioclimatic zone E has limits of 20-35. In 2010 and 2012, the weather conditions led the area to experience temperatures from warmer climatic zones (SWI ca. 36 and 35 respectively). For the 1961-1990 period, the mean annual air temperature was −1.4 • C and the annual precipitation was 750 mm (Cappelen, 2013). The sunlight hours between May and September range from 14 to 21 h. Outcalt's frost number (Nelson and Outcalt, 1987) indicates that discontinuous permafrost should be present, although no permafrost has been found. Nonetheless, thin lenses of ice may remain until late summer.

Measurements
We have used eddy covariance (EC) data on NEE, measured during the snow-free period from 2008 to 2015. Measurements typically started around the end of the snowmelt (ca. May-June) and extended until the freeze-in period (between September and October). Once the snow melts, the growing season (i.e. the part of the year when the weather conditions allow plant growth) has been reported as the most relevant period defining both spatial (Lund et al., 2010;Mbufong et al., 2014) and temporal (Aurela et al., 2004;Groendahl et al., 2007;Lund et al., 2012) CO 2 variability. The EC measurements were conducted at the EddyFen station ( Fig. 1b and c), located in a wet lowland, 40 m a.s.l. The EC tower is equipped with a closed-path infrared CO 2 and H 2 O gas analyser LI-7000 (LI-COR Inc, USA) and a 3-D sonic anemometer Gill R3-50 (Gill Instruments Ltd, UK). The anemometer was installed at a height of 2.2 m, while the air intake was attached 2.0 m above terrain on the steel stand. Adjacent to the EddyFen station, an independent system ( Fig. 1b and c) measures round-the-clock net CO 2 fluxes using an automatic chamber (AC) method based on Goulden and Crill (1997).
The transparent chambers, each covering a known surface area of 60 cm by 60 cm, with a height of 30 cm, can be opened and closed by the computer in succession for 10 min every hour. When the chamber closes, a CO 2 analyser (SBA-4, PP Systems, UK) monitors both the CO 2 concentration by a close loop of tubing (further information about the set up can be found in Mastepanov et al. (2013). Nearly 20 m from the EddyFen station, the automated SoilFen ( Fig. 1b and c) station provides environmental variables such as air and surface temperature (Vaisala HMP45C), soil temperature at different depths (Campbell scientific 10ST) and relative humidity (Vaisala HMP45C). Two kilometres from these stations, an automatic weather station provides complementary ancillary data such as short-and long-wave radiation (with a CNR1 instrument), photosynthetic active radiation (with a Kipp & Zonen PAR Lite instrument), precipitation (using an Ott Pluvio instrument) and snow depth (with a Campbell Scientific SR 50). The water table depth data were monitored using a piezometer located next to each of the six autochambers. Finally, a robust daily estimate of the timing of snowmelt was analysed at a pixel level from a timelapse camera (HP e427) located at 500 m a.s.l. (Westergaard-Nielsen et al., 2013).

Data handling 2.3.1 Data collection and pre-processing
Data collection from the EddyFen station was performed using Edisol software (Moncrieff et al., 1997). Raw data files were processed using EdiRe software (version 1.5.0.32, R. Clement, University of Edinburgh) calculating the CO 2 fluxes on a half-hourly basis. The flux processing integrated despiking (Højstrup, 1993), 2-D rotation, time lag removal by covariance optimization, block averaging, frequency response correction (Moore, 1986) and Webb-Pearman- Leuning correction (Webb et al., 1980). For more information, see Westergaard-Nielsen et al. (2013). Ancillary data (air temperature, soil temperature, incoming short-wave radiation, relative humidity, PAR and precipitation) were temporally resampled using R (R Development Core Team, 2015). Time-series-related packages such as zoo (Zeileis and Grothendieck, 2005), xts (Ryan and Ulrich, 2014) and lubridate (Grolemund and Wickham, 2011) were used to get the ancillary data aligned with the flux data on a half-hourly basis.

Generating robust and complete flux time series
Before the CO 2 flux time series were analysed, we applied three different processing techniques (u * filtering, gap filling and partitioning) to (1) filter the NEE data for quality, (2) fill the NEE gaps and (3) separate NEE into GPP and R eco . The identification of periods with insufficient turbulence conditions (indicated by low friction velocity u * ) is important for avoiding biases and uncertainties in EC fluxes. To control the data quality, the u * thresholds were bootstrapped by identifying conditions with inadequate wind turbulence according to the method described in (Papale et al., 2006). We subsetted the data to similar environmental conditions, aside from friction velocity: 8 years and 7 temperature classes. Within each year/temperature subset the u * threshold (5, 50 and 95 % of bootstrap) was estimated at 1000 samples per year. We used the subsequent gap filling and partitioning based on these different subsets to propagate the uncertainty of u * threshold estimation across NEE, GPP and R eco .
Our gap-filling method was similar to Falge et al. (2001), using the marginal distribution sampling (MDS) algorithm, re-adapted from Reichstein et al. (2005) in REddyProc (Reichstein et al., 2016). MDS takes into account similar meteorological data available with different window sizes (Moffat et al., 2007). Parallel to this approach, we also gap-filled the original EC NEE data with an independent AC NEE data set (2010)(2011)(2012)(2013). AC data were collected simultaneously with EC data, and so we can used them as a cross check. The EC NEE was predicted from AC NEE based on linear regression models. The subsequent product was gap-filled using the MDS algorithm (REddyProc).
We separated NEE into its two main components (GPP and R eco ) using two approaches: (1) the REddyProc partitioning tool (Reichstein and Moffat, 2014) and (2) a light-response curve (LRC) approach Lund et al., 2012). A brief description of each flux partitioning method is provided in the Supplement (Eq. S1). After the flux partitioning comparison, we used ReddyProc-based GPP and R eco estimates on further analyses.

Flux uncertainties
In order to estimate the NEE gap-filling uncertainty, we assessed three different sources of uncertainty. First, we addressed the 95 % confidence interval of the EC prediction based on AC data. Second, we inferred the random uncertainty of filled half-hourly values from the spread of variables with otherwise very similar environmental conditions. REd-dyProc uses the gap filling to also estimate an observation uncertainty for the measured NEE, by temporarily introducing artificial gaps (T. Wutzler and M. Migliavacca (BGC-Jena), personal communication). Finally, we assessed the effect of uncertainty in the estimate of the u * threshold. In the u * -NEE relationship we want to exclude the probably false low fluxes (absolute NEE values) at low u * . When choosing a lower u * threshold, the associated lower flux will contribute to the gap filling and the annual sums. Therefore, there is a tendency of a lower absolute NEE associated with lower u * . The difference between the 5 and 95 % of bootstrap provides a means of the uncertainties based on the u * filters. We summed and propagated all these sources of uncertainties over time. The GPP and R eco uncertainties include the bias from the oneto-one flux comparison obtained from each model. The micrometeorological sign convection used in this study present uptake fluxes (GPP) as negative, while the released fluxes (R eco ) are shown as positive.

Identifying environmental forcing
Snow-and phenology-related variables such as the end of the snowmelt period and the start, end and length of the growing season are important components that shape the Arctic CO 2 dynamics. In this study we defined the end of the snowmelt period as the day of year when more than 80 % of the surface of the fen was considered snow free; the threshold was chosen in agreement with suggestions previously reported in Hinkler et al. (2002) and Westergaard-Nielsen et al. (2015). For the start, end and length of the growing season (GS start , GS end , GS length ); the GS start and the GS end were defined as the first and last days on which the consecutive 3-day NEE average was negative (i.e. CO 2 uptake) and positive (i.e. CO 2 release) respectively (Aurela et al., 2004), while GS length is the number of days between GS start and GS end ).
A random forest machine-learning algorithm (Breiman, 2001;Pedregosa et al., 2011) was utilized in a data-mining exercise to identify how the environmental controls affect the variability of NEE, GPP and R eco . Random forest calculates the relative importance of explanatory variables over the response variables. Here, we use photosynthetic active radiation (PAR), air temperature (T air ), precipitation (Prec) and vapour pressure deficit (VPD) to explain the response of C fluxes (NEE, GPP and R eco ) to climate variability. Each decision tree in the forest is trained on different random subset of the same training data set. The random forest is a classifier that groups explanatory variables and, in each fi- nal cluster, a multiple linear regression is built to reproduce fluxes as function of driving factors. This approach has been used to extrapolate maps of biomass (Baccini et al., 2012;Exbrayat and Williams, 2015). This version of random forest sums the relative importance of each variable from 0 up to 100 %, which correspond to the fraction of decision in which a variable is involved to cluster the data. We applied random forests to assess the relative importance of PAR, T air , Prec and VPD at different temporal scales (hourly, daily, weekly and monthly), aggregating them at the timescale indicated and lumping all the years together. (Table S1; Supplement). Moreover, we evaluated the diurnal, seasonal and annual pattern for each explanatory variable (data binned per hour; this is one random forest per hour of the day, day of the year and year respectively). To make sure that these results were not an artefact of the partitioning method that is based on a relationship between hourly R eco and T air , we performed the same analyses using daytime and night-time only hourly NEE as respective proxies for GPP and R eco . Based on these results (Table S2) we concluded that the approach was robust for the Kobbefjord site.

Interannual and seasonal variation of environmental and phenological variables
The annual mean temperature documented from Nuuk The end of the snowmelt period and the growing season start and length presented high interannual variability (CVs were 9.5, 9.0 and 19.0 %). Kobbefjord became snow free on DOY 154 (3 June for non-leap years, SD = 15) on average. On average, the site switched from being a source of CO 2 to a sink (GS start ) on DOY 175 (24 June, SD = 20), and remained so (GS end ) until DOY 241 (29 July, SD = 8.4) (Table 1). The GS start and the GS length did not follow a consistent pattern among the analysed years, the growing season timing have fluctuated substantially. The high interannual variability of the GS start correlated with variations in temperature, end of snowmelt period and VPD (p < 0.05). The highest variability was observed during 2009-2012. The 2010's GS length was nearly twice as long as in 2011. Indeed, GS start in 2011 differs only by 26 days from the GS end in 2010.  Green represents C uptake while the orange-dark-red denotes C release. The solid lines represent the end of the snowmelt period, while the area within the dashed lines represents the period between the start and the end of the growing season.

Data processing and quality
The NEE gap filling and subsequent partitioning obtained from different approaches exposed inconsistencies in performance and specific uncertainties in the seasonal C budget calculation. During the eight study snow-free periods, data gaps made up 46.5 % of the record from the EddyFen station due to unfavourable micro-meteorological conditions, instrument failures, maintenance and calibration (Jensen and Christensen, 2014) but also due to the rejection of lowquality flux measurements or too low u * . In 2014 a major instrument failure forced the station to stop measurements in the middle of the season. In 2010 and 2012 there were two more interruptions in the measurements (data gaps of > 20 days), although the problems could be solved before the end of the season. Such prolonged gaps led to unreliable gap-filled NEE estimates. The REddyProc MDS algorithm tended to fill these large gaps with high peaks of respiration at noontime, coercing C uptake underestimation. For this reason, an independent AC NEE data set (2010)(2011)(2012)(2013) was tested to gap-fill EC data (Figs. 3 and S2). The R 2 obtained from the EC-AC correlations was always > 0.70 (2010: R 2 = 0.80, p < 0.001; 2011: R 2 = 0.72, p < 0.001; 2012: R 2 = 0.80, p < 0.001; 2013: R 2 = 0.84, p < 0.001). By using AC data, the proportion of missing data was reduced to 28 %, and we found that the random uncertainty from the combination of AC and MDS algorithm decreased by 5 % on average. By using the u * filtering and the AC data together with EC, there was an increase in ∼ 6 % in terms of C sink strength. Moreover, the propagated uncertainty in NEE never exceeded ±1.8 g C m −2 , mainly because the error related to u * filtering was low. Further, we hypothesized that different flux partitioning approaches would lead to different estimates of GPP and R eco . However, the results suggest a relatively good agreement (Fig. 4). There was a higher degree of agreement with regard to GPP (R 2 = 0.83) compared with R eco (R 2 = 0.30). LRC tended to estimate 12 and 15 % larger GPP and R eco respectively compared to REddyProc.

Interannual and seasonal variation of CO 2 ecosystem fluxes
Overall, land-atmosphere CO 2 exchange measured for the snow-free periods of 2008-2015, omitting 2011, acted as a sink of CO 2 , taking up −30 g C m −2 on average (range −17 to −41 g C m −2 ) ( Fig. 5; Table 2). The cumulative NEE showed a characteristic pattern during the measurement period (Fig. 5), with an initial loss of carbon in early spring right after snowmelt (also observed in Fig. 3), followed by an intense C uptake as assimilation exceeded respiratory losses, triggered by increases in temperature, PAR and vegetation growth. This transition point matched the growing season start, when NEE switched from positive values (a net C source) to negative values (a net C sink). Eventually, the ecosystem turned again into a net C source, defining the growing season end. Even with high interannual variability in terms of the end of snowmelt time and growing season start/length (Table 1), the results do not show a marked meteorological effect on the NEE. The ranges in annual GPP (−182 to −316 g C m −2 ) and R eco (144-  279 g C m −2 ) ( Table 2) were > 5 fold larger and more variable (CVs are 3.6 and 4.1 % respectively) than for NEE (0.7 %). There was a tendency towards larger GPP and R eco during warmer and wetter years (Fig. S3), but there were no warmer and drier years during the study period. The strongest growing season CO 2 uptake occurred in 2012 (NEE = −74.2 g C m −2 ; GS length = 78 days), followed by 2010 (NEE = −70.0 g C m −2 ; GS length = 85 days) (Tables 1  and 2). A lengthening of the growing season did not increase the net carbon uptake in this study. In other words, an earlier end of the snowmelt resulting in a longer growing season length did not lead to a stronger carbon sink.
The anomalous year, 2011, constituted a relatively strong source of CO 2 (41 g C m −2 ) and was associated with a major pest outbreak, which reduced GPP more strongly than R eco . Data on the larvae of the moth Eurois occulta, collected from pitfall traps in the surrounding Salixand Empetrumdominated plots, showed a strong peak at the beginning of the 2011 growing season (Lund et al., 2017), coinciding with high NEE and very low GPP (Fig. 4). In 2011 up to 2078 larvae were observed, while in other years only 14 (2008), 82 (2009), 186 (2010), 0 (2012) and 8 (2013) were observed. It is likely that the reduced primary production in the wetland area was a partial response to the Eurois occulta outbreak.
The daily aggregated NEE-GPP relationships displayed consistent linear correlation (2008-2015: R 2 = 0.77, p < 0.001) across the assessed years (Fig. 6a). The linear cor-relations were weaker in 2010 and 2011. A hysteresis was detected in 2010 (i.e. long growing season with higher R eco in autumn than in spring), while strong C releases were observed in 2011 across June and July. The relation between GPP and R eco , which can be understood as the degree of coupling between inputs and outputs of C and therefore the degree of C sink strength, showed non-linear patterns (Fig. 6b). The curved behaviour is likely because GPP increased more than R eco during early growing season, except for in 2011. Moreover, R eco lagged behind GPP due to (1) the vegetation green-up in the first part of the growing season and (2) the higher respiration rates due to increased biomass in the second part. The years with clearer hysteresis coincide with the years with positive temperature anomalies (i.e. 2010, 2012 and 2013) of the 2008-2015 series. It is worth mentioning the different directions (clockwise or anticlockwise) in the hysteresis observed in these years between June, July and August. The data suggest that the clockwise 2012 hysteresis was due to greater gross C cycling (GPP and R eco ) in June and July favoured by warmer conditions, while in 2010 (anticlockwise hysteresis), the higher gross C fluxes were measured in August with warmer and wetter conditions (Fig. S4).

Environmental forcing
The varied importance of meteorological variables (such as PAR, T air , VPD and precipitation) obtained from random forest at different temporal scales (hourly, daily, weekly and Where applicable: ± sum of the autochamber, random and u * filtering uncertainties. * incomplete growing season data set. monthly) showed differences in behaviour depending on the time aggregation utilized (Fig. 7). PAR dominated NEE and GPP while T air correlated the most with R eco in hourly averages, whereas T air became increasingly important at longer temporal aggregations for all the fluxes (Fig. 7). VPD and precipitation were not as important as the other variables while the use of water table depth in the analysis was discarded due to its very low impact on CO 2 fluxes. In general, NEE and GPP showed similar distributions of importance, reinforcing the linear relationships found between NEE and GPP (Fig. 6). The standard deviation of the variables' importance (across 1000 decision trees) tended to increase at coarser time aggregations.
Changes of environmental forcing (PAR, T air and VPD) across diurnal, seasonal and annual timescales reveal patterns of functional responses to C fluxes. The diurnal cycle analyses on hourly data showed the changes in importance between day-and night-time (Fig. 8). NEE and GPP had two predominant variables (T air and PAR) determining the variability at daytime. PAR was important at dawn (06:00 WGST) and dusk (20:00 WGST), while T air was more important at other times. This performance indicates a threshold response to PAR, and a more continuous response to temperature. On the other hand, R eco was mainly driven by T air at both night-time and daytime. VPD and PAR had a negligible impact on R eco . The seasonal pattern importance showed Figure 7. Importance of environmental variables PAR (yellow), T air (orange), Prec (pink) and VPD (green), explaining variability in NEE, GPP and R eco (partitioned by REddyproc) at different temporal aggregations (hourly, daily, weekly and monthly) when all the years were lumped together. Thick bars and error bars represent the mean ± standard deviation of the importance across 1000 decision trees.
PAR dominating NEE and GPP from early June to early October (Fig. 8), while T air and VPD became more important before and after the snow-free conditions. In terms of CO 2 emission (R eco ) the pattern is less clear and noisier, although T air appeared to be the most important variable. Finally, the annual pattern exposes a performance in line with previous results; i.e. PAR dominated NEE and GPP while R eco was more sensitive to variations of T air . Interestingly, the random forest analysis revealed a decrease in PAR's importance in 2011, the same year in which the sharp decrease in C sink strength was exposed.

Data processing and quality
The NEE gap filling and subsequent partitioning into GPP and R eco are needed to understand the CO 2 flux responses to the environmental forcing. However, these procedures expose unavoidable uncertainties in the seasonal C budget calculation (Table 2) and partial inconsistencies between approaches (Fig. 4). In this study, we used an MDS gap-filling technique, an enhancement to the standard look-up table.
Both methods have shown a good overall performance compared to other procedures such as non-linear techniques or semi-parametric models but slightly inferior to artificial neural network (Moffat et al., 2007). However, the MDS gap filling alone introduced NEE estimates out of range across the two extensive gaps in 2010 and 2012 (Fig. S2). Quantifying the uncertainty introduced by measurement gaps is complex (Falge et al., 2001;Moffat et al., 2007;Papale et al., 2006). One possibility would be a sensitivity analysis of time series with artificially introduced gaps (Dragomir et al., 2012;Pirk et al., 2017). However, the choice of gap length and position is difficult and would render uncertainty to the uncertainty assessment itself. Instead, we used the EC prediction based on independent autochamber (AC) measurements between 2010 and 2013. The agreement between EC and AC was always R 2 > 0.72 and p < 0.001, and the 95 % confidence intervals of the predictions were reported together with the resulting uncertainties (Table 2). Although the AC data itself incorporated a new source of uncertainty in the calculations, we consider this method to be less weak than an unreliable gap-filling estimate. We used the AC as platform with which to decrease the gap length and the total random uncertainty (Aurela et al., 2002) before the MDS algorithm was applied. AC was used together with MDS and was never used as an independent gap-filling procedure.
The NEE partitioning obtained from REddyProc and LRC suggests a relatively good agreement in model performance. The one-to-one comparison between different approaches found a better agreement with regard to GPP compared to R eco . In this analysis, REddyProc produced smoother R eco estimates compared to the noisier GPP estimates, whereas the LRC results were the other way around. This is mainly because measurement noise goes into GPP for the REd-dyProc method, and into R eco for the LRC method. REd- Figure 8. Diurnal, seasonal and annual importance of environmental variables PAR (yellow), T air (orange), and VPD (green), explaining variability in NEE, GPP and R eco . Thick lines and shading represent the mean ± standard deviation of the importance across 1000 decision trees. dyProc retrieves positive GPP values, whereas the LRC method results in negative R eco values. Both scenarios are not fully convincing, although it is not straightforward as to how they should be treated. Removing all positive GPP/negative R eco values would risk removing only one side of the extremes. Besides night-time-based (REddyProc) and daytimebased (LRC) partitioning approaches, several implementations have been proposed to improve the algorithm's performance. Lasslop et al. (2010) has modified the hyperbolic LRC to account for the temperature sensitivity of respiration and the VPD limitation of photosynthesis. Further, Runkle et al. (2013) proposed a time-sensitive multi-bulk fluxpartitioning model, where the NEE time series was analysed in 1-week increments as the combination of a temperaturedependent R eco flux and a PAR-dependent flux (GPP). However, it remains uncertain as to under which circumstances each partitioning approach is more appropriate, especially in the boundaries between low-and high-Arctic due to the lack of dark night during polar days (when light is not a limiting factor for plant growth). Since there are few methods with an unclear precision, an evaluation study on the effect of using different partitioning approaches along latitudinal gradients would be very beneficial to assessing the suitability for each method.

Interannual and seasonal variation of CO 2 ecosystem fluxes
The balance between the two major gross fluxes in terrestrial ecosystems, photosynthetic inputs (GPP) and respiration outputs (R eco ), displayed larger temporal variability than did NEE. These results suggest that both GPP and R eco were strongly coupled and sensitive to meteorological conditions such as insolation and temperature (Figs. 7 and 8). Interestingly, the tendency to warmer and wetter conditions led to greater rates of C cycling associated with larger GPP and R eco (Fig. S3). This result does not entirely coincide with Peichl et al. (2014), even though they performed a similar analysis for a Swedish boreal fen. This finding points towards the complexity in the response of wetland ecosystems towards changing environmental conditions. The response is dependent on many things, such as hydrological settings, and these differ between sites. In this study, larger rates of C uptake (GPP) were linked to larger rates of C release (R eco ), with the exception of the anomalous year 2011. The relative insensitivity of NEE to meteorological conditions during the snow-free period could be the result of the correlated response of ranked cumulative GPP and R eco (Fig. 5) Wohlfahrt et al., 2008). This site likely receives more precipitation than many other tundra ecosystems and has no permafrost; thus the NEE response to climate could be less variable. However, as Kobbefjord is located in a coastal area, it is not surprising that it receives high pre-cipitation, and other ecosystems such as coastal blanket bogs often receive even more precipitation without a clear impact of drought effect on the NEE sensitivity . Furthermore, permafrost adds another layer of complexity to the C dynamics (Christensen et al., 2004;Koven et al., 2011;Schuur et al., 2015). Although some studies showed similarities of CO 2 fluxes in various northern wetland ecosystems with and without permafrost , permafrost has a strong influence on the hydrology of peatlands (Åkerman and Johansson, 2008), and therefore their topography and distribution of vegetation (Johansson et al., 2013). Especially in the context of climate warming permafrost thaw can cause large changes to the ecosystems. Further, this study agrees with Parmentier et al. (2011) andLund et al. (2012), who suggested that a longer growing season does not necessarily increase the net carbon uptake. Here a more negative NEE indicated a stronger C sink (i.e.) in 2012 compared to 2010. Parmentier et al. (2011) hypothesized that this behaviour is due to site-specific differences, such as meteorology and soil structure, and that changes in the carbon cycle with longer growing seasons will not be uniform around the Arctic. Thus, the effects of climate change on the tundra C balance of are not straightforward to infer. NEE measured at Kobbefjord from 2008 to 2015 indicates a consistent sink of CO 2 (within a range of −17 to −41 g C m −2 ) with exception of the year 2011 (+41 g C m −2 ) ( Table 2). The year 2011, associated with a major pest outbreak, reduced GPP more strongly than R eco (Fig. 5) and Kobbefjord turned into a strong C source within an episodic single growing season. The return to substantial cumulative CO 2 sink rates following the extreme year of 2011 shows the ability of the ecosystem to recover from the disturbance (Lund et al., 2017). Indeed, the ecosystem not only shifted back from being a C source to a C sink, but it also changed rapidly from one year to the next. Thus we found evidence in Kobbefjord of ecosystem resilience to the meteorological variability, similar to other cases described in other northern sites (Peichl et al., 2014;Zona et al., 2014). Only a few reference sites have reported similar decreases in net C uptake, but in no case as large as the one observed here. Zona et al. (2014) described an effect of delayed responses to an unusual warm summer in Alaska. Their results suggested that vascular plants, which have enhanced their physiological activity during the warmer summer, might have difficulties readapting to cooler, but not atypical, conditions, which have provoked a significant decrease in GPP and R eco the following year. In their study, the ecosystem returned to be a fairly strong C sink after 2 years, suggesting strong ecosystem resilience. Moreover, Hanis et al., 2015 have reported comparable C sink-C source variations in a Canadian fen within the growing season due to changes in the water table depth. Drier and warmer than normal conditions have triggered an increase in C source strength. Finally, during an extensive outbreak of autumn and winter moths in a subarctic birch forest in Sweden, Heliasz et al. (2011) observed a similar decrease in net sink of C (most likely due to weaker GPP) across the growing season. However, the C source strength (NEE = 40.7 g C m −2 ) found in 2011 at Kobbefjord was higher compared to these other cases. To our knowledge, such abrupt disturbance concerning C sink strength in Arctic tundra has not been previously reported, excluding severely burned landscapes (Rocha and Shaver, 2011).
A combination of different factors could have led to the sharp change in C balance observed between 2010 and 2011, both physical and biological. The year 2010 had the warmest mean annual temperature (3.4 • C compared to the −0.4 • C mean annual temperature for 2008-2015) and the warmest mean wintertime temperature (−2.7 • C compared to the −6.79 • C mean for 2008-2015) (Fig. 2a). These climatic conditions generated the thinnest (maximum daily snow depth of 0.3 m compared to an average of 0.9 m) ( Table 1) and shortest-lasting snowpack. Consequently, 2010 had the longest growing season (85 days) and very high growing season C uptake (−70 g C m −2 ). Increases in temperature can lead to high respiration rates during early winter (Commane et al., 2017;Zona et al., 2016) but also during the following summer (Helfter et al., 2015;Lund et al., 2012), which is related to soil temperature and snow dynamics. Further, in Kobbefjord the year 2011 had one of the lowest mean annual temperatures and mean wintertime temperatures (−1.7 and −6.1 • C respectively), which created the thickest (maximum daily snow depth of 1.4 m) and the longest-lasting snowpack, leading to the shortest growing season for the study period (only 47 days). According to Lund et al. (2012), soils will be insulated from low temperatures when below thick snowpack, which acts as a lid and prevents R eco from being released to the atmosphere until the snowmelt period. Finally, an outbreak of larvae of the noctuid moth Eurois occulta occurred in 2011, overlapping the observed abrupt decrease in C sink strength. Although we cannot provide a quantification of change attributed to meteorological variations and biological disturbances, there is evidence showing that the moth outbreak could partially have decreased the C sink strength in Kobbefjord. In an undisturbed scenario, the meteorological conditions in 2015, colder and dryer than the mean 2008-2015 period (Fig. 2) but similar to 2011, would have stimulated similar behaviours in terms of C fluxes. However, the cumulative fluxes in 2015 (Fig. 5) followed analogous patterns compared to other years. This evidence agrees with the literature (Callaghan et al., 2012b;Lund et al., 2017) on the fact that tundra systems can fluctuate in sink strength influenced by factors such as episodic disturbances or species shifts, events which are very difficult to predict.

Environmental forcing
Our data indicate that the importance of the main environmental controls (radiation and temperature) for C fluxes did vary across diurnal, seasonal and annual cycles, but also between time aggregations. The hourly variability of NEE and GPP (Figs. 7 and 8) was mostly dependent on PAR because of the threshold nature on radiation control on GPP. Overall, the results indicate that environmental factors that can change rapidly such as PAR will have a high influence on short timescales (Stoy et al., 2014). The increased importance of PAR at 08:00 and 20:00 h WGST coincides with the sharp gradient in light at dawn and dusk (Fig. 8). The control of PAR on GPP is not a new finding itself, but the random forest approach helps to quantify its importance. There is no GPP at night, and therefore there will be a strong increase/decrease in GPP at dawn/dusk. The seasonal pattern also showed that radiation is the single main driver for NEE and GPP between early June and early October, supported by the longer daytime. Further, PAR appeared to be a limiting factor for annual NEE in 2011, increasing further the complexity around this anomalous year. These results agree with the literature (Groendahl et al., 2007;Stoy et al., 2014), suggesting that the uptake of CO 2 is partially controlled by radiation for the photosynthetic physiology at the leaf scale. Arctic plants are usually well adapted to environments with low light levels, reporting near-maximum rates ranging from 10 to 25 • C (Oechel and Billings, 1992;Shaver and Kummerow, 1992).
Photosynthesis is restricted by low temperature, so enzymatically driven processes such as carbon fixation are more sensitive to low temperature than the light-driven biophysical reactions (Chapin et al., 2011). In this paper the daily, weekly and monthly aggregated variability of C fluxes was primarily linked to T air . Moreover, the random forest analyses revealed a strong diurnal pattern with a marked contribution of T air to variations in NEE and GPP (both at night-time and between 08:00 and 18:00 h WGST). These results agree with Lindroth et al. (2007), who recognized T air as the key driver of NEE seasonal trends in northern peatlands. However, in this analysis both NEE and GPP had similar responses to common environmental forcing, contrary to the results in Reichstein et al. (2007). In order to circumvent the potential circularity conflicts based on the use of partitioning products, we filtered daytime NEE (true GPP) and night-time NEE (true R eco ), obtaining very similar results (Table S2). Further, our data also suggest that R eco is often dominated by air temperature. The patterns observed here are in agreement with findings on plant respiration dynamics (Heskel et al., 2016;Lloyd and Taylor, 1994;Tjoelker et al., 2001).
In this study, environmental drivers related to water availability such as VPD and precipitation were not found to be as influential as other assessed variables. We did not find significant relationships between CO 2 fluxes and the water table depth. Thus, there was no apparent water limitation on carbon dynamics during the 8-year period. However, the complex interactions based on changes in temperature and soil moisture particularly over full annual cycles and for sites with permafrost, should be further explored. Our results contrast with Strachan et al. (2015), who described water table depth as an important driver regulating the CO 2 balance, and others, who found that CO 2 emissions increase during dry years due to increased decomposition rates and a reduction in GPP Lund et al., 2007;Oechel et al., 1993;Peichl et al., 2014), whereas other sites act as sinks during relatively wet years (Lafleur et al., 1997). The fen in Kobbefjord is probably quite resistant to droughts since it is fed with water from the surroundings.

Conclusions
We have analysed eight snow-free periods in 8 consecutive years in a West Greenland tundra (64 • N) focusing on the net ecosystem exchange (NEE) of CO 2 and its photosynthetic inputs (GPP) and respiration outputs (R eco ). Here, the NEE gap filling exposed inherent uncertainties in the seasonal C budget calculation, but there were also inconsistencies between the flux partitioning approaches used. We find that Kobbefjord acted as a consistent sink of CO 2 during the years 2008-2015, except 2011, which was associated with a major pest outbreak. The results do not show a marked meteorological effect on the net C uptake. However, the relative insensitivity of NEE during the snow-free period was driven by the correlated, balancing responses of GPP and R eco , both more variable than NEE and sensitive to temperature and insolation. In this paper we show a tendency towards larger GPP and R eco during wetter and warmer years. The anomalous year 2011, affected by a biological disturbance, constituted a relatively strong source of CO 2 and reduced GPP more strongly than R eco . A novel analysis assessing the changes of environmental forcing across diurnal, seasonal and annual timescales unmasked patterns of functional responses to C fluxes.
Despite the fact that we analysed an 8-year data set, the results do not provide a complete picture due to the lack of year-round data . The snow season should be taken into account for a comprehensive understanding of complete C budget (Aurela et al., 2002;Commane et al., 2017;Zona et al., 2016) and the delayed effect of wintertime-based variables such as snow depth or snow cover on the C fluxes. Because some studies have suggested that GPP and R eco increase with observed changes in climate and NEE trends remain unclear , it is challenging to produce strong evidence while the data remains scarce and fragmented. Hence, there is a need for increased efforts in monitoring Arctic ecosystem changes over the full annual cycle Grøndahl et al., 2008). Future work is also required with C flux modelling in order to explore process-based insights of C exchange balance in the Arctic tundra and the interactions of photosynthesis and R eco with changes in C stocks.
Data availability. Measurement data from the Greenland Ecosystem Monitoring (GEM) programme are freely available from the GEM database (http://data.g-e-m.dk). Post-processed data and scripts are available from the authors upon request (elb@bios.au.dk).
Author contributions. EL-B, ML, MW, TRC and MPT designed the experiment. Data preparation and analysis were primarily performed by EL-B with contribution from ML (eddy covariance data processing, data quality control and LRC partitioning), MW and TRC (experimental set up), BUH (data gathering from Nuuk Ecological Research Operations, GeoBasis), AW-N (daily estimate of the timing of snowmelt) and J-FE (random forest approach). EL-B prepared the manuscript with active contributions from all coauthors.
Competing interests. The authors declare that they have no conflict of interest.