Scaling and balancing carbon dioxide fluxes in a heterogeneous tundra ecosystem of the Lena River Delta

The current assessments of the carbon turnover in the Arctic tundra are subject to large uncertainties. This problem can (inter alia) be ascribed to both the general shortage of flux data from the vast and sparsely inhabited Arctic region, as well as the typically high spatiotemporal variability of carbon fluxes in tundra ecosystems. Addressing these challenges, carbon dioxide fluxes on an active flood plain situated in the Siberian Lena River Delta were studied during two growing seasons with the eddy covariance method. The footprint exhibited a heterogeneous surface, which generated mixed flux signals that could be partitioned in such a way that both respiratory loss and photosynthetic gain were obtained for each of two vegetation classes. This downscaling of the observed fluxes revealed a differing seasonality in the net uptake of bushes (−0.89 µmol m⁻² s⁻¹) and sedges (−0.38 µmol m⁻² s⁻¹) in 2014. That discrepancy, which was concealed in the net signal, resulted from a comparatively warm spring in conjunction with an early snowmelt and a varying canopy structure. Thus, the representativeness of footprints may adversely be affected in response to prolonged unusual weather conditions. In 2015, when air temperatures on average corresponded to climatological means, both vegetation-class-specific flux rates were of similar magnitude (−0.69 µmol m⁻² s⁻¹). A comprehensive set of measures (e.g. phenocam) corroborated the reliability of the partitioned fluxes and hence confirmed the utility of flux decomposition for enhanced flux data analysis. This scrutiny encompassed insights into both the phenological dynamic of individual vegetation classes and their respective functional flux to flux driver relationships with the aid of ecophysiologically interpretable parameters. For comparison with other sites, the decomposed fluxes were employed in a vegetation class area-weighted upscaling that was based on a classified high-resolution orthomosaic of the flood plain. In this way, robust budgets that take the heterogeneous surface characteristics into account were estimated. In relation to the average sink strength of various Arctic flux sites, the flood plain constitutes a distinctly stronger carbon dioxide sink. Roughly 42 % of this net uptake, however, was on average offset by methane emissions lowering the sink strength for greenhouse gases. With growing concern about rising greenhouse gas emissions in high-latitude regions, providing robust carbon budgets from tundra ecosystems is critical in view of accelerating permafrost thaw, which can impact the global climate for centuries.

1 Introduction

Permafrost underlies between 12.8 % and 17.8 % of the land area in the Northern Hemisphere (Zhang et al., 2000). Large parts of this area coincide with the Arctic tundra, which is situated north of the boreal treeline and covers roughly 8 % of the global land surface (McGuire et al., 2012). As a consequence of the long-term carbon sink function, the underlying permafrost forms a carbon stock of global relevance: approximately 1300 Gt of soil organic carbon are stored in the circumpolar permafrost region (Hugelius et al., 2014). However, large fractions of this carbon pool may be remobilised in response to a warming climate, making the tundra a key ecosystem for climate change (Schuur et al., 2008).

The Arctic north of 60^∘ N latitude has warmed at a rate of 1.36 ^∘C per century since 1875, i.e. roughly twice as fast as the global average (Masson-Delmotte et al., 2013). And this rapid warming trend is projected to continue (Collins et al., 2013). Due to ambiguous model results and their large confidence intervals, it currently remains unclear whether the permafrost areas maintain their sink function or convert into a carbon source in the future (Heimann and Reichstein, 2008; Schuur et al., 2015). These uncertainties do not only arise from the limited knowledge of the physical thawing rates, the fraction of released carbon after thawing and the timescales of release but also from the general shortage of flux data in Arctic ecosystems (Ciais et al., 2013). The scarce data availability particularly applies to the extensive Siberian tundra, which covers around 3 million km², i.e. more than half of northern high-latitude tundra ecosystems (Chapin et al., 2005; Sachs et al., 2010). The low density of flux observation sites is due to both harsh environmental conditions as well as challenging logistics in these remote and sparsely inhabited areas that are often without line power. Consequently, current estimates of the tundra's sink strength for carbon dioxide are associated with large uncertainties: $-\mathrm{103}\pm \mathrm{193}$ Tg C yr⁻¹ (McGuire et al., 2012). The same issue applies to estimates that indicate a shift to a source for carbon dioxide: 462±378 Tg C yr⁻¹ (Belshe et al., 2013). The refinement of these macroscale estimates and the reduction of their uncertainties can be achieved via providing both more flux budgets (in particular from the Siberian tundra) and more reliable information on the variation in habitats (e.g. bogs, fens,) plus their associated surface heterogeneities (e.g. tussocks, hummocks).

Tundra ecosystems are frequently characterised by a pronounced vegetation patchiness with sharply defined boundaries between different vegetation classes (Shaver et al., 2007). Besides vegetation, surface classifications can also be based on differences in soil moisture, snow cover, permafrost features or combinations of them (Fox et al., 2008; Virkkala et al., 2018). The consequently high spatial variability in carbon fluxes complicates the estimation of robust carbon budgets that are accurate and precise. The omission of accounting for the spatial distribution of different surface types is likely to lead to incorrect budgets (Oechel et al., 1998). Therefore, an improved understanding of the effects of surface heterogeneity on these budgets, e.g. through a better characterisation of both spatial flux variability as well as associated key factors such as vegetation composition and structure, is necessary (Kade et al., 2012; Kwon et al., 2006). For quantifying vegetation properties, NDVI (normalised difference vegetation index), LAI (leaf area index) and foliar nitrogen content have been found suitable (Marushchak et al., 2013; Shaver et al., 2013). All of these predictors can be estimated by remote sensing, thereby neglecting the patchy nature of tundra ecosystems, but also offering the potential for macroscale modelling of both carbon dioxide budgets plus their prospective alterations through climate change. Such assessments are based on biome-level monitoring of the global warming-induced impacts on Arctic vegetation such as growing season prolongation as well as expansion of plant's growing range and size, e.g. the enhanced growth of shrubs and their northward migration into typical graminoid tundra ecosystems (Jia et al., 2009; Sweet et al., 2015). On the other side, microscale observations are crucial in order to reflect the individual biogeochemical dynamics in the mosaic of vegetation patches. The direct appraisal of the vegetation's responses to global warming through field surveys on the plant community level involves aspects such as enhanced primary productivity, deeper rooting depths, as well as augmented ground shading and snow accumulation trough taller growth forms (Myers-Smith et al., 2011; Sitch et al., 2007).

Chamber measurements operate on the microscale (10⁻²–10² m²) and form the common approach to differentiate the carbon dioxide exchange of multiple microforms (i.e. individual land cover components such as bare soil, water bodies, vascular plants, etc.) with the atmosphere (McGuire et al., 2012). Despite their widespread application, however, chamber measurements are associated with several drawbacks such as (i) a disturbance of the studied system due to collar insertion, (ii) a subjectivity in the selection of chamber locations, which is particularly momentous, if an upscaling of fluxes is intended, (iii) a lacking acquisition of a pronounced temporal flux variability on account of a usual confinement to discrete sampling, (iv) a limited spatial representativeness due to both the small sampled size and only a few replicate sites as a result of a high labour intensity, and (v) a decoupling of the sampled surface from the atmosphere that causes a modification of the environmental conditions in the headspace (Fox et al., 2008; Kade et al., 2012; Kutzbach et al., 2007a; Livingston et al., 2006; Riederer et al., 2014; Sachs et al., 2008; Wagner and Reicosky, 1992).

Alternatively, the non-intrusive as well as directly and continuously measuring eddy covariance technique circumvents most of these downsides. This method, however, operates on the mesoscale (10⁴–10⁶ m²) and yields turbulent fluxes that integrate across multiple microforms (Aubinet et al., 2012). The size and location of the sampled surface constantly shifts according to wind direction, wind speed, atmospheric stability, crosswind velocity and surface roughness (Detto et al., 2006). In the presence of a heterogeneous landscape (i.e. an irregular pattern of individual land cover components), the temporal variability in the observed fluxes is not only a result of the varying uptake/release rates of the individual microforms but also an outcome of the varying fractions of microforms in the sampled area. In addition, the footprint budgets may lack representativeness since the fractional composition of microforms within the footprint is likely to deviate from the microform distribution in the area of interest. In such an environment, budgets strongly depend on tower location, sensor height and wind field conditions, and are thus likely to exhibit a sensor location bias (Schmid and Lloyd, 1999). Moreover, heterogeneous flux signals also complicate the determination of model parameters, e.g. the light response curve of a specific vegetation type, if the microforms in the footprint exhibit strongly deviating characteristics (Lasslop et al., 2010). Despite these challenges in signal interpretation, a heterogeneous surface also provides the opportunity to both conduct a concurrent sampling of multiple microforms and study their respective carbon dioxide fluxes utilising only one eddy covariance instrumentation (Forbrich et al., 2011; Morin et al., 2017). Exploiting this potentially valuable information source requires the partitioning of the integrated flux into its microform-specific fluxes. Such a successful flux decomposition routine benefits from the advantages of eddy covariance measurements on the mesoscale whilst resolving the pronounced variability on the microscale. These decomposed fluxes in turn enable, in conjunction with a precise determination of the microforms' spatial coverages in the area of interest, the estimation of robust carbon dioxide budgets for a heterogeneous surface.

Addressing the problems of balancing carbon fluxes in a Siberian tundra ecosystem with both a heterogeneous surface and an unknown greenhouse gas sink/source strength, the objectives of this study are as follows: (i) elucidating the heterogeneity of the landscape with geospatial data, (ii) analysing the spatiotemporal variability of carbon dioxide fluxes with the eddy covariance technique, (iii) explaining this variability with a model-based approach, (iv) estimating robust carbon dioxide budgets that account for the heterogeneity of the area, and (v) combining these budgets with previously estimated methane budgets in order to determine the sink/source strength for greenhouse gases.

2 Material and methodology

2.1 Site description

The Lena River Delta, one of the largest deltas in the world, is located within the zone of continuous permafrost in northern Siberia (Fig. 1). One of its numerous islands is Samoylov Island (72^∘22^′ N, 126^∘28^′ E), which covers an area of 4.8 km² and features two geomorphological units: the Late Holocene river terrace in the eastern part and the active flood plain in the western part. The carbon dioxide exchange on the river terrace, which is characterised by ice-wedge polygonal tundra with sedges and mosses, has been repeatedly studied (Eckhardt et al., 2019; Kutzbach et al., 2007b; Runkle et al., 2013). In contrast to the river terrace, the flood plain has to date received scarce attention in terms of greenhouse gas fluxes, although active flood plain levels represent roughly 40 % of the soil-covered area of the Lena River Delta (Zubrzycki et al., 2013). Aside from the period of the annual spring flood, whose associated inundation is very variable in magnitude and duration, the flood plain on Samoylov Island stretches over an area between 1 km² (spring) and 2 km² (autumn). More importantly, the surface of the flood plain exhibits, in opposition to the river terrace, a distinct heterogeneity on the mesoscale.

Figure 1 Location of the Lena River Delta in northern Siberia indicated by the square. The dots point out sites that were utilised for the pan-Arctic comparison of carbon dioxide budgets (Table 3). The classification of the Arctic zones was based on vegetation occurrence (modified from AMAP, 1998). Accordingly, the treeline delimits the (terrestrial) Arctic; i.e. it corresponds with the boundary between sub-Arctic and low Arctic.

The central delta region is situated in a continental Arctic climate, which is characterised by very low temperatures and a low annual precipitation. In the distant town of Tiksi, located around 120 km southeast of Samoylov Island, a mean annual air temperature of −12.8 ^∘C was measured during 1936–2016 and a mean annual precipitation of 329 mm was gauged during 1956–2016 (AARI, 2017). Additional information on this study site can be found in Rößger et al. (2019a).

2.2 Experimental setup and data recording

An eddy covariance system was installed in the southern part of the flood plain, and the measurements covered two periods: 18 June to 2 October 2014 (107 d) and 9 June to 24 September 2015 (108 d).

The flux tower was equipped with a sonic anemometer (CSAT3, Campbell Scientific, UK) and a gas analyser for water vapour and carbon dioxide (LI-7500A, LI-COR Biosciences, USA). Both instruments were mounted at a height of 2.83 m and sampled with a frequency of 20 Hz. In addition, another eddy covariance system with the same instrumentation has already been deployed at a central position on the adjacent river terrace (Holl et al., 2019).

Supplementary measurements on the flood plain involved gathering data of both air temperature (HMP45, Campbell Scientific, UK) and photosynthetic photon flux density (SKP215, Skye Instruments, UK). These environmental variables were recorded on a logger (CR1000, Campbell Scientific, UK) in quarter-hourly intervals. Furthermore, a time-lapse camera (TLC200, Brinno, Taiwan) was mounted on the flux tower, pointing northeast, for monitoring the phenology during spring 2014 with the same interval of a quarter of an hour.

An acquisition of data on the footprint's surface structure was intended by employing ground-based LAI measurements, as this quantity is widely applied for characterising plant canopies. However, on account of the low height of lichens, mosses and sedges, the measurements with an upwards-pointing sensor (LAI-2200C, LI-COR Biosciences, USA) yielded poor quality results. Therefore, the measurements were terminated after a few test surveys.

2.3 Flux processing

The flux computation was carried out with the software EddyPro version 6.0.0 (LI-COR Biosciences, 2016) for 30 min flux intervals and followed the standard procedure. Detailed information on the executed (i) raw data processing (spike removal, coordinate rotation, block averaging, time lag compensation), the implemented (ii) flux correction scheme (density correction, spectral correction in low- and high-frequency ranges, flux error estimation), and the conducted (iii) quality assessment routine (stationarity test, integral turbulence characteristics test, skewness and kurtosis examination, energy flux quality verification, signal strength control, percentile removal) is provided in Rößger et al. (2019a).

For the footprint modelling, an analytical model for non-neutral stratification was employed (Kormann and Meixner, 2001). This model is based on a stationary gradient diffusion formulation with height-independent crosswind dispersion (Leclerc and Foken, 2014). When applying the solution of the resulting two-dimensional advection–diffusion equation for solving the power law profiles of both eddy diffusivity and mean wind velocity, it yielded a source weight function for each flux interval.

2.4 Surface structure

For studying the impact of the heterogeneous surface on the flux variability, the entire flood plain was mapped in August 2014 by employing helicopter-based visible aerial imagery. The resulting geo-referenced orthomosaic exhibited a resolution of 8.5 cm and hence provided very detailed spatial information, which was sufficient to resolve the pronounced heterogeneity of the surface. Based on maximum likelihood classification tools, the vegetation was classified employing a supervised classification routine on the orthomosaic (Fig. 2). In this process, four different land cover classes were utilised, two of which represent the vegetation.

Figure 2 Vegetation map of the flood plain on Samoylov Island obtained through supervised classification of a high-resolution orthomosaic. The flux tower was situated in the centre of the footprint isolines, which indicate the averaged area from which 10 %–90 % (increments of 10 %) of the flux originated during the measurement periods of 2014 and 2015 (footprint climatology). The small inset illustrates Samoylov Island being composed of flood plain (grey) and river terrace (white) plus the location of their respective flux towers.

Vegetation class 1 (“shrubs”) refers to sites, which were densely vegetated by large dwarf shrubs of the willow family such as Salix pulchra, Salix lanata, Salix hastata, Salix glauca, growing to a maximum height of around 1 m. This shrubby vegetation was located on a sandy ridge aligned in the north–south axis. The elevated area enabled, in conjunction with a spatially averaged maximum thaw depth of 0.93 m, a good drainage. Since the groundwater table remained at depths around 50 cm, the surface was mostly dry, forming favourable growing conditions for willow shrubs and a sparse cover of thin moss.

Vegetation class 2 (“sedges”) represents areas, which were dominated by sedges including Carex aquatilis, Carex chordorrhiza, Carex concolor, as well as species of Eriophorum and Equisetum. Also, small willow shrubs growing to a height of about 0.3 m were occasionally to be found. This predominantly graminoid vegetation was located in depressions around the dry ridge and exhibited a mean active layer depth of 0.69 m. Accordingly, the soil moisture conditions alternated between moist surfaces and wet patches with water levels up to 40 cm. The ample moisture attracted many mosses forming a dense cover of thick moss.

The two other classes, which do not occur in the 90 % contribution footprint around the flux tower, denote a large area of bare sand along the waterfront and some small water bodies mainly situated in the northern part of the flood plain (Fig. 2). The former class was not considered in the budget estimation since its carbon dioxide flux rates were (in comparison to the two vegetation classes) assumed negligible. The latter class was appended to vegetation class 2 as the few small water bodies surrounded by sedges were presumed to have similar flux rates. Further information on the classification routine is given in Rößger et al. (2019a).

2.5 Flux modelling

The model structure was based on the computation of the two components of the carbon dioxide flux.

$$\begin{array}{}\text{(1)}& F_{{\mathrm{CO}}_{\mathrm{2}}}=\mathrm{NEE}=\mathrm{TER}-\mathrm{GPP}\end{array}$$

$F_{{\mathrm{CO}}_{\mathrm{2}}}$ is the net carbon dioxide flux observed at the flux tower and is equal to the net ecosystem exchange (NEE). Its two components describe, respectively, the total ecosystem respiration (TER) and the gross primary productivity (GPP), both of which can be modelled simultaneously (Runkle et al., 2013).

$$\begin{array}{}\text{(2)}& \mathrm{NEE}=R_{\mathrm{base}}\cdot Q_{\mathrm{10}}^{\left(\frac{T_{\mathrm{air}}-T_{\mathrm{ref}}}{\mathit{\gamma }}\right)}-{\displaystyle \frac{P_{max}\cdot \mathit{\alpha }\cdot \mathrm{PPFD}}{P_{max}+\mathit{\alpha }\cdot \mathrm{PPFD}}}\end{array}$$

R_base denotes the basal respiration at the reference temperature (T_ref), which was set to 15 ^∘C, and a scaling factor (γ) was held constant at 10 ^∘C (Mahecha et al., 2010). Q₁₀ indicates the temperature sensitivity; i.e. this parameter is a value by which respiration multiplies/divides, when the temperature rises/drops by 10 ^∘C. P_max refers to the maximum photosynthetic potential and quantifies the theoretical maximum of photosynthesis at infinite irradiance. α represents the light sensitivity and states, as the initial quantum efficiency, the slope of the light response curve at irradiance of zero. All of the four (physiologically interpretable) parameters are best-fit parameters, which were estimated via non-linear ordinary least-squares regression utilising both air temperature (T_air) and photosynthetic photon flux density (PPFD) as explanatory variables. In order to take the heterogeneous surface structure into account, footprint information was included, forming the final model employed for carbon dioxide flux modelling.

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{NEE}=& \sum _{i=\mathrm{1}}^{i=\mathrm{2}}{\mathrm{\Omega }}_i\\ & \cdot \left(R_{\mathrm{base},i}\cdot Q_{\mathrm{10},i}^{\left(\frac{T_{\mathrm{air}}-T_{\mathrm{ref}}}{\mathit{\gamma }}\right)}-{\displaystyle \frac{P_{max,i}\cdot {\mathit{\alpha }}_i\cdot \mathrm{PPFD}}{P_{max,i}+{\mathit{\alpha }}_i\cdot \mathrm{PPFD}}}\right)\end{array}\end{array}$$

Another explanatory variable is the relative contribution of each vegetation class to the flux (Ω) that weights the two computed vegetation-specific fluxes. This variable was obtained through (i) computing the source weight function for a flux interval, (ii) spatially discretising this continuous function with a resolution of 1 m, (iii) adjusting the vegetation map to a resolution of 1 m, (iv) assigning each value of the source weight function to its spatially corresponding vegetation class, and (v) summing the values in each vegetation class.

The fitting procedure, in which only half-hourly quality-controlled flux data were employed, required the estimation of a large number of fitting parameters: R_base, Q₁₀, P_max, α for each vegetation class (Fig. 3). In order to avoid overparameterisation and equifinality problems, the model structure was gradually simplified along four different steps. These alterations in the parameterisation enabled the desired estimation of (i) reasonable seasonal courses of the fitting parameters, i.e. courses that displayed a predominantly smooth evolution with elevated values during the growing season and low values in the shoulder seasons, and (ii) meaningful and significant values for the fitting parameters, i.e. values that were within an acceptable range and their 95 % confidence interval did not overlap zero. Achieving both objectives provided the possibility to interpret the fitting parameters ecophysiologically.

Figure 3 Schematic overview of the model calibration along four steps, in which different parameterisations were applied to obtain meaningful estimates for the eight fitting parameters (R_base, Q₁₀, P_max, α for each vegetation class). The values in the boxes (e.g. 3-2-p model) denote the number of parameters to be fitted for each vegetation class.

In each of the four parameterisation steps, the respectively parameterised model was recalibrated for every day, applying a moving window with fixed/flexible window sizes and a step size of 1 d. In the initial step (step 1), which served the computation of representative Q₁₀ values, all of the eight fitting parameters were estimated in the model (4-4-p). Through its output, which encompassed eight best-fit time series, a representative Q₁₀ value was obtained for each vegetation class by determining the median out of the best-fit Q₁₀ values that fulfilled two requirements: statistical significance and an associated coefficient of determination (between observed NEE and modelled NEE) of 0.75 or greater. These two Q₁₀ values were held constant in the further fitting procedure. In the subsequent step (step 2), the simplified model (3-3-p) was run with six fitting parameters to be estimated, and a Gaussian bell curve was fitted to the time series of significant best-fit α values for each vegetation class. By adding/subtracting 30 % of the function values to/from these two replacement functions, a pair of encompassing threshold functions was respectively appended. These intervals around the replacement functions formed a range inside which best-fit α values were accepted. In the following step (step 3), the model (3-3-p) was run with the same parameterisation of the previous step. The model output was checked for α values inside the acceptable interval as well as significant values for R_base, P_max and α. If these criteria were satisfied, the accordingly modelled NEE was approved, and the fitting procedure proceeded to the next day. Alternatively, several models (3-2-p/2-3-p/2-2-p) were run employing α value(s) from the replacement function(s) for one or both vegetation classes, depending on which vegetation class insignificant and/or implausible best-fit parameters were in. The output was tested again, and in case of significant best-fit parameters, the modelled NEE was accepted, and the fitting procedure continued with the next day. Lastly, a replacement function for P_max was created by fitting a Gaussian bell curve to the time series of significant best-fit P_max values in both vegetation classes. In the final step (step 4), two greatly simplified models (2-1-p/1-2-p) were run with only three fitting parameters to be estimated, as well as a P_max value from the replacement function. Here, if not before, all fitting parameters have taken on meaningful and significant values, which ensured the computation of reliable NEE values. In addition to this brief explanation, a detailed description of the entire fitting process is attached in the Appendix.

Since the model was designed to simultaneously compute the component fluxes in both vegetation classes, it provided the capability for the decomposition of the observed fluxes into their separate flux contributions by the two vegetation classes. The reliability of this downscaling, however, was dependent on the restrictive acceptance of meaningful and significant values for the fitting parameters. The temporal integration of these partitioned fluxes and the subsequent projection of the resulting budgets on their corresponding areas on the flood plain formed the upscaling. The summation of both vegetation class budgets finally yielded a robust budget of the entire flood plain, which was designated as the area of interest. This budget, as opposed to the directly estimated footprint budget, did not exhibit a sensor location bias and hence allowed an unbiased appraisal of both the interannual variability and the sink/source strength. For the sake of comparability of the budgets between the years, the carbon dioxide budgets were calculated for the comparison period of 18 June to 24 September, where data were available in both years. The utilised methane budgets, being necessary for estimating greenhouse gas sink/source strengths, have been obtained for the same period and in a similar fashion to the carbon dioxide fluxes, i.e. an initial downscaling of the observed net fluxes and a subsequent upscaling of the decomposed fluxes on the flood plain (Rößger et al., 2019a).

Figure 4 Annual course of air temperature on Samoylov Island for the years 2014 and 2015, as well as the recent 20-year baseline (Boike et al., 2013, 2019). In each box plot, the central mark denotes the monthly median, and the bottom and top edges indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points excluding outliers. During the warm season, when flux data were available (June to September), 2014 was mostly warmer than 2015.

3 Results

3.1 Meteorological conditions

The mean air temperatures during the measurement periods in 2014 and 2015 amounted to 7.7 and 7.1 ^∘C, respectively. Furthermore, the respective precipitation sums totalled 92.3 and 130.4 mm. The assessment of these values was based on their comparison with long-term averages that were obtained for Samoylov Island between 1998 and 2018 (Fig. 4). The measurement period in 2014 was on average distinctively warmer and slightly drier, while the measurement period in 2015 featured the same mean temperature as the baseline but considerably more rain. The largest differences in air temperature between both measurement periods occurred in spring. Accordingly, the snowmelt in 2014 took place in a prolonged manner during mid-May already, whereas in 2015, the snowmelt was completed within a couple of days in early June, as usual.

3.2 Dynamics of observed fluxes

The carbon dioxide fluxes exhibited both a diurnal and a seasonal course with the following mean fluxes that were obtained by averaging half-hourly flux data for the subseasons in both 2014 and 2015 (Fig. 5). Between the snowmelt and the vegetative phase, the mean carbon dioxide fluxes remained slightly positive, indicating a prevalent respiration, while the vegetation largely remained dormant (0.26 µmol m⁻² s⁻¹). With the onset of the growing season in late June, stalks and foliage began to develop, and the uptake of carbon dioxide during daytime outweighed the release of carbon dioxide during nighttime (−1.06 µmol m⁻² s⁻¹). The intensity of this oscillation increased towards the onset of the reproduction phase in mid-July, where flowers and seeds developed. During this phase, the most negative fluxes occurred featuring a relatively constant magnitude (−1.77 µmol m⁻² s⁻¹). With the onset of the ripening phase in early August, bushes and sedges attained full maturity, and the flux amplitude of the diurnal cycle began to be progressively attenuated (−0.78 µmol m⁻² s⁻¹). During the nights of this period, the most positive fluxes occurred. Towards late August, the respiration exceeded photosynthesis again, indicating the onset of the senescence phase, which was associated with both leaf colouration and leaf drop (0.39 µmol m⁻² s⁻¹). After the end of the growing season in early September, when abscission was completed, the dominance of respiration continued to grow, leading to more positive mean carbon dioxide fluxes (0.55 µmol m⁻² s⁻¹).

Figure 5 Time series of observed carbon dioxide fluxes (after conducting the quality assessment) and modelled fluxes. During the growing season, which is indicated by an elevated variability between late June and early September, the daytime uptake directly followed the diurnal cycle of PPFD, while the nighttime release was dependent on air temperature.

3.3 Model calibration and performance

While the Q₁₀ values were optimised at constant values of 1.42 for vegetation class 1 and 1.48 for vegetation class 2, the other fitting parameters (R_base, P_max and α) displayed a seasonal course for each vegetation class in 2014 and 2015 (Figs. 6 and 7). The temporal evolution of α values could be well approximated with replacement functions, whose application reduced the noise not only in the seasonal courses of α but also in the seasonal courses of both R_base and P_max. In contrast to the replacement functions of α, which were created for both vegetation classes in both years, a replacement function for P_max was created only for vegetation class 1 in 2015 and for vegetation class 2 in 2014.

Figure 6 Time series of fitting parameters in 2014 for vegetation class 1 (index 1 and green confidence intervals) and vegetation class 2 (index 2 and yellow confidence intervals). The circles represent acceptable fits, while the respective reasons for reparameterisation such as insignificant values or values out of valid range are indicated by plus signs and squares. The triangles denote the fitting parameter(s), which caused a refit in the corresponding vegetation class.

Figure 7 Time series of fitting parameters in 2015 with the same symbols and colours as utilised in the previous figure.

The slightly simplified 3-3-p model, which was run at the start of step 3, yielded meaningful and significant values for the fitting parameters in 49 % of the modelled days including 2014 and 2015 (Fig. 3). In the further course of step 3, these goals were achieved by the gradually simplified 3-2-p/2-3-p/2-2-p models in 47 % of the modelled days. During the remaining 4 %, the greatly simplified 2-1-p/1-2-p models of step 4 were deployed. While the 3-3-p model was mainly employed during the summer season, the 3-2-p/2-3-p/2-2-p models were applied throughout the measurement periods with a focus on the shoulder seasons. The 2-1-p/1-2-p models were solely deployed during the shoulder seasons and more often during spring than during autumn. Hence, larger fluxes during the growing season could be more easily modelled in comparison to the remaining time, when lower fluxes associated with a less favourable signal-to-noise ratio prevailed.

R_base was the fitting parameter that could be estimated most confidently, as this parameter accounted for only 19 % of the insignificant values obtained during the fitting procedure in 2014 and 2015. While P_max caused 31 % of the insignificance values, α appeared to be the least certain fitting parameter representing the remaining 50 %. Furthermore, P_max featured most of the significant differences in best-fit values between both vegetation classes, i.e. the confidence intervals of vegetation classes 1 and 2 rarely overlapped, whereas best-fit α values exhibited the fewest significant differences.

On account of both the recalibration for each day as well as the coinciding variabilities of explanatory variables and the explained variable, the model was able to reproduce the observed fluxes very well (Fig. 5). This performance was expressed by a coefficient of determination (R²) of 0.88 for 2014 and 0.95 for 2015. Furthermore, the mean absolute error (MAE) amounted to 0.49 and 0.35 µmol m⁻² s⁻¹ for 2014 and 2015, respectively, while the root mean square error (RMSE) amounted to 0.75 and 0.52 µmol m⁻² s⁻¹. During the summer season, the model performed better in comparison to the shoulder seasons, where autumn in turn displayed a slightly better performance than spring.

3.4 Downscaling and upscaling of fluxes

The estimation of vegetation-class-specific parameter sets allowed the decomposition of the observed net fluxes. This downscaling yielded fluxes of NEE plus their component fluxes (TER and GPP) accounting for both vegetation classes in both years (Fig. 8). For the comparison period in 2014, the mean NEE amounted to −0.89 and −0.38 µmol m⁻² s⁻¹ for vegetation class 1 and vegetation class 2, respectively, and for the comparison period in 2015, −0.71 and −0.69 µmol m⁻² s⁻¹ (Table 1). In contrast to the similar mean net uptakes in 2015, the mean net uptakes in 2014 distinctly differed from each other. This discrepancy originated from the first half of the growing season (mid-June to early August), when the net uptake of vegetation class 1 was considerably larger relative to vegetation class 2. During the second half of the growing season (early August to late September), both net uptakes were rather similar again. Furthermore, the differences in the net uptakes between both years were governed by changes in GPP rather than in TER. In vegetation class 1, NEE in 2014 was only slightly greater in comparison to 2015, which can be attributed to a greater TER and a distinctly greater GPP. And in vegetation class 2, NEE in 2014 was smaller compared to 2015, which can be ascribed to a smaller TER and a clearly smaller GPP.

Figure 8 Time series of decomposed fluxes with 95 % confidence intervals for both vegetation classes. The width of the confidence intervals varied depending on both the flux magnitude and the number of fitting parameters in the chosen model. The decomposition revealed a distinct difference in the net uptake between both vegetation classes during the first half of the growing season in 2014, while the flux dynamics of both vegetation classes were rather similar during the remaining time and in 2015.

Table 1 Outcome of downscaling and upscaling carbon dioxide fluxes for the comparison period (18 June to 24 September) in both 2014 and 2015. The mean downscaled fluxes (± standard deviation) refer to the individual fluxes of both vegetation classes (Fig. 8). Their aggregation yielded cumulative fluxes, whose projection on their corresponding areas on the flood plain in turn returned upscaled fluxes (± combination of cumulative flux error and classification error). The uncertainty metrics derived from the comparison of the applied vegetation map, which was obtained through supervised classification of aerial imagery, with another vegetation map that was acquired via ground-based surveys (Rößger et al., 2019a).

The aggregation of the decomposed fluxes over the comparison period yielded individual budgets, whose multiplication with the corresponding fractional coverages on the flood plain formed the upscaling (Table 1). The subsequent summation of both vegetation-class-specific net uptakes returned the net uptake of the entire flood plain for the comparison period: $-\mathrm{4.42}\pm \mathrm{0.49}$ Mmol in 2014 and $-\mathrm{6.17}\pm \mathrm{0.66}$ Mmol in 2015. The stated uncertainties were obtained by means of standard error propagation techniques including both cumulative flux error and classification error, where the former was an order of magnitude smaller than the latter. Dividing these budgets by the total area of the flood plain yielded mean flood plain budgets of $-\mathrm{4.22}\pm \mathrm{0.47}$ and $-\mathrm{5.89}\pm \mathrm{0.63}$ mol m⁻² (Table 2). These budgets consider the surface heterogeneity; i.e. they are corrected for the sensor location bias, plus they contain an areal reference and thus enable an appropriate comparison with other sites.

Table 2 Comparison of the sink/source strengths between flood plain and river terrace for the comparison periods in 2014 and 2015 (Holl et al., 2019; Rößger et al., 2019a). Accounting for methane's radiative efficiency as a potent greenhouse gas, the methane budgets were converted to carbon dioxide equivalents with a factor of 34, which corresponds to methane's global warming potential based on a time horizon of 100 years including climate carbon feedbacks (Myhre et al., 2013). The flood plain budgets are given for each vegetation class and for the total area. These budgets are the result of a scaling procedure, which included fairly large classification errors that caused distinctly greater uncertainties in comparison to the river terrace budgets, which derived from a representative footprint and hence did not undergo any scaling processes. In comparison to the flood plain, the polygonal tundra on the river terrace took up less carbon dioxide per square metre, but also released less methane, resulting in a similar (2014) and weaker (2015) sink strength for greenhouse gases.

3.5 Greenhouse gas balances

The evaluation of the flood plain's sink/source strength for greenhouse gases required the corresponding methane emission budgets and their conversion to carbon dioxide equivalents (Rößger et al., 2019a). Despite methane's minor percentage of roughly 3 % in the entire greenhouse gas exchange (specified in molar units), its carbon dioxide equivalents diminished the greenhouse gas sink strength (given by the carbon dioxide net uptake) by half in 2014 and by one-third in 2015. Accordingly, the greenhouse gas balances specify that the flood plain formed a moderate sink of $-\mathrm{2.21}\pm \mathrm{0.61}$ mol CO₂ eq. m⁻² and a stronger sink of $-\mathrm{3.81}\pm \mathrm{0.74}$ mol CO₂ eq. m⁻² during the warm season in 2014 and 2015, respectively (Table 2). The lower sink strength in 2014 was a result of a reduced carbon dioxide net uptake rather than an augmented methane efflux. And this reduced carbon dioxide net uptake in turn was caused by a lowered net uptake in vegetation class 2 that effectively counteracted the elevated early season net uptake in vegetation class 1. This class constituted a stronger greenhouse gas sink than vegetation class 2 in both years, which is mainly due to the fact that methane emissions were only present in vegetation class 2. Since these emissions hardly changed between the years along with the negligible methane release in vegetation class 1, the interannual variability in the greenhouse gas sink strength was governed by the carbon dioxide net uptake.

These balances are the first greenhouse gas budgets of a flood plain in the Lena River Delta. Based on these budgets, the sink strength of the adjacent river terrace, where another eddy covariance system has been in operation for many years, could finally be put in context within the domain of the Lena River Delta (Table 2). In 2014 and 2015, the flood plain sequestered per square metre roughly 20 % and 60 % more carbon dioxide, respectively, but it also emitted approximately 70 % more methane. Hence, the flood plain constituted a sink for greenhouse gases that resembled (2014) or was 1.5 times (2015) the sink strength of the polygonal tundra on the river terrace.

4 Discussion

4.1 Assessment of the flux decomposition model

The partitioning of carbon dioxide fluxes was conducted during the Arctic summer, when fully dark conditions during the nights were absent. Consequently, a partitioning approach that is based on fitting parameters to nighttime respiration followed by extrapolating these fits to daytime, and subsequently subtracting the estimated daytime respiration flux from the observed net flux to obtain the photosynthesis flux is confronted with elevated uncertainties (Reichstein et al., 2005). The partitioning approach of the present study avoids this problem since the parameter fitting employs the entire data set. However, the model may have a shortcoming in the small number of environmental driving parameters, which may oversimplify the complex biogeochemical processes involved in the carbon dioxide exchange between soils, plants and the atmosphere. While the entire temperature sensitivity of the modelled NEE is manifested through changes in TER, the effect of temperature on the biochemical reactions in GPP is neglected (Haraguchi and Yamada, 2011). At the same time, no correlation between air temperature and model residuals could be detected (homoscedasticity), which indicates that the temperature-induced variability was sufficiently considered. The confounding effect of a high vapour pressure deficit (VPD), which tends to take place in the afternoon, leading to a limited photosynthetic activity, was not taken into account (Lasslop et al., 2010). However, only very few days with low humidity (VPD >10 hPa) occurred, and the typically asymmetric diurnal cycle of NEE could not be found on these days. A missing linkage of the model with potential flux limitations through a low soil moisture is deemed appropriate given the constantly high moisture availability in the permafrost-affected soils at the study site (Gao et al., 2017; Minkkinen et al., 2018). The diverse effect of direct and diffuse solar radiation on photosynthetic efficiency was also not taken into consideration (Williams et al., 2014). This effect plays a tangential role for the low sedges but adds uncertainty to the light response curves calculated for the larger shrubs. Further uncertainty may also be appended by a potential inaccuracy in both surface classification and footprint model. While the former is deemed appropriate due to extensive ground truthing (in the form of a comparison between the classification results and comprehensive field surveys), the latter is difficult to assess. However, the employed footprint model is a widely applied tool within the flux community, and it constitutes a suitable model for this study site in a flat tundra landscape with low roughness lengths (Thomas Foken, personal communication, 2015). More importantly, the flux decomposition method, as carried out in the present study, may approach methodical limits, if the surface classes in the footprint are too uniformly distributed and/or their individual flux rates are too similar. Whether the assignment of flux rates from a mixed signal to individual surface classes is still possible under these circumstances may be an objective of further studies at other sites.

4.2 Validation of the decomposed fluxes

The flux decomposition yielded insights into the flux dynamics of both investigated vegetation classes. The validity of these dynamics and hence the reliability of the employed model is examined utilising four approaches.

Firstly, it has been demonstrated that the photosynthetic cycle of a canopy during a growing season is linked to its seasonal changes in greenness (Peichl et al., 2014; Sonnentag et al., 2012). The evolution of canopy greenness can be examined by determining the green chromatic coordinates (g_cc) of a target area in images obtained by digital repeat photography (Richardson, 2012). Employing the images from the time-lapse camera on the flux tower, this method yielded g_cc values for vegetation class 1 with a central tendency that is significantly greater than the one of the g_cc values for vegetation class 2 (P<0.05). These differences in greenness substantiate the most prominent result of the flux decomposition: the greater photosynthesis of vegetation class 1 at the onset of the growing season in 2014 (Fig. 9).

Figure 9 Outcome of the phenocam approach, i.e. determining the green chromatic coordinates (g_cc) of two target areas in time-lapse images that were taken between 18 June and 6 July 2014. The g_cc values were computed for both vegetation classes and depict the fraction of the green colour in relation to the three primary colours in the RGB colour space. While the scatter plot (a) displays daily means of photosynthesis versus their corresponding g_cc values, the box plots (b) visualise the difference in the distribution of g_cc values between both vegetation classes. The significantly greater greenness in vegetation class 1 was associated with larger photosynthetic rates, whereas vegetation class 2 was characterised by a less green canopy and thus a lower photosynthetic activity (P<0.05).

Secondly, during periods with a certain wind direction and atmospheric stability, the fetches of some observed fluxes were strongly dominated by only one vegetation class as opposed to the commonly mixed signals. Thus, observed fluxes that are accompanied with a large contribution of one vegetation class (Ω>0.7) were compared to fluxes that were modelled for the same vegetation class. The choice of an Ω of 70 % rested in the desire to identify a trade-off between both applying many fluxes for a broad statistical basis (low Ω) and utilising many fluxes without a mixed fetch for an accurate evaluation (large Ω). Both observed and modelled fluxes match very well as indicated by a mean R² of 0.88 and a mean RMSE of 0.82 µmol m⁻² s⁻¹. When MAE is applied as an intuitive error metric, the decomposed fluxes are associated with a mean error of roughly 0.56 µmol m⁻² s⁻¹. The frequent similarity of the vegetation-class-specific flux rates, however, reduces the effectivity of this validation test. Therefore, the observed fluxes governed by one vegetation class were also compared to fluxes modelled for the other class. This counter-check caused a rise in mean RMSE and MAE by 89 % and 99 %, respectively, thus lending further credibility to the modelled flux rates. It can be assumed that this rise would be far greater if the flux rates of both vegetation classes were less similar.

Thirdly, closed chamber measurements have been carried out with an opaque chamber during mid-June 2014 in vegetation class 2 east of the flux tower (Benjamin Runkle and Alex Sabrekov, personal communication, 2016). Similar to the respiration modelled for this class, a mean carbon dioxide flux with a standard deviation of 2.1±0.9 µmol m⁻² s⁻¹ was observed. This mean, however, is based on five individual discontinuous chamber measurements and thus conclusive to only a limited extent, since taking the spatial variability into account is crucial, when fluxes are scaled between eddy covariance and chamber measurements (Oechel et al., 1998). Chamber measurements have also been conducted in vegetation class 1 but only to the exclusion of shrubs due to their chamber-incompatible size.

Fourthly, the discussion of the obtained fitting parameters and their comparison with values estimated at other sites gives further confidence in the validity of the decomposed fluxes (Figs. 6 and 7):

• The estimated R_base values follow a temperature-driven seasonal cycle, in which R_base,2 is mostly lower than R_base,1. A smaller autotrophic respiration can be attributed to the lesser biomass of the sedges, and a smaller heterotrophic respiration can be ascribed to both increased soil moisture and decreased soil temperature, which in turn hamper microbial activity in the depressions (Hobbie et al., 2000; Walz et al., 2017). For comparison with values found at other sites, a mean peak season TER was computed for vegetation class 1 (2.8 µmol m⁻² s⁻¹) and vegetation class 2 (2.3 µmol m⁻² s⁻¹). While the latter respiratory rate corresponds to the mean mid-growing season TER of 2.2 µmol m⁻² s⁻¹, which was estimated for northern peatlands, the former rate is greater (Frolking et al., 1998; Laurila et al., 2001). The comparatively large respiration in vegetation class 1 is likely due to both the large willow shrubs (fostering autotrophic respiration) and the large active layer depth (facilitating heterotrophic respiration).

• The estimated Q₁₀ values of 1.42 and 1.48 are well within the range of 1.3≲Q₁₀≲1.5, which was retrieved across different ecosystems and climates (Mahecha et al., 2010). Furthermore, the fact that Q_10,1 was lower than Q_10,2 is in accordance with a concept, which suggests a correlation between a lower/greater soil temperature sensitivity and a drier/wetter tundra (Olefeldt et al., 2013).

• The estimated P_max values also follow a seasonal course reflecting the growth and senescence of the canopy. The reason for $P_{max,\mathrm{1}}$ being greater than $P_{max,\mathrm{2}}$ is the larger biomass of the bushes relative to the sedges. The values agree well to the maximum assimilation rates of approximately 15.9 and 11.1 µmol m⁻² s⁻¹ that are, respectively, found for Salix pulchra and Carex aquatilis during the peak of the Arctic growing season (Oberbauer and Oechel, 1989; Tieszen, 1975). Given a mean mid-growing season P_max of 8.6 µmol m⁻² s⁻¹ for northern peatlands, $P_{max,\mathrm{2}}$ (8.9 µmol m⁻² s⁻¹) constitutes a representative uptake capacity, whereas $P_{max,\mathrm{1}}$ (12.3 µmol m⁻² s⁻¹) suggests a comparatively large potential for sequestering carbon dioxide (Frolking et al., 1998; Laurila et al., 2001). Another aspect that indicates the reliability of the estimated P_max values is their relationship with the NDVI as seen at many other tundra ecosystems (Mbufong et al., 2014; Shaver et al., 2007). Regarding both growing seasons, the footprint's NDVI was greater in 2015, suggesting a more active vegetation than in 2014 (ORNL, 2017). Similarly, the P_max values of both vegetation classes, in particular the values of the more abundant vegetation class 2, were greater during 2015. Satellite records for tundra landscapes are, however, often confounded by various effects that are particularly profound in high-latitude regions (Stow et al., 2004). Therefore, satellite-derived NDVI values of tundra ecosystems may need to be double-checked with optical sampling in the field, if they are applied to resolve interannual differences (Gamon et al., 2013).

• The estimated α values amount to 0.042 (α₁) and 0.04 (α₂), and are thus greater than the mean mid-growing season α of northern peatlands amounting to 0.023 (Frolking et al., 1998; Laurila et al., 2001). The high light sensitivity indicates an efficient physiology enabling a considerable photosynthetic activity at low irradiance levels. The ratio between the α values of both vegetation classes is similar to the ratio of quantum yields that were compiled from tundra landscapes in Alaska and Sweden during the Arctic peak season: 0.038 for Salix spp. and 0.03 for Carex spp. (Shaver et al., 2007).

4.3 Interpretation of diurnal, seasonal and interannual flux variabilities

On the diurnal scale, the temporal variability in the carbon dioxide fluxes was controlled by meteorological conditions. For comparison among climate-relevant trace gases at this site, the methane fluxes exhibited a larger temporal variability, which was rather governed by the spatial variability, i.e. the constantly varying source area composition in the fetch (Rößger et al., 2019a). On the interannual scale, the carbon dioxide fluxes displayed, in contrast to methane fluxes, a larger variability, which was driven by typical abiotic factors such as snowmelt timing and initial growing season temperatures (Aurela et al., 2004; Groendahl et al., 2007). Due to the inhomogeneous surface, however, biotic factors such as canopy structure and distribution also provide explanatory power.

In 2015, the rapid snowmelt coincided with the spring flood, thus enabling both vegetation classes to commence their canopy development concurrently in early June. The growing season was initialised in mid-June by mosses, which are much more abundant in vegetation class 2 (Fig. 8). Mosses are, in contrast to vascular plants, able to start assimilating right after snowmelt since their photosynthetically active tissue can be maintained over winter (Oechel, 1976). From this point until late September/early October, mosses formed a basal net uptake. Considerable moss activity until late autumn has also been observed on the nearby river terrace (Eckhardt et al., 2019; Kutzbach et al., 2007b). Furthermore, mosses can account for distinctly more than half of total photosynthesis, as demonstrated for graminoid areas with high moss cover (Douma et al., 2007; Sommerkorn et al., 1999). However, it is possible that mosses did not fully photosynthesise throughout the growing season due to their tendency to lower their photosynthetic capacity under high irradiance (Murray et al., 1993). This light stress depends on cloudiness, solar altitude, moss structure and shadowing by vascular plants, altogether promoting a late-season activity of mosses while other plants were already dormant (Zona et al., 2011). On top of the basal moss activity, the shrubs of vegetation class 1 exhibited a larger net uptake until the growing season peak around late July/early August, after which the sedges of vegetation class 2 dominated the carbon dioxide exchange. The fact that Carex spp. started growing earlier than Salix spp. has also been observed at other sites; however, considerable variation exists in the timing of phenological events both among and within species (Chapin et al., 1992; Wielgolaski et al., 1975).

In 2014, air temperatures were higher than the monthly long-term means throughout the measurement period (Fig. 4). During the early and slow snowmelt in mid-May, the low sedges and mosses remained buried in the depressions with accumulated snow longer than the large bushes on the elevated ridge with less snow. Thereby, the willow twigs were exposed to daytime temperatures above freezing leading to the development of catkins in late May already. Hence, vegetation class 1 was more advanced in its phenology than vegetation class 2 at the onset of the growing season. The consequence was the substantially larger net uptake of the shrubs until the seasonal peak in early August. Apparently, the shrubs largely benefitted from elevated early growing season temperatures, an effect that has also been found favourable for shrub encroachment in the Arctic (Myers-Smith et al., 2011). Incidentally, shrubs have been growing on Samoylov Island only since the 1960s (Eva-Maria Pfeiffer, personal communication, 2017). Besides the delayed phenological development, the low carbon sequestration of vegetation class 2 during that period can also be attributed to a soil moisture deficit-induced decline in net assimilation of mosses, as they are prone to desiccation due to both missing roots and the absent ability to actively regulate their internal water content (Turetsky et al., 2012). Similar to the other year, vegetation class 2 dominated the net uptake after the growing season peak, in particular during late August, which can likely be ascribed to enhanced moss activity.

Table 3 Putting the carbon dioxide budgets of the flood plain in perspective with budgets estimated by means of the eddy covariance method at other Arctic sites. A negative/positive difference in the budgets indicates a carbon dioxide sink strength of the comparison site being greater/lower than the one of the flood plain on Samoylov Island. For an appropriate comparison, the budgets of this study's site were recomputed for the day of the year (DOY) periods of the other sites. In this process, representative 2015 flux data were basically utilised, unless the budget of the other site derived from 2014. If the budget period of the other site exceeded the period of this study’s budget, the (few) missing days in the start/end were filled with the daily sum of the first/last day in the budget period. The numbers next to the bars refer to the numbers in the circumpolar map (Fig. 1). Compared to various Arctic sites, the flood plain on Samoylov Island constitutes a distinctly strong carbon dioxide sink.

4.4 Appraisal of the budgets' representativeness

The spatial representativeness of the observed fluxes can be assessed with the sensor location bias (Schmid and Lloyd, 1999). If the flux rates of the considered surface classes are similar, as in 2015, the deviation between the respective surface class compositions in footprint and area of interest plays a minor role (Table 1). In 2014, the sensor location bias came into effect as the flux decomposition unveiled a varying seasonality between both vegetation classes that was concealed in the net signal. In this case, a quantitative comparison of the flux budgets with other sites lacks validity due to a potentially non-representative surface class composition; i.e. the comparison of the flood plain's greenhouse gas budgets with the budgets of the river terrace must remain restricted to Samoylov Island and cannot be extended on the Lena River Delta (Table 2). The revealing outcome of the flux decomposition proves its utility for an enhanced interpretation of eddy covariance data by gaining insights into the phenological dynamic of individual vegetation classes. It also demonstrates that climatologically unusual conditions can adversely affect the representativeness of the footprint, resulting in the potential need to regularly examine the representativeness of apparently homogeneous footprints, in particular during prolonged unusual weather conditions as biased budgets may otherwise be estimated.

The temporal representativeness of the obtained budgets may thus be constrained on the interannual scale. As the air temperatures in 2015 better correspond to long-term means than in 2014, the 2015 budgets are better suited for an intersite comparison (Table 3). Moreover, the obtained budgets also possess a confined validity on the annual scale since they only cover a period that is similar to the growing season. Outside this period, no uptake of carbon dioxide occurs, implying a lower year-round sink strength for greenhouse gases. This assumption is based on the accumulating evidence that the release of carbon dioxide and methane is not negligible during the very cold winter – in contrast to the traditional view of a wintertime inactivity in Arctic ecosystems (van der Molen et al., 2007). For instance, at multiple sites in Alaska, the cold season release of carbon was found to equal 1–2 times the warm season net uptake (Euskirchen et al., 2012; Oechel et al., 2014; Zona et al., 2016).

4.5 Comparison of the budgets with other Arctic sites

Across various Arctic flux sites, the flood plain of Samoylov Island exhibits a carbon dioxide sink strength that is distinctly greater than the average (Fig. 1 and Table 3). This aspect appears noteworthy when local conditions are taken into consideration: the mean net radiation during the growing season is lower than for most Arctic sites, and the underlying permafrost displays one of the lowest ground temperatures in the world (Boike et al., 2013; Obu et al., 2018; Romanovsky et al., 2010). The diminishing effects of these climate factors are counterbalanced by the deposition of nutrients in the course of spring flooding (van Huissteden et al., 2005). Among the three great Siberian rivers draining into the Arctic Ocean (Ob, Yenisei, Lena), the Lena River ranks first in terms of total suspended matter (Cauwet and Sidorov, 1996). A large portion of this matter is transported during the annual spring flood, thereby regularly mitigating the nutrient limitation that affects many Arctic ecosystems (Beermann et al., 2014; Fedorova et al., 2015).

More specifically, the net uptake of the flood plain on Samoylov Island is distinctly weaker compared to flood plains of the Siberian rivers Kolyma and Indigirka (Kittler et al., 2017; Parmentier et al., 2011). Other Siberian sites encompass Seida and Lavrentiya, which exhibit a similar and stronger net uptake, respectively (Marushchak et al., 2013; Zamolodchikov et al., 2003). Furthermore, the flood plain's net uptake is considerably stronger than budgets of high Arctic sites in Svalbard, Greenland and Canada (Lafleur et al., 2012; López-Blanco et al., 2017; Lüers et al., 2014; Lund et al., 2012). In comparison with sites in either the low Arctic or sub-Arctic, no general conclusions can be drawn, which is likely due to the ubiquitously high spatiotemporal flux variability in the Arctic region. Also, no uniform picture emerges in comparison to Scandinavian peatlands (Aurela et al., 2002, 2009; Fox et al., 2008). Compared with sites in the northern part of the north slope of Alaska, the flood plain exhibits a substantially stronger net uptake (Oechel et al., 2014; Raz-Yaseef et al., 2017); in the southern part, however, similar net uptakes seem to prevail (Euskirchen et al., 2016).

5 Conclusion

The core of the present study is the advanced scaling options of the demonstrated flux decomposition methodology, i.e. fitting a set of area-weighted, surface class-specific flux models to an observed flux. In this way, two major advantages could be gained. Firstly, downscaling net flux signals from the mesoscale to the microscale yielded flux rates for homogeneous landscape units, therefore generating valuable insights into seasonal variability and functional flux to flux driver relationships of major tundra vegetation types. Moreover, these unbiased flux rates offer the possibility to aid the calibration of macroscale models or the validation of their subgrid variability. Secondly, upscaling the decomposed flux rates to a larger area circumvented the sensor location bias of the study site and thus yielded defensible flux budgets, which take the pronounced surface heterogeneity into account. Moreover, the values estimated for the fitting parameters (in particular P_max) provide the opportunity to contribute to the estimation of carbon dioxide budgets on the macroscale (e.g. pan-Arctic) based on their relationships with remote-sensing-derived parameters such as NDVI.

While the aggregated seasonal flux rates of both pre-defined classes (bushes and sedges) were mostly similar, the flux decomposition revealed a varying seasonality that was hidden in the net signal during a comparatively warm spring period. Accordingly, a seasonal difference between locally observed and regionally estimated fluxes can emerge in response to climatologically unusual conditions. This aspect may gain importance against the projected rise in weather extremes in the course of climate change. Beyond such anomalous situations, the flux decomposition may also be important in a general context, as footprints are frequently assumed homogeneous, but surfaces are seldom entirely homogeneous (depending on the desired scale and the examined greenhouse gas). In this context, the flux decomposition methodology can be adopted in other tundra ecosystems as well as regions outside periglacial environments and hence may be supportive in the fields of landscape ecology, experimental agronomy, catchment hydrology and biogeochemical modelling.

Appendix A

The modelling of carbon dioxide fluxes was based on a fitting procedure, which was designed to enable an ecophysiological interpretation of eight fitting parameters: R_base, Q₁₀, P_max, α for each vegetation class. Thus, their best-fit values had to be in a meaningful range and statistically significant (P<0.05). Moreover, their values had to be arranged in reasonable seasonal courses with greater values during the growing season and lower values during the shoulder seasons. This overall target was achieved via four steps (Fig. 3).

Step 1: Prior to calibrating the model, the following requirements needed to be satisfied for every window: flux samples were available for at least 80 % of the period, the air temperature spread did not fall below 12 ^∘C, and the mean air temperature did not drop below −10 ^∘C. Imposing these requirements was meant to ensure robust and representative fits. The model was first fitted to the observed fluxes of each day, utilising a moving window with a fixed size of 14 d. The choice of a suitable window size was based on identifying an optimum between two conflicting demands: the window size ought to have been as small as possible to capture most of the variability in fluxes, whereas the window size ought to have been as large as possible to obtain less noisy time series of preferably significant values for the eight fitting parameters. Running the model with varying window sizes and counting the number of significant values for each model run yielded the following: increasing the window size caused the number of significant values to rise, soon to level off, and eventually at a window size of 14 d, to remain at similar values.

The purpose of this step involved the fixation of Q₁₀ in order to prevent overparameterisation, and moreover, alterations in temperature sensitivity were thought to be less plausible and hence expected to be negligible. This assumption was confirmed by the time series of estimated Q₁₀ values displaying an implausible variability, whereas the other fitting parameters presented a rather seasonal course. Based on the deliberation of negligible alterations in temperature sensitivity during both years, the model was run for 2014 and 2015 together during this step, whereas the model was respectively run for 2014 and 2015 during the next steps. The application of a larger period provided more data points, with the aid of which Q₁₀ could be fixed at a more representative value. Two final Q₁₀ values were determined for each vegetation class by calculating the median out of all estimated best-fit Q₁₀ values, which met the following two requirements: statistical significance and an associated coefficient of determination of at least 0.75 between observed NEE and modelled NEE.

Step 2: The model was run with Q₁₀ being fixed throughout the measurement periods in 2014 and 2015, applying a moving window with a fixed size of 14 d and a step size of 1 d again. The requirements checked before fitting, as laid down in the previous step, had to be met again except the requirement of a sufficient air temperature spread.

The aim of this step comprised the creation of two replacement functions for α after six best-fit parameters were estimated. The necessity for replacement functions arose through large peaks in the time series of α. These peaks tended to occur at the onset of the growing season and were hence deemed spurious. Large α values would have promoted photosynthesis, which was of rather minor magnitude at that time of the year. In order to reproduce the low observed NEE, the erroneously elevated GPP would have been counteracted by a mistakenly enhanced TER utilising a large R_base. This outcome would have exacerbated the potential problem of equifinality and thus also complicated the interpretation of the fitting parameters.

To prevent this adverse circumstance, two replacement functions, one for each vegetation class, were calculated by fitting a Gaussian bell curve to the time series of significant α values. In addition, two threshold functions were computed for each replacement function by adding/subtracting 30 % of the function values to/from the replacement function. Hence, the threshold functions formed an interval around the replacement function within which the estimated α values were accepted during the further procedure. The threshold of 30 % was visually selected since this value generated an interval outside of which only the peaks were situated; i.e. spurious and meaningful α values could be reliably separated.

Step 3: The model was initially run with the same parameterisation as in the previous run but employing a moving window with a flexible size for every day. The application of a flexible window allowed a closer reproduction of the variability in the observed data through adjusting its size. However, since small windows were in conjunction with a small amount of flux samples, which increased the risk of estimating insignificant parameters, every fit required a minimum of 240 flux samples, which equals 5 d with 48 fluxes per day. Based on this setting, the model was run, and the estimated parameters were checked for significance. If one best-fit parameter was insignificant, the window size was increased by 1 d, and the model was run again. This procedure was repeated until a maximum window size of 20 d, if all fitting parameters were not significantly estimated before utilising a preferably smaller window size.

The objective of this step included the bulk of model calibration within the fitting procedure. Hence, after the initial model run of this step, its output was inspected in two respects: the significance of the remaining fitting parameters (R_base, P_max, α) and the location of α values (inside or outside the acceptance interval). In case of all six fitting parameters being significant and the two fitted α values being situated between the respective thresholds, the estimated NEE was accepted and appended to the modelled time series. If one criterion/both criteria was/were not fulfilled, another model with less parameters was run employing, again, a moving window with a flexible size for every day. This simplification comprised the application of α values adopted from the previously defined replacement function. The model choice depended on the vegetation class, where the criterion/criteria was/were not satisfied. Hence, α values from the replacement function were employed for either one or both vegetation classes. For instance, if the fitting parameters of only one vegetation class were insignificant, only this vegetation class was refitted applying a replacement function α value whilst reutilising the retained significant fitting parameters of the other vegetation class. Subsequently, the significance of the refitted parameters was examined. If all parameters were significant, the correspondingly estimated NEE was added to the modelled time series. Any remaining insignificance values were otherwise addressed in the next step.

In preparation for the next step, two replacement functions for P_max were created. This fitting parameter was chosen over R_base since P_max featured more insignificant values than R_base. Once again, a Gaussian bell curve was fitted to the time series of significant P_max values of each vegetation class.

Step 4: Towards the end of the procedure, a greatly simplified model was run for each day, applying a moving window with a fixed size. This size corresponded to the average of all window sizes found during the previous step.

The goal of this step encompassed the remaining model calibrations for a complete time series of modelled NEE. To achieve this target, the model included only three best-fit parameters: R_base twice and P_max once. The second P_max value for the other vegetation class was adopted from its previously calculated replacement function. This confined parameterisation was, given a constant amount of observed flux samples, associated with an elevated number of degrees of freedom, which in turn allowed a more precise estimation of the remaining fitting parameters; i.e. their confidence intervals were smaller. In this way, all best-fit parameters were significant and could be utilised for reliable modelling of NEE.