*Article* **Hydrological Connectivity in a Permafrost Tundra Landscape near Vorkuta, North-European Arctic Russia**

**Nikita Tananaev 1,\* , Vladislav Isaev <sup>2</sup> , Dmitry Sergeev <sup>3</sup> , Pavel Kotov <sup>2</sup> and Oleg Komarov <sup>2</sup>**

	- tpomed@yandex.ru (V.I.); kotovpi@mail.ru (P.K.); tpomed@rambler.ru (O.K.)

**Abstract:** Hydrochemical and geophysical data collected during a hydrological survey in September 2017, reveal patterns of small-scale hydrological connectivity in a small water track catchment in the north-European Arctic. The stable isotopic composition of water in different compartments was used as a tracer of hydrological processes and connectivity at the water track catchment scale. Elevated tundra patches underlain by sandy loams were disconnected from the stream and stored precipitation water from previous months in saturated soil horizons with low hydraulic conductivity. At the catchment surface and in the water track thalweg, some circular hollows, from 0.2 to 0.4 m in diameter, acted as evaporative basins with low deuterium excess (*d*-excess) values, from 2‰ to 4‰. Observed evaporative loss suggests that these hollows were disconnected from the surface and shallow subsurface runoff. Other hollows were connected to shallow subsurface runoff, yielding *d*-excess values between 12‰ and 14‰, close to summer precipitation. 'Connected' hollows yielded a 50% higher dissolved organic carbon (DOC) content, 17.5 ± 5.3 mg/L, than the 'disconnected' hollows, 11.8 ± 1.7 mg/L. Permafrost distribution across the landscape is continuous but highly variable. Open taliks exist under fens and hummocky depressions, as revealed by electric resistivity tomography surveys. Isotopic evidence supports upward subpermafrost groundwater migration through open taliks under water tracks and fens/bogs/depressions and its supply to streams via shallow subsurface compartment. Temporal variability of isotopic composition and DOC in water track and a major river system, the Vorkuta River, evidence the widespread occurrence of the described processes in the large river basin. Water tracks effectively drain the tundra terrain and maintain xeric vegetation over the elevated intertrack tundra patches.

**Keywords:** permafrost hydrology; Russian Arctic; water tracks; hydrological connectivity; stable water isotopes; dissolved organic carbon; electrical resistivity tomography; taliks

### **1. Introduction**

Hydrologic connectivity is a complex concept referring to water transfer in the landscape, or between landscapes, or, more generally, within or between the water cycle units, and its (dis)continuity along the major water transport pathways acting on the watershed [1–3]. It includes both lateral water transfer along slopes, including channelized runoff, and vertical water transfer between surface and subsurface compartments, or between different groundwater aquifers at different depths [4]. Connectivity exists between larger domains, e.g., surface runoff and groundwater flow, landscape elements and fluxes—*structural connectivity*, and between processes—*functional connectivity* [4,5].

Permafrost significantly alters the water cycling through the affected landscapes compared to that in temperate catchments [6–8]. In continuous permafrost, water transport is mostly confined to the active layer, which rarely exceeds 3 m depth, and to vertical and lateral talik zones. Water migration in soils is partly driven by processes related to phase transition in soils [9,10]. The hydrological system structure is simplified, and the

**Citation:** Tananaev, N.; Isaev, V.; Sergeev, D.; Kotov, P.; Komarov, O. Hydrological Connectivity in a Permafrost Tundra Landscape near Vorkuta, North-European Arctic Russia. *Hydrology* **2021**, *8*, 106. https://doi.org/10.3390/ hydrology8030106

Academic Editors: Il-Moon Chung, Sun Woo Chang, Yeonsang Hwang and Yeonjoo Kim

Received: 19 June 2021 Accepted: 17 July 2021 Published: 22 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

existing connections between compartments are exposed to observation [11,12]. Open talik zones are frequently detected under the largest lakes using geophysical methods [13–15] and may serve as pathways for both upward and downward water migration, connecting surface water to intra- and subpermafrost aquifers [16,17]. In discontinuous permafrost, with a deeper active layer (down to 5 m), persistence of residual thaw layers, permafrost fragmentation, and abundance of talik zones, the potential for water exchange between the compartments is significantly higher [18–20].

Stable water isotopes, <sup>2</sup>H (deuterium) and <sup>18</sup>O, are widely used to track water sources and hydrologic connectivity across the compartments and ecosystem classes, including boreal and permafrost-affected catchments [21–24]. The isotopic evaporation signal allows the tracing of connectivity between wetlands and perennial streams and the modeling of surficial wetland runoff contribution to streams during summer [25]. The isotope mass balance method reveals the specifics of the permafrost thaw cycle in continuous permafrost [26,27] and is successfully used in studies of both contemporary lakes and lacustrine paleoenvironments [28].

In discontinuous permafrost, taliks of different kinds, i.e., vertical open and closed taliks, connected to lateral intra-permafrost taliks and residual thaw layers, are responsible for conveying water from the slopes toward the streams [29,30]. In the Northern Yenisey region, isotopically heavier water, originating from late summer precipitation and thermokarst lakes subject to evaporation, was found to contribute significantly to the winter runoff through the residual thaw layer, an interface between the seasonally freezing layer and the top of the permafrost [31]. In the Ob River basin, a strong evaporation signal persists in most river samples in late autumn and around the spring freshet peak dates, evidencing subsurface connections between lakes and rivers of the region [32]. Lake-to-river connectivity is also maintained through sub-lacustrine taliks, both open and closed, developing even under shallow thermokarst lakes [33]. Geophysical techniques, notably electrical resistivity tomography, are useful in describing the complex frozen ground configuration in discontinuous permafrost [34,35].

The climate change presently occurring in the Arctic, followed by the deepening of the active layer, may lead to the rebuilding of existing connectivity patterns through changes in the saturated zone boundary, to an increase in non-frozen soil volume where water migration is possible, and to an increase in groundwater discharge on hillslopes [36], affecting future dissolved organic carbon (DOC) and other constituent fluxes [37]. Trends in regional climate and hydrology may also imply changes in fluvial activity [38], though the latter is showing only minor signs in the first-order fluvial network of the region. The potential effects of climate change and permafrost degradation on hydrological connectivity and water and material fluxes, including biogeochemical cycling, are still poorly understood. Recent reviews acknowledge important knowledge gaps in the subsurface hydrology of permafrost regions, including existing water and carbon transport pathways, their relation to frozen grounds, and alterations in subsurface routing resulting from permafrost degradation [39,40]. Ultimately, permafrost degradation is expected to alter hydrological connectivity in affected catchments, resulting in water flow redistribution between the surface and compartments and changes in seasonal water discharge [20,39].

This study was conceived to address these gaps and to better understand hydrological connectivity and water and DOC transport in a discontinuous permafrost environment at a small scale. We present new data on water stable isotope composition and DOC concentrations from several subarctic streams and water bodies in minor tundra water track catchments near Vorkuta, north-European Russia. These data are used to trace water origin in these water objects under late summer conditions, close to the maximum thaw period and to evaluate microscale hydrologic connectivity in the landscape. Geophysical survey data are used to support the discussion on surface water interaction with groundwater. The study region, with its mild and humid subarctic climate, may serve as a model region for other permafrost regions in transition under observed climate change.

#### **2. Study Area**

Fieldwork was performed in September 2017 in north-European Russia, on the margin of the Bolshaya Zemlya tundra region, about 30 km to the south-west from Vorkuta, Komi Republic, the closest large city (Figure 1a). The major sampling effort was concentrated around Khanovey key study site (N 67◦17.1930 , E 63◦39.2520 ; Figure 1b), an abandoned settlement for railroad workers on the right bank of the Vorkuta River, where the seasonal permafrost research station is maintained by the Department of Geocryology, Moscow State University [41]. The studied location occupies a typical periglacial landscape at the southern margin of the Bolshaya Zemlya tundra, a hilly terrain with gently rolling slopes dissected by hummocky depressions and first-order stream valleys. Permian bedrock, exposed locally in the Vorkuta River bluffs, is overlaid with Quaternary deposits with thickness from 10 to 15 m, mostly loams and loamy clays, with variable ice content and cryostructure [41]. Surficial loamy clays are highly thixotropic, and easily lose structural integrity on stress. Sandy deposits were described in shallow excavations in the adjacent areas but were never encountered at the watershed in question. A topsoil organic layer is omnipresent and has a thickness of 0.3–0.5 m. *Hydrology* **2021**, *8*, x FOR PEER REVIEW 4 of 16

**Figure 1.** Geographical location of the study site, (**a**) in north-European Russia, (**b**) at the Vorkuta River valley slope; (**c**) reference orthophoto image of the studied water track section, north is up, T1 and T2 denote the electric resistivity tomography profiles. **Figure 1.** Geographical location of the study site, (**a**) in north-European Russia, (**b**) at the Vorkuta River valley slope; (**c**) reference orthophoto image of the studied water track section, north is up, T1 and T2 denote the electric resistivity tomography profiles.

The meso-scale topography is dominated by *uvals*, a local name for smoothed hilly chains. These hilly chains have a submeridional orientation and elevation between 170 and 200 m a.s.l. and are divided by the valleys of the Usa River and its major right tributaries, the Vorkuta River and Seida River. The *uvals'* surface is an undulating plain, hosting numerous lakes, peatlands, and a network of overwetted depressions, with mires presumably of thermokarst origin. The hillslopes descending toward the major rivers are transformed by the joint action of fluvial processes, thermal erosion, and linear thermokarst, with the widespread occurrence of water tracks. This network evolves from chaotic in the interfluve zone to a linearly shaped fluvial network hosting intermittent streams (Figure 2), then

**Figure 2.** Typical midslope landscape of the study region, with minor (first-order) water tracks on the background, a second-order water track feature crossing the image left to right, all easily detectable by its contrasting foliage color, elevated tundra patches, and intertracks with *Betula nana* L.

The studied territory occupies the interfluve area between the Vorkuta River and its right tributary, the Lyok-Vorkuta River, and the valley slope inclined toward the east, to the Vorkuta River valley (Figure 1b). This slope is dissected by three first-order water track valleys, one of which was studied in detail in its middle and lower reach, downstream from the railroad crossing (Figure 1c). The water track valley is oriented west to east. Its headwaters are connected, via a network of wet depressions and mires, to the headwaters of all major neighboring water tracks, so that no interfluve exists between the water track systems draining in different directions. For this reason, the basin area of the

and lichen patches.

(**c**)

again becoming poorly organized toward the foot slope with fen-like features. Overwide valleys surrounding the remaining elevated tundra patches resemble those of organic-rich and wide water track classes described by Trochim et al. [42]. **Figure 1.** Geographical location of the study site, (**a**) in north-European Russia, (**b**) at the Vorkuta River valley slope; (**c**) reference orthophoto image of the studied water track section, north is up, T1 and T2 denote the electric resistivity tomography profiles.

*Hydrology* **2021**, *8*, x FOR PEER REVIEW 4 of 16

(**a**) (**b**)

**Figure 2.** Typical midslope landscape of the study region, with minor (first-order) water tracks on the background, a second-order water track feature crossing the image left to right, all easily detectable by its contrasting foliage color, elevated tundra patches, and intertracks with *Betula nana* L. and lichen patches. **Figure 2.** Typical midslope landscape of the study region, with minor (first-order) water tracks on the background, a second-order water track feature crossing the image left to right, all easily detectable by its contrasting foliage color, elevated tundra patches, and intertracks with *Betula nana* L. and lichen patches.

The studied territory occupies the interfluve area between the Vorkuta River and its right tributary, the Lyok-Vorkuta River, and the valley slope inclined toward the east, to the Vorkuta River valley (Figure 1b). This slope is dissected by three first-order water track valleys, one of which was studied in detail in its middle and lower reach, downstream from the railroad crossing (Figure 1c). The water track valley is oriented west to east. Its headwaters are connected, via a network of wet depressions and mires, to the headwaters of all major neighboring water tracks, so that no interfluve exists between the water track systems draining in different directions. For this reason, the basin area of the Interfluves are scarcely vegetated because of the snow cover removal by heavy winds, active in the region during winter, and subsequently lower ground temperatures, hence only lichens are omnipresent at these surfaces, associated with *Arctous alpinus* (L.) Nied. and crowberry (*Empetrum nigrum* L.). Locally, moss-dominated mires with *Sphagnum* spp. occur in topographical depressions in the interfluve belt. Gentle slopes are covered by creeping willow (*Salix arctica* Pall.); dwarf birch (*Betula nana* L.); and northern Labrador tea (*Rhododendron tomentosum* Harmaja); and small deciduous shrubs and plants: blueberry (*Vaccinium cyanococcus* Rydb.), blackberry (*Vaccinium myrtillus* L.), lingonberry (*Vaccinium vitis-idaea* L.), bearberry (*Arctostaphylos uva-ursi* Spreng.), and crowberry. Water track valleys are willow dominated, mainly *Salix phylicifolia* L., associated with *Equisetum arvense* L. and *Carex* spp.

The studied territory occupies the interfluve area between the Vorkuta River and its right tributary, the Lyok-Vorkuta River, and the valley slope inclined toward the east, to the Vorkuta River valley (Figure 1b). This slope is dissected by three first-order water track valleys, one of which was studied in detail in its middle and lower reach, downstream from the railroad crossing (Figure 1c). The water track valley is oriented west to east. Its headwaters are connected, via a network of wet depressions and mires, to the headwaters of all major neighboring water tracks, so that no interfluve exists between the water track systems draining in different directions. For this reason, the basin area of the studied water track, is estimated, with significant uncertainty, to be around 0.901 <sup>±</sup> 0.055 km<sup>2</sup> . This uncertainty is not related to the digital elevation model (DEM) resolution but reflects the fluvial network structure, as seen in Figure 1c. No clear line separating the two neighboring water track catchments can be easily drawn, because the flow direction can hardly be determined in the interconnected polygonal network in intertrack spaces.

The regional climate is subarctic, summer is short and cool, and winter is long and cold, lasting over eight months from October to May (Table 1). At the same time, the period without negative daily temperatures is only 70 days in an average year. Mean annual daily temperature observed at Vorkuta meteorological station (N 67◦29.520 , E 63◦58.530 ) is −5.6 ◦C. Precipitation is approximately 530 mm, of which from 50% to 70% falls as snow, which can occur at any month of the year.


**Table 1.** Mean monthly air temperature (T, ◦C) and precipitation (P, mm), observed at Vorkuta meteorological station (1927–2019).

Permafrost is continuous, with thickness varying from 50 to 100 m, and mean annual ground temperatures between −0.5 and −1.0 ◦C, at zero annual amplitude depth, around 12 m. A residual thaw layer occurs annually between the base of the seasonally freezing layer, ca. 2 to 3 m, and the top of permafrost at depths of 3.5 to 4.0 m. The prevailing permafrost cryostructure is massive; within the active layer, distinct traces of melted segregation ice lenses can be found–small lenticular unconformities parallel to the ground surface, which were filled by ice during winter and melted later in summer. Important cryogenic processes include thermokarst and fluvial thermal erosion. Frost boils are common cryogenic features, mostly occurring in a narrow belt surrounding the water track valleys.

#### **3. Materials and Methods**

Field observations were performed from the 5–19 September 2017. Water samples for the analyses of stable water isotope, <sup>2</sup>H and <sup>18</sup>O, were collected regularly, once in 2–3 days, from the Vorkuta River and the stream at the water track thalweg, draining into the river near the base camp at the Khanovey station, near its mouth (Figure 1c). Multiple samples were taken along the water track thalweg in its lower reach, downstream from the railroad crossing. Several samples were taken from natural hollows, circular depressions from 0.2 to 0.4 m in diameter, occurring on the ground surface in the water track valley and on slopes. Several soil pits, 40 to 90 cm deep, were dug at various locations in the water track valley and in the open tundra to sample shallow subsurface groundwater. Rainwater was sampled in Vorkuta, from an intense rain shower occurring on 11 September 2017. Subpermafrost groundwater was sampled from an artesian well, No.39-B, near the Khanovey railway station, feeding from a regional aquifer at a depth of around 80 m. Water samples were collected in 15 mL conical centrifuge tubes, sealed with Parafilm©, and stored at 4 ◦C before they were transported to the lab.

Stable water isotope samples (*n* = 35) were analyzed by multiflow-isotope ratio mass spectrometry (MF-IRMS, with instruments from Elementar, Germany) at SHIVA Isotopic Platform, EcoLab, Toulouse, France, in December 2017. The internal standard was Vienna Standard Mean Ocean Water (VSMOW2), and resulting values were expressed in δ notation relative to this standard [43]. The analytical precision of the method is ±0.1‰ for δ <sup>18</sup>O, and ±1.0‰ for δ <sup>2</sup>H. Each water sample was measured in duplicate and averaged, so each δ <sup>18</sup>O/δ <sup>2</sup>H value presented in the paper is a mean value. Deuterium excess (*d*-excess, or *d*ex, ‰) was calculated as *d*ex = δ <sup>2</sup><sup>H</sup> <sup>−</sup> <sup>8</sup>·<sup>δ</sup> <sup>18</sup>O.

Dissolved organic carbon (DOC) samples (*n* = 33) were collected in 20 mL LDPE bottles, pre-washed with weak sulfuric acid and rinsed with MilliQ water. All samples were acidified in the field directly after collection with two drops of 30% H2SO<sup>4</sup> to suppress biological activity and stored at 4 ◦C until transported to the lab. The analyses were carried out at VNIRO (All-Russian Research Institute of Fisheries and Oceanography, Moscow), on a Shimadzu TOC-Vcph analyzer (Shimadzu, Japan), and are precise to ±0.1 mg/L.

Electric resistivity tomography (ERT) surveys were performed with 'SKALA-64' ERT station (NEMPHIS, Russia), using a combined three-electrode protocol (AMN-MNB) with a remote electrode installed at a distance of 600–800 m from the profile, which was considered as infinity. The ERT surveys were done at currents between 35 and 70 mA, and survey data were treated with Res2dInv (Geotomo Software, Malaysia, https://www.geotomosoft. com/, accessed on 13 July 2021) and X2ipi (Aleksey Bobachev, Moscow State University, http://x2ipi.ru/en, accessed on 13 July 2021) software.

Three shallow boreholes around the study area were instrumented with ground temperature sensors and dataloggers: one borehole was equipped with a Hobo 2-sensor datalogger at 0.5 and 2.0 m depth, and two boreholes had Hobo 4-sensor dataloggers at 0.5, 1.0, 2.0, and 5.0 m depths. The ground temperature readings are accurate to ±0.1 ◦C. Ground temperatures from −0.2 to −0.5 ◦C at 5 m depth were observed in two of them, and around +0.8 ◦C in the third borehole, in the water track thalweg.

#### **4. Results**

#### *4.1. Water Stable Isotopes*

The closest stations of the IAEA Global Network for Isotopes in Precipitation (GNIP) [44] network, providing baseline data for the isotopic composition of meteoric waters, are located in Pechora and Salekhard, several hundreds of kilometers from the studied location (Table 2). Their data show the effect of different vapor sources and the transformation of isotopic composition of regional precipitation as air masses cross the Polar Ural Mountains. At the Pechora site, more than 800 km to the SW from Vorkuta (N 65◦07.300 , E 57◦08.980 ), the local meteoric water line (LMWL) is close to the global one (GMWL). In Salekhard (N 66◦32.200 , E 66◦37.940 ), ca. 300 km to the SE from Vorkuta on the eastern side of the Ural Mountains in the Ob' River estuary, the LMWL plots below the GMWL with an intercept *b* = 1.83 (±1.79).

**Table 2.** Isotopic composition of the major water sources in the studied region and at closest GNIP stations.


 *n* is number of samples; where separated by a slash, first value refers to number of δ O samples, second value, to the number of δ H and *d*-excess values.

The local meteoric water line (LMWL), plotted using all data except the most enriched samples subject to evaporative loss, is close to the global one (GMWL) and plots slightly above GMWL (Figure 3). The LMWL equation is

$$
\delta^2 \mathcal{H} = 7.65 \times \delta^{18} \mathcal{O} + 9.8 \text{ (\%\)},
\tag{1}
$$

The single precipitation event during the field campaign occurred on 10 September 2017 and was sampled in Vorkuta city (N 67◦28.920 , E 64◦01.480 ), about 30 km to the northeast from the Khanovey study site. It shows a more depleted isotopic signature compared to long-term September averages for Pechora (Table 2), and a relatively high *d*-excess value, evidencing significant kinetic fractionation and distant moisture sources.

The local evaporation line (Figure 3) is drawn through a rain sample point at the bottom left of the plot, and across the field of points corresponding to samples with a high degree of evaporative loss. An evaporation effect is clear in a thermokarst lake sample that has a negative *d*-excess value and isotopically enriched composition, as well as in several samples from microtopographical hollows, where the *d*-excess was low positive (Table 2). The slope of the local evaporation line is around 5.0, when connecting rainfall to the highly enriched thermokarst lake point, and around 6.2 when only hollows and soil pits are considered; the averaged evaporation line, shown on Figure 3, yields a slope of 5.7. slope of 5.7.

slightly above GMWL (Figure 3). The LMWL equation is

**Figure 3.** Local meteoric water line for water samples, collected in the Khanovey study area, in relation to the global MWL and the potential evaporation line. **Figure 3.** Local meteoric water line for water samples, collected in the Khanovey study area, in relation to the global MWL and the potential evaporation line.

The local meteoric water line (LMWL), plotted using all data except the most enriched samples subject to evaporative loss, is close to the global one (GMWL) and plots

The single precipitation event during the field campaign occurred on 10 September 2017 and was sampled in Vorkuta city (N 67°28.92′, E 64°01.48′), about 30 km to the north-east from the Khanovey study site. It shows a more depleted isotopic signature compared to long-term September averages for Pechora (Table 2), and a relatively high *d*-excess value, evidencing significant kinetic fractionation and distant moisture sources. The local evaporation line (Figure 3) is drawn through a rain sample point at the bottom left of the plot, and across the field of points corresponding to samples with a high degree of evaporative loss. An evaporation effect is clear in a thermokarst lake sample that has a negative *d*-excess value and isotopically enriched composition, as well as in several samples from microtopographical hollows, where the *d*-excess was low positive (Table 2). The slope of the local evaporation line is around 5.0, when connecting rainfall to the highly enriched thermokarst lake point, and around 6.2 when only hollows and soil pits are considered; the averaged evaporation line, shown on Figure 3, yields a

δ2H = 7.65 × δ18O + 9.8 (‰), (1)

The isotopic composition of sampled rivers and streams departs significantly from the GMWL toward higher *d*-excess and δ18O (Table 2 and Figure 3) and plots uniformly above the other end member points on the *d*-excess–δ18O diagram (Figure 4), evidencing various degrees of water fractionation in surface and shallow subsurface compartments [23]. It is presumably closer to the average isotopic composition of July and August rains, The isotopic composition of sampled rivers and streams departs significantly from the GMWL toward higher *d*-excess and δ <sup>18</sup>O (Table 2 and Figure 3) and plots uniformly above the other end member points on the *d*-excess–δ <sup>18</sup>O diagram (Figure 4), evidencing various degrees of water fractionation in surface and shallow subsurface compartments [23]. It is presumably closer to the average isotopic composition of July and August rains, in accordance with the lower relative humidity of these summer months. *Hydrology* **2021**, *8*, x FOR PEER REVIEW 8 of 16

in accordance with the lower relative humidity of these summer months.

**Figure 4.** The *d*-excess–δ18O diagram of sampled water bodies in the Khanovey study site. **Figure 4.** The *d*-excess–δ <sup>18</sup>O diagram of sampled water bodies in the Khanovey study site.

DOC content was in hollows, while for other compartments it was more stable.

**Figure 5.** The dissolved organic carbon (DOC) concentrations across different water bodies in the

The ERT surveys show a highly diverse distribution of high- and low-resistivity grounds in the studied sections (Figure 6), in most cases correlating closely with the hy-

Below the point cloud comprised mostly of surface waters, two separate end members are plotted in the opposite corners of the *d*-excess–δ18O diagram (Figure 4). In the top left corner, a single sample of deep subpermafrost groundwater is plotted, depleted in heavy 18O isotope, with high *d*-excess and isotopic signature consistent with that of a confined groundwater aquifer [45]. In the bottom right corner, samples that have undergone substantial evaporative loss are plotted, as discussed above, including the lake and several hollows (Figure 4). Below the point cloud comprised mostly of surface waters, two separate end members are plotted in the opposite corners of the *d*-excess–δ <sup>18</sup>O diagram (Figure 4). In the top left corner, a single sample of deep subpermafrost groundwater is plotted, depleted in heavy <sup>18</sup>O isotope, with high *d*-excess and isotopic signature consistent with that of a confined groundwater aquifer [45]. In the bottom right corner, samples that have undergone substantial evaporative loss are plotted, as discussed above, including the lake and several hollows (Figure 4).

The DOC concentration measured at the Khanovey study site, is relatively low, averaging 10.4 ± 5 mg/L across the dataset. However, it is highly variable across water body types (Figure 5). In general, the highest DOC content is observed in standing water, i.e., bogs and hollows, whilst it is significantly lower in streams and rivers. The most variable

Khanovey region.

*4.3. Electric Resistivity Tomography (ERT)* 

drographic network's features.

#### *4.2. Dissolved Organic Carbon 4.2. Dissolved Organic Carbon*

and several hollows (Figure 4).

**Figure 4.** The *d*-excess–δ

The DOC concentration measured at the Khanovey study site, is relatively low, averaging 10.4 ± 5 mg/L across the dataset. However, it is highly variable across water body types (Figure 5). In general, the highest DOC content is observed in standing water, i.e., bogs and hollows, whilst it is significantly lower in streams and rivers. The most variable DOC content was in hollows, while for other compartments it was more stable. The DOC concentration measured at the Khanovey study site, is relatively low, averaging 10.4 ± 5 mg/L across the dataset. However, it is highly variable across water body types (Figure 5). In general, the highest DOC content is observed in standing water, i.e., bogs and hollows, whilst it is significantly lower in streams and rivers. The most variable DOC content was in hollows, while for other compartments it was more stable.

<sup>18</sup>O diagram of sampled water bodies in the Khanovey study site.

Below the point cloud comprised mostly of surface waters, two separate end members are plotted in the opposite corners of the *d*-excess–δ18O diagram (Figure 4). In the top left corner, a single sample of deep subpermafrost groundwater is plotted, depleted in heavy <sup>18</sup>O isotope, with high *d*-excess and isotopic signature consistent with that of a confined groundwater aquifer [45]. In the bottom right corner, samples that have undergone substantial evaporative loss are plotted, as discussed above, including the lake

*Hydrology* **2021**, *8*, x FOR PEER REVIEW 8 of 16

**Figure 5.** The dissolved organic carbon (DOC) concentrations across different water bodies in the Khanovey region. **Figure 5.** The dissolved organic carbon (DOC) concentrations across different water bodies in the Khanovey region.

#### *4.3. Electric Resistivity Tomography (ERT) 4.3. Electric Resistivity Tomography (ERT)*

The ERT surveys show a highly diverse distribution of high- and low-resistivity grounds in the studied sections (Figure 6), in most cases correlating closely with the hydrographic network's features. The ERT surveys show a highly diverse distribution of high- and low-resistivity grounds in the studied sections (Figure 6), in most cases correlating closely with the hydrographic network's features. *Hydrology* **2021**, *8*, x FOR PEER REVIEW 9 of 16

**Figure 6.** Electrical resistivity tomography profiles across the transects T1 (top) and T2 (bottom), Scheme 1. c for geographical reference. Vertical axis is altitude, in m a.s.l.; the lateral distance from the starting point, in m, is given along the profile. **Figure 6.** Electrical resistivity tomography profiles across the transects T1 (**top**) and T2 (**bottom**), Figure 1c for geographical reference. Vertical axis is altitude, in m a.s.l.; the lateral distance from the starting point, in m, is given along the profile.

Geophysical evidence suggests a relatively thick and steady permafrost, exceeding 30 m, persisting under raised non-dissected tundra patches and peat plateaus; a shallow and thin permafrost under a recent mire on the left side of T1 transect (see Figures 1c and 6); and open taliks under bogs, fens, or hummocky depressions. In all cases, the talik walls are subvertical, with thaw bulbs slightly expanding downward. This underscores the vulnerability of contemporary permafrost and also suggests significant groundwater Geophysical evidence suggests a relatively thick and steady permafrost, exceeding 30 m, persisting under raised non-dissected tundra patches and peat plateaus; a shallow and thin permafrost under a recent mire on the left side of T1 transect (see Figures 1c and 6); and open taliks under bogs, fens, or hummocky depressions. In all cases, the talik walls are subvertical, with thaw bulbs slightly expanding downward. This underscores the vulnerability of contemporary permafrost and also suggests significant groundwater circulation in the talik zones [29].

The available data on water stable isotopes and their variability, presented in Table 2 and Figure 7, allow a generalized description of hydrological connectivity at a scale of a

Elevated tundra patches underlain by sandy loams were detached from the hydrological system at the time of observations. The upper part of their soil profile, above 0.6–0.9 m, was unsaturated. The saturated zone groundwater was sampled in a soil pit and had δ18O = −10.55‰, δ2H = −81.52‰, and *d*ex = 2.9‰ (Figure 4), showing signs of evaporative loss and representing the isotopically heavier precipitation of preceding summer months. Other soil pits were opened in the water track and secondary drainage network thalwegs and contained water, intermediate between sampled September rainfall and stream/river water, with δ18O = −13.3‰…−13.4‰, δ2H = −89.9‰…−95.2‰, and *d*ex = 9.9‰…13.1‰. One soil pit sample appeared close to the subpermafrost groundwater, potentially evidencing upward groundwater migration and discharge through the

talik zones. This will be discussed in Section 5.4 in more detail.

*5.1. Hydrological Connectivity at a Catchment Scale* 

circulation in the talik zones [29].

minor water track catchment.

**5. Discussion** 

#### **5. Discussion**

#### *5.1. Hydrological Connectivity at a Catchment Scale*

The available data on water stable isotopes and their variability, presented in Table 2 and Figure 7, allow a generalized description of hydrological connectivity at a scale of a minor water track catchment. *Hydrology* **2021**, *8*, x FOR PEER REVIEW 10 of 16

**Figure 7.** Boxplots of stable water isotope and *d*-excess values in the water track catchment in the Khanovey study site; median, 25%, and 75% probability of exceedance shown as bars, highest and lowest as whiskers, and outliers as separate points (see also Table 2). **Figure 7.** Boxplots of stable water isotope and *d*-excess values in the water track catchment in the Khanovey study site; median, 25%, and 75% probability of exceedance shown as bars, highest and lowest as whiskers, and outliers as separate points (see also Table 2).

On the tundra surface, multiple hollows were sampled, which were expected to show signs of evaporative loss. Biological fractionation from photosynthetic plant and algae activity and preferential evaporation of lighter oxygen isotope 16O could potentially result in enrichment in heavy oxygen isotope, but this process is not expected to play an important role in late autumn, when the plants are in the end of their annual lifecycle. Surprisingly, only three out of six sampled hollows were evaporative basins, with mean δ18O = −10.6‰ ± 0.6‰ and *d*ex = 3.23‰ ± 0.39‰. The other three hollows contained water that was isotopically similar to stream and river water, with δ18O = −12.4‰ ± 0.4‰, δ2H = −87.1‰ ± 5.1‰, and *d*ex = 11.7‰ ± 3‰, assuming their direct connection to shallow subsurface groundwater and the water track stream. This connection can be assured by the transmissivity feedback effect [46], occurring widely in the organic topsoil outside the Elevated tundra patches underlain by sandy loams were detached from the hydrological system at the time of observations. The upper part of their soil profile, above 0.6–0.9 m, was unsaturated. The saturated zone groundwater was sampled in a soil pit and had δ <sup>18</sup>O = <sup>−</sup>10.55‰, <sup>δ</sup> <sup>2</sup>H = <sup>−</sup>81.52‰, and *<sup>d</sup>*ex = 2.9‰ (Figure 4), showing signs of evaporative loss and representing the isotopically heavier precipitation of preceding summer months. Other soil pits were opened in the water track and secondary drainage network thalwegs and contained water, intermediate between sampled September rainfall and stream/river water, with δ <sup>18</sup>O = <sup>−</sup>13.3‰ . . . <sup>−</sup>13.4‰, <sup>δ</sup> <sup>2</sup>H = <sup>−</sup>89.9‰ . . . <sup>−</sup>95.2‰, and *d*ex = 9.9‰ . . . 13.1‰. One soil pit sample appeared close to the subpermafrost groundwater, potentially evidencing upward groundwater migration and discharge through the talik zones. This will be discussed in Section 5.4 in more detail.

elevated tundra patches on better drained soils. Heavily silted and clayey soils, observed locally in the catchment, are saturated and highly thixotropic. At rest, a stiffened thixotropic layer is expected to show barrier functions for water migration [47], barring water infiltration, leading to quick saturation of the overlying soil, and initiating rapid water drainage through the organic topsoil. Multiple open-surface bogs in the drainage network depressions exist in the region, and one of them was sampled for water stable isotopes. This water object shows an isotopic signature intermediate between the hollows and the closest water track stream, which is consistent with the local buffer role of such bogs on the microscale slopes [48], transforming pluvial runoff and isotopic signal and conveying it to streamflow. On the tundra surface, multiple hollows were sampled, which were expected to show signs of evaporative loss. Biological fractionation from photosynthetic plant and algae activity and preferential evaporation of lighter oxygen isotope <sup>16</sup>O could potentially result in enrichment in heavy oxygen isotope, but this process is not expected to play an important role in late autumn, when the plants are in the end of their annual lifecycle. Surprisingly, only three out of six sampled hollows were evaporative basins, with mean δ <sup>18</sup>O = <sup>−</sup>10.6‰ <sup>±</sup> 0.6‰ and *<sup>d</sup>*ex = 3.23‰ <sup>±</sup> 0.39‰. The other three hollows contained water that was isotopically similar to stream and river water, with δ <sup>18</sup>O = <sup>−</sup>12.4‰ <sup>±</sup> 0.4‰, δ <sup>2</sup>H = <sup>−</sup>87.1‰ <sup>±</sup> 5.1‰, and *<sup>d</sup>*ex = 11.7‰ <sup>±</sup> 3‰, assuming their direct connection to shallow subsurface groundwater and the water track stream. This connection can be assured by the transmissivity feedback effect [46], occurring widely in the organic topsoil

follow the coevolution of their isotopic composition (Figure 8). In the major stream, the Vorkuta River, it was gradually shifting toward lighter δ18O and δ2H values, reflecting the autumn flow recession and increasing groundwater input and also a potential change in rainwater isotopic signature (Figure 8a,b). No similar variation in δ18O values was ob-

*5.2. Temporal Evolution of Water Isotopic Composition*

outside the elevated tundra patches on better drained soils. Heavily silted and clayey soils, observed locally in the catchment, are saturated and highly thixotropic. At rest, a stiffened thixotropic layer is expected to show barrier functions for water migration [47], barring water infiltration, leading to quick saturation of the overlying soil, and initiating rapid water drainage through the organic topsoil. Deuterium excess values were highly variable with time in the Vorkuta River, jumping from 12‰ to 18‰ over the time span of several days (Figure 8b). Runoff inputs from sources adjacent to the observation point were responsible for the first peak, which coincided with the same *d*-excess peak in the water track stream (Figure 8d), while more distant sources could potentially contribute to the second peak. The Vorkuta River has a large basin, around 4000 km2 at the sampling site, and smoother variations in water

served in the water track stream, but at the same time, a progressive depletion in δ2H was

Multiple open-surface bogs in the drainage network depressions exist in the region, and one of them was sampled for water stable isotopes. This water object shows an isotopic signature intermediate between the hollows and the closest water track stream, which is consistent with the local buffer role of such bogs on the microscale slopes [48], transforming pluvial runoff and isotopic signal and conveying it to streamflow. chemistry are expected. However, the water isotopic signal of the river could be reset by human activities upstream, i.e., the dam of the Heat Production Station No. 2 upstream of Vorkuta, where the flow is almost ceasing during low-flow periods. In this case, our samples reflect the changes in the isotopic composition of major Vorkuta River tributaries, including Ayach-Yaga River and Yun'-Yaga River, as well as numerous minor trib-

#### *5.2. Temporal Evolution of Water Isotopic Composition* utaries and water tracks. The low temporal resolution of the survey may add to the lack of smooth variations on the graph.

limb of the stormflow hydrograph [49].

recorded (Figure 8c).

*Hydrology* **2021**, *8*, x FOR PEER REVIEW 11 of 16

Paired samples were taken in the Vorkuta River and at the water track mouth to follow the coevolution of their isotopic composition (Figure 8). In the major stream, the Vorkuta River, it was gradually shifting toward lighter δ <sup>18</sup>O and δ <sup>2</sup>H values, reflecting the autumn flow recession and increasing groundwater input and also a potential change in rainwater isotopic signature (Figure 8a,b). No similar variation in δ <sup>18</sup>O values was observed in the water track stream, but at the same time, a progressive depletion in δ <sup>2</sup>H was recorded (Figure 8c). A DOC peaks follows the *d*-excess peak with a three-day lag in the water track stream (Figure 8d). The origin of these peaks is unclear, as is their causal relation. We may hypothesize that the *d*-excess peak at a small catchment scale could have been produced by rain events, when first, pre-event water yielding high *d*-excess values is pushed from the water storage within the catchment. The DOC-rich event waters from the subsurface compartment reach the fluvial network several days later, inducing hysteretic behavior of DOC, when its values are lower on the rising limb and higher on the falling

**Figure 8.** Temporal evolution of isotopic composition (**a**,**c**), *d*-excess and dissolved organic carbon (DOC) concentrations (**b**,**d**) in the Vorkuta River (**a**,**b**) and the water track (**c**,**d**). Numbers denote changes in variables shown (1) on the left axis and (2) on the right axis of each graph. **Figure 8.** Temporal evolution of isotopic composition (**a**,**c**), *d*-excess and dissolved organic carbon (DOC) concentrations (**b**,**d**) in the Vorkuta River (**a**,**b**) and the water track (**c**,**d**). Numbers denote changes in variables shown (1) on the left axis and (2) on the right axis of each graph.

Deuterium excess values were highly variable with time in the Vorkuta River, jumping from 12‰ to 18‰ over the time span of several days (Figure 8b). Runoff inputs from sources adjacent to the observation point were responsible for the first peak, which coincided with the same *d*-excess peak in the water track stream (Figure 8d), while more

distant sources could potentially contribute to the second peak. The Vorkuta River has a large basin, around 4000 km<sup>2</sup> at the sampling site, and smoother variations in water chemistry are expected. However, the water isotopic signal of the river could be reset by human activities upstream, i.e., the dam of the Heat Production Station No. 2 upstream of Vorkuta, where the flow is almost ceasing during low-flow periods. In this case, our samples reflect the changes in the isotopic composition of major Vorkuta River tributaries, including Ayach-Yaga River and Yun'-Yaga River, as well as numerous minor tributaries and water tracks. The low temporal resolution of the survey may add to the lack of smooth variations on the graph.

A DOC peaks follows the *d*-excess peak with a three-day lag in the water track stream (Figure 8d). The origin of these peaks is unclear, as is their causal relation. We may hypothesize that the *d*-excess peak at a small catchment scale could have been produced by rain events, when first, pre-event water yielding high *d*-excess values is pushed from the water storage within the catchment. The DOC-rich event waters from the subsurface compartment reach the fluvial network several days later, inducing hysteretic behavior of DOC, when its values are lower on the rising limb and higher on the falling limb of the stormflow hydrograph [49].

This effect is significantly better expressed and occurs twice in the larger Vorkuta River runoff most probably because its large catchment effectively integrates hydrological signals from numerous rain events in its different parts. It is important to note that the *d*-excess variability during the peak is comparable in the small water track and in the larger Vorkuta River, while the DOC increase is significantly higher in the smaller stream.

#### *5.3. DOC Export from the Catchment*

The DOC concentrations are relatively low and only slightly vary within each water object type, except the small water track stream (Figure 5). Its variability in pits and hollows follows the variability detected in the isotopic signature. The hollows subject to evaporative loss and disconnected from subsurface flow (Figure 4), also yielded lower DOC values, 11.8 ± 1.7 mg/L on average. The pits and hollows connected to fast subsurface flow had DOC concentrations about 50% higher than disconnected ones (17.5 ± 5.3 mg/L), showing the important role they play in the DOC lateral fluxes. Lower DOC content in disconnected hollows may reflect multiple processes, from photodegradation to aerobic microbial degradation, with complex interaction between them [50].

The water track stream discharge is estimated to vary between 2 and 4 L/s, this estimate is based on occasional hydrologic observations from previous years, performed around the same dates in early September and in comparable weather conditions. Using the average DOC concentration in the water track, 9.45 mg/L, and an average daily discharge of around 3 L/s, the daily DOC export from the water track catchment for an average day in late autumn equals 2.4 kg/km<sup>2</sup> .

#### *5.4. Subpermafrost Groundwater Input*

An 80 m deep artesian well sampled during the course of the study intercepts subpermafrost groundwater with a distinct isotopic signature (Table 2, Figure 4). Isotopically similar water, highly depleted in <sup>18</sup>O (δ <sup>18</sup>O = <sup>−</sup>15.8‰) and with *<sup>d</sup>*-excess above 15‰, was sampled from a soil pit near the borehole located in a water track thalweg, between 48 and 72 m from the start point of the ERT T2 profile (see Figures 1c and 6 for reference). This soil pit exposed sandy loams down to 0.9 m and, when opened, remained almost dry for 2 to 3 h, until the water level at the pit slowly settled at a 0.6 m depth from the surface. The water source in this soil pit, located in the topographical depression can be either from the surrounding elevated tundra patches or from the subsurface compartment, including intra- and subpermafrost groundwater. The ERT data suggests the existence of a thin relict permafrost layer, underlain by non-frozen ground and interrupted by subvertical talik zones. We suppose that water observed in this soil pit originates from the deep

groundwater aquifer, close in isotopic composition to the groundwater from the artesian well (see Table 2).

The elevated tundra patches are shown to yield significantly isotopically heavier water, and we cannot propose any viable mixing mechanism or end member that could change the isotopic composition of this water in a way that brings it close to the discussed soil pit water. The occasional evaporation of the water sample during storage would lead to heavier isotopic composition of this sample, and not a depleted one, so this possibility can also be ruled out. From this, we can conclude that, even if the single sample is not an entirely convincing evidence, when combined with geophysical data, it is highly likely that there exists an important connection between shallow subsurface groundwater in water track valleys with deeper groundwater aquifers through upward water migration.

#### *5.5. Water Tracks: Current and Future Development*

Water track are widespread periglacial features of the Russian Arctic, and they dominate the slope topography in studied region of north-European Arctic Russia. Unlike Western Siberia and north-eastern Russia, the water track drainage network is highly developed, with deeply incised water track valleys overgrown by willows, and intertrack areas with a polygonally shaped secondary drainage network resembling thermokarst-related patterns of permafrost degradation.

Hydrological connectivity in the water track landscape controls surface and subsurface water fluxes and is, as such, an important ecohydrological factor of the vegetation community structure. Linear thermokarst features of the Bylot Island, Nunavut, were recently reported to control tundra landscapes by promoting the transition from wet to mesic tundra vegetation in areas adjacent to thermo-erosional gullies [51]. At the Khanovey study site, and generally in the north-European Arctic Russia, water tracks appear to play a comparable role, sustaining drier habitats on the tops of elevated tundra patches and moist habitats in the water track thalwegs.

Contemporary climate in the Khanovey region is capable of maintaining wet tundra habitats, as evidenced by the abundance of fens in the water track channels on the interfluves of major rivers of the studied region, at slope summits, shoulders, and partially on backslopes. At these positions, slope steepness is insufficient to rapidly convey water downslope even through water tracks, which leads to the persistent waterlogging of such locations. At footslope positions with steeper slopes, water tracks are increasingly efficient to drain the adjacent tundra. Surface subsidence accompanied by permafrost degradation (linear thermokarst) increases the local height difference between the elevated tundra patches and the water track thalwegs and enhances drainage of these patches. As a result, the intertrack elevated tundra hosts mostly xeric communities, with dwarf shrubs and ericaceous species, i.e., red bearberry, crowberry, Labrador tea, and dwarf birch. Water tracks lower the groundwater table and drain the surrounding tundra effectively enough to maintain xeric habitats.

The water track network in the study site and adjacent territories is still developing, both vertically and laterally. The studied water track valley (Figure 1c) is supposedly freshly incised, because its longitudinal profile is not in equilibrium and hosts several waterfalls, up to ca. 1.5 m high. The evolution of the drainage network at the intertrack surface continues. The dominant local slope is directed toward the water track valley, rather than toward the base level of the Vorkuta River. Because of this, we observe the development of new linear depressions on the left side of the water track valley, with depth between 0.15 and 0.30 m, associated with terrain highly disturbed by hummocks.

Climate change is expected to alter the functioning of the water track system of the region, but the direction of future change is unclear. Permafrost degradation is expected to promote the gradual lateral thawing of elevated tundra patches now frozen [52]. Visual inspection of satellite imagery shows that there is significant difference between water track networks on the south- and north-facing slopes; hence, increasing insolation and air temperature in the future climate are expected to produce a significant geomorphic response.

#### **6. Conclusions**

The hydrological snapshot of a small tundra catchment in North-European Russia provides several insights into the hydrological connectivity within its limits. Elevated tundra patches appear to be disconnected from the stream, while the slopes and the riparian zone contribute actively, though locally, to the stream runoff. Minor hollows are found to be either connected to the shallow subsurface runoff or disconnected from it. The difference in connectivity is traceable via *d*-excess; the disconnected hollows act as evaporative basins and have lower *d*-excess values. The connected hollows also serve as important DOC flux conveyors, with DOC concentrations up to 50% higher than in disconnected hollows.

**Author Contributions:** Conceptualization and methodology, N.T.; field investigation, N.T., V.I., D.S., P.K. and O.K.; writing—original draft preparation, N.T.; writing—review and editing, N.T., D.S. and P.K.; visualization, N.T. and O.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was conducted with support from RuNoCore: Russian–Norwegian researchbased education in cold region engineering (#CPRU-2017/10015), funded by the Norwegian Center for International Cooperation in Education (SIU). The manuscript was prepared under the implementation of the state task and the research plan of the IEG RAS, Research Topic AAAA-A19-119021190077- 6, and the Melnikov Permafrost Institute SB RAS, Research Topic AAAA-A20-120111690008-9.

**Data Availability Statement:** The analyzed dataset is available through Figshare, an open repository, and accessible by DOI:10.6084/m9.figshare.14502003.

**Acknowledgments:** The authors acknowledge the support from the Northern Railroad services in arranging the field surveys.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**

