*Article* **Hydrochemistry of Medium-Size Pristine Rivers in Boreal and Subarctic Zone: Disentangling Effect of Landscape Parameters across a Permafrost, Climate, and Vegetation Gradient**

**Oleg S. Pokrovsky 1,\*, Artem G. Lim 2, Ivan V. Krickov 2, Mikhail A. Korets 3, Liudmila S. Shirokova 1,4 and Sergey N. Vorobyev <sup>2</sup>**


**Abstract:** We studied two medium size pristine rivers (Taz and Ket) of boreal and subarctic zone, western Siberia, for a better understanding of the environmental factors controlling major and trace element transport in riverine systems. Our main objective was to test the impact of climate and land cover parameters (permafrost, vegetation, water coverage, soil organic carbon, and lithology) on carbon, major and trace element concentration in the main stem and tributaries of each river separately and when considering them together, across contrasting climate/permafrost zones. In the permafrost-bearing Taz River (main stem and 17 tributaries), sizable control of vegetation on element concentration was revealed. In particular, light coniferous and broadleaf mixed forest controlled DOC, and some nutrients (NO2, NO3, Mn, Fe, Mo, Cd, Ba), deciduous needle-leaf forest positively correlated with macronutrients (PO4, Ptot, Si, Mg, P, Ca) and Sr, and dark needle-leaf forest impacted Ntot, Al, and Rb. Organic C stock in the upper 30–100 cm soil positively correlated with Be, Mn, Co, Mo, Cd, Sb, and Bi. In the Ket River basin (large right tributary of the Ob River) and its 26 tributaries, we revealed a correlation between the phytomass stock at the watershed and alkaline-earth metals and U concentration in the river water. This control was weakly pronounced during high-water period (spring flood) and mostly occurred during summer low water period. Pairwise correlations between elements in both river systems demonstrated two group of solutes—(1) positively correlated with DIC (Si, alkalis (Li, Na), alkaline-earth metals (Mg, Ca, Sr, Ba), and U), this link originated from groundwater feeding of the river when the labile elements were leached from soluble minerals such as carbonates; and (2) elements positively correlated with DOC (trivalent, tetravalent, and other hydrolysates, Se and Cs). This group reflected mobilization from upper silicate mineral soil profile and plant litter, which was strongly facilitated by element colloidal status, notably for low-mobile geochemical tracers. The observed DOC vs DIC control on riverine transport of low-soluble and highly mobile elements, respectively, is also consistent with former observations in both river and lake waters of the WSL as well as in soil waters and permafrost ice. A principal component analysis demonstrated three main factors potentially controlling the major and TE concentrations. The first factor, responsible for 26% of overall variation, included aluminum and other low mobile trivalent and tetravalent hydrolysates, Be, Cr, Nb, and elements strongly complexed with DOM such as Cu and Se. This factor presumably reflected the presence of organo-mineral colloids, and it was positively affected by the proportion of forest and organic C in soils of the watershed. The second factor (14% variation) likely represented a combined effect of productive litter in larch forest growing on carbonate-rich rocks and groundwater feeding of the rivers and acted on labile Na, Mg, Si, Ca, P, and Fe(II), but also DOC, micronutrients (Zn, Rb, Ba), and phytomass at the watershed. Via applying a substituting space for time approach for south-north gradient of studied river basins, we predict that

**Citation:** Pokrovsky, O.S.; Lim, A.G.; Krickov, I.V.; Korets, M.A.; Shirokova, L.S.; Vorobyev, S.N. Hydrochemistry of Medium-Size Pristine Rivers in Boreal and Subarctic Zone: Disentangling Effect of Landscape Parameters across a Permafrost, Climate, and Vegetation Gradient. *Water* **2022**, *14*, 2250. https:// doi.org/10.3390/w14142250

Academic Editor: Maurizio Barbieri

Received: 29 June 2022 Accepted: 14 July 2022 Published: 18 July 2022

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

**Copyright:** © 2022 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/).

climate warming in northern rivers may double or triple the concentration of DIC, Ca, Sr, U, but also increase the concentration of DOC, POC, and nutrients.

**Keywords:** metals; carbon; nutrients; trace elements; landscape; permafrost; river; watershed; boreal

#### **1. Introduction**

Climate change, which is mostly pronounced in high latitudes, strongly impacts the chemistry of rivers and streams and can bring yet unknown consequences on carbon, nutrient, and metal export from land to ocean thus enhancing the retroactive link to climate change drivers [1–4]. There are numerous case studies of climate change impact on surface and groundwater sustainability and hydrochemistry, such as those performed in temperate/subtropical regions [5,6]. However, huge northern unpopulated territories of the world remain poorly covered by a coupled hydrochemical/hydrological approach. This is especially true for pristine permafrost-bearing regions located at high latitudes and containing sizable amount of carbon in the form of peat, organic litter, and vegetation. An example is Western Siberian Lowland (WSL), located in the gradient of climate and permafrost zones within essentially the same lithological background, minimal variations in river runoff and relief and moderate to negligible anthropogenic activity. The interest of the WSL is that it contains huge peat resources and presents rather shallow, essentially discontinuous to sporadic/isolated permafrost, highly vulnerable to thawing [7–10]. For this relatively large territory (2 million km2), extensive studies of small [11–19] and large [20] river dissolved, colloidal and particulate load, chemical composition of soil ice and suprapermafrost waters [21–24], and gaseous regime of rivers and lakes [25–27] have been performed. However, the majority of these studies except that of the Ob River [20] dealt with single site sampling of a given river, without addressing the spatial variability of riverine solutes within a watershed. This is a clear shortcomings of current state of knowledge of western Siberian rivers, because without sufficient spatial coverage, one cannot test the impact of various landscape factors on water hydrochemistry within the same river basin.

In the present study, we used a coupled hydrochemical/landscape (land cover) approach that allows revealing the main environmental factors controlling the hydrochemical composition of the river water. Up to now, such an approach has been efficiently used only for single-point sampling of small rivers [13–19] and one large Siberian River [20] but has not been implemented at the scale of the whole watershed of medium-size rivers. Toward better understanding of land cover control on riverine C, nutrient and metal transport, we chose two medium size rivers of western Siberia, Taz (Swatershed = 150,000 km2) and Ket (Swatershed = 95,000 km2) located within permafrost-bearing forest-tundra/tundra and permafrost-free taiga biomes, respectively. Both rivers drain through similar sedimentary deposits overlayed by peatland and forest/tundra; they are virtually pristine (population density < 1 people km<sup>−</sup>2) and have no industrial or agricultural activity on their watershed. Within the on-going drastic climate change in Siberia, the southern river (Ket) can be considered as an extreme scenario of long-term transformation of the northern River (Taz) in case of complete disappearance of discontinuous permafrost and northward migration of the taiga forest. Therefore, our primary objective in this work was to test the impact of climate and land cover parameters (permafrost, vegetation, water coverage, soil organic carbon, and lithology), on carbon, major and trace element concentration in the main stem and tributaries of each river separately and when considering them together, across contrasting climate/permafrost zones. A flowchart of methodology and approach used in the present study is illustrated in Figure 1. We anticipate that, via employing similar approach for two sufficiently contrasting and yet representative river basin of the WSL, we can provide essential information for coupled land—River C, nutrient and metal fluxes for climate modeling of the region and to foresee future, long-term changes in much large territory

of strongly understudied permafrost-affected Eurasian lowlands, which includes North Siberia, Anabar, Kolyma and Yana-Indigirka lowlands, with an overall territory more than 2 million km2.

**Figure 1.** A flow chart of methodological approach used in the present study.

#### **2. Study Site and Methods**

#### *2.1. Ket River Basin*

We sampled the Ket River main stem and its 26 tributaries during the peak of the spring flood and the end of summer baseflow. The catchment area of sampled watersheds ranged from 94,000 km<sup>2</sup> (Ket's mouth) to 20 km<sup>2</sup> (smallest tributary). The studied watersheds presented rather similar lithology, climate, and vegetation, given its generally west—east orientation, Figure 2; Ref. [28]). The Ket River basin (right tributary of the Ob River) is poorly accessible and can be considered as essentially pristine comprising about 50% forest and 40% of wetlands while having virtually negligible agricultural and forestry activity. The river basin is poorly populated (0.27 person km<sup>−</sup>2) and, compared to left tributaries of the Ob River, lacks road infrastructure due to absence of hydrocarbon exploration activity. We therefore consider this river as a model one for medium size rivers of the boreal zone of western Siberia Lowland covered by forests and bogs. As such, the results on river water hydrochemistry obtained from this watershed can be extrapolated to much larger territory, comprising several million km<sup>2</sup> of permafrost-free taiga forest and bog biome, extending over 3000 km between the left tributaries of the Yenisei River in the east and Finland in the west. The mean annual air temperatures (MAAT) of the Ket River basin is −0.75 ± 0.15 ◦<sup>C</sup> and the mean annual precipitation is 520 ± 20 mm y−1. The lithology is represented by silts and sands with carbonate concretions overlayed by quaternary deposits (loesses, fluvial, glacial, and lacustrine deposits). The dominant soils are podzols in forest areas and histosols in peat bog regions.

**Figure 2.** Map of the two studied river basins (Taz in the north and Ket in the south). The sampling points of the main stem and tributaries are shown by red and black circles, respectively.

#### *2.2. Taz River Basin*

The Taz River sampled in this work included the main stem and its 17 tributaries whose catchment area ranged from 149,000 km2 at the Taz's mouth to 25 km2 of smallest sampled tributary. Alike the Ket River basin, all sampled catchments exhibited similar lithology (clays, silts and sands overlayed by quaternary loesses, fluvial, glacial, and lacustrine deposits), but more contrasting climate and vegetation due to its north to south orientation (Figure 2). The MAAT ranges from −4.6 ◦C in the headwaters (Tolka village) to −5.4 ◦C in the low reaches (Tazovskiy town). The mean annual precipitation is 500 mm y−<sup>1</sup> in the central part of the basin (Krasnoselkup) and 600 mm y−<sup>1</sup> in the low reaches of the Taz River. Given its landscape and climate features, the Taz River basin includes: (1) the upper (southern) part ("headwaters") which comprises 400–800 km upstream of the river mouth, where the permafrost is sporadic to discontinuous and the dominant vegetation is forest-tundra and taiga, and (2) the low reaches (northern) part, located 0–400 km upstream of the mouth where the permafrost is continuous to discontinuous and the dominant biome is tundra and forest-tundra.

#### *2.3. Sampling*

For the Ket River, the peak of annual discharge in 2019 occurred in the end of May whereas in August, the discharge was three to five times lower. From 18 May to 28 May 2019, and from 30 August to 2 September 2019, we started the boat trip in the middle course of the Ket River (Beliy Yar), and moved, first, 475 km upstream the Ket River till its most headwaters, and then moved 834 km downstream till the river mouth. We stopped each 50 km along the Ket River and sampled for major hydrochemical parameters, suspended matter, and bacteria. We also moved several km upstream of selected tributaries to sample for river hydrochemistry.

In the Taz River, the peak of annual discharge in 2019 occurred in the middle of June (5600 m<sup>3</sup> s−1; in August, the discharge was 5 times lower). Therefore, the month of July (average discharge is 2300 m<sup>3</sup> s<sup>−</sup>1; range 1920–3370 m3 s<sup>−</sup>1) can be considered as the end of spring flood period. From 12 July to 16 July 2019 we moved 800 km downstream the Taz River from its most headwaters (Tolka village) to the low reaches (Tazovskiy town). We stopped each ~50 to 100 km of the boat route and sampled for hydrochemical parameters, river suspended matter and total bacterial number of the main stem. For sampling the tributaries, we moved 500–1500 m upstream of the confluence zone.

#### *2.4. Analyses*

The dissolved oxygen (CellOx 325; accuracy of ±5%), specific conductivity (TetraCon 325; ± 1.5%), and water temperature (±0.2 ◦C) were measured in situ at 20 cm depth using a WTW 3320 Multimeter. The pH was measured using portable Hanna instrument via combined Schott glass electrode calibrated with NIST buffer solutions (4.01, 6.86, and 9.18 at 25 ◦C), with an uncertainty of 0.01 pH units. The river water was sampled in pre-cleaned polypropylene bottle from 20–30 cm depth in the middle of the river and immediately filtered through disposable single-use sterile Sartorius filter units (0.45 μm pore size). The first 20 mL of filtrate was discarded. The DOC and dissolved inorganic carbon (DIC) were determined by a Shimadzu TOC-VSCN Analyzer (Kyoto, Japan) with an uncertainty of 3% and a detection limit of 0.1 mg/L. Blanks of Milli-Q water passed through the filters demonstrated negligible release of DOC from the filter material. Specific ultraviolet absorbance (SUVA254) was measured via ultraviolet absorbance at 254 nm using a 10-mm quartz cuvette on a Bruker CARY-50 UV-VIS spectrophotometer and divided by the concentration of DOC (L mg−<sup>1</sup> m<sup>−</sup>1). The concentrations of C and N in suspended material (particulate organic carbon and nitrogen (POC and PON, respectively)) were determined via filtration of 1 to 2 L of freshly collected river water (at the river bank or in the boat) with pre-weighted GFF filters (47 mm, 0.8 μm) and Nalgene 250-mL polystyrene filtration units using a Mityvac® manual vacuum pump. Particulate C and N were measured using catalytic combustion with Cu-O at 900 ◦C with an uncertainty of ≤0.5% using Thermo Flash 2000 CN Analyzer at EcoLab, Toulouse. The samples were analyzed before and after 1:1 HCl treatment to distinguish between total and inorganic C; however, the ratio of Corganic:Ccarbonate in the river suspended matter (RSM) was always above 20 and the

contribution of carbonate C to total C in the RSM was equal in average 0.3 ± 0.3% (2 s.d., *n* = 30).

Total microbial cell concentration was measured after sample fixation in glutaraldehyde, by flow cytometry (Guava® EasyCyteTM systems, Merck, Darmstadt, Germany). Cells were stained using 1 μL of a 10 times diluted SYBR GREEN solution (10000×, Merck), added to 250 μL of each sample before analysis. Particles were identified as cells based on green fluorescence and forward scatter [29].

All analytical approaches used in this study for major and trace element analyses followed methods developed for western Siberian organic-rich surface waters [17,30,31]. The samples were preserved via refrigeration 1 month prior to analysis. Major anion (Cl, SO4 <sup>2</sup>−) concentrations were measured by ion chromatography (HPLC, Dionex ICS 2000) with an uncertainty of 2%. International certified samples ION, PERADE, and RAIN were used for validation of the analyses. Major cations (Ca, Mg, Na, K), Si, and ~40 trace elements were determined with an Agilent iCap Triple Quad (TQ) ICP MS using both argon and helium modes to diminish interferences. About 3 μg L−<sup>1</sup> of In and Re were added as internal standards along with three various external standards. Detection limits of TE were determined as 3× the blank concentration. The typical uncertainty for elemental concentration measurements ranged from 5–10% at 1–1000 μg/L to 10–20% at 0.001–0.1 μg/L. The Milli-Q field blanks were collected and processed to monitor for any potential contamination introduced by our sampling and handling procedures. The SLRS-6 (Riverine Water Reference Material for Trace Metals certified by the National Research Council of Canada) was used to check the accuracy and reproducibility of analyses [32,33]. Only those elements that exhibited good agreement between replicated measurements of SLRS-6 and the certified values (relative difference < 15%) are reported in this study.

#### *2.5. Landscape Parameters and Water Surface Area of the Ket and Taz River Basin*

The physio-geographical characteristics of the Ket and Taz tributaries and several sampling points of the main stem of both rivers were determined by applying available digital elevation model (DEM GMTED2010), soil, vegetation, and lithological maps. The landscape parameters were typified using TerraNorte Database of Land Cover of Russia (Ref. [34]; http://terranorte.iki.rssi.ru, accessed on 1 February 2020). This included various types of forest (evergreen, deciduous, needleleaf/broadleaf), grassland, tundra, wetlands, water bodies, and riparian zones. The climate parameters of the watershed were obtained from CRU grids data (1950-2016) [35] and NCSCD data (Ref [36]; doi:10.5879/ecds/00000001), respectively, whereas the biomass and soil OC content were obtained from BIOMASAR2 [37] and NCSCD databases. The lithology layer was taken from GIS version of Geological map of the Russian Federation (scale 1: 5,000,000, http://www.geolkarta.ru/, accessed on 1 February 2020).

#### *2.6. Data Analysis*

Element concentrations for all datasets were tested for normality using the Shapiro– Wilk test. In case of the data were not normally distributed, we used non-parametric statistics. Comparisons of major and trace element concentration in the main stem and the tributaries during two sampling seasons (Ket) and summer baseflow (Taz) were conducted using a non-parametric Mann Whitney test at a significance level of 0.05. For comparison of unpaired data, a non-parametric H-criterion Kruskal–Wallis test was used to reveal the differences between different seasons and between the main stem and tributaries. The Pearson rank order correlation coefficient (*p* < 0.05) was used to determine the relationship between major and trace element concentrations and main landscape parameters of several points of the main stem and all tributaries, as well as other potential drivers of TE concentrations in the river water such as pH, O2, water temperature, specific conductivity, DOC, SUVA (aromaticity), DIC, Fe, Al, particulate carbon and nitrogen, and total bacterial number.

#### **3. Results and Discussion**

#### *3.1. Spatial and Seasonal Variation of Elements in the Ket River and Control of River Hydrochemistry by Landscape Parameters*

All the primary data on dissolved (<0.45 μm) element concentration in both rivers together with relevant landscape parameters of the sampling points are available from the Mendeley Repository [28]. The spatial variations of element concentration in the main stem and tributaries of the Ket River were generally lower than the differences between the two seasons. This is illustrated by a plot of several major and trace elements along the full length of the river basin (Figure 3A–E), and further confirmed by a Mann–Whitney U test (Table S1 of the Supplementary Materials) which shows that the maximal differences in element concentrations are observed between seasons in both the main stem and the tributaries. The differences between the Ket River main stem and the tributaries were always lower or not pronounced for the same season. This allowed to calculate the mean concentrations of elements across the full length of the main stem and among all tributaries for each seasons (Figure S1 and Table S2 of the Supplementary Materials). Analyses of the differences in element concentration between the flood and the baseflow of the Ket River allowed distinguishing three group of solutes as illustrated in Figure 4. The first groups comprised elements having a factor of 2 to 10 higher concentrations in both main stem of the Ket River and tributaries during spring flood compared to summer baseflow, and included DOC, low soluble low mobile "lithogenic" hydrolysates (Be, Al, Ti, Cr, Ga, Y, Zr, Nb, REE, Hf, Bi, Th), some divalent heavy metals (Ni, Cu, Cd), oxyanions (Sb), Cs, Tl, and SO4. The second group typified highly mobile elements exhibiting lower concentrations during spring flood compared to the baseflow and included DIC, POC, alkaline and alkaline earth metals (Li, Na, Mg, Ca, Sr, Ba), labile nutrients (Si, P, Mn), oxyanions (As, Mo), Fe and U(VI). Finally, a group of elements demonstrated rather similar (within 30%) concentrations during both seasons and included Cl, micronutrients (V, Co, Zn, Se, Rb), Ge, W, and Cd.

A pairwise (Pearson) correlation of major and trace element concentration in the main stem and tributaries of the Ket River with main land cover parameters of the watershed was tested for both seasons, but notable correlations were observed only for the summer baseflow (Table S3). These correlations were detected only for alkaline-earth elements (Mg, Ca, Sr) and U, whose concentrations positively correlated (RPearson ≥ 0.5; *p* < 0.05) with phytomass stock on the watersheds (Figure 5). However, the observed correlations do not necessarily indicate a direct control but may be the consequence of the fact that carbonate-bearing loesses are most favorable substrates for productive forests compared to clay, sand and peat –rich soils [20]. The carbonate minerals of these substrates are known to act as sizable sources of alkaline-earth metals and uranium [the latter is carried in the form of highly labile uranyl-carbonate complexes] in shallow subsurface and groundwater feeding the river during summer baseflow [16,17]. Other elements did not demonstrate any sizable correlations to the landscape factors. The most likely reason for paucity of such a control is highly homogeneous coverage of the river basins by forest, bogs, and riparian zones with essentially the same climate, runoff, vegetation, soil, and total phytomass stock, represented by similar tree species of the boreal taiga of permafrost-free zone of the WSL.

Overall, the dominance of the season over space in solute control in the river basin let us to consider, for further analysis, only the same season (summer baseflow) for both rivers. This allowed better understanding the main driving factors (landscape parameters) of dissolved major and trace element variation along the full length of the main stem of both rivers as well as within each river tributaries.

**Figure 3.** Al (**A**,**B**), Fe (**C**,**D**), Mn (**E**,**F**), Sr (**G**,**H**) and Mo (**I**,**G**) concentration in the Ket River main stem (open circles) and tributaries (solid diamonds) during spring flood (**A**,**C**,**E**,**G**,**I**) and summer baseflow (**B**,**D**,**F**,**H**,**J**).

**Figure 4.** Mean ratio of element concentrations between spring flood period and summer baseflow in the Ket River main stem and tributaries.

**Figure 5.** Example of positive correlations between Mg (**A**), Ca (**B**), Sr (**C**), and U (**D**) and landscape factors of the Ket River in summer.

*3.2. Major and Trace Element Spatial Variation over the River Main Stem and among Tributaries of the Taz River Basin and Land Cover Control*

The variations of major and trace element concentration along the Taz River were also rather low as illustrated in a plot of some major and trace element concentrations over the river distance, from the headwaters to the mouth (Figure 6). The concentration of DOC, Ptot, NO3, NO2, Al, Fe, Cd, Ba, and Bi increased southward, with an increase in mean annual temperature. In contrast, in the northward direction, with an increase in tundra and discontinuous permafrost coverage, the concentrations of Cl, SO4, Li, B, Na, K, V, Ni, Cu, Y, REE, and U increased, which may stem from a combination of sea-salts release from the underlying marine clays and silts and far-range (>300 km) atmospheric transfer from the Norilsk smelters (V, Ni, Cu). The variations among tributaries were generally larger but did not exhibit any systematic evolution between the upper and lower reaches of the Taz River basin (Table S4). Therefore, we can hypothesize that the primary factor controlling solute concentration in the tributaries is land cover (see Section 3.3 below). In the main stem and tributaries, the difference in dissolved element concentration between the upper (southern) and lower (northern) part of the river has not exceeded 20–30% which was often within the standard deviation of the mean values. The only exception is Mn concentrations in the tributaries were two times higher in the south compared to the north.

**Figure 6.** SO4 <sup>2</sup><sup>−</sup> (**A**), Al (**B**), DIC (**C**), Fe (**D**), Ca (**E**) Mn (**F**), Mg (**G**), and Sr (**H**) concentration in the Taz River main stem and tributaries during summer baseflow.

Pairwise correlations of major and trace element concentration in the Taz River basin with main physico-geographical, geocryological, and climatic features of the watershed revealed several potential drivers (Table S5). The tundra coverage of the watershed, corresponding to northward directions and proximity to the sea exhibited strong positive

(RPearson > 0.60, *p* < 0.01) correlations with Li, B, Na, K, Cl, SO4, and U. These elements likely originate from sea salts of the former marine clays dominating the bedrocks of river catchments in the northern part of WSL [17] but also the deposition of marine aerosols in the form of snow [38]. The latter is also known to be enriched in Ni and Cu, reflecting the proximity of the low Taz reaches to the Norilsk Cu-Ni smelters in the northern part of WSL. This can explain positive (RPearson > 0.50, *p* < 0.05) correlations between the tundra coverage and the concentrations of Ni and Cu in the rivers of the Taz basin. Presumably, the atmospheric deposition of metal-rich aerosols can be transferred, via plant litter decay and surface runoff, to the river main stem and tributaries in the northern part of the basin.

The vegetation also exhibited notable impact on elements concentrations in the river water. Thus, presence of larch trees (deciduous needle-leaf forest coverage) positively correlated with concentrations of DOC and macronutrients (Si, Mg, Ca, Sr) and Cs. Light coniferous and broadleaf mixed forest positively impacted concentrations of DOC, macro- (PO4, NO3, NH4) and micro-nutrients (Mn, Fe, Co, Mo, Ba, Cd). Finally, the organic carbon content in upper 0–30 and 0–100 cm of soil positively correlated with some micro-nutrients (Mn, Co, Mo) but also Be, Sc, Cd, Sb, and Bi. Noteworthy is that neither the watershed surface area nor the permafrost coverage, which are the two potential drivers of element geochemistry in surface waters i.e., [39–42], correlated with major and trace elements of the Taz River basin. Furthermore, we were not able to detect any control of other landscape features such as riparian zone, wetlands, and recent burns.

#### *3.3. Common Features of Spatial Distribution of Major and Trace Element in Two River Basins during Summer*

Analysis of pairwise correlations described in Sections 3.1 and 3.2 for two river basins individually revealed a number of common features in terms of apparent control of landscape parameters on element concentration in river waters. This allowed correlating the hydrochemical composition of the water column with major environmental factors that are likely to operate on both river basins as illustrated for some elements in Figure 7. The first factor is forest biomass in general and, in particular, the coverage of the watershed by larch trees (deciduous light needle leaf forest), which are most efficient for recycling the nutrients such as Li, B, Si, Zn, Rb, and Ba [43]. The second environmental parameter is bog coverage of the watershed, which was positively correlated with SUVA and low mobile lithogenic elements and divalent transition metals whose transport is facilitated by high DOC concentration (V, Ni, Cu, Y, heavy REE). The impact of bogs on retention of these elements was demonstrated in North European boreal rivers [44,45]. Noteworthy that the bog presence facilitated only the transport of heavy REE and not the light ones. In the WSL surface waters, the former are known to form much stronger complexes with allochthonous DOM from bogs and peat waters and less prone to be carried as organo-ferric colloids [30,46–48].

A small number of elements were associated with total bacterial number (SUVA, P, V, Fe, and As) and presumably marked shallow subsurface discharge of Fe(II)-rich waters from the adjacent peatlands and biotically driven formation of large size Fe hydroxide colloids, stabilized by allochthonous (aromatic) organic carbon. Strong co-precipitation of P, V, and As with Fe hydroxides is known from laboratory experiments with organicrich (peatland) waters in the presence of soil and aquatic bacteria [49–52]. The area of riparian zones of the river positively correlated with DIC, Sr, Ba, U, and particulate organic nitrogen. While the first four elements could mark an enhanced discharge of deep and subsurface groundwaters in the riparian zones, notably of larger rivers, the correlation with PON could illustrate the generation of N-rich particles in the sediments of highly productive floodplains [13,18]. Among common climate parameters of both river basins, the precipitation did not impact the hydrochemical composition. However, the mean annual air temperature encompassed the contrast between northern and southern rivers and correlated with labile components of the river water (S.C., DIC, Li, B, Si, Ca, Sc, Zn, Se,

Rb, Sb, Cs, Ba, and U). This was fully consistent with previous observations of large (Ob, Ref. [20]; Lena, Ref. [53]) and small [15–17] Siberian rivers.

**Figure 7.** Pairwise correlations of several major and trace elements with two most pronounced landscape parameters—Phytomass stock at the watershed and bog coverage, total bacterial count (TBC) and mean air temperature, for both Taz and Ket river main stem and tributaries.

Given that the season was an important driving factor of element concentrations in the Ket River (see Section 3.1), we attempted a PCA treatment of elementary dataset of both rivers during summer baseflow. This analysis revealed three main factors capable of explaining 26, 14, and 6.6% of total variability, respectively (Figure S2, Table S6). The first factor acted positively on aluminum and other low mobile hydrolysates such as Be, Al, Ti, Cr, Ga, Y, Zr, Nb, Y, REE, Hf, Th and elements which are known to be strongly complexed with DOM such as Cu and Se [54]. Presumably, it reflected mobilization of low soluble elements from lithogenic (silicate) minerals in the form of organo-mineral (essentially organo-aluminum) colloids as it is known from studies in soil porewaters of the WSL regions [24]. The second factor negatively acted on labile Na, Mg, Si, Ca, P, and Fe (probably as Fe(II)). The F2 was also positively linked to DOC, micronutrients (Zn, Rb, Ba), and phytomass at the watershed (primarily, light needle-leaf trees) as well as mean temperature and the proportion of younger (<25 Ma) rocks. In the WSL, these rocks contain carbonate concretions and partially carbonated loesses. The second factor thus likely represented a combined effect of productive litter in larch forest growing on carbonate-rich rocks and groundwater feeding of the rivers by soluble, highly mobile elements. Finally, the third factor, although exhibited limited explanation capacity, was strongly linked to pH, specific conductivity, DIC, Sr, and U and clearly marked a direct impact of bicarbonate-rich, slightly mineralized waters that reside in shallow (Taz) and deeper (Ket) subsurface reservoirs containing carbonate minerals and actively discharge into the river during baseflow.

Noteworthy is the similarity of the two main groups of major and trace elements, DOC- and DIC-related, based on pairwise correlations between element concentrations. The DIC concentration in both river basins significantly (*p* < 0.05) correlated with S.C., Si, alkalis (Li, Na), alkaline-earth metals (Mg, Ca, Sr, Ba), and U. This correlation marks an enhanced mobility of these elements in the form of ionic forms and neutral molecules (uranyl as carbonate/hydroxide complexes), essentially controlled by discharge of DIC and mobile element—rich groundwaters either at the river bank or in the hyporheic zone. These mechanisms of labile element mobilization are fairly well-known across the boreal zone [55]. The DOC correlated with trivalent (Sc, Al, Ga, Y, REE), tetravalent (Ti, Zr, Hf, Th), and other (Be, Cr, Nb) hydrolysates, Se and Cs. These elements are essentially present in the river water in the form of Al-organic colloids or organic complexes [14,54,56–58]. These two main groups of solutes were evidenced not only in large and small rivers and streams [14,16,17,20,59] but in lentic waters of the region, such as lakes and ponds [60,61], supra-permafrost waters [23], peat porewaters [22–24], and peat ice [21].

#### *3.4. Climate Change Consequences on Element Concentration in WSL Rivers (Vegetation and Permafrost/Lithology Control) Using a Substituting Space for Time Approach*

The main limitation of our approach is a restricted seasonal coverage in river sampling and lack of daily/weakly discharge measurements, which did not allow calculating elementary export fluxes (yields). As such, predictions of climate change impact on riverine export of this particular pristine basins can be performed solely in terms of river water hydrochemistry. Note however that an approach of element annual yields is well developed for other European Arctic [55] and small western Siberian rivers [13,59] including anthropogenically affected Ob River in its middle course [19]. Therefore, taking the advantage of highly pristine character of Taz and Ket River basins (in contrast to the neighboring Pur River basin and left tributaries of the Ob River, which are strongly impacted by gas and oil industry), one can use these river systems to approximate future changes in river hydrochemistry linked to on-going climate change and permafrost thaw. The validity of this approach is based on rather weak spatial variations within each basin discovered in the present study (notably over the course of the main stem but also among tributaries) in contrast to sizable variations in solute concentrations between the northern and southern river basin.

For this, we can employ a "substituting space for time" approach which postulates, in a broader context, that spatial phenomena which are observed today can be used to describe past and future events [62]. Such an approach has been successfully used in western Siberia since the pioneering work of Frey and Smith [12] to model possible future changes in small rivers [13,59], lakes [61], soil waters [23], and permafrost ice [21]. Indeed, with permafrost and forest boundary shift northward over next decades [9,63–65], the northern part of the Taz River (tundra and forest-tundra of continuous to discontinuous permafrost) can be entirely transformed into southern part-like territory of taiga and forest-tundra biome with discontinuous to sporadic permafrost, whereas the entire Ket River basin can be used as a surrogate for hydrochemistry for the southern part of the Taz River. The mean element concentrations in the main stem and tributaries in the two sub-basins of the Taz and Ket River are illustrated in Figure 8. It can be seen that only a few elements exhibited sizable contrast in river water concentration between the southern and northern part of the Taz River basin and the Ket river catchment. Therefore, applying a substituting space for time approach, we can anticipate that the concentration of DIC and mobile elements (Li, B, Mg, Ca, Sr, Rb, Sb, W, U) and, to a lesser degree, DOC, Si, and micronutrients (Mn, Fe, Zn, Mo) will increase in the low reaches of the Taz River.

**Figure 8.** Elements concentration in Upper Taz (0–400 km from the mouth, northern part), Low Taz (400–800 km from the mouth, southern part), and Ket River (main stem and tributaries). Box represents median, whiskers correspond to lower and upper quartile.

Overall, the main consequences of the climate change in western Siberia on middle size river basin hydrochemistry seems to be primarily linked to the changes in dominant vegetation and biomass at the river watershed. In the northern part of the WSL, the change in hydrological pathways due to permafrost thaw might enrich the river waters in soluble elements due to enhanced underground water feeding (DIC, alkaline-earth elements (Ca, Sr), oxyanions (Mo, Sb), and U). The thickening of the active layer and enhanced involvement of thawed peat layer down to underlying mineral horizons may increase the export of trivalent and tetravalent hydrolysates in the form of organo-ferric colloids. However, at the short-term scale, due to two counterbalanced sources and sinks (permafrost thaw and plant uptake), the overall impact of the climate change on inorganic solute export by rivers from the land to the ocean may be smaller than that traditionally viewed for organic carbon.

#### **4. Conclusions**

Toward a better understanding of land cover control on dissolved (<0.45 μm) major and trace elements in the river water of high latitude regions, we selected two mediumsize pristine rivers of the northern, permafrost-bearing and southern, permafrost-free region of the world's largest peatland, the Western Siberia Lowland (Taz and Ket River, respectively). We sampled the main stem and tributaries while encompassing a large gradient in river basin size, permafrost, and vegetation coverage at essentially similar lithological substrate, runoff, and relief. In the Ket River, the difference in major and trace solute concentration between two seasons was larger than the difference between the main stem and the tributaries. Given that the primary factor controlling the river solute concentration was the season, this allowed straightforward comparison of two river basins during the summer baseflow. Furthermore, sizable variations in solutes among different tributaries of each watershed (depending on land cover) allowed testing the impact of first-order environmental factors on river water hydrochemistry.

The similarity of landscape parameters (mainly, taiga vegetation, partially bogs and permafrost coverage, and, in a lesser degree, the lithology) controlled the major and trace element concentration pattern in two medium-size rivers studied in this work. A number of landscape parameters of the main stem and tributaries correlated with nutrients (B, Si, Zn, Rb, Ba with phytomass), low mobile lithogenic elements, and divalent transition metals whose transport is facilitated by high DOC (V, Ni, Cu, Y, heavy REE with bogs), or specific components reflecting bacterially controlled processing of DOM and groundwater discharge in the water column or within the hyporheic zone (Fe, P, SUVA, V, As). The mean annual air temperature encompassed the contrast between northern and southern rivers and correlated with labile elements of the river water (DIC, Li, B, Si, Ca, Sc, Zn, Se, Rb, Sb, Cs, Ba, and U), which was fully consistent with previous observations of large and small Siberian rivers.

Via quantitative comparison of element concentration along the south–north gradient of studied river basins and applying a substituting space for time approach, we infer that permafrost and vegetation shift northward may sizably (a factor of 2 to 3) enrich the river water in the north in highly mobile elements (DIC, Ca, Sr, U), and, in a lesser degree, in DOC, POC, and nutrients (Si, Mn, Zn, Fe, Mo), whereas the current concentrations of other macro (P) and micro-nutrients (V), and insoluble hydrolysates – geochemical tracers (Y, REE, Zr) in the low reaches of Taz might decrease.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/w14142250/s1, Figure S1: Seasonal mean ± SD concentration elements exhibiting higher concentrations during spring flood compared to summer baseflow in tributaries and main channel Ket' in spring flood (blue) and summer (early fall) baseflow (orange). Figure S2: Results of PCA of ~65 hydrochemical abd 16 land cover variables in ~ 60 sampling points of the main stem and tributaries of the Ket and Taz River basin, collected during summer baseflow. Table S1A. Mann-Whitney U Test comparison concentration in different season (flood vs. baseflow) in tributaries and main channel of the Ket River. Red color is for statistically significant differences. Table S1B. Mann-Whitney U

Test comparison concentration tributaries and main channel in different seasons. Table S2. Mean ± SD concentrations of elements and other measured parameters in the Ket River main stem and tributaries. Table S3. Pairwise Pearson correlations of major and trace element concentrations in the Ket River (main stem and tributaries) during summer baseflow and main hydrochemical parameters of the water column and land cover. Significant (*p* < 0.05) correlations are given in red and most significant (*p* < 0.01) are highlighted in pink. Table S4. Hydrochemical parameters of the Taz River basin (stem and tributaries), July 2019. Table S5. Pearson pairwise correlations of major and trace elements of the river water (Taz and tributaries) with landscape parameters, vegetation coverage and climate. Significant (*p* < 0.05) correlations are given in red and most significant (*p* < 0.01) are highlighted in pink. Table S6. Results of PCA of ~65 hydrochemical abd 16 land cover variables in ~ 60 sampling points of the main stem and tributaries of the Ket and Taz River basin, collected during summer baseflow. Significant (*p* < 0.05) correlations are given in red and most significant (*p* < 0.01) are highlighted in pink.

**Author Contributions:** O.S.P. designed the study and wrote the paper; A.G.L. and I.V.K. performed sampling, analysis, and their interpretation; S.N.V. and O.S.P. were responsible for the choice of sampling objects and statistical treatment. L.S.S. was in charge of DOC and bacterial analyses and their interpretation; M.A.K. provided the GIS data for river watersheds from available databases. All authors have read and agreed to the published version of the manuscript.

**Funding:** RSF grant 22-17-00253, RFBR grant 20-05-00729, the TSU Development Program "Priority-2030", and grant "Kolmogorov" of MES (Agreement No 075-15-2022-241).

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** We acknowledge support from RSF grant 22-17-00253, RFBR grant 20-05-00729, the TSU Development Program "Priority-2030", and grant "Kolmogorov" of MES (Agreement No 075-15-2022-241).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **How to Implement User-Friendly BLMs in the Absence of DOC Monitoring Data: A Case Study on Bulgarian Surface Waters**

**Tony Venelinov <sup>1</sup> and Stefan Tsakovski 2,\***


**Abstract:** The metal bioavailability concept is implemented in the Water Framework Directive (WFD) compliance assessment. The bioavailability assessment is usually performed by the application of user-friendly Biotic Ligand Models (BLMs), which require dissolved metal concentrations to be used with the "matching" data of the supporting physicochemical parameters of dissolved organic carbon (DOC), pH and Cadissolved. Many national surface water monitoring networks do not have sufficient matching data records, especially for DOC. In this study, different approaches for dealing with the missing DOC data are presented: substitution using historical data; the appropriate percentile of DOC concentrations; and combinations of the two. The applicability of the three following proposed substitution approaches is verified by comparison with the available matching data: (i) calculations from available TOC data; (ii) the 25th percentile of the joint Bulgarian monitoring network DOC data (measured and calculated by TOC); and (iii) the 25th percentile of the calculated DOC from the matching TOC data for the investigated surface water body (SWB). The application of user-friendly BLMs (BIO-MET, M-BAT and PNEC Pro) to 13 surface water bodies (3 reservoirs and 10 rivers) in the Bulgarian surface waters monitoring network outlines that the suitability of the substitution approaches decreases in order: DOC calculated by TOC > the use of the 25th percentile of the data for respective SWB > the use of the 25th percentile of the Bulgarian monitoring network data. Additionally, BIO-MET is the most appropriate tool for the bioavailability assessment of Cu, Zn and Pb in Bulgarian surface water bodies.

**Keywords:** BLM; TOC; DOC; Fe; historical data

#### **1. Introduction**

The determination of metal species and their bioavailability depends crucially on Dissolved Organic Matter (DOM). The presence of DOM is specified as a Dissolved Organic Carbon (DOC) concentration, contrary to Natural Organic Matter (NOM), which is a collective organic component of water. NOM is a complex mixture of leaf litter, organic acids, proteins and many more complex organic molecules. The Total Organic Matter (TOM) constitutes Particulate Organic Matter (POM) and DOM. The organic carbon results of unfiltered samples are reported as Total Organic Carbon (TOC) and samples filtered through a 0.45 μm filter are reported as DOC. Both results equal the mass of carbon present in the mixture of organic compounds in the raw water/filtrate [1]. Sometimes TOC is reported instead of DOC, but it should be stressed that the DOC concentration will always be less than the concentration of TOC. DOC generally accounts for about 50% of DOM in typical surface water bodies and is in the range between 1 and 15 mg/L. DOC mainly comprises humic substances of natural origin, which is formed as a result of the breakdown of plant and animal tissues by chemical and biological processes or from anthropogenic sources [2].

**Citation:** Venelinov, T.; Tsakovski, S. How to Implement User-Friendly BLMs in the Absence of DOC Monitoring Data: A Case Study on Bulgarian Surface Waters. *Water* **2022**, *14*, 246. https://doi.org/10.3390/ w14020246

Academic Editor: Liudmila S. Shirokova

Received: 28 November 2021 Accepted: 13 January 2022 Published: 15 January 2022

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

**Copyright:** © 2022 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/).

The bioavailability of Cu, Ni and Zn depends on DOC concentrations [3–5], with higher concentrations of DOC resulting in reduced toxicity. Evidence suggests that the pH also has a significant impact on metal bioavailability and toxicity in aquatic environments [6,7]. Major cations (Ca2+, Na+, K+ and Mg2+) can also protect aquatic organisms against free metals by competing on the biotic binding sites, e.g., fish gills [8]. The presence of dissolved Fe and Al can affect the binding of the abovementioned metals to DOC due to their very high affinity for complexation by DOC. In such conditions, the availability of binding sites for the less strongly bound metals is reduced, resulting in an increase in their bioavailability for aquatic organisms. Normally, bodies of water that are high in Fe and Al are associated with a lower water pH (6.5–7.5), under which conditions both metals form insoluble precipitates.

DOC is a composite parameter for a heterogeneous mixture of polyfunctional polymers and contains functional groups (ligands) that bind free metal ions (the most toxic inorganic metal fraction) to reduce the interaction between free metal ions and aquatic organisms [9]. Discharge derived DOC (such as sewage effluents) is likely to be of a different composition to natural DOC and comprises proteins, amino acids, polysaccharides and synthetic chelating agents, such as ethylenediamine tetra acetic acid (EDTA). The current evidence base suggests that this latter type of DOC has a much greater capacity for binding metals compared with similar concentrations of naturally sourced DOC [10]. In bodies of water with low pH, DOM–metal complexes are weaker, resulting in a release of bioavailable ions and, therefore, an increase in the toxicity of the metal to aquatic organisms [11]. Taking DOM binding into account, the effect of considering bioavailability in environmental quality standards (EQSbioavailable) [12] has been assessed by the use of user-friendly biotic ligand models (BLMs), which has been proposed as a tool for regulating agencies. In the Netherlands, the assessment of bioavailability is incorporated into policy documents [13]. Feasibility studies on the implementation of a bioavailability-based approach for metals have been undertaken by several member states [14–19]. Guidance on how bioavailability may be incorporated into compliance assessments, classifications and local risk assessments is included in the recent EU Technical Guidance for deriving Environmental Quality Standards (EQS) under the European Water Framework Directive (WFD) [20].

High-quality surface water monitoring data provide tools for the assessment of potential risks to aquatic species. Obtaining such reliable data should follow the principles of surface water chemical monitoring under the WFD [21], including sampling, preservation and analysis [22]. The analytical methods for every analyte of interest should comply with the minimum performance requirements as stated in the QA/QC Directive [23]. A common implementation strategy of the WFD to achieve a good status and prevent any deterioration in the status (ecological and chemical) of European water bodies led to the development of methods for deriving EQSs to be used for the status assessment of surface water bodies (SWBs). Following these requirements, as a general prerequisite for WFD monitoring, EQS compliance assessment and subsequent decision-making are reliably achieved. An EQS is a limit for environmental disturbances, in particular from an ambient concentration of pollutants and wastes, which determines the maximum allowable degradation of environmental media. A water body cannot achieve a good chemical or ecological status if the EQS for any WFD Priority or Specific Substance is exceeded.

The compliance assessment follows a tiered approach, which: (i) compares the annual average concentration, calculated from monitoring data with the EQSbioavailable; (ii) uses user-friendly tools (based on Biotic Ligand Models) for the calculation of local metal bioavailability for a comparison between the measured dissolved metal concentration at a site and the EQS; and (iii) considers the local background concentrations of metals as a part of the EQS in the risk assessment [8]. The compliance assessment, which takes into account the bioavailability and uses simplified and user-friendly bioavailability tools (such as BLMs), requires that the concentration of dissolved metals (Cu, Zn, Ni, Pb and Mn) preferably be used along with data for the supporting physicochemical parameters of DOC, pH and dissolved Ca [24]. Due to temporal and spatial variations, data for all of the required input parameters (dissolved metal concentrations, dissolved Ca, DOC and pH) should ideally be "matched" to increase the reliability of the results. Without these, the simplified tools will either not run or not run reliably. The term "matched" means that the required input parameters are all determined in the same sample, obtained from the site of interest, to enable the identification of potential risks. A few of the member states (including Bulgaria) are just starting the implementation of bioavailability-based EQS and some of the required data are not fully available, especially for DOC, which is not routinely measured in Europe. In such cases, two approaches are possible when dealing with the missing DOC data: the use of historical monitoring data (e.g., data for similar and/or neighbouring catchment types) and the use of substitute data. Such approaches should be carefully chosen due to the quality of the available data—the dataset is likely to have been collected with inadequate sampling protocols, insufficient analytical sensitivity and LOQ, for a different purpose than the implementation of an EQS compliance assessment, with different sampling frequency, etc. [25,26]. Different approaches to estimating DOC have been proposed, including: UV absorbance [27,28]; colour measurements [29]; dissolved iron [30]; fluorescence measurements [31]; and a modifying factor that is dependent on the optical properties of DOC [32]. Any of the proposed approaches to estimating DOC should be based on locally derived empirical relationships and proven to be valid for the studied territory. Such data can be used for feasibility and screening assessments, but these need to be assessed on a case-by-case basis.

The present study aims to propose a methodology to the national environmental bodies for the implementation of BLM in compliance assessments when DOC data are missing. The methodology comprises a selection of appropriate substitute approaches, which are applicable for the environmental authorities and their validation, with the matching data using the three most widely used BLMs: BIO-MET; M-BAT; and PNEC Pro. In this manner, the metal bioavailability could be assessed using the most suitable pair of substitute approaches and a user-friendly BLM.

#### **2. Materials and Methods**

Data for the concentrations of DOC, Fedissolved and TOC were obtained by the Bulgarian Environmental Agency surface water monitoring programme. The database returned over 27,000 measurement results for around 950 SWBs for the period 1 January 2010– 31 December 2020 (Figure 1).

The methods used by the different laboratories followed standardized procedures. For the determination of pH, EN ISO 10523:2012 was used; for Ca, EN ISO 14911:2002 and ISO 6058:2002 were used; for DOC and TOC, EN 1484:2001 was used, including one result for TOC that was obtained by EN ISO 6468:2006. Cu, Zn, Ni and Pb were measured using EN ISO 17294-2:2016. The basic statistics of the dataset are presented in Table 1. No outlier tests were performed and where the halved limit of quantitation (LOQ) data was presented (for either DOC, Fedissolved or TOC), the whole dataset was omitted. Where the calculated DOC/TOC ratio exceeded one, the concentrations were not used in the subsequent calculations.

**Table 1.** The basic statistics of the monitoring data for the Bulgarian surface waters (2010–2020).


**Figure 1.** Data availability map.

User-friendly tools (based on the BLMs) were used in the metal bioavailability calculations. For the application of BIO-MET (https://bio-met.net/, last accessed 25 October 2021), M-BAT (https://www.wfduk.org, last accessed 25 October 2021) and PNEC Pro (http://www.pnec-pro.com/, last accessed 25 October 2021), data for the concentrations of pH, DOC, TOC and Cadissolved were obtained by the open data portal of the Bulgarian Environmental Agency surface water monitoring programme (https://data.egov.bg/data/ view/9ee2ca78-051d-4ec8-b8f8-1a1e93ee637e, last accessed 25 October 2021). A total of 13 water bodies—3 reservoirs (Dospat Dam, Pchelina Dam and Iskar Dam) and 10 rivers (Strumeshnitsa, Mesta, Yantra, Dospat, Dragovishtitsa, Negovanka, Eleshnitsa, Struma, Rusenski Lom and Iskar) were used for the comparison during the implementation of the BLMs, and for which the matching datasets for DOC, TOC, pH, Ca, Cu, Zn, Ni and Pb were available for 2020.

#### **3. Results**

#### *3.1. Substitutes for Missing DOC*

The Bulgarian Environmental Agency surface water monitoring programme revealed a total of 1829 measurement results for DOC for all monitored surface water bodies since 2010 and 13,234 measurement results for dissolved Fe. The "matching" of the data indicated 805 results for DOC and Fedissolved. Unfortunately, nearly half of the results could not be used for correlation between the parameters since most of the Fedissolved concentrations were reported as LOQ/2. The analysis showed a poor correlation (R2 = 0.38, n = 445) when "matched" data was used (Figure 2).

**Figure 2.** The relationship between the measured concentrations of DOC and dissolved Fe in Bulgarian surface waters.

The Bulgarian Environmental Agency surface water monitoring programme contained 11,945 results for TOC concentrations in the 2010–2020 period. The "matching" data indicated 736 results for DOC and TOC. The analysis showed a significant correlation (R<sup>2</sup> = 0.95, n = 736) when matched data was used (Figure 3) and a conversion factor DOC = 0.81 × TOC could be applied for missing DOC data.

**Figure 3.** The relationship between the measured concentrations of DOC and TOC in Bulgarian surface waters.

#### *3.2. Substitutes in the Absence of Matching Data*

The UK [9] and France [18] have adopted an approach to use the 25th percentile of the DOC concentrations for estimating DOC when the matching data are unavailable. In this paper, the calculated DOC concentrations using the TOC data (DOC = 0.81 × TOC) was used as a substitute. Furthermore, the 25th percentile of the DOC concentrations (2.4 mg/L) from the Bulgarian DOC data obtained in 2020 (n = 1711) was used and compared with the results of the matching data. Additionally, the 25th percentile of the calculated DOC concentrations using the measured TOC concentration for the respective water body was used and compared with the results of the matching data.

#### *3.3. Application of User-Friendly Models (BIO-MET, M-BAT and PNEC Pro) in Compliance Assessment*

Following the tiered approach, the user-friendly BLMs were implemented on 13 water bodies (3 reservoirs and 10 rivers) with the available matching datasets for DOC, TOC, pH, Ca, Cu, Zn, Ni and Pb for 2020. The average annual concentrations that were calculated from the monitoring data were compared with the EQSbioavailable (1 μg/L for Cu, 4 μg/L for Ni, 8 μg/L for Zn and 1.2 μg/L for Pb), which was derived from the WFD [21] and national ordinance [33].

The results show 13 exceedances for Cu, 6 exceedances for Zn and 4 exceedances for Pb (Table 2). No exceedances of the average annual concentrations were observed for Ni, compared with the EQS.


**Table 2.** The exceedances for Cu, Zn, Ni and Pb in the studied Bulgarian surface water bodies.

When the average annual concentrations exceeded the EQS, the samples moved to the second tier where the Biotic Ligand Model-based tools (BIO-MET, M-BAT and PNEC Pro) were used for the calculation of local metal bioavailability and a comparison between the measured dissolved metal concentration at a site and the EQS was made. The DOC concentration played an important role in determining the bioavailable dissolved metal concentration.

For the calculations, the following data were used: the measured DOC concentrations; the calculated DOC from the matching TOC data (DOC = 0.81 × TOC); the 25th percentile (2.4 mg/L) of the DOC and calculated TOC for 2020 of all surface water bodies (n = 1711); and the 25th percentile of the calculated DOC from the matching TOC data for 2020 for each of the water bodies studied. An indication of whether there were exceedances was expressed as the Risk Characterization Ratio (RCR). If the RCR was lower than one, the site was considered as having "passed" the assessment and no further investigation was necessary. If the value of the RCR was greater than one, this indicated an exceedance of the EQS and led to a progression to Tier 3.

#### 3.3.1. BLM Comparisons for Copper

The user-friendly BLM models were implemented for the 13 water bodies where the annual concentrations of Cu for 2020 exceeded the generic EQSbioavailable (1 μg/L). First, the BLM models were performed with the matching data (n = 63) for dissolved Cu, pH, dissolved Ca and DOC, and after that, the DOC values were substituted with: the calculated DOC from the available TOC data (DOC = 0.81 × TOC); the 25th percentile of the DOC data from the Bulgarian monitoring network results from 2020 (P25\_BG); and by the 25th percentile of DOC data that was calculated from the TOC for the respective water body (P25\_SWB), as described earlier. The calculated RCRs for all of the approaches are presented in the Supplementary Material (Figures S1–S3).

The comparisons between the Cu bioavailability results obtained by the BLM models are presented in Figures 4–6. The correlations between the BIO-MET calculated RCR using matching and substituted data were very good (Figure 4). For the RCR values using percentiles approaches: P25\_BG and P25\_SWB were overestimated with slopes of 2.7 and 1.8, respectively. These environmental conservative estimates were more pronounced when the 25th percentile of the data from the Bulgarian monitoring network was used. The results from the M-BAT model followed a similar pattern to those of the BIO-MET (Figure 5), but the R2 values dropped to 0.54 for P25\_BG and 0.69 for P25\_SWB. The comparison between the RCR values obtained by the PNEC Pro model revealed lower correlations and slopes that were significantly different from one for all pairs (Figure 6). It should be mentioned that the comparison was restricted by the validated applicability domain of the PNEC Pro model, which decreased the number of calculated RCR values.

**Figure 4.** The comparison of the BIO-MET calculated RCR for Cu between the matching data and substituted data (P25\_BG—the 25th percentile of the available DOC data and/or the DOC data calculated by the TOC from the Bulgarian monitoring network results during 2020; P25\_SWB—the 25th percentile of the DOC data calculated by the TOC for the respective water body; TOC—the matching DOC data calculated by TOC).

**Figure 5.** The comparison of the M-BAT calculated RCR for Cu between the matching data and substituted data (P25\_BG—the 25th percentile of the available DOC data and/or the DOC data calculated by the TOC from the Bulgarian monitoring network results during 2020; P25\_SWB—the 25th percentile of the DOC data calculated by the TOC for the respective water body; TOC—the matching DOC data calculated by TOC).

**Figure 6.** The comparison of the PNEC Pro calculated RCR for Cu between the matching data and substituted data (P25\_BG—the 25th percentile of the available DOC data and/or the DOC data calculated by the TOC from the Bulgarian monitoring network results during 2020; P25\_SWB—the 25th percentile of the DOC data calculated by the TOC for the respective water body; TOC—the matching DOC data calculated by TOC).

#### 3.3.2. Comparisons for Zn

The user-friendly BLM models were implemented for the six water bodies (two reservoirs and four rivers) where the annual concentrations of Zn for 2020 exceeded the generic EQSbioavailable (8 μg/L) and the matching data (n = 34) for dissolved Zn, pH, dissolved Ca and DOC were available. The calculated RCRs are presented in the Supplementary Material (Figures S4–S6).

The comparison between the Zn bioavailability results obtained by the three BLM models using the matching and substituted DOC data is presented in Figures 7–9. The correlations between the calculated RCRs using the matching and substituted data for all BLMs were very good (R<sup>2</sup> > 0.90). The slopes varied between 1 and 1.5, decreasing in the order: P25\_BG > P25\_SWB > TOC.

**Figure 7.** The comparison of the BIO-MET calculated RCR for Zn between the matching data and substituted data (P25\_BG—the 25th percentile of the available DOC data and/or the DOC data calculated by the TOC from the Bulgarian monitoring network results during 2020; P25\_SWB—the 25th percentile of the DOC data calculated by the TOC for the respective water body; TOC—the matching DOC data calculated by TOC).

**Figure 8.** The comparison of the M-BAT calculated RCR for Zn between the matching data and substituted data (P25\_BG—the 25th percentile of the available DOC data and/or the DOC data calculated by the TOC from the Bulgarian monitoring network results during 2020; P25\_SWB—the 25th percentile of the DOC data calculated by the TOC for the respective water body; TOC—the matching DOC data calculated by TOC).

**Figure 9.** The comparison of the PNEC Pro calculated RCR for Zn between the matching data and substituted data (P25\_BG—the 25th percentile of the available DOC data and/or the DOC data calculated by the TOC from the Bulgarian monitoring network results during 2020; P25\_SWB—the 25th percentile of the DOC data calculated by the TOC for the respective water body; TOC—the matching DOC data calculated by TOC).

#### 3.3.3. Comparisons for Lead

The user-friendly BLM models were implemented for the four water bodies (two reservoirs and two rivers) where the annual concentrations of Pb for 2020 exceeded the EQSbioavailable (1.2 μg/L) and the matching data (n = 22) for dissolved Pb, pH, dissolved Ca and DOC were available. The calculated RCRs are presented in the Supplementary Material (Figures S7–S9).

The comparison between the Pb bioavailability results obtained by the three BLM models using the matching and substituted DOC data are presented in Figures 10–12. The performed linear regression analysis only showed a good agreement between the RCRs

obtained using the matching data and the DOC calculated by TOC. The best result achieved by the percentile approaches was for the BIO-MET calculated RCR using P25\_SWB as the substituted data.

**Figure 10.** The comparison of the BIO-MET calculated RCR for Pb between the matching data and substituted data (P25\_BG—the 25th percentile of the available DOC data and/or the DOC data calculated by the TOC from the Bulgarian monitoring network results during 2020; P25\_SWB—the 25th percentile of the DOC data calculated by the TOC for the respective water body; TOC—the matching DOC data calculated by TOC).

**Figure 11.** The comparison of the M-BAT calculated RCR for Pb between the matching data and substituted data (P25\_BG—the 25th percentile of the available DOC data and/or the DOC data calculated by the TOC from the Bulgarian monitoring network results during 2020; P25\_SWB—the 25th percentile of the DOC data calculated by the TOC for the respective water body; TOC—the matching DOC data calculated by TOC).

**Figure 12.** The comparison of the PNEC Pro calculated RCR for Pb between the matching data and substituted data (P25\_BG—the 25th percentile of the available DOC data and/or the DOC data calculated by the TOC from the Bulgarian monitoring network results during 2020; P25\_SWB—the 25th percentile of the DOC data calculated by the TOC for the respective water body; TOC—the matching DOC data calculated by TOC).

#### **4. Discussion**

The first step in the implementation of BLMs in the case of missing DOC monitoring data is the possible substitution of the available historical data from the Bulgarian Environmental Agency surface water monitoring programme. Contrary to the practice used in Great Britain [29], the correlation obtained between the DOC and Fedissolved concentrations (Figure 2) shows that it cannot be used for the substitution of missing DOC data. The result obtained from the comparison between the available DOC and TOC data reveals that a conversion factor of DOC = 0.81 × TOC can be applied when DOC data are missing (Figure 3). This outcome is in agreement with the approach used in Finnish [34] and Swedish [11] networks, where TOC is traditionally measured and the applied conversion factors are 0.94 and 0.90, respectively. Similar to the Bulgarian monitoring data, a relationship (DOC = 0.83×TOC) is used for BLM purposes in the State of Oregon, USA [35]. Similar conversion factors were reported for rainfall and runoff samples in the USA (DOC = 0.84 × TOC) [36] and in Polish groundwater samples (DOC = 0.93×TOC) [37].

The second step of the proposed scheme for the implementation of BLMs deals with cases where no matching data are available, which is caused by the different monitoring frequencies used for metals and DOC. To overcome this problem, the environmental authorities of England and Wales [7] and France [18] and the EPA of the USA [38] propose the use of the 25th percentile of measured DOC concentrations as the input for the userfriendly BLMs. The proposed approaches in this study for dealing with the absence of matching data includes: (i) the substitution of missing DOC data using available TOC data (DOC = 0.81 × TOC) and (ii) the use of the 25th percentile of the data consisting of measured and substituted DOC values. This allows for the significant broadening of the implementation of BLMs in the Bulgarian monitoring network (Figure 1).

The comparison of the following suggested substitution approaches was performed for 13 water bodies where all data records were available.: (i) DOC calculated by TOC; (ii) 25th percentile of the joint (measured and calculated) Bulgarian monitoring network DOC data; and (iii) 25th percentile of the respective SWB with matching data.

The application of BLMs for assessing the Cu bioavailability showed the closest agreement between the BIO-MET and M-BAT results (Figures 4 and 5). While both methods confirmed the applicability of the proposed substitution approaches, those based on the 25th percentile provided more conservative estimates. It should be mentioned that the correlations were better between the BIO-MET calculated RCRs using matching data and data substituted by the 25th percentile. The PNEC Pro results were unsatisfactory for all BMLs and were an indication that none of the substitution approaches could be used for this user-friendly BLM (Figure 6).

The comparison between the Zn bioavailability results obtained by all BLMs using matching and all types of substituted data confirmed the applicability of all of the proposed approaches (Figures 7–9). Again, the most conservative estimates are provided by the 25th percentile of Bulgarian monitoring network data, followed by the 25th percentile of SWB data and the DOC data calculated by TOC.

The application of BLMs for the assessment of Pb bioavailability revealed a clear separation between the substituted data calculated by TOC on one side and the data substituted by the 25th percentile on the other (Figures 10–12). The comparison of the calculated RCRs by all BLMs between the matching data and the DOC data calculated by TOC showed very good correlations and slopes that were close to one. The comparison between the matching data and the 25th percentile substituted data yielded lower R<sup>2</sup> values and the slopes of the PNEC Pro comparison differed significantly from one. When using the BLM calculations in Tier 2, the results showed a great reduction in the exceedances of the samples (Table 3). For Cu, when BIO-MET was applied, 5 out of 63 samples were found to be exceedances when the measured DOC was used, 4 exceedances were found when TOC value was used (DOC = 0.81 × TOC), 8 exceedances were found when the 25th percentile of the DOC concentration was used and 7 exceedances were found when the 25th percentile of the converted TOC for the respective water body was used, in contrast to the direct comparison of the EQS to the measured copper values (57 out of 63).

Similar results were obtained using M-BAT (7, 6, 7 and 9 exceedances out of 63, respectively) and PNEC Pro (9, 8, 10 and 5 exceedances out of 63, respectively).

The results showed a reduction in the exceeding samples for BIO-MET (2, 5, 7 and 6, respectively), M-BAT (6, 6, 7 and 7, respectively) and PNEC Pro (6, 5, 7 and 6, respectively) compared with the direct comparison of the EQS to the zinc values (14 out of 34 exceeding values).

The results from the application of the BLM tools showed the greatest reduction in the exceeding samples for Pb—only one RCR value was found to be above one—when PNEC Pro was used with the 25th percentile of the Bulgarian DOC data. The direct comparison with the EQS (Tier 1) showed 8 exceedances out of 22 samples.

Generally, the user-friendly BLM tools returned a lower number of RCR exceedances when the measured DOC concentration was used compared with the other approaches, with BIO-MET returning the lowest number of exceedances.


**Table 3.** The exceedances for Cu, Zn, Ni and Pb in the studies of Bulgarian surface water bodies after the application of the BLMs.

#### **5. Conclusions**

The metal bioavailability assessment by user-friendly BLMs requires inputs for the dissolved metal concentrations and the supporting physicochemical parameters, such as DOC, pH and dissolved Ca. The main obstacle in front of the environmental authorities for the implementation of BLMs is the lack of matching data, especially for DOC. In cases where the required data are not fully available, alternative approaches are needed to fill the missing inputs.

In this study, various alternative approaches were tested to support the environmental decision-making bodies. The proposed methodology included three approaches for the substitution of missing DOC data: (i) calculation from available TOC data; (ii) the 25th percentile of the joint (measured and calculated by TOC) Bulgarian monitoring network DOC data; and (iii) the 25th percentile of the joint DOC data for the investigated SWB. For the next step, the proposed substitution approaches were validated by a comparison with the available matching data from the Bulgarian monitoring surface waters network. The comparisons were performed for the three most widely used BLMs: BIO-MET; M-BAT; and PNEC Pro. The described methodology would allow the environmental authorities to

estimate the bioavailability of certain metals and then choose the best possible substitution approach by the implementation of the most appropriate BLM, which would cover as many SWBs as possible without the matching data being available.

The application of the abovementioned methodology to the Bulgarian surface waters monitoring data outlined that the suitability of the substitution approaches decreases in the following order: DOC calculated by TOC > use of the 25th percentile of the data for the respective SWB > use of the 25th percentile of the Bulgarian monitoring network data. The results obtained by the implementation of BIO-MET rendered it the most appropriate tool for the bioavailability assessment of Cu, Zn and Pb in Bulgarian surface waters.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/w14020246/s1, Figure S1: The calculated RCR for Cu using the BIO-MET and measured DOC, the TOC, the 25th percentile of the Bulgarian DOC data and the 25th percentile of the DOC = 0.81 × TOC for the respective water body; Figure S2: The calculated RCR for Cu using the M-BAT and measured DOC, the TOC, the 25th percentile of the Bulgarian DOC data and the 25th percentile of the DOC = 0.81 × TOC for the respective water body; Figure S3: The calculated RCR for Cu using the PNEC Pro and measured DOC, the TOC, the 25th percentile of the Bulgarian DOC data and the 25th percentile of the DOC = 0.81 × TOC for the respective water body; Figure S4: The calculated RCR for Zn using the BIO-MET and measured DOC, the TOC, the 25th percentile of the Bulgarian DOC data and the 25th percentile of the DOC = 0.81 × TOC for the respective water body; Figure S5: The calculated RCR for Zn using the M-BAT and measured DOC, the TOC, the 25th percentile of the Bulgarian DOC data and the 25th percentile of the DOC = 0.81 × TOC for the respective water body; Figure S6: The calculated RCR for Zn using the PNEC Pro and measured DOC, the TOC, the 25th percentile of the Bulgarian DOC data and the 25th percentile of the DOC = 0.81 × TOC for the respective water body; Figure S7: The calculated RCR for Pb using the BIO-MET and measured DOC, the TOC, the 25th percentile of the Bulgarian DOC data and the 25th percentile of the DOC = 0.81 × TOC for the respective water body; Figure S8: The calculated RCR for Pb using the M-BAT and measured DOC, the TOC, the 25th percentile of the Bulgarian DOC data and the 25th percentile of the DOC = 0.81 × TOC for the respective water body; Figure S9: The calculated RCR for Pb using the PNEC Pro and measured DOC, the TOC, the 25th percentile of the Bulgarian DOC data and the 25th percentile of the DOC = 0.81 × TOC for the respective water body.

**Author Contributions:** Conceptualization, T.V. and S.T.; methodology, T.V. and S.T.; software, T.V.; formal analysis, T.V. and S.T.; investigation, T.V.; resources, T.V.; data curation, T.V.; writing—original draft preparation, T.V.; writing—review and editing, S.T.; visualization, T.V.; supervision, S.T.; project administration, T.V. and S.T.; funding acquisition, T.V. and S.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Education and Science, Bulgarian National Science Fund, Grant DN 19/15 and by the National Science Program "Environmental Protection and Reduction of Risks of Adverse Events and Natural Disasters", approved by the Resolution of the Council of Ministers № 577/17.08.2018 and supported by the Ministry of Education and Science (MES) of Bulgaria (Agreement № Д01-363/17.12.2020).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The authors gratefully acknowledge the financial support of the Ministry of Education and Science, Bulgarian National Science Fund (Grant DN 19/15) and the National Science Program "Environmental Protection and Reduction of Risks of Adverse Events and Natural Disasters" (Agreement № Д01-363/17.12.2020).

**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**

