Water track distribution and effects on carbon dioxide flux in an eastern Siberian upland tundra landscape

Shrub expansion in tundra ecosystems may act as a positive feedback to climate warming, the strength of which depends on its spatial extent. Recent studies have shown that shrub expansion is more likely to occur in areas with high soil moisture and nutrient availability, conditions typically found in sub-surface water channels known as water tracks. Water tracks are 5–15 m wide channels of subsurface water drainage in permafrost landscapes and are characterized by deeper seasonal thaw depth, warmer soil temperatures, and higher soil moisture and nutrient content relative to adjacent tundra. Consequently, enhanced vegetation productivity, and dominance by tall deciduous shrubs, are typical in water tracks. Quantifying the distribution of water tracks may inform investigations of the extent of shrub expansion and associated impacts on tundra ecosystem carbon cycling. Here, we quantify the distribution of water tracks and their contribution to growing season CO2 dynamics for a Siberian tundra landscape using satellite observations, meteorological data, and field measurements. We find that water tracks occupy 7.4% of the 448 km2 study area, and account for a slightly larger proportion of growing season carbon uptake relative to surrounding tundra. For areas inside water tracks dominated by shrubs, field observations revealed higher shrub biomass and higher ecosystem respiration and gross primary productivity relative to adjacent upland tundra. Conversely, a comparison of graminoid-dominated areas in water tracks and inter-track tundra revealed that water track locations dominated by graminoids had lower shrub biomass yet increased net uptake of CO2. Our results show water tracks are an important component of this landscape. Their distribution will influence ecosystem structural and functional responses to climate, and is therefore of importance for modeling.


Introduction
Climate warming is relaxing constraints on vegetation productivity in arctic tundra ecosystems imposed by low temperatures and short growing seasons (ACIA 2005, Myers-Smith et al 2011, IPCC 2014. Observational evidence indicates that warming and associated consequences, such as permafrost thaw, stimulate vegetation productivity (Walker et al 2006a, Elmendorf et al 2012a, 2012b, Natali et al 2012. Remote sensing observations and proxy data indicate widespread increases in vegetation productivity in recent decades (Forbes et al 2010, Beck and Goetz 2011, Guay et al 2014 often accompanied with an expansion of shrub cover (Myers-Smith et al 2011, Frost and Epstein 2014, IPCC 2014. The potential for widespread shrub expansion is of particular importance because accompanying changes in albedo (Loranty et al 2011) will act as a positive feedback to climate warming (Chapin et al 2005, Bonfils et al 2012. However, shrub expansion is heterogeneous in space (Tape et al 2012), and mounting evidence suggests that shrub responses to climate are strongest in areas with relatively high soil moisture (Elmendorf et al 2012b, Myers-Smith et al 2015, Swanson 2015 and nutrient contents (Tape et al 2012).
Water tracks are 5-15 m wide bands of high productivity vegetation resulting from subsurface drainage and associated increases in soil temperature, moisture, and nutrient availability driven by topography and variability in active layer depth (Chapin et al 1988, Hastings et al 1989. Water tracks contain many of the conditions thought to be necessary for shrub expansion, and have been identified as a component of the landscape capable of supporting tall shrub cover (Swanson 2015). Several studies have examined the carbon (C) cycle effects of water tracks and upslope water additions, finding increased rates of primary production and ecosystem respiration (ER) that indicate water tracks may be important for understanding the tundra C balance at landscape and regional scales (Chapin et al 1988, Oberbauer et al 1991, 1989. Despite the known effects of water tracks on vegetation composition, C and nutrient cycling, and their likely importance for shrub expansion, they have only been mapped across a limited region on Alaska's North Slope (Walker and Maier 2007) and their distribution across Siberian landscapes remains unknown. Quantifying the areal extent of water tracks will help improve our understanding of landscape scale variability in the drivers of shrub expansion and C cycling in tundra ecosystems. This will be especially important given recent findings indicating that tall shrub distribution is limited by environmental factors, including soil wetness and thaw depth, and that tall shrubs grow along water tracks (Swanson 2015). Moreover, climate induced changes in arctic precipitation patterns (ACIA 2005, IPCC 2014) may lead to divergent ecosystem C cycle responses in these areas of complex terrain (e.g. Riveros-Iregui et al 2012). As a result, the goals of this study were to quantify the areal extent of water tracks and associated variability in vegetation biomass and C flux in a northeast Siberian tundra landscape. To achieve these objectives we used geospatial data to investigate water track distribution, field observations to assess environmental conditions and ecosystem-scale C flux patterns, and a simple ecosystem model alongside leaf area derived from remotely sensed data to assess landscape-scale patterns in C fluxes.

Site description
The study was conducted in northeastern Siberia in the Sakha Republic, Russia, along the eastern bank of the Kolyma River (figure 1). Field observations were collected at Krutaya Drisva (69.34°N, 161.51°E), approximately 32 km south of the river's mouth at the Arctic Ocean, approximately 30 m above sea level, on a 2.5% northwest-facing slope. Average recorded January, July, and annual temperature and annual precipitation at the Ambarchik Bay meteorological observatory ∼40 km from the site are −29.1±4. 1°C, 7.6±2.2°C, −11.7±2.4°C, and153±75.0 mm, respectively, from 1950 to 2009 (Global Historical Climatology Network: http://ncdc.noaa.go v/ ghcnm/v 3.php). The mean temperature during the field study period from 9 July and 18 July 2014 was 14.6°C, and mean relative humidity was 77.7%. The site contained tussock tundra and a water track (figure 2). This area is part of a region of tundra with rolling topography-which contributes to increased drainage-down slope from a nearby mountain range. The tundra within the site contained a mix of shrubs (Betula nana and Salix spp.), tussock forming sedge (Eriophorum vaginatum) and other sedge vegetation (Carex spp). The water track matched the description of Chapin et al (1988) and was dominated by large patches of B. nana and Salix spp. alternating with equally large patches of Carex spp.

Geospatial datasets
Geospatial analysis was carried out using one composite raster image and a digital surface model (DSM). These data were acquired from the Polar Geospatial Center (PGC) at the University of Minnesota. The image was orthorectified and corrected to top-atmosphere reflectance by the PGC. This data encompassed areas of tundra and boreal forest along the eastern bank of the Kolyma River, including the study site at Krutaya Drisva. The image was captured by World-View 2 on 15th August 2011 and covers approximately 1165 km 2 . It consists of eight spectral bands: RED (630-690 nm), red edge (705-745 nm), coastal blue (coastal: 400-450 nm), blue (450-510 nm), green (510-580 nm), yellow (585-628 nm), near infrared 1 (NIR1: 770-895 nm) and near infrared 2 (NIR2: 860-1040 nm) with 1.85 m nadir resolution. The DSM was generated from stereo-pairs of high-resolution imagery and covered approximately 3000 km 2 of the region. It contained a single band showing elevation at 2 m resolution with a vertical resolution of <1 m not accounting for variations in vegetation canopy height. The only other peak growing season image available from the region was a QuickBird image mosaic acquired 4-11 July 2011.
The scope of our geospatial analysis was limited to the area of these images that encompassed tundra. We considered tundra to be high latitude locations with low growing vegetation and a lack of trees (Walker et al 2005). Tree line was determined by eye and a polygon encompassing the 'tundra' was digitized manually. This polygon was then used to clip an area of approximately 448 km 2 from the raster files for use in geospatial analysis.

Water track mapping
A number of potential methods exist for automatically classifying water tracks and surrounding tundra including supervised, unsupervised and object oriented spectral classifications as well as topographic analysis including flow path mapping. Given that water tracks are visible to the naked eye in highresolution imagery, manual digitization would be the most effective means of classifying this data. However manual digitization is less feasible on a large scale and a spectral classification approach was chosen. In order to determine the best method for generating the final map of water tracks a variety of different techniques mentioned above were used to classify a small sample image-approximately 32 km 2 . These classifications were compared to a map of water tracks outlined using heads up digitizing. Maximum likelihood supervised classification of all the bands present in the WorldView 2 image was deemed the best approach and was thus used in the final analysis followed by use of the boundary clean tool in ArcMap.
The entire data set was then classified and water track area was calculated (see figure 3 for overview). Training points were selected in water tracks, intertrack tussock tundra, inter-track shrub tundra and a variety of other landforms for a total of 15 training samples. These classifications were then aggregated into the broad categories of water track and inter track tundra using the reclassify tool. The normalized difference vegetation index (NDVI) was calculated on a per cell basis from the image using the equation below (equation (1) NDVI compares red and near infrared radiation reflected, which can be used to determine the presence and photosynthetic capacity of vegetation (Tucker 1979). As an indicator of photosynthetic capacity NDVI is inherently tied to leaf area and gross primary production (GPP) (Williams et al 2006. Surface water was mapped using a supervised classification and the normalized difference water index (NDWI) (equation (2) NDWI compares costal blue and NIR2 to detect the presence of standing water in areas greater than a single pixel (Wolf 2012). Representative areas of surface water and tundra were selected from within the World View 2 image and a supervised classification was run using these points. The total area of water tracks, tundra and standing water in the study area was then calculated. The zonal statistics tool in ArcGIS and the maps produced above were used to calculate the average of leaf areas indices (LAI) and NDVI for water tracks and tundra in the entire World View 2 image.
The accuracy of our classification was assessed using 500 m diameter buffers around fifteen points randomly generated across the study area. In each location well-developed water tracks (Walker and Maier 2007) were hand digitized using the World View 2 image. The output of the supervised classification was then compared to these hand digitized samples. Our initial classification was tailored to avoid errors of commission at the expense of increased omission.

Landscape-scale CO 2 flux modeling
In order to quantify the effects of water tracks on landscape-scale CO 2 dynamics we used the model of tundra net ecosystem CO 2 exchange (NEE) developed by Shaver et al (2007). We chose this modeling approach because it enabled a more robust assessment of water track CO 2 fluxes by covering a wider range of vegetation and environmental conditions than would have been possible using field observations alone (described below). The model relies on the principle of functional convergence of NEE, which means that a single parameterization of an NEE model applies across all tundra plant functional types. The model requires LAI, air temperature (airT), and incident photosynthetically active radiation (I) as inputs, and has been successfully applied for numerous tundra vegetation types at a range of spatial and temporal scales (Shaver et  The World View 2 NDVI map described above was used to produce a map of LAI with the equation and constants described by Street et al (2007) (equations (1) and (3)). We used the parameterization covering all vegetation types, as we did not have data We then combined this map of LAI with 434 records of half hourly meteorological data (described below) in order to model NEE, GPP, and ER using the Raster package (Hijmans and van Etten 2012) in R (R Development Core Team 2015). NEE was modeled as in Shaver et al (2013Shaver et al ( , 2007 We parameterized the model using the pan-Arctic parameterization from Shaver et al (2007Shaver et al ( , 2013: the light-saturated photosynthetic rate per unit leaf area or P maxL (14.747 μmol m −2 leaf s −1 ); the Beer's law extinction coefficient or k (0.5 m −2 ground m −2 leaf area); the initial slope of the light response curve or E 0 (0.041 μmol CO 2 fixed μmol −1 photons absorbed); the basal respiration rate or R 0 (1.177 μmol CO 2 m 2 leaf s −1 ); the constant β (0.046°C −1 ); and a constant representing respiration in deep soil horizons or R x (0.803 μmol CO 2 ). Air temperature (airT) and the top of the canopy photon flux (I) was obtained from half hourly meteorological data collected in the field, as described in the following section. Negative NEE values represent a net ecosystem sink of CO 2 .
The modeling operations created 1302 raster layers at 1.84 m resolution spanning the course of the observational period. Parallel computing was implemented in R using the spatial.tool and foreach packages (Greenberg, 2014, Revolution Analytics and Weston 2014). These rasters were aggregated to calculate average GPP, NEE and ER over the ∼9 day observational period. The output of the classification was combined with these to calculate average GPP, NEE and ER for inter-track tundra and water tracks.

Observations of CO 2 flux and environmental conditions
All field data were collected between 9 July and 18 July 2014. We established 12 plots inside a water track and 12 in the surrounding tundra. These 0.36 m 2 plots were spaced ∼10 m apart. The tundra and water track plots were established in areas dominated by deciduous shrubs (n=6) or graminoids (n=6).
We measured NEE under light and dark conditions over the course of three days, resulting in a total of 84 pairs of light and dark measurements. CO 2 fluxes were measured using the static chamber technique (e.g., Shaver et al 2007). For each measurement a 0.3 m tall 0.36 m 2 Plexiglas chamber was placed over the ground and vegetation, and air was circulated between the chamber and a Licor infrared gas analyzer (LI840) to measure the CO 2 concentration every second for approximately 1 min. The chamber was flushed with ambient air before each measurement, and two fans inside the chamber ensured air mixing. After initial quality screening NEE and ER were calculated from light and dark fluxes respectively using the rate of change in CO 2 concentration over the measurement period (R Development Core Team 2015) and the ideal gas law (equation (1), Shaver et al 2007). We then calculated GPP as the difference between NEE and ER.
Each flux measurement was accompanied by measurements of soil temperature (5 cm depth), air temperature, and internal chamber temperature using a Fisher Scientific Traceable Dual-Channel Thermometer, and radiometric surface temperature using an Apogee Infrared Radiometer. Soil moisture was measured using a Decagon Devices GS3 Ruggedized Soil Moisture, Temperature and Electrical Conductivity Sensor (5 cm depth). NDVI was measured using a Decagon Devices Spectral Reflectance Sensor mounted at the end of a 1 m long PVC pipe and held approximately 1.5 m above the plot.
On three separate days, thaw depth was measured at two corners of each plot using a metal thaw depth probe. At the conclusion of the fieldwork the basal diameters of all B. nana and Salix spp. present in each plot were measured using calipers and converted to biomass using allometric equations developed for this region (Berner et al 2015). The environmental data from all 24 plots was averaged by plot and analyzed using both single factor and nested analysis of variance (ANOVAs) in R (R Development Core Team 2015). To compare differences between water-tracks and inter-track tundra, we grouped data by plot type (i.e. water track or inter-track tundra) across vegetation types using a single factor ANOVA. We also examined interactive effects of vegetation type (i.e. shrub or graminoid) and plot type by nesting vegetation within plot type in a nested ANOVA.
A meteorological station (Onset Corp, Bourne MA), which was located ∼500 m from the water track, measured and recorded air temperature, barometric pressure, PAR, relative humidity and moisture content at the half hourly intervals from 9 July to 18 July 2014. These data were used in conjunction with measured fluxes and as model inputs.

Distribution of water tracks
The total areas of tundra, water tracks and surface water in our study area were 408.3 km 2 , 33.4 km 2 (7.4%) and 6.8 km 2 , respectively (figure 4(a)). Water tracks had a mean NDVI and LAI of 0.61 and 0.65, and inter-track tundra had a mean NDVI and LAI of 0.56 and 0.47, respectively (figures 4(b) and (c)).
Our classification had overall accuracy 90.1%, user accuracy of 92.9% in inter-track tundra and user accuracy of 62.2% in water tracks (see supplemental table). Our initial classification avoided errors of commission at the expense of increased omission. This conservative classification approach ensured that other tundra landforms with dense shrub cover, such as floodplains, were not falsely classified as water tracks leading to false increases in water track area, NDVI and LAI. Twenty percent of errors of commission occurred within 20 m of the digitized water track polygons suggesting that our vector digitization presented too discrete and binary a boundary between water tracks and surrounding inter track tundra leading to the omission of the fringes of water tracks and poorly developed water tracks successfully detected by the per-pixle classifier.

Field observations
Soil moisture was higher in the water tracks (P<0.001) compared to inter-track tundra, and significant differences in soil moisture existed between plots of different vegetation and landform types (P=0.001; table 1). Average soil temperature at 5 cm was also higher in the water track plots compared to tundra plots (P<0.001; table 1); however, ground surface temperature was higher in inter-track tundra compared to water track plots (P<0.001; table 1). Average thaw depth was deeper inside the water track than in the surrounding tundra (P=0.001; table 1). Thaw depth was also deeper in water track shrub plots compared to inter-track shrub plots, but did not differ between graminoids plots inside and outside of the water track (P=0.05). NDVI was higher in water track plots when compared to inter-track tundra (P<0.001; table 1) and was also higher in shrub plots of the same landform type (P<0.001; table 1). This confirms the trend observed in the satellite imagery.
Total shrub biomass was significantly greater in water track shrub plots than in inter-track shrub plots (P<0.001; table 2). This trend was driven by higher B. nana biomass in water track shrub plots (P<0.001) in spite of lower Salix spp. biomass when   At the landscape scale, average modeled GPP and ER were higher inside the water track and NEE was more negative relative to inter-track tundra (table 3).
Average plot-level ER was significantly higher in the water track when compared to inter-track tundra (P<0.001; table 4). Average plot level GPP was also higher and there was slightly more CO 2 uptake (i.e., more negative NEE) observed.

Discussion
Water tracks composed 33.4 km 2 or 7.4% of the 448.63 km 2 of tundra analyzed, and they were greater net sinks of CO 2 during the observation period, based on modeled results. Although differences in observed C fluxes between water track and inter-track tundra were not statistically significant, the magnitude and direction of observed differences were similar to modeled differences which encompassed a much larger range of vegetation and environmental conditions. Observations of greater soil moisture, plant biomass, thaw depth, NDVI and C uptake in water tracks relative to inter-track tundra are consistent with patterns observed in other studies and posit a mechanism for increased productivity in water tracks (Chapin et al 1988, Hastings et al 1989, Oberbauer et al 1991. Other authors have attributed these increases in productivity to increased N availability resulting from deeper thaw depths, warmer soil and possibly the liberation and movement of nutrients by flowing water (Chapin et al 1988, Hastings et al 1989, Oberbauer et al 1989, Harms and Jones 2012, Harms et al 2014.

Water track mapping
Water tracks constitute a large area within the study region (7.4% of the total or 7.6% if standing water is excluded). A number of methods using different combinations of bands and supervised and unsupervised classification schemes were investigated for use in classifying water tracks in this study. The classification method, which worked best in small-scale tests (maximum likelihood supervised classification of all the bands present in the WorldView 2 image followed by use of the boundary clean tool in ArcMap), was used in the final analysis and captured the greatest number of bands at the highest resolution, thus providing the most information about a given pixel. However as a result of our conservative classification approach it consistently underestimated the area of water tracks when compared to digitization by eye and is more than likely an under estimate. Twenty percent of errors of commission occurred within 20 m of the digitized water track area suggesting that our digitization may have been too discrete and failed to account for poorly-developed water tracks (e.g. Walker and Maier 2007). Judging water track presence or absence based upon a vector boundary may have lead to the omission of the fringe of water tracks or poorly developed water tracks. The classifiers pixel based approach better captures the detailed boundary between water tracks and tundra potentially leading to apparent commission.
Given that water tracks result from subsurface flow and permafrost topography, inputs of water vary across time and space, and water track formation has been observed at decadal time scales these figures of water track area may be specific to the landscape and time period investigated (Ostendorf and Reynolds 1993, Osterkamp et al 2009).
Differences in remotely sensed NDVI between water tracks and inter-track tundra agree well with field observations showing greater vegetation biomass  and productivity from this (tables 1 and 2) and previous studies (Chapin et al 1988, Hastings et al 1989, Oberbauer et al 1989. NDVI is heavily correlated albeit nonlinearly with vegetation biomass in arctic ecosystems (Jia et al 2003(Jia et al , 2006. These spatial patterns in NDVI are likely the result of increased thaw depth, soil temperature and soil moisture, which in turn affect vegetation structure and function in ways similar to those observed by other authors (Chapin et al 1988, Hastings et al 1989, Oberbauer et al 1989. Given their large area, previous observations of variations in soil C flux and differences in vegetation composition and soil conditions, water tracks are relevant to the study of landscape scale tundra C cycling and ecosystems responses to climate warming including shrub expansion (Oberbauer et al 1991, Swanson 2015).

Environmental conditions
Physical environmental conditions within water tracks differed from adjacent inter-track tundra in several ways. As expected, soils were wetter in the water tracks, and the wetter soils may in turn be affecting other variables (table 1). The higher thermal conductivity associated with wetter soils along with flowing subsurface water and feedbacks between shrubs, soil temperature and permafrost is likely contributing to higher soil temperatures and deeper thaw depths in water tracks-although measurements were taken prior to the period of maximum permafrost thaw and active layer temperature (Chapin et al 1988, Hinzman et al 1991, O'Donnell et al 2009, Blok et al 2010, Shiklomanov et al 2010, Foster 2015. Vegetation type on the other hand does not appear to be affecting soil temperature in any significant way-albeit over a short measurement period. The observed lower surface temperature in water tracks-contrary to the trend in soil temperature-is likely due to increased latent heat flux associated with increased soil moisture and increased transpiration by vegetation biomass as well as shading of the soil surface by shrubs ( Tundra ecosystems are generally nutrient limited, and increased thaw depth coupled with soil temperature and subsurface flow may directly affect nitrogen and phosphorus availability and explain many of the changes in vegetation structure observed in water tracks (table 2; Chapin et al 1988, Hastings et al 1989, Giblin et al 1991, Oberbauer et al 1991, Naito and Cairns 2011

Carbon flux
Variability in ecosystem C fluxes inferred from modeling are consistent with observed differences in NDVI and biomass (tables 1-3) and those observed in other studies (Boelman et al 2003, Raynolds et al 2006, Sweet et al 2015. Both modeled GPP and ER were higher leading to more negative NEE or greater net uptake of C in the water tracks relative to inter-track tundra. Increases in ER resulting from increased biomass, deeper thaw depths, and warmer soils were complemented by increases in GPP (table 3). According to the model water tracks account for 9.9% of total GPP, 9.0% of total ER and 12.1% of NEE during the middle of the growing season, slightly higher than would be expected, given their areal distribution of 7.4% of tundra land area. GPP also increased slightly in water track graminoid plots despite decreased shrub biomass. These estimates only account for C flux during a fraction of the peak of the growing season, and so the relative contribution of water tracks to growing season or annual landscape-scale C fluxes may be different. Phenological differences in vegetation communities during early and late season could affect C fluxes integrated over the entire growing season (Zeng et al 2013, Sweet et al 2014. Differences may also exist outside of the growing season, where C fluxes are governed largely by co-variation in snow cover and soil temperature (Webb et al 2016). Increased snow depth resulting from topographic relief and increased shrub cover has been documented in water tracks and other shrub dominated locations (McFadden et al 2001, Pomeroy et al 2006. Despite greater snow depths the rate of snowmelt is often accelerated in shrub tundra decreasing the magnitude of these effects (Pomeroy et al 2006). Variations in snow depth and snow melt have the potential to impact phenology, growing season start and length and thus vegetation structure and function ( Though not all differences were significant, observed fluxes show differences between water tracks and inter-track tundra similar to those observed in the modeled data. There is a significant difference in ER between water track and inter-track tundra (P<0.001; table 3) and this increase in ER is offset by an increase in NEE leading to greater net uptake of C in the water track plots. Absolute values for mean fluxes differ between observed and modeled results; however, this should be expected given the limited temporal and spatial scale over which the observed fluxes were collected and variations in LAI corresponding to vegetation type resulting from the NDVI-LAI modeling parameters used. The model incorporates a much larger range of vegetation and meteorological conditions, and performs well in predicting broad trends in ecosystem CO 2 fluxes for water track and intra-track tundra that are consistent with observations.

Conclusion
The results presented here show that water tracks constitute a substantial area within tundra ecosystems for this landscape in northeastern Siberia. Water tracks have soil conditions capable of supporting tall shrub expansion, (Naito and Cairns 2011, Swanson 2015) and therefore knowing their spatial distribution may help to refine predictions of vegetation change in the Arctic. The proportional contribution of water tracks to landscape-level C balance is slightly higher than their areal extent. However, differences in vegetation, and moisture and nutrient availability may lead to divergent C cycle responses relative to upland tundra under continued climate warming. Thus, knowledge of water track distribution may also improved understanding of current and future spatial variability in tundra C cycling at the landscape scale. The distribution and extent of water tracks is likely controlled by a combination of topographic, permafrost driven and climatic factors. Future work is necessary to determine how the spatial extent of water tracks varies across the tundra biome.