A long-term (2002 to 2017) record of closed-path and open-path eddy covariance CO2 net ecosystem exchange fluxes from the Siberian Arctic

Abstract

Ground-based observations of land–atmosphere fluxes are necessary to progressively improve global climate models. Observed data can be used for model evaluation and to develop or tune process models. In arctic permafrost regions, climate–carbon feedbacks are amplified. Therefore, increased efforts to better represent these regions in global climate models have been made in recent years. We present a multi-annual time series of land–atmosphere carbon dioxide fluxes measured in situ with the eddy covariance technique in the Siberian Arctic (72${}^{\circ }$22${}^{\prime }$ N, 126${}^{\circ }$30${}^{\prime }$ E). The site is part of the international network of eddy covariance flux observation stations (FLUXNET; site ID: Ru-Sam). The data set includes consistently processed fluxes based on concentration measurements of closed-path and open-path gas analyzers. With parallel records from both sensor types, we were able to apply a site-specific correction to open-path fluxes. This correction is necessary due to a deterioration of data, caused by heat generated by the electronics of open-path gas analyzers. Parameterizing this correction for subperiods of distinct sensor setups yielded good agreement between open- and closed-path fluxes. We compiled a long-term (2002 to 2017) carbon dioxide flux time series that we additionally gap-filled with a standardized approach. The data set was uploaded to the Pangaea database and can be accessed through 10.1594/PANGAEA.892751.

Introduction

The release of the Arctic's belowground carbon (C) pools to the atmosphere can potentially act as a positive feedback on climate change. Organic material that is now stored in the permanently frozen soil and largely inaccessible for microbial decomposition might become available under a warming climate resulting in an increased release of greenhouse gases from Arctic regions . At the same time, the Arctic vegetation responds to ongoing warming with a greening trend , probably enhancing summer carbon assimilation. Although the importance of permafrost carbon pools for a potential amplification of climate change has been widely recognized e.g.,, the earth system models analyzed for the Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC) did not include permafrost carbon emissions.

While efforts to include permafrost dynamics in global climate models have been made recently e.g.,, models can be improved by using ground-based flux measurements for calibration and validation. assessed the carbon balance of the Arctic tundra combining ground-based observations and process and atmospheric inversion models. The authors found that the uncertainty with which a carbon balance can be quantified is still very large, with upper and lower uncertainty bounds indicating the Arctic tundra as a sink for carbon at one and as a C source at the other bound. conclude that reducing uncertainties in regional estimates based on observational data relies on high-quality ground-based measurements that should be placed strategically, e.g., along hydrological or vegetation gradients. In situ gas flux measurements from the Arctic are, however, still scarce. Moreover, the available data are biased towards Alaska; observations from the Eurasian Arctic are even more scarce . To be able to distinguish climate-change-related flux responses from interannual variability, long-term data sets are essential as recently argued by .

Within the scope of this publication, we aimed at creating a high-quality, long-term ${\mathrm{CO}}_{\mathrm{2}}$ flux data set from a polygonal tundra site in the Russian Arctic. We had the opportunity to analyze a 16-year record of eddy covariance data that includes periods with simultaneous measurements from two different (closed-path and open-path) ${\mathrm{CO}}_{\mathrm{2}}$ gas analyzer types. Our objective was to consistently process the data while following standardized quality control methods to allow for comparability between the different years of our record and with other data sets. We additionally aimed at cross-calibrating open-path and closed-path ${\mathrm{CO}}_{\mathrm{2}}$ fluxes and at gap filling the data set by employing the method of that is widely used in the FLUXNET community.

Site description

The investigation site is located on Samoylov Island in the southern central part of the Lena River delta at 72${}^{\circ }$22${}^{\prime }$ N, 126${}^{\circ }$30${}^{\prime }$ E (see Fig. ). The fan-shaped delta covers an area of roughly 30 000 ${\mathrm{km}}^{\mathrm{2}}$ and is characterized by a network of channels and more than 1500 islands . Being the largest delta in the Arctic and one of the largest worldwide , it lies in the continuous permafrost zone with permafrost depths of about 500 to 600 m . Mean annual permafrost temperatures range around $-\mathrm{9}$ ${}^{\circ }$C at 10 m depth , making the Lena River delta one of the coldest permafrost regions on earth. inferred an annual mean soil temperature of $-\mathrm{8.6}$ ${}^{\circ }$C at 10.7 $\mathrm{m}$ depth from a 2006 to 2011 time series of temperature measurements in a borehole on Samoylov Island. Based on long-term hydrological observations in the delta area, found an increase in discharge as well as in sediment flux indicating recently intensified thawing of ice complex sediments in the region.

Location of Samoylov Island (center of b) in the Lena River delta (a). Map data from OpenStreetMap contributors, under Open Database License.

[Figure omitted. See PDF]

divides the delta area into three main geomorphological units. The oldest, ice-rich river terrace consists of fine-grained sediments with high organic content. It developed as an eroded Pleistocene plane characterized by polygonal ground and thermokarst processes. The second largest unit consists of Late Pleistocene to Early Holocene sandy sediments with low ice content and covers 23 % of the northwestern part . Samoylov Island is part of the third unit, the Mid- to Late Holocene river terrace , which makes up about two-thirds of the delta .

The island itself consists of two morphological units, an annually flooded, modern floodplain (1.49 ${\mathrm{km}}^{\mathrm{2}}$) in the west and a Late Holocene river terrace (2.85 ${\mathrm{km}}^{\mathrm{2}}$) in the east, which lies 10 to 16 m a.s.l. and is not flooded regularly . The data presented here were collected with eddy covariance systems installed on the elevated river terrace. In contrast to the modern floodplain, the river terrace's surface is patterned due to frost action that formed a wet polygonal tundra landscape consisting of mostly low-centered and some high-centered ice-wedge polygons as well as thermokarst lakes and channels. Due to the underlying permafrost and thereby hampered drainage, water-saturated soils or ponds form in the polygon centers, whereas on the rims, which can be elevated by up to 50 cm above the centers, a drier, moderately moist water regime prevails . Accordingly, the vegetation community in the wetter centers is dominated by hydrophytic sedges (Carex aquatilis, Carex chordorrhiza, Carex rariflora) and mosses (e.g., Limprichtia revolvens, Meesia longiseta, Aulacomnium turgidum). Mesophytic dwarf shrubs (e.g., Dryas octopetala, Salix glauca), forbs (e.g., Astragalus frigidus) and mosses (e.g., Hylocomium splendens, Timmia austriaca) dominate on the rims . Maximum summer leaf coverage was estimated by to be 0.3 for vascular plants and 0.95 for mosses and lichens at both polygon centers and rims. The river terrace as a whole is composed of polygon rims with a coverage of 60 % to 65 % and of depressed surfaces (including vegetated and water-filled polygon centers as well as lakes and channels) that cover the remaining 35 % to 40 % of area .

An arctic–continental climate with low mean annual temperatures prevails in the Lena River delta. Although precipitation is low as well, the climate can be considered humid as evaporation rates are low due to low ambient temperatures, and relative humidity is high . Based on long-term (1998 to 2017) in situ measurements on Samoylov Island, inferred an annual mean air temperature of $-\mathrm{12.3}$ ${}^{\circ }$C, the coldest and warmest months being February and July with mean temperatures of $-\mathrm{32.7}$ and 9.5 ${}^{\circ }$C, respectively. For the period from 1998 to 2011, estimated total annual precipitation to be composed of $\mathrm{124}\pm \mathrm{57}$ mm summer rainfall and $\mathrm{65}\pm \mathrm{35}$ mm snowfall. Interannual variability in rainfall was, however, very high, with a maximum of 199 $\mathrm{mm}$ and a minimum of 48 $\mathrm{mm}$. Snowmelt usually starts in mid-May and lasts until early June. Snow accumulation typically commences between late September and early October. Between 1998 and 2011, the snow season lasted on average $\mathrm{224}\pm \mathrm{18}$ days and the snow-free period $\mathrm{138}\pm \mathrm{18}$ days. Snow depth was reported by , averaging 0.3 $\mathrm{m}$ between 2002 and 2017 with a maximum of 0.8 $\mathrm{m}$ in 2017. Beginning in early to mid-June, the soil starts to thaw from the top, forming the so-called active layer. report a mean active layer depth in August of 49 $\mathrm{cm}$ with a maximum of 79 $\mathrm{cm}$ between 1998 and 2011. The closest WMO (World Meteorological Organization) weather station is located on the continent, around 110 $\mathrm{km}$ southeast of Samoylov Island in the city of Tiksi (WMO ID 21824). Between 1936 and 2017 the mean air temperature reported from Tiksi is $-\mathrm{12.74}$ ${}^{\circ }$C, mean annual precipitation amounts to 304.5 $\mathrm{mm}$ . While the mean air temperature in Tiksi is very similar to the 20-year mean from Samoylov Island, average annual precipitation appears to be much higher in Tiksi than in the delta region. explain this divergence with the fact that Tiksi is located on the coast of the Laptev Sea and surrounded by mountains.

The soils of Samoylov Island were classified as Gelisols by based on work by according to the US Soil Taxonomy . At a subgroup level, typical soils of the river terrace are Glacic Aquiturbels, which developed on the polygon rims and are characterized by the translocation of soil material due to freeze–thaw processes (cryoturbation). In the wetter polygon centers Typic Historthels formed. On the more sand-rich active floodplain, Typic Aquorthels and Typic Psammorthels dominate. According to the FAO World Reference Base for Soil Resources , the diverse soils of Samoylov Island belong to the reference soil group of Cryosols. estimated the soil organic carbon (SOC) stocks for the upper meter of the island's two major landscape units to be $\mathrm{29}\pm \mathrm{10}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ for the river terrace and $\mathrm{14}\pm \mathrm{7}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ for the active floodplain.

Eddy covariance (EC) tower positions on the river terrace of Samoylov Island and surface class distribution according to . Photographic image of the entire island (top right corner) from .

[Figure omitted. See PDF]

Methods

Instrumentation

We used the eddy covariance (EC) technique to determine half-hourly gas and energy fluxes. The EC method requires high-frequency (typically $>\mathrm{10}$ $\mathrm{Hz}$) raw gas concentration and three-dimensional wind velocity measurements. A comprehensive description of the EC approach is given, for example, by . We recorded carbon dioxide (${\mathrm{CO}}_{\mathrm{2}}$) and water vapor concentrations as well as three-dimensional wind velocity with changing instrumentation on three different tower structures, all located on the central river terrace of Samoylov Island between 2002 and 2017 (see Fig. ). We deployed open-path (OP) as well as closed-path (CP) gas analyzers, at times simultaneously. Models, manufacturers and years of deployment are given in Table . Between the different setups, CP intake tube lengths varied from 5 to 8 $\mathrm{m}$. OP analyzers were always installed inclined by about 10${}^{\circ }$ from the vertical, as suggested in the analyzer manuals. Raw data were recorded at 20 Hz except for the periods 22 August 2009 to 19 July 2010 (10 $\mathrm{Hz}$) and 31 August 2012 to 17 May 2013 (5 $\mathrm{Hz}$). Until 29 April 2014, all raw data were recorded on a CR3000 data logger (Campbell Scientific, UK). From then on, CP analyzer and anemometer data were logged on a CR3000, whereas OP analyzer and anemometer data were recorded on a LI-7550 data logger (LI-COR Biosciences, USA). Although data coverage is biased towards the growing season, the data set contains considerably more shoulder season and winter fluxes in its second half from 2010 to 2017 (see Table ). The availability of year-round ancillary meteorological data, also increasing, resulted in gap-filled flux time series covering each half hour of the two years 2010 and 2016 (see Fig. ).

List of deployed instrument types. All infrared gas analyzers were manufactured by LI-COR Biosciences (USA), R3 sonic anemometers were built by Gill Instruments Ltd. (UK) and CSAT3 anemometers by Campbell Scientific Ltd. (UK).

	Gas analyzer		Anemometer		Data coverage
Year	Closed path	Open path	Model	Height, m	Date range	Days
2002	LI-7000	n/a	R3	3.65	12 July to 3 September	53
2003	LI-7000	n/a	R3	3.65	19 July to 22 October	95
2004	LI-7000	n/a	R3	3.65	28 May to 20 July	53
2005	LI-7000	n/a	R3	4	17 July to 1 September	46
2006	LI-7000	n/a	R3	4	5 June to 19 September	106
2007	n/a	LI-7500	CSAT3	2.4	11 July to 23 August	36
2008	n/a	LI-7500	CSAT3	2.4	22 April to 26 September	157
2009 I	n/a	LI-7500	CSAT3	2.4	10 April to 14 June	65
2009 II	n/a	LI-7500	CSAT3	4.15	15 July to 29 December	167
2010	LI-7000	LI-7500	CSAT3	4.15	1 January to 31 December	359
2011	LI-7000	LI-7500	CSAT3	4.15	1 January to 22 August	233
2012	n/a	LI-7500	CSAT3	4.15	13 July to 10 November	120
2013	LI-7000	LI-7500A	CSAT3	4.15	4 May to 5 November	185
2014	LI-7000	LI-7500A	CSAT3	4.15	21 February to 29 October	250
2015	LI-7000	LI-7500A	CSAT3	4.15	6 May to 31 December	239
2016	LI-7000	LI-7500A	CSAT3	4.15	1 January to 19 November	323
2017	LI-7000	LI-7500A	CSAT3	4.15	1 January to 30 September	272

Flux processing

Prior considerations

Due to the contrasting designs of OP and CP analyzers, these sensor types have distinct signal response characteristics that we considered during data processing. The most apparent constructional difference between OP and CP gas analyzers is the presence or absence of a housing for the measurement cell that contains the optical path. In a CP instrument, the measurement cell is housed, whereas the optical path of an OP analyzer is exposed to the atmosphere. CP systems are typically more bulky and installed at the base of an EC tower, from where tubing leads to an intake close to the anemometer. Sample air is drawn into the cell with a pump. OP sensors are commonly installed in close proximity to the anemometer and do not require a pump, which greatly reduces the power consumption of OP instruments compared to CP setups. Due to the tubing acting as a low-pass filter, the response to high-frequency concentration variations is systematically attenuated in CP setups as opposed to OP systems . Moreover, the severity of frequency dampening can vary nonlinearly with environmental conditions, especially with relative humidity .

Infrared gas analyzers typically measure gas densities and report the number of molecules per volume of air. To be able to refer the mass of a gas to the mass of air, gas densities are transformed to mixing ratios using air density. However, as the optical path of an OP gas analyzer is exposed to the varying temperature, pressure and humidity conditions of the atmosphere, air density in the measurement cell fluctuates mainly due to thermal expansion/contraction and water dilution/concentration. This effect, which leads to faulty concentration readings of OP instruments and thereby to incorrect flux estimates, has first been described by . The authors proposed two flux correction terms to compensate for these density fluctuation effects that are referred to as Webb–Pearman–Leuning (WPL) terms and have since been verified experimentally and theoretically and are routinely applied in OP EC studies. Especially at times of low gas fluxes, WPL terms can become orders of magnitude larger than raw gas fluxes . CP analyzers have the advantage of controlled temperature and pressure conditions in the measurement cell, allowing for the sample-wise calculation of mixing ratios rather than molar densities and thereby avoiding the need to apply air density fluctuation correction terms after raw flux calculation.

Major drawbacks of OP instruments, especially in harsh environments, are (1) their downtime during adverse weather conditions (e.g., precipitation) and (2) flux biases due to sensor self-heating . The OP self-heating effect was first recognized due to apparent off-season ${\mathrm{CO}}_{\mathrm{2}}$ uptake in flux time series obtained with LI-7500 (LI-COR Biosciences, USA) OP gas analyzers. However, recently found that this effect is not limited to cold conditions but extends throughout all seasons. The necessary corrections can be substantial but decrease greatly when the sensor is not mounted vertically but inclined instead as shown by and .

Processing steps

We performed separate flux processing steps on OP and CP data sets and computed half-hourly fluxes using the software EddyPro (LI-COR Biosciences, USA). An overview of the processing steps is given in Table . We detected and removed raw data spikes according to , with a maximum of 1 % accepted spikes and a maximum of three samples as consecutive outliers. We applied an angle of attack correction, i.e., compensation for flow distortion induced by the anemometer frame , on wind velocity data collected with the R3 (Gill Instruments Ltd., UK) anemometer. The majority of the wind velocity records come, however, from a CSAT3 (Campbell Scientific, UK) instrument, for which this correction is not necessary. Coordinate rotation to align the anemometer $x$ axis to the current mean streamlines was calculated as double rotation according to . For OP fluxes, we compensated for air density fluctuations due to thermal expansion/contraction and water dilution/concentration following . Because simultaneous water vapor concentration, cell temperature and cell pressure measurements from inside the CP analyzer were available, ${\mathrm{CO}}_{\mathrm{2}}$ concentrations from this sensor could be converted directly into mixing ratios, i.e., concentrations referring to dry air of constant temperature , making further corrections for density fluctuations unnecessary. We compensated for CP time lags by using the automatic time lag optimization option in EddyPro. For this procedure, prior to processing the complete data set, time lags were determined for a subperiod of raw data by covariance maximization . A searching window around the median of the time lags found (nominal time lag, $T_{\mathrm{nom}}$) is defined by $T_{\mathrm{nom}}\pm \mathrm{3.5}\times \mathrm{MAD}$, where $\mathrm{MAD}$ is the median absolute deviation of the time lags found. When processing the complete data set, EddyPro performed a covariance maximization of vertical wind velocity and the scalar of interest for each half hour and then checked whether the time lag found fell within the searching window defined before. If not, $T_{\mathrm{nom}}$ was used as the time lag. Water vapor concentration time series were binned in 10 relative humidity (RH) classes, and the procedure was applied to each class, resulting in 10 different nominal time lags. ${\mathrm{CO}}_{\mathrm{2}}$ concentrations were not binned in humidity classes. We computed CP time lag statistics annually and within a year if pump speeds or instrumental setups varied. OP time lags were determined by covariance maximization within a searching window of $-\mathrm{10}$ to 10 s. We evaluated OP time lags statistics, binned in classes of wind direction sectors, later on in the course of quality filtering.

Eddy covariance flux processing steps. Partly differing processing was applied to raw data from closed- and open-path analyzers. OP and CP fluxes were computed consistently for the whole period from 2002 to 2017. Setup-dependent statistics (for time lags and in situ spectral correction methods) were evaluated annually or if tower position, CP pump speed or any other analyzer metadata changed.

Processing step	Method
	Closed-path data	Open-path data
Spike detection	raw data spike removal
and removal
Angle of attack	from 2002 to 2006 during Gill	n/a, sensor was not deployed
correction	anemometer deployment	between 2002 and 2006
Axis rotation	Double rotation
Detrending	linear
Correction for air	sample-wise conversion of raw	application of WPL terms
density fluctuations	data to mixing ratios	to fluxes
Time lag compensation	covariance maximization with	covariance maximization
	nominal time lag from statistics
Spectral corrections for
high-pass filtering	analytic
low-pass filtering	in situ/analytic	analytic
instrument separation		n/a
EddyPro version	$\ge \mathrm{6.0.0}$

Spectral attenuation in the high- and the low-frequency spectral range was compensated for according to the following methods. Low-frequency signal loss due to the finite averaging time used for flux calculations (30 min) and due to linear raw data detrending was corrected for following the method of for both OP and CP fluxes. High-frequency signal loss of OP fluxes due to path and volume averaging of the sonic anemometer and the gas analyzers as well as due to the separation between the two instruments were corrected for with the analytical approach of . High-frequency signal loss of CP fluxes due to spectral attenuation by the intake tube and volume averaging in the measurement cell were corrected for using the in situ method of . For each measurement period with a unique instrumental setup and CP pump speed, we determined the cutoff frequency of a first-order low-pass filter from ensemble means of 30 min power spectra of ${\mathrm{CO}}_{\mathrm{2}}$ concentration and sonic temperature time series data. The spectral correction factor was then parameterized as a function of the cutoff frequency found and the mean wind speed for stable and unstable atmospheric conditions as described by . Before using them for ensemble spectra estimations, the 30 min power spectra were quality-filtered by applying the scheme of and by omitting half hours that were assigned quality class 2 according to . High-frequency noise was removed from the ensemble means of ${\mathrm{CO}}_{\mathrm{2}}$ concentration power spectra before the determination of the cutoff frequency where it was deemed necessary. High-frequency signal losses due to crosswind and vertical separation of the sample air tube intake and the anemometer were corrected for according to .

Quality filtering

We set EddyPro to calculate quality flags according to that represent flux quality in three classes (0, 1 and 2), with 0 denoting the highest- and 2 denoting the lowest-quality class. This quality evaluation is based on tests for stationarity and developed turbulence and thereby indicates whether general EC assumptions about atmospheric conditions were met during a flux calculation period. Flux quality assessment was largely based on the scheme of . In the data set available for download, we included one column for each analyzer type containing this quality flag. Additionally, we applied six further screening steps and flagged fluxes of low quality. If a flagged flux was not already assigned to class 2 according to Mauder and Foken (2004), we set the quality flag to 2. In our opinion, fluxes of quality class 2 should be omitted from further analysis. They are included in the reported data set for the sake of completeness. We performed the six additional flagging steps in the following sequence. An overview of these filtering steps including the number of flagged values is given in Table .

Additional quality flagging steps after flux processing. Flagged fluxes were assigned to quality class 2 if not in this class already according to the quality assessment. As CP time lag detection quality had been addressed earlier during flux processing in EddyPro, it was not screened at this stage.

	Applied to		No. of flagged fluxes
Step	OP fluxes	CP fluxes	OP	CP
1: raw data skewness/kurtosis	yes	yes	23 769 (23 %)	12 043 (18 %)
2: instrument signal strength	yes	no	6951 (7 %)	n/a
3: time lag detection quality	yes	no	20 277 (20 %)	n/a
4: absolute concentration limits	yes	yes	223 (0.2 %)	2261 (3 %)
5: exclusion of outliers when simul-	yes	n/a	346 (0.3 %)	n/a
taneous CP fluxes close to zero
6: absolute flux limits	yes	yes	634 (0.6 %)	102 (0.6 %)

In step 1, skewness and kurtosis were computed with EddyPro for the half-hourly high-frequency raw data time series of ${\mathrm{CO}}_{\mathrm{2}}$ concentration, vertical wind speed and sonic temperature. If any of these statistics was outside certain intervals (skewness: $[-\mathrm{2},\mathrm{2}]$; kurtosis: $[\mathrm{1},\mathrm{8}]$; equivalent to the hard flag defined by ), ${\mathrm{CO}}_{\mathrm{2}}$ flux values were flagged.

In step 2, OP fluxes were additionally filtered for an instrument signal strength indication ($\mathrm{AGC}$) recorded from the LI-7500 sensor. Along with a software upgrade, this diagnostic value was renamed $\mathrm{RSSI}$, and its definition was changed. We therefore recalculated the $\mathrm{AGC}$ values for sensors not running on firmware version 6.6 and above (before July 2013). According to the old $\mathrm{AGC}$ definition in the LI-7500 manual, typical clean window values range between 55 % and 65 %. As dirt accumulates on the windows (or anywhere in the optical path), the $\mathrm{AGC}$ value will increase up to 100 %. The new $\mathrm{RSSI}$ value takes 100 % for clean windows and decreases as windows get dirtier. In order to obtain one consistent diagnostic variable for the cleanness of the optical path, $\mathrm{AGC}$ was converted to the $\mathrm{RSSI}$ range. $\mathrm{AGC}$ values smaller than 44 were set to 44, then AGC values were mapped to the RSSI range as follows. $$\mathrm{RSSI}\left(\mathrm{AGC}\right)=\mathrm{188}-\mathrm{2}\cdot \mathrm{AGC}.$$ We flagged OP ${\mathrm{CO}}_{\mathrm{2}}$ flux values when $\mathrm{RSSI}\le \mathrm{60}$.

As quality control of the half-hourly time lag detection results was not applied during OP flux processing in EddyPro, we additionally screened OP time lags to identify low-quality flux values in step 3. We divided the time lag data set into subsets of different instrumental setup and binned the time lags of these subsets in 36 10${}^{\circ }$ wind direction sectors. We used the 25th and 75th percentiles per class as filter thresholds. We flagged OP flux values with associated time lags outside the range spanned by these thresholds. Because we computed CP fluxes in EddyPro considering and compensating for low time lag detection quality, we did not perform this type of filtering step on CP fluxes.

In step 4, we flagged CP as well as OP fluxes when 30 min average concentration measurements were larger than 450 ppm or smaller than 300 ppm. ${\mathrm{CO}}_{\mathrm{2}}$ concentrations outside this range indicate dirty OP gas analyzer optics or technical problems of the CP air sampling system (sudden pump speed changes due to brownouts, blocked filters, etc.).

To filter dubious, large OP fluxes that coincided with reasonable CP fluxes, we selected all OP fluxes when simultaneously measured CP values ranged between $-\mathrm{2}$ and 2 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$. Step 5 only affected OP data from this subset. We calculated the 99th and 1st percentile of this group and flagged fluxes from it when they lay outside this percentile range.

In step 6, we flagged remaining outliers in both the CP and OP data sets by using the 0.1st and 99.9th percentile ($-\mathrm{3.5423}$ and 3.3473 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) of the CP time series after the concentration limits filter as absolute limits to define an acceptable range of OP and CP flux values.

Open-path self-heating correction

To account for self-heating errors induced by the LI-7500 sensor electronics, we corrected OP fluxes as described by . The authors use WPL-corrected fluxes and add a correction term that accounts for self-heating effects of vertically installed instruments. In their approach, use a scaling factor $\mathit{\xi }$, taking values between 0 and 1, to trim the correction for inclined analyzer setups. With simultaneously available CP fluxes, we were able to estimate this scaling factor specifically for our site and periods of unique instrumental setups. As suggested by , we optimized this parameter with a nonlinear least squares method in Matlab (v. 9.2). We determined $\mathit{\xi }$ for periods of different instrumental setups and separately for night (incoming shortwave radiation $<\mathrm{20}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) and day (incoming shortwave radiation $\ge \mathrm{20}$ ${\mathrm{Wm}}^{-\mathrm{2}}$) conditions using the following equation:

$$F_c=F_{c,\mathrm{WPL}}+\mathit{\xi }{\displaystyle \frac{(T_{\mathrm{s}}-T_{\mathrm{a}}){\mathit{\rho }}_c}{r_{\mathrm{a}}{T}_{\mathrm{a}}}},$$ where $F_c$ ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is the true ${\mathrm{CO}}_{\mathrm{2}}$ flux, $F_{c,\mathrm{WPL}}$ ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is the WPL-corrected OP ${\mathrm{CO}}_{\mathrm{2}}$ flux, $T_{\mathrm{s}}$ ($\mathrm{K}$) is the instrument surface temperature, $T_{\mathrm{a}}$ ($\mathrm{K}$) the ambient air temperature, $r_{\mathrm{a}}$ ($\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$) the aerodynamic resistance and ${\mathit{\rho }}_c$ ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) the ambient ${\mathrm{CO}}_{\mathrm{2}}$ density. Prior to $\mathit{\xi }$ optimization, we also estimated the instrument surface temperature $T_{\mathrm{s}}$ following the parameterization of separately for nighttime and daytime:

$$\begin{array}{cc}{\displaystyle }& {\displaystyle }{T}_{\mathrm{s},\mathrm{day}}=\mathrm{0.93}(T_{\mathrm{a}}-T_{\mathrm{0}})+\mathrm{3.17}+T_{\mathrm{0}}\mathrm{and}\\ & {\displaystyle }& {\displaystyle }{T}_{\mathrm{s},\mathrm{night}}=\mathrm{1.05}(T_{\mathrm{a}}-T_{\mathrm{0}})+\mathrm{1.52}+T_{\mathrm{0}},\end{array}$$ with $T_{\mathrm{s},\mathrm{day}}$ ($\mathrm{K}$) and $T_{\mathrm{s},\mathrm{night}}$ ($\mathrm{K}$) as instrument surface temperature estimates and $T_{\mathrm{0}}$ set to 273.15 $\mathrm{K}$. We determined the scaling factor as a parameter of Eq. (2), being the modified approach from . For function fitting, we assumed CP fluxes of quality classes 0 and 1 as true fluxes. We used WPL-corrected OP quality class 0 fluxes and the surface temperature estimates described above as independent variables. Before parameter optimization, we quality-screened the correction term (expression to the right-hand side of $\mathit{\xi }$ in Eq. 2) and removed spikes ranging within the uppermost or lowest percent of its distribution. Throughout all years, $\mathit{\xi }$ is larger during the daytime than at nighttime but generally small, adding mostly below 1 % of the full correction term to the uncorrected flux (see Table ). In four of the seven available years with simultaneous CP and OP fluxes, nighttime $\mathit{\xi }$ optimization converged to values below zero. Before applying the correction models to these periods, we set nighttime $\mathit{\xi }$ estimates to the median of the years yielding parameter values that, including their 95 % confidence bounds, ranged above zero. We used this value and the median of all daytime model optimizations to calculate corrected OP fluxes at times without parallel CP measurements. We did not correct OP fluxes when radiation measurements or correction term estimates were not available. Correlation between CP and OP fluxes improved throughout all quality classes by applying the self-heating correction (see Table ), while fluxes indicating net ${\mathrm{CO}}_{\mathrm{2}}$ uptake were affected more strongly than fluxes above zero (see Fig. ).

Effect of the self-heating correction on the correlation between open-path (OP) and closed-path (CP) fluxes (a). Correlations were quantified using Spearman's rank correlation coefficient $r_{\mathrm{s}}$ and Pearson's correlation coefficient $r$. Only quality class 0 is shown. Negative fluxes are affected more strongly by the correction than positive fluxes (b).

[Figure omitted. See PDF]

Estimates of scaling factor $\mathit{\xi }\pm \mathrm{95}$ % confidence intervals used for open-path flux correction. $\mathit{\xi }$ describes the portion of the self-heating correction term, given by for vertically installed instruments, that is needed to correct OP fluxes determined with inclined gas analyzers. The scaling factor was optimized as a parameter of a nonlinear function where CP data were regarded as true fluxes. It was therefore determined for years when parallel CP and OP measurements were available. In the case of an optimization converging to unreasonable values (below 0), we used the median of the remaining $\mathit{\xi }$ estimates.

Year	Daytime $\mathit{\xi }$	Nighttime $\mathit{\xi }$
2010	$\mathrm{0.0076}\pm \mathrm{0.0012}$	$\mathrm{0.0071}\pm \mathrm{0.0013}$
2011	$\mathrm{0.0116}\pm \mathrm{0.0009}$	$\mathrm{0.0068}\pm \mathrm{0.0015}$
2013	$\mathrm{0.0150}\pm \mathrm{0.0007}$	$\mathrm{0.0104}\pm \mathrm{0.0009}$
2014	$\mathrm{0.0094}\pm \mathrm{0.0006}$	0.0071
2015	$\mathrm{0.0050}\pm \mathrm{0.0010}$	0.0071
2016	$\mathrm{0.0051}\pm \mathrm{0.0005}$	0.0071
2017	$\mathrm{0.0069}\pm \mathrm{0.0005}$	0.0071

Spearman's rank correlation coefficient $r_{\mathrm{s}}$ and Pearson's correlation coefficient $r$ between closed-path (CP) and open-path (OP) fluxes with and without the applied self-heating correction. The agreement between CP and OP fluxes increases throughout all quality classes after OP correction.

		Quality class 0	Quality classes 0, 1	Quality classes 0, 1, 2
$r_{\mathrm{s}}$	OP uncorrected	0.896	0.866	0.508
	OP corrected	0.907	0.871	0.512
$r$	OP uncorrected	0.894	0.871	0.042
	OP corrected	0.904	0.877	0.055

Carbon dioxide flux gap filling

We used the CP and the corrected OP fluxes (see Fig. ) to compile a ${\mathrm{CO}}_{\mathrm{2}}$ flux time series. We aimed at keeping as many measured data points as possible, while omitting records with large uncertainty. We accepted all CP values of quality classes 0 and 1. At time steps where no CP fluxes were available, we selected OP values of the same quality classes. The resulting time series contains 75 921 data points. Additionally, we filled the remaining gaps in the time series using the marginal distribution sampling (MDS) method as first presented by . This method employs two types of model value calculations. The environmental variables global radiation, air temperature and water vapor pressure deficit are binned in classes and combined in a lookup table (LUT). In the case of a gap, flux values related to similar environmental conditions can be looked up and used for averaging and gap filling. The setup of different LUTs for fixed time periods was first described by . This process can be refined by the use of moving time windows around gaps, as applied by . The second model type implemented in the MDS algorithm exploits the commonly high autocorrelation of gas flux time series. The mean diurnal variation (MDV) technique was also first described by and uses the average of available gas flux measurements from adjacent days at the same hour of day to fill a flux gap. The MDS method has found wide application, as it has, for example, been the standard technique within the processing pipeline of the FLUXNET2015 data set, which includes over 1500 site years of data. The algorithm of combines a screening procedure of the available data for similar environmental conditions (lookup table steps) and the use of an MDV method (diurnal cycle steps) if a gap could not be filled within the lookup table steps. Both techniques include moving windows with variable sizes that are increased until a solution can be found. Large gaps are skipped. To run the gap-filling algorithm, we used the REddyProc routine that is accessible through a web-based service hosted by the Department of Biogeochemical Integration at Max Planck Institute Jena. The $R$ routine that is executed on this server is a further-developed and extended version of the approach and is described by . We did not use the friction velocity filter or the flux partitioning capabilities of the REddyProc online tool. Gap filling resulted in 131 908 data points. The provided data set includes quality flags for each gap-filled value that depend on the method used and time window size, as defined by . These flags take values between 0 and 3, with 0 denoting measurement data, 1 indicating the most reliable and three least reliable gap-filled fluxes. To assess the overall quality of the gap-filling result, the MDS algorithm, in a stepwise manner, treats single available values as gaps and fills them according to the described scheme. Pearson's correlation coefficient between our compiled ${\mathrm{CO}}_{\mathrm{2}}$ flux time series and the MDS quality assessment run, where these values were treated as artificial gaps, is 0.92, with a root mean squared error of 0.31 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$.

Multi-annual carbon dioxide flux time series compiled from fluxes measured with closed-path and open-path sensors on Samoylov Island's river terrace. Fluxes of quality class 2 are not shown. Self-heating errors in the OP data set have been corrected for. Additionally, the result from gap filling this time series with the MDS method is shown. The number of values given for the gap-filled time series include measured fluxes.

[Figure omitted. See PDF]

Flux uncertainty estimation

Flux uncertainty can be regarded as a combination of a systematic and a random part. While the attempt should be made to remove systematic biases, random errors cannot be corrected for . However, statistical methods exist to estimate the uncertainty in a flux measurement due to random errors. We used three different approaches from the literature to quantify random uncertainty and addressed fluxes with a suspected large bias by correcting for it during processing or by filtering in the course of quality assessment.

Normalized mean contributions of the surface classes defined by to the eddy covariance footprint. Values were averaged over each subperiod and normalized to sum up to 1. Additionally, the average non-normalized sum of all surface class contributions is given as the column “median image contribution”. These values indicate how sufficient the classified area is to describe the EC footprint. Non-normalized half-hourly contributions of the single classes are given in the provided data set.

Year	Tundra		Water		Median image
	Dry	Wet	Overgrown	Open	contribution
2002	0.71	0.17	0.07	0.05	0.88
2003	0.70	0.17	0.07	0.05	0.87
2004	0.71	0.16	0.07	0.06	0.88
2005	0.71	0.17	0.07	0.05	0.87
2006	0.70	0.17	0.07	0.06	0.86
2007	0.54	0.37	0.06	0.02	0.73
2008	0.53	0.34	0.09	0.04	0.77
2009 I	0.54	0.32	0.08	0.06	0.72
2009 II	0.64	0.19	0.09	0.08	0.71
2010	0.65	0.18	0.09	0.08	0.73
2011	0.67	0.18	0.08	0.07	0.79
2012	0.67	0.18	0.08	0.07	0.80
2013	0.69	0.17	0.08	0.06	0.83
2014	0.66	0.18	0.08	0.07	0.77
2015	0.66	0.18	0.08	0.08	0.78
2016	0.65	0.18	0.09	0.08	0.74
2017	0.67	0.18	0.08	0.07	0.82

Description of columns included in the data set file.

Column name	Unit/format	Description
Date/time (Local)	yyyy-mm-ddTHH:MM	Time stamp referring to end of 30 min flux calculation period in local time (UTC $+\mathrm{9}$ h).
Date/time (UTC)	yyyy-mm-ddTHH:MM	Time stamp referring to end of 30 min flux calculation period in UTC.
CP ${\mathrm{CO}}_{\mathrm{2}}$ flux	$\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$	Closed-path ${\mathrm{CO}}_{\mathrm{2}}$ flux
QC CP ${\mathrm{CO}}_{\mathrm{2}}$ flux	dimensionless	Closed-path ${\mathrm{CO}}_{\mathrm{2}}$ flux quality classes 0, 1 and 2
CP ${\mathrm{CO}}_{\mathrm{2}}$ flux rand unc	$\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$	Closed-path ${\mathrm{CO}}_{\mathrm{2}}$ flux random uncertainty estimate
OP ${\mathrm{CO}}_{\mathrm{2}}$ flux	$\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$	Open-path ${\mathrm{CO}}_{\mathrm{2}}$ flux
OP corr ${\mathrm{CO}}_{\mathrm{2}}$ flux	$\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$	Corrected open-path ${\mathrm{CO}}_{\mathrm{2}}$ flux
QC OP ${\mathrm{CO}}_{\mathrm{2}}$ flux	dimensionless	Open-path ${\mathrm{CO}}_{\mathrm{2}}$ flux quality classes 0, 1 and 2
OP ${\mathrm{CO}}_{\mathrm{2}}$ flux rand unc	$\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$	Open-path ${\mathrm{CO}}_{\mathrm{2}}$ flux random uncertainty estimate
${\mathrm{CO}}_{\mathrm{2}}$ flux comp	$\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$	Time series compiled of open- and closed-path quality class 0 and 1 fluxes
${\mathrm{CO}}_{\mathrm{2}}$ flux gf	$\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$	Gap-filled ${\mathrm{CO}}_{\mathrm{2}}$ flux time series
QC ${\mathrm{CO}}_{\mathrm{2}}$ flux gf	dimensionless	Quality flag of gap-filled fluxes, between 0 and 3
${\mathrm{CO}}_{\mathrm{2}}$ flux gf std	$\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$	Standard deviation of gap-filled flux estimates, calculated from the data used for averaging
FP CC dry	dimensionless	Contribution of surface class “dry tundra” to the eddy covariance footprint
FP CC wet	dimensionless	Contribution of surface class “wet tundra” to the eddy covariance footprint
FP CC ove	dimensionless	Contribution of surface class “overgrown water” to the eddy covariance footprint
FP CC wat	dimensionless	Contribution of surface class “open water” to the eddy covariance footprint

Most importantly, systematic errors are introduced when underlying EC assumptions are not met. Using the method of that combines an assessment of well-developed turbulence and steady-state conditions, we identified biased fluxes and flagged them. Other sources of systematic errors that we addressed include, for example, the angle of attack correction of faulty sonic anemometer readings, filtering for low instrument signal strength, the OP self-heating correction, and compensations for high-frequency loss and air density fluctuations (see Sect. 3.2.2, 3.3 and 3.4). Although we are confident that we applied corrections for systematic errors both rigorously and carefully enough, biases were certainly not always removed efficiently. The quality flags included in the data set, reflect a level of confidence based on the assessment of general EC assumptions and our six additional quality filtering steps (see Sect. 3.3).

To be able to include a random uncertainty estimate for each individual OP and CP flux in the provided data set, we set EddyPro to calculate random uncertainty estimates following . The authors developed a method that aims at quantifying flux uncertainty associated with turbulence sampling errors. These errors can contribute largely to the total random error as they refer to the insufficient sampling of large eddies with high spectral energy. Due to the stochastic nature of turbulence, this type of error is random. To estimate its magnitude, the so-called integral turbulence timescale (ITS) is first determined by expressing the covariance of vertical wind velocity and gas concentration as a function of a lag time between these two time series. The ITS is then given by integrating the cross-correlation function theoretically from 0 to infinity, in practice, however, until an upper lag time limit is reached. The upper limit can be defined in three different ways in EddyPro. We used the definition of the normalized cross-correlation function reaching a value of $\mathrm{1}/e=\mathrm{0.369}$ to determine an upper lag time limit used for integration. While the normalized cross correlation should reach zero with increasing lag time in theory, in practice it sometimes does not. The setting we used on the one hand provides the least conservative estimate of the ITS but on the other hand offers computational efficiency and makes sure that an upper limit for integration can be reliably found. With the ITS, a flux uncertainty can be determined by calculating the variance of an EC flux or, as put it, by calculating the variance of the covariance. This ensemble variance would approach zero with the averaging time approaching infinity. In the data set available for download, a random uncertainty estimate calculated with the method of is given for each OP and CP flux (see Table ). Random uncertainties based on ITS estimation observations increase with absolute fluxes with mean values of 0.16 and 0.05 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for OP and CP fluxes (see Fig. ). OP random uncertainty estimates are generally larger and more scattered with respect to the corresponding flux values.

Random uncertainty estimates for all closed-path (CP) and open-path (OP) ${\mathrm{CO}}_{\mathrm{2}}$ fluxes calculated using estimates of the integral turbulence timescale (ITS), the successive-observations approach and results from gap filling (GF), and the paired-observations approach during periods with simultaneous OP and CP records.

[Figure omitted. See PDF]

As the random uncertainty estimate described above specifically addresses the turbulence sampling error, other sources of random flux errors such as the noise introduced by the different components of the measurement system are neglected. With simultaneous measurements from two sensors, we could additionally estimate random errors for the measurement system as a whole during times when the data sets from both sensors overlapped. We followed the paired-observations approach as presented by and calculated a random error estimate $\mathit{ϵ}$ as $$\mathit{ϵ}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}\cdot (F_{\mathrm{CP}}-F_{\mathrm{OP}}),$$ with the closed-path and open-path ${\mathrm{CO}}_{\mathrm{2}}$ fluxes $F_{\mathrm{CP}}$ and $F_{\mathrm{OP}}$ of quality classes 0 and 1 in $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$. The distribution of $\mathit{ϵ}$ estimates is shown in Fig. . The $\mathit{ϵ}$ values calculated with OP fluxes corrected for the self-heating error have a mean close to zero and are distributed more symmetrically than the $\mathit{ϵ}$ values calculated with uncorrected OP fluxes. The mean of this distribution is shifted from its mode as well as from zero, indicating a much stronger systematic component within the measurement error. This result increases our confidence that the OP self-heating correction we applied was successful in removing a systematic bias from the data. Further following , we used the $\mathit{ϵ}$ system error data set from the overlap period to generate flux uncertainty estimates for bins of increasing OP flux ranges. We sorted the $\mathit{ϵ}$ values into 20 corresponding flux bins between $-\mathrm{2}$ and 2 $\mathrm{\mu }$mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ and calculated an uncertainty estimate for each bin $\mathit{\sigma }(\mathit{ϵ}{)}_i$ as $$\mathit{\sigma }(\mathit{ϵ}{)}_i=\sqrt{\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{N_j}}\sum _{j=\mathrm{0}}^{N_j}|{\mathit{ϵ}}_{i,j}-{\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_i|.$$ Results show (see Fig. ) a similar data range and pattern of uncertainty estimates in relation to associated fluxes like the half-hourly values calculated following .

As a third method of random uncertainty estimation, we simplified the successive-observations approach from by using results of the quality run performed during MDS gap filling (see Sect. 3.5). We selected the time steps when a flux observation and an MDS value that was estimated using a 1-day window and the MDV technique were available. We used the standard deviation of the fluxes measured at the same hour of day within a 1-day window, as an uncertainty estimate of the observed flux. Results are shown in Fig.  and also increase with rising absolute fluxes in the same ranges as random uncertainties due to turbulence sampling error or measurement system error do.

We included the results obtained with ITS estimation in the uploaded data set considering the similarity between the uncertainty–flux relations calculated with independent methods as well as due to the advantage of a distinct uncertainty estimate for each sensor and time step.

Distributions of the measurement system errors $\mathit{ϵ}$ estimated using the paired-observations approach for differences between closed path and corrected (a) as well as uncorrected (b) open-path (OP) fluxes.

[Figure omitted. See PDF]

Footprint modeling

In order to quantify the cumulative contribution of distinct surface classes to the EC source area, we evaluated the two-dimensional analytical footprint formulation described by in combination with a 0.14 $\mathrm{m}$ $\times \mathrm{0.14}$ $\mathrm{m}$ resolution surface classification of Samoylov Island's central river terrace provided by . The authors divide the surface into four classes based on hydrology and vegetation communities, as illustrated in Fig. . presented an analytical solution to the crosswind-distributed advection–diffusion equation described by and . Using the analytical model of , the authors solved the power-law profiles of horizontal wind speed and eddy diffusivity by relating them to the Monin–Obukhov similarity theory, including the stability dependence of the exponents in the power laws at a certain height. We implemented the equations given in as a Matlab (v. 9.2) function and added a quality filter, omitting calculations when friction velocity was larger than 0.9 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ or smaller than 0 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, wind speed was below zero or above 20 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, the crosswind standard deviation was below zero or above 3 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, or Monin–Obukhov length was smaller than ${\mathrm{10}}^{-\mathrm{3}}$ $\mathrm{m}$ or larger than ${\mathrm{10}}^{\mathrm{4}}$ $\mathrm{m}$. Prior to half-hourly footprint calculations, we additionally determined roughness length statistics for annual subsets of data and binned them in 2 ${}^{\circ }$ wind direction classes. The medians of these classes were used in the subsequent half-hourly footprint estimation, depending on the mean wind direction during these 30 min. We evaluated the footprint model at the same resolution that was used by to classify the surface (i.e., 0.14 $\mathrm{m}$ $\times \mathrm{0.14}$ $\mathrm{m}$). We could thereafter assign a probability of being the EC source area to each classified pixel and sum up the probabilities of all pixels belonging to the same surface class to estimate the contribution of each class. This process of combining an EC source area estimation with a land cover classification is similar to what has been applied and described in more detail by .

Mean surface class composition of the eddy covariance footprint during 17 subperiods of four different tower setups at three locations on Samoylov Island.

[Figure omitted. See PDF]

Discussion

Although we did our best to ensure the consistency and appropriateness of the data processing workflow for the presented net ecosystem exchange (NEE) time series, due to technical and logistical constraints during 16 years of field work, disparities in the experimental setup exist which may challenge its integrity. The EC tower was relocated twice, and the measurement height was changed three times (see Fig.  and Table ). These changes of tower location and measurement height affected the source area and hence the surface types sampled during flux measurements. Most notably, between July 2007 and June 2009, the EC tower was placed about 650 $\mathrm{m}$ southwest of its original position at the center of Samoylov Island, in an area with an increased coverage of the surface class “wet tundra”. This is revealed by the footprint analysis (Fig. ). While the EC footprint is dominated by the surface class “dry tundra” throughout the time series, during subperiods 2007, 2008 and 2009 I the contributions of “wet tundra” to the measured flux are significantly higher. To check the effect of the shifts in tower location and measurement height on cumulative ${\mathrm{CO}}_{\mathrm{2}}$-C fluxes, we calculated flux sums for a period when flux time series without gaps were available in most years. The overlapping period covers days of year 200 to 234, i.e., part of the growing season in all years except for 2004 (see Fig. ). Interannual variability of cumulative C fluxes in years with constant tower location (and measurement height) appears to be large and driven by a more complex set of variables than shifts in surface class contributions only. Flux sums from the periods when EC tower relocation led to a significant shift in EC footprint composition are well within the range of the distribution of cumulated fluxes from years with a more homogeneous EC fetch area. We therefore assume that, at least with respect to budget calculations, the presented long-term time series is not disrupted and can be regarded as representative of a polygonal tundra site dominated by “dry tundra”. For a more in depth analysis of flux dynamics, footprint information should and can be considered by users of the data set. Recently, a comparison between surface class level NEE models based on chamber measurements with EC fluxes, using the half-hourly footprint information provided in this data set for scaling, yielded good agreement between the results obtained with both methods . We regard the availability of half-hourly footprint information in the presented NEE data set an attribute that sets it apart from other studies and holds possibilities for comprehensive analyses.

Apart from the changes in anemometer height, other deviations of the general instrument setup occurred due to limitations in data storage during two winter periods when the acquisition frequency was reduced to 5 and 10 $\mathrm{Hz}$. demonstrated in a field experiment that fluxes calculated from raw data recorded at frequencies below 20 $\mathrm{Hz}$ compare well with fluxes derived from high-frequency raw data. Differences arise as an increase in random noise and not as a systematic bias. High-frequency noise removal before ensemble spectra estimation in EddyPro is effective in limiting the effect of increased noise on the quality of transfer function estimation in the process of spectral correction. Overall spectral correction in EddyPro is expressed as a spectral correction factor (SCF) which comprises the effect of all applied compensations for high- and low-frequency loss. Raw fluxes are multiplied with the respective SCFs during processing. We compared the SCF distributions of the two abovementioned winter periods with statistics of the remaining parts of the time series when data were recorded at 20 $\mathrm{Hz}$. SCF deviations between the different acquisition frequencies are minor (see Fig. ), which implies that systematic differences between fluxes calculated from raw data of different temporal resolutions are in fact small; random uncertainties increase, however.

Comparison of cumulative ${\mathrm{CO}}_{\mathrm{2}}$ flux sums of different years during the same day-of-year range.

[Figure omitted. See PDF]

Spectral correction factor statistics for periods with different acquisition frequencies.

[Figure omitted. See PDF]

Scientific overview

While results on methane exchange fluxes and the soils' methane production and oxidation potential are more prominent in the publication record e.g.,, the literature on ${\mathrm{CO}}_{\mathrm{2}}$ flux time series recorded with the same measurement system presented in this publication is available for distinct years. Flux processing has, however, been streamlined only now. The length of the time series, the addition of detailed footprint information, the site-specific correction of OP fluxes, and the coherent processing and quality filtering distinguishes the data set at hand from past publications like the contribution made to the FLUXNET2015 data set .

Ongoing analysis of the long-term data set (Kutzbach, unpublished) inter alia confirms what has been found in the past . The polygonal tundra of Samoylov Island appears to be a robust growing season ${\mathrm{CO}}_{\mathrm{2}}$-C sink, whereas this sink strength can vary so much interannually that prolonged low-level respiratory ${\mathrm{CO}}_{\mathrm{2}}$-C loss during the cold season can offset ${\mathrm{CO}}_{\mathrm{2}}$-C uptake during the vegetation period. Reduced summer uptake has been observed for both the coldest and warmest summers. found that with frequent early season heat spells, the temperature-induced increase in respiratory release can exceed the rise in photosynthetic uptake. Recently, all data from this publication have been contributed to the Arctic Data Center's chamber and EC synthesis project “Reconciling historical and contemporary trends in terrestrial carbon exchange of the northern permafrost-zone” that aims at identifying seasonal and interannual C flux dynamics and its drivers based on a newly established pan-Arctic database.

In context with the improvement of earth system models (ESMs), carbon dioxide fluxes from Samoylov Island can be especially of use due to the site's comparably high moss cover. Using data from Samoylov, found that current ESMs miss an observed early season ${\mathrm{CO}}_{\mathrm{2}}$ uptake peak suspected to be connected to the earlier onset of moss photosynthesis in comparison with vascular plants. Although there have been advances and, e.g., developed a dynamic moss model for JSBACH , noted that the simulated ${\mathrm{CO}}_{\mathrm{2}}$ uptake and release terms combining vascular vegetation and moss carbon fluxes did not agree with observational data. The fact that the Samoylov Island NEE data set has now been extended and its quality has been greatly improved holds the opportunity to estimate the performance of updated ESM versions that are set up to represent carbon fluxes in the moss layer better.

The data set was uploaded to the Pangaea database and can be accessed through 10.1594/PANGAEA.892751. The included columns are given in Table . Ancillary long-term time series of meteorological and soil variables from Samoylov Island are available from and can be accessed through 10.1594/PANGAEA.891142.

Conclusions

We are confident that the presented carbon dioxide land–atmosphere flux data set is of high quality and is likely to be of value to the scientific community. We screened the data carefully and applied filtering rules to identify erroneous data, taking into account sensor diagnostics, time lag statistics and the presence of atmospheric conditions that allow for a robust application of the EC method. We followed standardized processing and quality control/assurance routines to allow for comparability between different years from our site as well as with flux time series from other tundra environments. With OP measurements being paralleled by CP measurements in 7 years, we had the opportunity to correct for self-heating errors in our OP measurements with a site-specifically scaled correction term rather than using default correction methods e. g.. We could therefore address different sensor setups with different correction terms and thereby improve our OP data set, as the self-heating effect has distinct impacts on sensors installed at different inclinations. We quantified the contribution of certain soil and vegetation community types to each half-hourly EC footprint, taking into account varying roughness lengths throughout different years and wind direction sectors. We estimated the cumulative probability of being the EC source area for the four main surface classes on Samoylov Island's river terrace by using a land cover classification and by computing an analytical EC footprint model. Multi-annual results show (see Table ) that on average the combination of different surface classes within the EC footprint is representative of the surface composition of the whole river terrace that developed as a polygonal tundra landscape. According to , the river terrace is composed of 65 % “dry tundra”, 19 % “wet tundra” and 16 % ponds (sum of “open water” and “overgrown”). On average, the surface class compositions within the EC footprint are very similar to these values. Deviations arise, however, in the years between 2007 and 2009, when the tower location was shifted from the center towards the southwestern cliff of Samoylov Island. Nevertheless, the contributions of each surface class to the EC footprint are not only available on average, as presented in Table , but half-hourly in the uploaded data set, ensuring that EC source area deviations are quantifiable by a potential user. A total of 16 years of consistently processed and quality-controlled carbon dioxide fluxes from a polygonal tundra landscape typical of Arctic lowlands are a valuable addition to the already existing database of ${\mathrm{CO}}_{\mathrm{2}}$ net ecosystem exchange observations from the Arctic, especially because of the site's location in Northern Siberia, from where only limited data are available up to now. Furthermore, analysis of this NEE time series is not limited to the gas flux data only. An extensive data stream of meteorological and soil variables between 2002 and 2017 has recently been published by . The authors made their records publicly accessible on the two long-term repositories Pangaea (10.1594/PANGAEA.891142) and Zenodo (https://zenodo.org/record/2223709, last access: 1 February 2019). The fact of ancillary ecosystem variables available in parallel enables a potential user to put the gas flux dynamics reported in this publication into context with the variability of other ecosystem properties and potential flux drivers. We regard this type of analysis as vital to understanding interannual variability of gas fluxes and are working on it ourselves (Kutzbach, unpublished).

DB, JB, MG, IF, LK and EP conceptualized and administered the research activity planning and execution and acquired funds for it. LB, DH, LK, ML, BR, PS, TS and CW conducted the investigation. DH and CW analyzed the data; DH created the visualizations. DH wrote the original draft; DH, LK, BR, TS and CW reviewed and edited the original draft.

The authors declare that they have no conflict of interest.

Acknowledgements

Without the dedicated work of many scientists, logistics experts and engineers over the years, we would not have been able to present this long-term eddy covariance NEE data set. We want to thank Niko Bornemann, Tim Eckhardt, Mauel Helbig, Lars Heling, Oliver Kaufmann, Zoé Rehder, Norman Rößger, Norman Rüggen, Günter Stoof and Waldemar Schneider for their commitment, diligence and ingenuity. We thank Jakob Sievers for providing us with a starting point for the Matlab implementation of the footprint model and Norman Rößger for sharing his analysis of the long-term meteorological data from Tiksi with us. This work was supported through the Cluster of Excellence CliSAP (EXC177), Universität Hamburg, funded through the German Science Foundation (DFG), by the European Commission through the project PAGE21 (FP7-ENV-2011, 282700), and by the German Ministry of Education and Research (BMBF) through the projects CarboPerm (03G0836A) and KoPf (03F0764A). Edited by: David Carlson Reviewed by: two anonymous referees