REE concentrations obtained from all the prepared samples as well as the supplemental data from [1] are shown in Table 1. Total REE concentrations (∑REE) for surface water samples range from 0.44 to 7.2 ppb with an average of 1.58 ppb. ∑REE concentrations from WWTP samples range from 0.27 to 0.72 ppb with an average of 0.48 ppb. Total rare earth elements and yttrium (∑REY) concentrations for surface water samples range from 0.57 to 8.4 ppb with an average of 2.0 ppb. ∑REY concentrations for WWTP samples range from 0.40 to 0.73 ppb with an average of 0.52 ppb. The supplemental data did not include Y concentrations, so it is not included in ∑REY. However, Y is often included with REEs as it behaves similarly. Overall, WWTP samples have less dissolved ∑REE and ∑REY than surface water samples.

**Table 1.** Concentrations of rare earth elements in the sampled Triangle Area waters and WWTP influents and effluents. For WWTP samples, -I indicates the influent sample and -E indicates the effluent samples from the named WWTP.

1 MC: Morgan Creek, NC: Northeast Creek, NR: Neuse River (where 3 and 4 refer respectively to the upstream and downstream samples from [1]), JL: Jordan Lake, CR: Cape Fear River, MF: Mason Farm WWTP, T: Triangle WWTP.

> Surface water samples from 2017 show enrichment of heavy REEs (HREE; Gd to Lu) in relation to lighter REEs (LREE; La to Eu) with LaSN/YbSN values ranging from 0.35 to 0.98 except for one sample—MF-E (Table 2). LaSN/YbSN is a ratio used to show relative enrichment in the lanthanides; HREE enrichment is indicated when the ratio is <1. The enrichments from the samples agree with those of the previous study of REEs in the Neuse River and is representative of REE enrichment in the geology underlying the streams of this region [1]. WWTP samples also show HREE enrichment with LaSN/YbSN ranging from 0.048 to 0.44. This enrichment is, on average, smaller than that of the sampled surface waters. There is one surface water sample (JL-2) that shows a slight LREE enrichment with a LaSN/YbSN value of 1.04.

**Table 2.** Select quantities for samples including LaSN/YbSN used to indicate fractionation between LREE and HREE, Gd anomaly (GdSN/GdSN\* as quantified in Equation (1)), ppt of Gd estimated to be of anthropogenic origin (Gdanth), and percent of anthropogenic Gd of total Gd in sample (%Gdanth). For WWTP samples, -I indicates the influent sample and -E indicates the effluent samples from the named WWTP.

1 MC: Morgan Creek, NC: Northeast Creek, NR: Neuse River (where 3 and 4 refer respectively to the upstream and downstream samples from [1]), JL: Jordan Lake, CR: Cape Fear River, MF: Mason Farm WWTP, T: Triangle WWTP.

#### *3.2. PAAS-Normalization and Anthropogenic Gd*

PAAS-normalized REE distributions of samples show a relatively flat pattern (Figure 2). The main exception to this is the positive Gd anomaly present in WWTP samples, samples downstream of WWTPs, and the Jordan Lake samples. The Neuse River samples show the most nonlinear patterns beyond the Gd anomaly with a small negative Ce anomaly in each sample. Both negative Ce anomalies and occasional positive Eu anomalies are common in pristine waters unaffected by anthropogenic inputs [26,27]. Unlike other studies on REEs as microcontaminants, there does not seem to be any contamination from anthropogenic La or Sm [7].

The Gd anomaly (GdSN / GdSN\*), approximated Gd concentration (Gdanth), and percent of total Gd estimated to be from anthropogenic origin (%Gdanth) are also shown in Table 2. The highest Gd anomalies are found in the WWTP samples ranging from 58.0 to 481 with an average of 301. Surface water sample Gd anomalies range from 1.2 to 12.0 with an average of 1.4 in the samples upstream of WWTPs and 5.2 in the downstream samples. Jordan Lake samples have anomalies of 4.0 and 4.5 with the downstream of the lake having a lower anomaly of 1.73. Estimated Gdanth concentrations range from 3.9 to 689 ppt. The highest Gdanth concentrations were all from WWTP samples with a range of 0.39 to 0.69 ppb. The %Gdanth in these samples were estimated to be 98.1% to 99.8%. The upstream samples had the smallest percentages with estimates from 5.9% to 26.9% Jordan Lake samples contained %Gdanth values of 72.5% and 75.5% with a decrease to 36.3% at the downstream.

**Figure 2.** PAAS-normalized REE distributions for: (**A**) Morgan Creek and Mason Farm WWTP samples; (**B**) Northeast Creek and Durham WWTP sample; (**C**) Neuse River samples; and (**D**) Jordan Lake samples.

#### **4. Discussion**

#### *4.1. Anthropogenic Gd and Its Sources*

The occurrence of a small positive Gd anomaly (no higher than 1.4) is common to natural seawater due to the "tetrad effect" from Gd's half-filled 4f electron shell and REE scavenging by organic surface ligands which possess a negative Gd anomaly [28,29]. Therefore, to account for a margin of error, anomalies less than 1.4 will be considered to indicate pristine waters. In this study, only two of the Gd anomalies have a Gd anomaly of less than 1.4. These samples are the upstream samples for Northeast Creek and Morgan Creek. The Neuse River upstream shows anomaly values of 1.5 and 1.8, indicating anthropogenic Gd. This is likely because the Neuse River receives wastewater effluent from other WWTPs on streams that feed into it (such as the Smith Creek WWTP) prior to reaching the WWTP location assessed in this study. All other locations had an anomaly of at least 1.7, indicating their anomalies are all due to anthropogenic inputs. The presence of anthropogenic Gd in Jordan Lake (Gd anomalies of 4.0 and 4.5) and its downstream (Gd anomaly of 1.7) at locations roughly 14, 25, and 32 kilometers respectively downstream of the Mason Farm WWTP display the persistence of anthropogenic Gd in the environment and its efficacy as a wastewater tracer.

As surface water samples have the highest ∑REY, WWTP effluent is likely not a major contributor of overall REY concentrations in the studied streams. However, there is a large increase in estimated %Gdanth between upstream and downstream samples (14.7 to 90.9% for Morgan Creek; 5.9 to 64.5% for Northeast Creek, and 26.9 to 39.5% and 37.3 to 70.9% for the Neuse River). This indicates that WWTP effluent is a large contributor of anthropogenic Gd to the downstream of these rivers. Influent and effluent samples all had an estimated %Gdanth of at least 98.1%.

The most likely source of anthropogenic Gd in WWTP samples are Gd-based contrast agents (GBCAs). GBCAs are the most common MRI contrast agent and have been widely used since approval by the FDDA in 1988. Gd3+ as a free ion can compete with Ca2+ making it toxic to biological systems, thus it is chelated to an organic ligand to make it safe for injection [30]. Though recent studies have shown that a portion of the injected Gd can accumulate in the brain, bones, and tissues, the majority is urinated out of the body within 12 h of exposure [23,31]. Effluent released into surface waters from WWTPs contain all

excreted Gd as conventional wastewater treatment processes do not remove the element. Patients are typically injected with low doses of Gd, but with the number of MRI scans conducted each year and the increase in scans per year, longitudinal studies have shown significant increases in the amount of anthropogenic Gd in waters of developed regions over the years [17,32].

#### *4.2. Health and Ecotoxicological Effects*

The first evidence of health effects from usage of GBCAs were noticed 18 years after their introduction. While GBCAs are not considered to be harmful, repeated injection can lead to noticeable Gd retention in the body and cause nephrogenic systemic fibrosis (NSF) [31,33,34]. With regards to Gd as a potential anthropogenic contaminant in natural waters, not much is known. Natural Gd in surface water is partially removed from the dissolved pool of water due to its tendency to coagulate into a colloid in low salinity waters [16]. As GBCAs injected into the body must be highly stable, anthropogenic Gd in the environment should likely be nonreactive in comparison and stay in the dissolved pool. For example, mussels in locations downstream of WWTPs have been shown to incorporate anthropogenic Sm and La into their shells, but not anthropogenic Gd, suggesting that it is not bioavailable in comparison to geogenic Gd [35]. This distinction along with GBCAs' relatively long half-life of several weeks explains how anthropogenic Gd can be used as a tracer for wastewater [36].

There are few studies that investigate the health effects of Gd alone. Most studies focus on other REEs such as Sm or La or investigate REE patterns. Some recent studies reported that REEs can accumulate in humans due to ingestion of food with high REE concentrations in contaminated soil [37,38]. Thus, drinking water with high REE concentrations could also be a potential source of REE accumulation. While anthropogenic Gd has not been bioavailable in streams, LREEs have been found to bioconcentrate in carps and specific REEs have been shown to be toxic in high enough concentrations to various organisms [39–41]. However, our sample with the highest ∑REE is magnitudes below one of the lowest REE concentration limits (1.39ppm of La) where ecotoxicological effects were found [40]. While the waters sampled in this study do not approach levels that cause known ecotoxicological effects, at least one other study has noted ∑REY levels above such limits [10]. Thus, especially with Jordan Lake as a crucial source of drinking water, Gd and other REY should be monitored as a future micropollutant in the region.

#### **5. Conclusions**

Few studies exist that show the presence of anthropogenic Gd in North American waters. This study analyzed dissolved Gd concentrations in rivers of North Carolina's Triangle Area with attention to WWTPs. Samples upstream and downstream of WWTPs were obtained as well as samples of effluent and influent of two of the three selected WWTPs. We found Gd anomalies in all samples but two of the upstream samples from pristine streams. A large change in the percentage of anthropogenic Gd (%Gdanth) exists from upstream to downstream, demonstrating WWTP effluents large contribution of anthropogenic Gd to natural waters. The presence of anthropogenic Gd anomaly in samples taken from over 30 kilometers downstream of a WWTP indicates anthropogenic Gd's persistence in the environment and its utility as a wastewater tracer. Despite that there are no current health hazards or ecotoxicological effects at the measured concentrations in this region, Gd anomaly and total REE concentrations should be monitored because of their ability to have toxic health effects.

**Author Contributions:** J.M.Z. and C.C. analyzed data and wrote the manuscript with the support and supervision of X.-M.L., Q.W. collected samples and processed them. J.M.Z. wrote the manuscript with inputs from all the authors. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the University of North Carolina at Chapel Hill to X.-M. Liu.

**Data Availability Statement:** Not Applicable.

**Acknowledgments:** The authors thank field sampling help from WWTP staff in the Triangle Area, North Carolina. X.-M.L. and C.C. would like to acknowledge funding support from the Department of Geological Sciences, University of North Carolina at Chapel Hill. J.Z. would also like to acknowledge the support from UNC's Institute for the Environment IDEA undergraduate research program.

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

#### **References**

## *Article* **Modeling Management and Climate Change Impacts on Water Pollution by Heavy Metals in the Nizhnekamskoe Reservoir Watershed**

**Yury Motovilov \* and Tatiana Fashchevskaya**

Laboratory of Regional Hydrology, Water Problems Institute of Russian Academy of Sciences, 119333 Moscow, Russia<EMAIL_ADDRESS>
**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** The semi-distributed physically based ECOMAG-HM model was applied to simulate the cycling of the heavy metals (HM) Cu, Zn, and Mn, and to identify spatial and temporal patterns of heavy metal pollution in water bodies of a large river catchment of the Nizhnekamskoe reservoir (NKR) in Russia. The main river of the catchment is the Belaya River, one of the most polluted rivers in the Southern Urals. The model was tested against long-term data on hydrological and hydrochemical monitoring of water bodies. It is shown that the pollution of rivers is formed mainly due to diffuse wash-off of metals into rivers from the soil-ground layer. Numerical experiments to assess the impact of water economic activities on river pollution were carried out by modeling scenarios of changes in the amount of metal discharged with wastewater, a disaster with a salvo discharge of pollutants, and the exclusion of anthropogenic impact on the catchment to assess selfpurification of the basin. Modeling of chemical runoff in accordance with the delta-change climatic scenario showed that significant changes in water quality characteristics should not be expected in the near future up to 2050.

**Keywords:** ECOMAG-HM model; river basin; heavy metals; pollution modeling; climate and management change impacts

#### **1. Introduction**

Heavy metals (HMs) are one of the most dangerous pollutants of the environment. The main anthropogenic sources of HM are various enterprises of ferrous and non-ferrous metallurgy, mining enterprises, chemical enterprises, fuel installations, cement plants, electroplating industries, and transport. HM released into the environment are not decomposed, but rather only redistributed between its individual components. HMs accumulate relatively easily in soils, but their removal is slow. For example, the period of semi-removal of zinc from the soil is up to 500 years, while copper takes up to 1500 years [1]. In micro quantities, HMs are necessary for the processes of vital activity of organisms. Excessive amounts of HM are easily accumulated by the organs and tissues of hydrobionts and humans, adversely affecting their health and increasing the environmental risks of morbidity [2,3]. One of the main factors of the temporary detoxification of natural waters from HM ions is their adsorption by suspended particles, deposition in the form of poorly soluble inorganic compounds, and burial in bottom sediments.

HMs are one of the most common pollutants in the surface waters of Russia [4]. Their content in water bodies often exceeds the maximum allowed concentrations (MAC) established by the Russian Federation (Table 1), not only for fisheries and water use, but also for drinking water. Cases of high and extremely high pollution of HMs are recorded in individual water bodies, at which their concentrations reach values of 30–50 MAC or more.

**Citation:** Motovilov, Y.; Fashchevskaya, T. Modeling Management and Climate Change Impacts on Water Pollution by Heavy Metals in the Nizhnekamskoe Reservoir Watershed. *Water* **2021**, *13*, 3214. https://doi.org/10.3390/ w13223214

Academic Editors: Guilin Han and Zhifang Xu

Received: 27 September 2021 Accepted: 8 November 2021 Published: 12 November 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/).

In recent decades, a significant decrease in the amount of wastewater discharged through point sources and, accordingly, the amount of pollutants, has been observed in most river catchments of the Russian Federation [5,6]. However, the expected decrease in the HM content in many water bodies is not marked [4,7]. This discrepancy indicates the need to deepen the understanding of the influence of various processes occurring in river basins on the formation of river water quality. Such investigations are complicated by a number of circumstances. Firstly, the existing network of state hydrochemical monitoring of water bodies in Russia has large spatial and temporal sparsity, which makes it difficult to reliably identify the regularities of the hydrochemical regime and the formation of the quality of river water under the influence of changing natural and anthropogenic factors. Secondly, due to the imperfection of the accounting system of water user enterprises that provide statistical reports on wastewater discharges, it is quite difficult to correctly assess the contribution of point sources to water pollution [8–10]. Third, monitoring of diffuse sources of pollution in the catchment area is not carried out. The best indicators of diffuse pollution are medium and, especially, small rivers, but hydrochemical monitoring points are mainly located only in large rivers in places of organized wastewater discharge: Individual large industrial enterprises, thermal power plants, nuclear power plants, as well as in the locations of cities and large settlements. Prospects for studying the patterns of HM content in river basins, predicting its changes with possible changes in climate and economic activity, and planning water protection measures in low-light conditions with observational data are associated with the construction of mathematical models describing these patterns.

Modeling the water quality of river catchments has a relatively long history. Intensive development of methods for assessing river pollution associated with the use of fertilizers in agriculture began in the 1960s. These methods were based mainly on empirical relationships for calculating pollutant washouts from agricultural fields using correlations with climatic and physico-geographic parameters [11].

One of the first conceptual field-scale models of water quality was the CREAMS model (Chemicals Runoff and Erosion from Agricultural Management Systems) [12]. The model includes a description of processes in hydrology, erosion, plant, and chemistry sub-models and is used in agricultural management practices. The main limitations of the CREAMS are related to the small size of the simulated area bound by the agricultural field. The ideology of the CREAMS model was the impetus for the creation of many modifications of conceptual models, describing water quality formation in areas with nonpoint spatial pollution (reviews can be found in [13–16]).

Since the 2000s, semi-distributed models of the river runoff and water quality formation have been used to evaluate the effects of alternative management decisions on water resources and non-point source pollution in river basins. Among others, the Soil and Water Assessment Tool (SWAT) is widely used for the simulation of river discharge and

pollution loads at different scales: From the continental scale [17], through the regional scale [18], to the scale of individual river catchments [19,20]. There is evidence in using semi-distributed models for the evaluation of measures to reduce contaminant loads and improve the ecological status of water bodies in Germany (SWIM, Soil and Water Integrated Model [21]; MONERIS, Modeling Nutrient Emissions in River Systems [22]; METALPOL model [23]) and in Sweden (HYPE, Hydrological Predictions for the Environment [24]) in the framework of the implementation of the European Water Framework Directive. In semi-distributed models, the watershed is divided into several subbasins—Hydrological Response Units (HRUs). Each HRU is divided into several layers, and the vertical and horizontal fluxes of water, sediment, and nutrients loadings that move into the river network are modeled. Previous authors [24] drew comparisons between several semi-distributed water quality models based on their complexity, input data constraints, and performance of the models. Overviews of the water quality formation models can be found in [25–28].

Most of the semi-distributed models of water quality formation in river basins were used mainly for biogenic elements (nitrogen and phosphorus). Much less research is devoted to modeling the formation, transport, and transformation of heavy metals in rivers, which may create serious environmental problems in the region. Such models are usually developed for mining areas [29] with widespread waste rock dumps, tailing dumps [19], and contaminated soil [30], as well as in areas with a wide spread of "dirty" industries [31]. For example, in [19], the load of HM transported by the highly contaminated mine drainage river has been assessed using a combination of the SWAT model (for hydrological simulation) and an empirical equation statistically relating local discharge with the concentration of elements. A statistical approach using a multilayer perceptron neural network model was also used in [29] to forecast the content of HM pollutants in stormwater sediments on the basis of atmospheric data and catchment physico-geographic characteristics. In [31], a detailed heavy metal transport and transformation module was developed and combined with the SWAT model for the purpose of simulating the fate and transport of metals at the watershed scale. This modification of the SWAT-HM model was calibrated and validated to simulate Zn and Cd dynamics in the watershed, which has been impacted by mining activities for decades. The authors noted that the developed model has a high potential for application in environmental risk analysis and pollution control to provide scientific support for pollution control and remediation decisions.

One physically based model of runoff and water quality formation in river basins developed in Russia is the ECOMAG-HM model (Ecological Model for Applied Geophysics– Heavy Metals) [30,32]. The model is well adapted to the description of processes in river basins with a snow-type water regime, as well as the structure and composition of hydrometeorological information in Russia. The hydrological module of the ECOMAG model has been repeatedly tested on many large river basins and is used in the practice of hydrological calculations, forecasts, and water resource management in the Russian Federation [33]. The hydrochemical module was applied to assess possible changes in river water pollution under various scenarios of the Pechenganikel mining complex in the north-west of the Kola Peninsula, which creates environmental problems in the region [30]. Besides, the model was verified using hydrochemical monitoring data in a large river basin to study the regime of HM concentrations in the watercourses of the Nizhnekamskoe reservoir (NKR) catchment in Russia (for copper [32] and zinc [34]) [35].

The main purpose of the study is to assess the ECOMAG-HM model's potential to study spatio-temporal patterns of HM runoff formation in a large river basin (with the example of the NKR basin) to quantify scenarios of the water protection measures under possible changes in water economic activity and climate. The specific objectives of the study were (a) to evaluate the capabilities of the model for assessing HM pollution of river waters in various sections of the river network to solve the problem of designing a hydrochemical monitoring network; (b) to assess the contributions of diffuse and point anthropogenic sources to HM pollution in the NKR watershed; (c) to assess changes in water quality of water bodies under various scenarios of anthropogenic load on the river basin; and (d) to assess changes in the hydrological and hydrochemical regime of rivers under a scenario of possible climatic changes. In this paper, in the section of the testing results of the model, more attention is paid to the study of the spatio-temporal patterns of manganese in river waters. The results of the copper and zinc simulations in the river system are used for comparison with the manganese regime, as well as to illustrate the HM regimes under various scenarios of changes in water economic activity and climate.

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

#### *2.1. Study Area: The Nizhnekamskoe Reservoir Watershed*

The Nizhnekamskoe reservoir was created in the valley of the Kama River (the largest tributary of the Volga River) in 1979. The area of the NKR watershed between the Nizhnekamsk (Naberezhnye Chelny city) and Votkinsk (Tchaikovsky city) hydropower plants is 186,000 km2, most of it (142,000 km2) occupied by the basin of the Belaya River River (Figure 1). About 2/3 of the catchment area in the western and central parts is flat territories, and the eastern part is the Ural folded mountain region. The average height of the basin is 392 m, and the highest is 1654 m (Yaman-Tau Mountain). Forest occupies about 50% of the territory.

**Figure 1.** Map of the Nizhnekamskoe Reservoir (NKR) watershed and location of monitoring posts: Hydrological gauges (circles) and hydrochemical monitoring points (triangles). The inset shows the location of the NKR watershed (pink) in the Volga basin (gray) and Eurasian continent.

The climate of the territory is continental. There are a number of transitions from the climate of semi-arid steppe regions, where the annual precipitation ranges from 300–400 mm and the average annual air temperature is about 3 ◦C, to more humid areas (north-eastern and eastern mountain-forest), where the annual precipitation exceeds 600 mm and the average annual air temperature is below 1 ◦C. The rivers are mainly fed by snow. More than 60% of the annual runoff passes during the spring flood. The average annual lateral inflow of water to the NKR is 36.5 km3, of which 26.1 is accounted for by the flow of the Belaya River. To compensate for the shortage of water resources, which often occurs in low-water years, about 400 reservoirs and ponds with a volume of more than 100,000 m3, as well as many smaller ponds, operate in the basin. The largest reservoirs are Pavlovskoye on the Ufa River, Nugushskoye on the Nugush River, and Yumaguzinskoye on the Belaya River.

A feature of the geological structure of the territory is the wide distribution of ore deposits and the significant concentration of ore elements in rocks. More than 3000 deposits and manifestations of various types of minerals have been discovered in the NKR basin [36]. The soils of the catchment area (chernozems, sod-podzolic, gray forest) are characterized by a high content of humus and heavy mechanical composition. Well-drained mountain soils

are common in the east of the region. Soils inherit the chemical composition of soil-forming rocks. In the soil, HMs are included in the composition of humic substances, forming strong complexes with humic and fulvic acids, and are not removed from it for hundreds and thousands of years [1]. The heavy mechanical composition of soils contributes to the increased accumulation of HM. About 60% of the studied territory is erosively dangerous land as a result of water and wind erosion. Eroded soils enriched with trace elements contribute to the entry of HMs into water bodies with sediments [36,37].

The industrial development of the region began almost 300 years ago, was associated with the development of mineral deposits, and took place without taking into account environmental restrictions. "Ancient" anthropogenic-transformed mining landscapes, modern industrial enterprises for the extraction and processing of mineral resources, large settlements, and their infrastructure objects are sources of additional metal supply to the catchment area. In the eastern part of the catchment, such sources are mining enterprises, while in the western and central parts, they are enterprises of oil production, oil refining, chemistry and petrochemistry, metallurgy, mechanical engineering and energy, and objects of storage of production and consumption waste.

As mentioned above, in recent decades, the volume of pollutants entering water bodies from controlled point sources has been decreasing for reasons that are not fully established: Perhaps due to the decrease in water consumption [5,6] or due to the decrease in the number of monitored water users providing statistical environmental reports of industry, settlements, and communal services on the 2-TP (vodkhoz) form [8]. In [38], an assessment of the correctness of the information contained in the reporting forms of 2-TP (vodkhoz) was evaluated. According to the results of the assessment, a large number of cases were revealed when the amount of metal actually contained in wastewater significantly exceeded the statistical data, and in some cases, the excess was 1–2 orders of magnitude. The revealed differences can be explained by the shortcomings of the existing system for monitoring the composition of wastewater at enterprises. In the absence of automatic continuous monitoring of the HM content in wastewater at enterprises, the assessment of the HM discharge is carried out either on the basis of individual episodic water samples, or indirectly by the number of products produced. Such approaches inevitably introduce significant errors in determining the actual HM discharges from point sources [9,10].

#### *2.2. ECOMAG-HM Model*

The semi-distributed physically based ECOMAG-HM model [30,32] was developed to simulate the spatio-temporal dynamics of hydrological and hydrochemical cycle components in large river basins. The model operates with a daily time step and consists of two main blocks: The hydrological submodel of runoff formation and the hydrochemical submodel of pollutant migration and transformation. The first submodel describes the main processes of a land hydrological cycle: Vertical fluxes (infiltration of rain and snowmelt water into the soil, percolation of water into deeper soil horizons and groundwater zone, evapotranspiration), changes in water and energy contents in different components of the geosystem (snow cover formation and melting, soil freezing and melting, surface water, soil moisture, groundwater level), and lateral fluxes (formation of surface, subsurface, groundwater flow), which form the river runoff.

The hydrochemical submodel describes the processes of migration and transformation of conservative pollutants in catchments. It accounts for pollutant accumulation on land surfaces, dissolution by melt and rainwater, penetration of dissolved forms of pollutants into soil, and interactions with soil solution and soil solid phases. Pollutants are transported by surface, subsurface, and groundwater flows to the river network and form the lateral diffuse inflow into rivers. The submodel of the transfer and transformation of pollutants in the river system include the lateral diffuse inflow and load of pollutants from point anthropogenic sources to the rivers and routed through the river network to the outlet of the catchment, taking into account the exchange of pollutants between the river water and riverbed.

The Information-Modeling Complex (IMC) ECOMAG was developed for easy use and adaptation of the model for river basins under various projects of information support in water resources management, forecasting water and hydrochemical regimes, and research. The IMC computer technology contains the calculation module of the mathematical model, databases, and instruments for information and technological support of this module. Thematic digital maps (a digital elevation model, hydrographic network, soils, and landscapes) are applied for automated division of the river basin into elementary watersheds (model cells, analog HRU in SWAT) and a schematization of the river network using the specific GIS tool Ecomag Extension. Data on hydrometeorological, hydrochemical, and water management monitoring as well as instruments for geoinformation processing of these spatial data are involved in the IMC ECOMAG for calibration of the model parameters, verification, testing of the model, and calculations for various projects.

The concept, algorithms, equations, and results of the application of the hydrological module in operative water management and forecast practice were described in many studies [33] and applied to river basins of different scales (including the largest in the northern hemisphere: The Volga, Lena, Selenga, Amur, Mackenzie, etc.) located in different geographic zones with different runoff formation conditions and hydrological regimes of water bodies. The geography of application of the hydrochemical module of the ECOMAG-HM model for environmental problems is more moderate. A description of the hydrochemical module and results of validation of the model were published in [30,32,34,35].

#### *2.3. Data Preparation and Model Setup*

For calculations on the model, boundary conditions were set in the form of daily spatial meteorological fields, concentrations of heavy metals in atmospheric precipitation, and pressure groundwater feeding the upper groundwater zone. Daily data on 56 meteorological stations located in the NKR watershed for the 1979–2011 period were used to construct meteorological fields of air temperature, precipitation, and air humidity deficit. The concentrations of metals in atmospheric precipitation were set by constant values based on the weighted average concentrations given in [39]. Furthermore, the concentrations of metals in the pressure groundwater were set as constant values, based on the data given in [40]. In addition, as information on the point anthropogenic sources of river water pollution, mean annual data on metal discharges with wastewater in 12 large settlements in the Belaya River were set on the basis of State statistical reporting forms 2-TP (vodkhoz) for the 2002–2007 period. Maps of heavy metals content in the arable soil layer on the territory of the Republic of Bashkortostan [41], which occupies a large part of the NKR catchment area, were used to set the initial conditions for the hydrochemical module of the ECOMAG-HM model.

To calibrate the model parameters and verify the hydrological module, data on daily river discharge at 5 hydrological gauges were used. The data of observations of HM concentrations at points of the Russian State Monitoring Survey on the Belaya River and its tributaries were used to calibrate the parameters and test the hydrochemical module: For copper and zinc, at 35 points for the 2004–2007 period; for manganese, at 26 points for the 2002–2007 period. As a rule, observations of the metal content in river waters were carried out once a month.

The spatial schematization of the catchment area and the river network in the NKR basin was produced using the Ecomag Extension tool on the basis of a digital elevation model. As a result, 503 elementary watersheds were identified in the catchment; their average area is about 400 km2.

The values of most of the model parameters were set from the IMC ECOMAG databases using available maps of land surface characteristics (relief, soil, vegetation) or were determined on the basis of literary data [42]. Some parameters of the model were corrected using the calibration procedure by comparisons with daily runoff hydrographs

and observations of HM concentrations in river waters at observation stations within the framework of hydrological and hydrochemical monitoring of water bodies. Numerical experiments have shown that the most sensitive parameters of the hydrological submodel are the saturated hydraulic conductivities of soils, the evaporation index, and the degree– day coefficient of snow melting. In the hydrochemical submodel, the constants of the sorption equilibrium of HM in soils and the parameter accounting for the exchange of metals between river water and bottom sediments were refined by calibration.

#### *2.4. Scenarios of the Water Economic Activity Changes Experiments*

In order to assess the impact of water economic activity changes on river water pollution, a series of numerical experiments was evaluated, in which the following scenarios were performed:

#### 2.4.1. Scenario 1. Change in the Amount of Wastewater Entering the River Waters

Water pollution was simulated under various scenarios of changes in the amount of metal discharged in wastewater at 0.1, 10, 20, 40, 60, 80, and 100 times in relation to the existing level (represented as 2-TP (vodkhoz) forms) in all localities. In fact, such changes may be due to an increase in industrial production, or in some cases, as a result of more objective and correct presentation of information in the 2-TP (vodkhoz) reports on the amounts of HM discharged by water users in some localities [9,10,38].

#### 2.4.2. Scenario 2. The Simulation of a Disaster

The simulation of an extreme situation that led to a significant amount of HM entering the watercourse was carried out. The reason for this situation may be accidents at liquid waste storage facilities, multiple discharges of wastewater as a result of unstable operation of enterprises, unauthorized discharges of contaminated wastewater without treatment, etc. The algorithm of events was as follows. As the input to the model, various amounts of metals were dumped at a certain point in the river network during the day. Then the concentrations of metals in the river network below the emergency point were calculated.

#### 2.4.3. Scenario 3. Exclusion of Anthropogenic Impact on the Catchment Area

The estimation of the time scale of self-cleaning of the catchment from HM in the scenario of complete exclusion of anthropogenic impact on the NKR basin was performed. It was assumed that HMs were not dumped in the wastewater, and the entry of metals into the atmosphere with industrial emissions was also excluded (the concentration of metals in precipitation was equal to zero). The concentration of metals in the pressure groundwater decreased to the lowest level in the range of typical values for the underground water of the studied region. Under these conditions, numerical experiments were carried out to estimate the dynamics of heavy metals in river waters for 400 years ahead. A series of historical data of meteorological observations at weather stations for 33 years, from 1979 to 2011, was used as a meteorological forcing. The final results of calculations of the river basin characteristics for 31 December 2011 were recorded at the starting point of 1 January 1979, and thus the calculations for the 33-year time series were repeated many times.

#### *2.5. Description, Evaluation, and Processing of Climate Scenario Data*

To construct a scenario of the future climate for the NKR basin, the trends in annual temperature, precipitation, and runoff for several meteorological and hydrological stations located in the watershed area over the period of 1956–2007 were analyzed. It was found that at the turn of the 1980s, there were changes in climate signals between the climate parameters for two periods. The average air temperature for the 1980–2007 period increased by about 1 ◦C compared to the previous period of 1956–1979. The increase in annual precipitation over the 1980–2007 period compared to the previous period averaged for the whole considered meteorological stations was about 10%. The difference in the mean annual runoff was about 10% in the gauges of the Belaya River during these periods (see examples in Figure 2). Modeling the water and chemical runoff for the future period, up to 2050, was carried out using the so-called "delta change" (DC) climatic scenario, assuming changes in climate parameters to be the same as for the previous observation periods. So, the obtained increments of the annual daily air temperature of 1 ◦C and daily precipitation multiplied to 1.1 were introduced to the corresponding historical daily series of climate parameters for the last observation period. It should be noted that similar trends in runoff changes for the region under consideration over the past decades and in the near future are given in the review [43].

**Figure 2.** Changes in mean annual air temperature (**a**) and precipitation (**b**) in Ufa station and the mean annual discharges (**c**) of the Belaya River at Birsk gauge for the 1956–2007 period (dotted lines are trends).

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

#### *3.1. ECOMAG-HM Testing of the Model*

#### 3.1.1. Model Verification

The hydrological submodel was calibrated and verified by comparing the observed and simulated runoff hydrographs at several hydrological stations in the NKR basin over the 1978–2013 period. The results of the tests showed satisfactory agreement between the observed and simulated runoff characteristics according to the Nash–Sutcliffe efficiency criterion for the daily hydrographs, and the PBIAS criterion for the annual runoff volumes [32,33]. Figure 3 shows these comparison results at four stations over the 2002–2007 period, which will be later used for the simulation of hydrochemical characteristics.

**Figure 3.** Observed (blue) and simulated (red) hydrographs of the Belaya River and inflow into NKR.

An additional successful test of the model for the rivers not covered by hydrometric observations was performed by comparing the map of the mean annual specific runoff estimated on the basis of averaging the simulated daily runoff in the nodes (centers of gravity) of elementary watersheds over many years with the map of specific runoff given in SNIP-2.01.14-83 [32,44]. The latter is used in Russia for the estimation of the hydrological characteristics in the design of hydraulic structures on rivers in the absence of hydrometric observations.

The hydrochemical submodel was tested by comparing observed and simulated concentrations of copper [32] and zinc [34] in different sections of rivers of the NKR basin. The rare frequency of observations at the hydrochemical monitoring points (usually up to 12 observations per year), as well as significant errors in determining HM concentrations in water samples (~50%), do not allow the full use of traditional statistical hydrological criteria for comparing the results of model calculations of intra-annual HM concentrations and episodic data of hydrochemical measurements [45]. Characteristics averaged over longer time periods can be determined more reliably based on such data.

Examples of the comparison of the intra-annual distribution of Mn concentrations (daily and averaged over quarter (three-month) periods) are shown in (Figure 4). As can be seen, seasonal changes in the Mn content are characterized by increased concentrations during high-water periods (snowmelt and rainfall runoff) and a decrease during the summer and winter low-flow periods. Correlation analysis showed that the maximal correlation coefficients, R, between the observed and simulated concentrations at all monitoring points in the basin, averaged by three months (quarters), were obtained for copper in the second quarter (the spring flood period, R = 0.55), for zinc in the second and third quarters (R = 0.50 and R = 0.69, respectively), and for manganese in all quarters (R varies in the range from 0.68 (fourth quarter) to 0.85 in the first quarter).

Spatial differences in the long-term mean annual HM concentrations in the Roshydromet hydrochemical monitoring points located in the order from the upper reaches of the Belaya River toward its inflow into the NKR differ significantly (Figure 5). Table 2 shows the correspondences between the observed and simulated HM concentrations and the range of changes in monitoring points. The correlation coefficients between simulated and measured concentrations are 0.58 for copper, 0.60 for zinc, and 0.80 for manganese.

Additional verification of the model was carried out by comparing the calculated and actual inflow of pollutants into the NKR. The model allows for calculating the balance of chemical components entering the river network and estimating the contribution of point (controlled wastewater discharges) and diffuse sources to river water pollution. As the information on point sources, yearly data on metal discharges with wastewater in 12 large localities were used. The actual inflow of metals from the catchment area to the NKR was estimated by multiplying the observed annual lateral inflow of water into the reservoir by the measured average annual concentration of metal at the last hydrochemical monitoring station on the Belaya River before its confluence with the NKR (Dyurtuli station) (Table 3).

**Table 2.** Observed and simulated HM concentrations and its ranges of changes in monitoring points (μg·L−1).

**Table 3.** Simulated components of HM load into the river network from the NKR watershed and the actual and simulated HM load into the NKR (t·year−1).

**Figure 4.** Observed (green) and simulated (red) Mn concentrations at certain hydrochemical monitoring points over the 2002–2007 period: On the left, the dynamics with daily resolution; on the right, the mean intra-annual concentrations by quarters.

**Figure 5.** Observed (green) and simulated (red) long-term mean annual heavy metal (HM) concentrations at the hydrochemical monitoring points located in the order from the upper reaches of the Belaya River toward its inflow into the NKR: (**a**) Cu; (**b**) Zn; (**c**) Mn. The arrow shows the direction of river flow.

Model balance calculations show that the total wash-off of metals in the rivers from the NKR catchment averaged over the period of simulation is formed by approximately 80–85% due to leaching from the soil-ground layer. Surface wash-off of metals into rivers usually does not exceed 20%.  A significant proportion of metals washed out from the catchment area accumulate with sediments in the river bottom: On average, ~46% copper, 62% zinc, and 30% manganese. These estimates are in good agreement with the results of HM retention in the Elbe River systems obtained using the METALPOL model [23]: For copper, on average, 51% with a range of variation from 35 to 60% for different tributaries; and for zinc, on average, 45% with a range of variation from 30 to 54%.

An example of the inter-annual dynamics of various sources of river water pollution and a comparison of the actual and simulated Mn load in the NKR is shown in Figure 6. The differences between the compared annual values can be largely explained by significant errors in determining the concentration of heavy metals in river water.

**Figure 6.** Inter-annual dynamics: (**a**) Simulated components of Mn load in the river network; (**b**) The actual and calculated Mn load in the NKR.

Inter-annual changes in the HM load in the river network and NKR are characterized by the increase in high-water years (2002, 2007), decrease in medium years (2003, 2005), and minimum values in low-water years (2004, 2006). Such inter-annual changes in the HM load occur when HMs enter the river network mainly from diffuse sources from the catchment area, which are affected by climatic factors. Thus, the results presented in this section indicate that, on the whole, the model adequately reproduces the main spatio-temporal patterns of HM runoff formation in the river basin and HM content in the river network.

#### 3.1.2. Mapping of the Simulated HM Characteristics

Regularities of water quality formation in river basins are associated with the study of genetic components of river and chemical runoff; therefore, among the mapped characteristics, we will consider the HM concentrations in river waters and the HM wash-off with various genetic components of river runoff. The algorithm for determining spatial distributions of components of HM runoff based on the developed model was as follows. Daily fields of weather characteristics in the catchment area were constructed using data of meteorological stations, and the ECOMAG-HM model was applied to simulate runoff and HM removal in the local river network by surface, soil, and ground flow for all elementary watersheds. By averaging the daily fields over a multi-year period, maps of the mean annual specific HM load in the river network by various genetic components were obtained for cooper [32] and zinc [34].

The analysis of these maps demonstrates the spatial differences in the predominance of various genetic components in the formation of hydrochemical metal runoff in the NKR basin. Over most of the catchment area, the total HM load into the river network was formed mainly due to subsurface (soil and ground) runoff. The HM wash-off by surface water is almost an order of magnitude smaller, except for the western province of the basin, where the wash-off by surface runoff was comparable or slightly higher than the HM leaching by subsurface runoff.

The map of the mean annual Cu, Zn, and Mn load in the river network shows the close correlation with the spatial distribution of the initial HM content in soils (Figure 7) due to the predominance of HM leaching from the underground zone in the water quality formation.

**Figure 7.** Initial HM contents in the soils, long-term mean annual HM wash-off, and long-term mean annual HM concentrations in the river network.

In turn, the spatial distribution of HM load affects the distribution of HM concentrations in the river network. A map of the simulated mean annual HM concentrations in river water was created (Figure 7) to analyze the spatial distribution of HM content in the NKR drainage basin, including the rivers not covered by hydrochemical observations. This map was obtained by averaging the simulated daily HM concentrations in the elements of the modelled river network over the period of simulation. The line thickness indicates the HM concentration in the river network in accordance with the legend.

As can be seen from Figure 7, small rivers that mainly flow near foci with a high content of HM in soils have higher levels of HM pollution. These areas with increased river pollution are delineated with shafts. Figure 7 clearly shows that a narrow area of increased concentrations of zinc in rivers is formed near the focus of increased zinc content in the soil in the area of the Ufa River catchment. The area with a high copper content in the rivers is noticeably larger and mainly covers the north-eastern part of the NKR basin. The distribution of the increased content of manganese in the soil in the catchment area is spotty; therefore, foci of increased concentrations of manganese in the river water are found everywhere, with the exception of the southwestern part of the basin. The maximum mean annual HM concentrations in these streams exceed the similar average for the other

part of the NKR basin by 2–6 times and reach 11 μg·L−1 for copper, 29 μg·L−1 for zinc, and 166 μg·L−1 for manganese. When small, polluted tributaries flow into larger rivers, HM concentrations in river waters decrease as a result of dilution by cleaner water in larger watercourses. A similar effect was obtained in the results of simulating HM concentrations in river waters near polluted areas in the Liuyang River basin in China impacted by mining activities using the SWAT-HM model [31]. As already noted, hydrochemical monitoring points are usually located on large and medium-sized rivers in the immediate vicinity of large settlements and industrial enterprises, so the existing monitoring network does not cover small rivers that are most susceptible to diffuse pollution.

The analysis of simulated results shows that the water quality in the river monitoring points does not meet the MAC standards for fishing use regarding copper (by 2–7 times) and manganese (by 4–19 times). The entire range of changes in zinc concentrations is within the MAC standard. The water quality regarding copper and zinc in all monitoring points is suitable for drinking use. The content of manganese in river water exceeds the MAC standards for drinking use only at several points in some seasons of the year.

#### *3.2. Modeling the Impact of Economic Activity on River Water Pollution* 3.2.1. Change in the Amount of Metals Discharged in Wastewater

**Intra-annual distribution.** Figure 8 shows the simulated dynamics of the copper content in the Belaya River at Ufa city with actual wastewater discharges, as well as with the 60-times increase in copper content in wastewater. At the current level of the copper content in wastewater, its maximum concentration in the river water is simulated during high water periods (snowmelt and rainfall runoff) as a result of the metal load, mainly from diffuse sources. When the multiplicity of the increase in the content of metals in wastewater relative to the existing level exceeds a certain critical value, the following pattern occurs: During periods of high water, the lowest concentrations of metals begin to be observed in contrast to the original case (Figure 8). This occurs as a result of diluting highly polluted industrial wastewater in rivers with significantly less polluted snow and rainwater. Nomograms can be constructed for hydrochemical monitoring stations, showing the "critical" levels of anthropogenic load (the multiplicity of increase in metals discharged from wastewater relative to the existing level), at which such changes occur based on the results of numerical experiments. In particular, in the example shown in Figure 9, the relationship between the river discharges and maximal concentrations of copper in the Belaya River at Ufa city is reversed, when the discharge of copper into the river increases by 20 or more times.

**Figure 8.** Simulated hydrograph (blue) and dynamics of the copper content in the Belaya River (Ufa city) with the actual content of copper in wastewater (1) and with an increase in the copper content by 60 times (60) over the 2004–2007 period.

**Figure 9.** Relationship between water discharges at the Belaya River (Ufa city) and the concentration of copper in river water with the actual content of copper in wastewater (1) and under scenarios of its change at 0.1, 10, 20, 40, 60, 80, 100 times.

Numerical experiments have shown (Figure 9) that with an increase in discharge of the Belaya River near the city of Ufa over 3500 m3·s−1, copper concentrations in river water asymptotically approach the value of ~6 μg·L−1 under various scenarios of wastewater loads. In other words, the contribution of point sources of pollution due to dilution is minimized during periods of high water, and the concentration of metals in river water is mainly determined by the maximum exchange capacity of the catchment leaching of HM from the watershed by surface and subsurface flow during intensive snowmelt or rainfall that forms the maximum runoff. This constant concentration can be named the equilibrium concentration of diffuse saturation of the river basin for this type of HM.

**Mean annual characteristics.** With an increase in the HM content in wastewater, the contribution of point sources to the river's pollution in the NKR basin increases (Figure 10a). Thus, with a 20-fold increase in HM content in wastewater, the mean annual contribution of point sources to river water pollution is ~20% for copper, 44% for zinc, and 2% for manganese. A 100-fold increase in HM load increases the contribution of point sources to river pollution to ~53%, 80%, and 8%. An increase in the content of HM in wastewater also leads to an increase in their inflow by river water to the reservoir (Figure 10b). Thus, a 100-fold increase in metals causes an increase in the inflow into the NKR of copper by ~2 times, zinc by ~4 times, and manganese by 1.1 times.

#### 3.2.2. Simulation of a Disaster with a Salvo Discharge of Pollutants

Figure 11 shows simulated changes in the maximum concentration of copper along the length of the river under various scenarios of copper salvo discharge (5000, 4000, 3000, 2000, 1000, 500, and 100 kg·day−1) into the Belaya River near Salavat city. As the contamination spot is removed from the site of the salvo discharge, the concentration of copper decreases intensively due to the dilution by the tributaries that flow into the Belaya River, and copper deposits on the river bottom along with sediments. As early as passing Ufa city (~300 km from the accident site), the maximum concentrations of copper are reduced by about two orders of magnitude.

**Figure 10.** Simulated anthropogenic impact on the NKR basin with changes in the metal content in the wastewater: (**a**) Contribution of point sources to the river water pollution, (**b**) wash-off of metals to the NKR.

**Figure 11.** Changes in the maximum concentration of copper at several hydrological stations along the Belaya River under various scenarios of copper dumping in the Salavat city (5000 kg, . . . , 100 kg).

3.2.3. Simulation of Self-Cleaning of the NKR Watershed under Exclusion of Anthropogenic Impact on the Catchment

The estimation of time scales of the catchment's self-purification from HM according to scenario 3 was carried out by simulating the dynamics of HM reserves with the complete exclusion of any anthropogenic impact on the basin over a 400-year period. Based on the maps of metal content in the soil [41], it was determined that in the NKR catchment, the initial reserves of copper are ~1000, zinc ~550, and manganese ~11,350 thousand tons (Figure 12a). The annual flow of copper, averaged over a 33-years period, from the catchment to the river network at the beginning of the calculation period was 147, zinc ~156, and manganese ~1745 t·year−1 (Figure 12b). After 200 years, calculations showed a decrease in the removal of copper, zinc, and manganese from the catchment area, respectively, to 140, 147, and 1700 t·year−1. After 400 years, the average annual removal of metals from the river network continued to decline, and amounted to 135, 140 t, and 1640 t year−1, respectively. Thus, in the absence of external impacts, the catchment area was slowly cleared and the content of copper in soils decreased to ~950, zinc to ~500, and manganese to ~ (Figure 12a) thousand tons. As a result of the reduction of metals reserves in the catchment area, their content in the river water decreased. Calculations have shown

that over a 400-year period, due to the self-cleaning of the NKR basin, the average annual concentrations of copper, zinc, and manganese in river water decreased slightly by 6–8%.

**Figure 12.** Changes in Mn reserves in soils (**a**), wash-off intensity, and Mn concentration in river waters (**b**) over a 400-year self-purification period.

#### *3.3. Sensitivity of the Water and Chemical Runoff to Climatic Changes*

Modeling the water and chemical runoff for the future period, up to 2050, was carried out using a DC climatic scenario, assuming changes in the climate parameters to be the same as the previous observation periods (increments of the daily air temperature of 1 ◦C and daily precipitation of 10%). Table 4 presents the results of numerical experiments for the simulation of river runoff characteristics under both historical data series and the DC scenario for the 2002–2007 period.

**Table 4.** Simulated runoff characteristics under historical data series (Hist.) and DC scenario (DC) for the 2002–2007 period.

The mean annual inflow into the NKR under the DC scenario increased by 7% with a range from 4 to 12% in different years compared to the runoff simulated under historical meteorological data. Increments of the maximal discharges of the spring runoff in different years are multidirectional and vary from −14 to 22%. Both increments of runoff and changes in maximal discharges are not related to the amount of annual runoff. They are determined by the conditions of runoff formation during spring in specific years, and in particular, by the factor of "friendly" spring. Let us analyze this in more detail on the example of flow hydrographs for 2002–2003 (Figure 13).

The winter and pre-spring period of 2002 was characterized by an unstable temperature regime with frequent thaws. The increase of 1 ◦C in the DC scenario caused earlier snowmelt and the formation of a much larger first wave of spring runoff than those simulated from historical meteorological data. The maximum discharge in the second wave of spring runoff decreased by 14% due to the fact that a significant part of the snow water equivalent was consumed during the first wave of spring runoff.

**Figure 13.** Simulated hydrographs of lateral inflow into the NKR under historical (blue) and DC scenarios (red) over the 2002–2003 period.

The winter and pre-spring season of 2003 was much more severe in terms of temperature conditions than in 2002. Therefore, the increase in temperature by 1 ◦C did not cause earlier snow melting. The increased maximal snow water equivalent due to both the addition of precipitation in the winter season and more intense snowmelt in the spring under the DC scenario caused a 22% increase in maximal discharge of spring runoff, with the peak 18 days earlier than under the historical scenario.

It is more difficult to generalize the patterns of the influence of climate change on the characteristics of water quality. On the whole, the mean annual Mn concentrations in the last downstream hydrochemical station on the Belaya River (Durtuli gauge) averaged over 7 years in the DC scenario increased by 5% compared to those modeled based on historical meteorological data (Table 5).

**Table 5.** Simulated water quality characteristics under historical data series and DC scenario for the 2002–2007 period.

The changes in annual Mn loads into the NKR in different years are proportional to increments in the mean annual runoff multiplied by the changes in annual concentrations at the last monitoring gauge on the Belaya River. The maximum Mn concentrations averaged over a 7-year period decreased slightly by 2%, and the range of changes in different years vary from 4 to −10%.

Let us compare these values with assessments of the climate change impact on the quantity and quality of water in the Elbe River basin performed using the process-oriented SWIM model [21]. Using climate change scenarios close to the one we employed, river runoff will increase by an average of 11% by the middle of the century. The predicted changes in the concentrations of nitrogen and phosphorus nutrients show multidirectional

changes. Uncertainties and high spatial variability of characteristics within the basin may be associated with the schematization of the Elbe basin into a number of separate tributary catchments and with the need to calibrate the model parameters for each of these sub-catchments. In the ECOMAG-HM model, the whole basin and spatial variability of characteristics are simulated with a single set of model parameters.

#### **4. Conclusions**

The semi-distributed, physically based ECOMAG-HM model developed for modeling the hydrochemical cycle of heavy metals at the catchment scale was calibrated and verified for the large watershed of the Nizhnekamskoe reservoir according to hydrometeorological, hydrochemical, and water management monitoring data. The model was performed by comparing the observed and simulated river runoff hydrographs and behavior of HM (Cu, Zn, and Mn) concentrations in monitoring points on rivers of the NKR basin. Comparing results showed that seasonal changes in the metal content are characterized by increased concentrations during the high-water period (snowmelt and rainfall runoff). Inter-annual changes in the observed and simulated HM load into the river network and NKR are characterized by its increase in high-water years and minimum values in low-water years. Such seasonal and inter-annual changes in the HM load occur when HMs mainly enter the river network from diffuse sources of pollution on the catchment area caused by climatic factors. Model balance calculations have shown that the mean annual wash-off of metals into the rivers from the NKR catchment is approximately 80–85% attributed to leaching from the soil-ground layer, and the surface wash-off of metals usually does not exceed 20%. The contribution of HMs entering the river network with wastewater discharges is small and does not exceed 4%. A significant part (from 30 to 62%) of metals washed into rivers accumulate with sediments at the river bottom.

The simulated map of the total HM wash-off into the river network closely correlates with the spatial distribution of the initial content of HM in soils due to the predominance of HM leaching from the soil-ground zone in the diffuse sources of river pollution. The initial HM content in soils affects the distribution of HM concentrations in the river network to a large degree. As a result, the most polluted small rivers flow near areas with a high content of HM in the soils. This leads to the important consequences: Without maps of the initial content of heavy metals in soils, it is impossible to simulate spatial differences in water pollution in the river network.

Numerical experiments have shown that with an increase in discharge of the Belaya River more than a certain critical value, the contribution of point sources of pollution is minimized, and the concentration of metals in river water is mainly determined by the maximum exchange capacity of the catchment leaching HM from the watershed by surface and subsurface flows during intensive snowmelt or rainfall that forms the maximum flow; for the city of Ufa, the critical discharge is about 3500 m3·s−1, above which the copper concentrations in river water asymptotically approach the value of ~6 μg·L−1 under various scenarios of wastewater loads. This constant concentration can be called the equilibrium concentration of diffuse saturation of the river basin for this type of HM.

The simulation of time scales of the catchment's self-purification from HMs under the scenario with complete exclusion of anthropogenic impacts on the basin has shown that over a 400-year period, the mean annual HM concentrations in the river water decreased by 6–8%. The low rate of reduction of heavy metal concentrations allows us to recommend their current annual values as background values in the NKR basin.

Modeling the river runoff for the future period up to 2050 under a DC climatic scenario showed that the mean annual inflow into the NKR increases by 7% with a range from 4 to 12% in different years. No significant changes in the characteristics of water quality for heavy metals should be expected in the near future. For example, averaged over 7 years, the mean annual Mn concentration increased by 5%, the maximum Mn concentration decreased slightly by 2%, and the annual Mn load into the NKR increased by about 10% compared to those modeled based on historical meteorological data.

The presented results show that the ECOMAG-HM model can be used as an information support tool for assessing water quality characteristics under possible climatic and anthropogenic changes in scientific and applied projects.

**Author Contributions:** Methodology, writing—review and editing, supervision, funding acquisition, Y.M.; data curation, formal analysis, visualization, project administration, T.F.; conceptualization, validation, investigation, resources, writing—original draft preparation, Y.M. and T.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was carried out under the Governmental Order to the Water Problems Institute, Russian Academy of Sciences, subject no. 0147-2019-0001.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**

## *Article* **Natural and Anthropogenic Controls of Groundwater Quality in Sri Lanka: Implications for Chronic Kidney Disease of Unknown Etiology (CKDu)**

**Sen Xu 1, Si-Liang Li 1,2,\*, Fujun Yue 1, Charitha Udeshani 1,3 and Rohana Chandrajith 3,4**

**Abstract:** Poor groundwater quality in household wells is hypothesized as being a potential contributor to chronic kidney disease of unknown etiology (CKDu) in Sri Lanka. However, the influencing factors of groundwater quality in Sri Lanka are rarely investigated at a national scale. Here, the spatial characteristics of groundwater geochemistry in Sri Lanka were described. The relationships of groundwater quality parameters with environmental factors, including lithology, land use, and climatic conditions, were further examined to identify the natural and anthropogenic controlling factors of groundwater quality in Sri Lanka. The results showed that groundwater geochemistry in Sri Lanka exhibited significant spatial heterogeneity. The high concentrations of NO3 − were found in the districts that have a higher percentage of agricultural lands, especially in the regions in the coastal zone. Higher hardness and fluoride in groundwater were mainly observed in the dry zone. The concentrations of trace elements such as Cd, Pb, Cu, and Cr of all the samples were lower than the World Health Organization guideline values, while some the samples had higher As and Al concentrations above the guideline values. Principal component analysis identified four components that explained 73.2% of the total data variance, and the first component with high loadings of NO3 −, hardness, As, and Cr suggested the effects of agricultural activities, while other components were primarily attributed to natural sources and processes. Further analyses found that water hardness, fluoride and As concentration had positive correlations with precipitation and negative correlations with air temperature. The concentration of NO3 − and water hardness were positively correlated with agricultural lands, while As concentration was positively correlated with unconsolidated sediments. The environmental factors can account for 58% of the spatial variation in the overall groundwater geochemistry indicated by the results of redundancy analysis. The groundwater quality data in this study cannot identify whether groundwater quality is related to the occurrence of CKDu. However, these findings identify the coupled controls of lithology, land use, and climate on groundwater quality in Sri Lanka. Future research should be effectively designed to clarify the synergistic effect of different chemical constituents on CKDu.

**Keywords:** groundwater; water quality; heavy metals; CKDu; dry zone; fluorosis; arsenic

#### **1. Introduction**

Chronic kidney disease (CKD) is a very common public health problem that can be observed in many parts of the world and a higher prevalence is reported from many countries including the USA, Australia, and Japan [1,2]. A CKD of unknown etiology (CKDu) has

**Citation:** Xu, S.; Li, S.-L.; Yue, F.; Udeshani, C.; Chandrajith, R. Natural and Anthropogenic Controls of Groundwater Quality in Sri Lanka: Implications for Chronic Kidney Disease of Unknown Etiology (CKDu). *Water* **2021**, *13*, 2724. https://doi.org/10.3390/w13192724

Academic Editor: Fernando António Leal Pacheco

Received: 27 August 2021 Accepted: 30 September 2021 Published: 1 October 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/).

been reported in some countries in recent years in which the etiology is not recognizable [3]. Such a disease with unknown etiology has been reported in certain parts of the world, especially in Africa, Central America, and Asia [4,5]. CKDu is also found in the rural dry zone regions of Sri Lanka, particularly in the North Central Province adjacent regions for more than two decades [6,7]. Due to its remarkable geographical distribution in some specific regions where groundwater is the main drinking water source, it is generally suspected that long-term exposure to various nephrotoxic elements through drinking groundwater is an important risk factor [7–9]. In the rural dry zone regions of Sri Lanka, groundwater is the primary source of water for both drinking and cooking, and CKDu is considered to be related to poor groundwater quality [10–12].

Groundwater is an indispensable limited resource, playing an important role in the social and economic development throughout the world. High temperature and limited precipitation in the rural dry zone regions of Sri Lanka pose challenges in the supply of adequate water to meet daily needs [10]. Groundwater provides domestic water and also contributes irrigation water supply for rural communities in Sri Lanka, since no or less treatment is often required [10,13]. However, the deterioration of groundwater quality and the depletion of water resources due to population growth and excessive use of pesticides and agrochemicals, pose threats to public health [13]. Previous studies found that the prevalence of CKDu is related to the recharge sources of groundwater, and is significantly higher with the groundwater stagnated as well as groundwater recharged from regional flow paths [8,14]. Therefore, it is suggested that the source, recharge mechanism, and flow pattern of groundwater, as well as geological conditions that would cause natural contamination of groundwater, may be the main causative factors for CKDu [14,15].

Anthropogenic activities such as agricultural and sanitation practices, and environmental factors such as lithology, land use, and climatic conditions impact on the groundwater quality in Sri Lanka [16–18]. Proposed factors related to etiology of CKDu in groundwater include but not be limited to fluoride [8,19–24], hardness [8,15,22], major ions [12,20], heavy metals and metalloids [18,25,26], and agrochemical residues [18,27]. Despite groundwater quality related to the occurrence of CKDu in Sri Lanka has been studied extensively, few studies have been carried out to analyze the natural and anthropogenic factors that control the groundwater geochemistry in Sri Lanka. Therefore, the objective of this paper is to investigate the spatial characteristics of groundwater geochemistry, and its relationships with environmental factors, and also to provide a comprehensive review regarding the relationships between groundwater quality and the occurrence of CKDu in the dry zone of Sri Lanka.

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

#### *2.1. Study Area*

Sri Lanka (6–10◦ N, 79–82◦ E) is an island located in the Indian Ocean, just 800 km north of the equator, close to southeastern coast of the Indian subcontinent, and includes 25 administrative districts with a total area of 65,610 km2. The length of the island is 440 km from north to south with a maximum width of 226 km. Although the area of the island is relatively small, it has a widely changing topography, geology, and climate within a short distance imparting some unique environmental features (Figure 1) [28]. The topographical configuration of the island is a highland region located in the center, surrounded by a vast lowland plain. Over 90% of the island is underlain by Precambrian rocks. According to a global lithological map [29], metamorphic rocks cover about 68.2% of the island, followed by acid plutonic rocks (14.4%) and unconsolidated sediments (12%). As per the year 2017, the percentages of forest, cropland, and impervious surface of the island were 60.4%, 30.9%, and 0.6%, respectively. The island is situated in the tropical, Indian ocean monsoonal climate region that contains three climate zones based on the amount of precipitation, named as the wet, dry and intermediate zones [13,16]. The western slopes of the central highlands have the highest intensity of precipitation that exceeds 4000 mm per annum according to the records of the average for years 1970–2000 [30]. The lowest precipitation

is recorded in the Mannar Island with values about 900 mm per annum. CKDu is found particularly in the dry zone, metamorphic terrain that occupies two-thirds of the country with the annual rainfall of about 1250 mm, while air temperature varies from 28 to 30 ◦C, with humidity of around 70% [16].

**Figure 1.** Maps showing (**a**) land use, (**b**) lithology, (**c**) elevation, (**d**) annual precipitation, (**e**) annual mean temperature, and (**f**) CKD prevalence rates across the most affected districts in Sri Lanka [10].

#### *2.2. Data Source*

Groundwater quality data including the concentrations of NO3 −, hardness, fluoride, heavy metals and metalloids (As, Cd, Fe, Mn, Cu, Cr, Pb, and Al) of groundwater samples collected from the wells in Sri Lanka were obtained from Groundwater Quality Atlas of Sri Lanka [31]. A total of 1304 samples were collected and analyzed from 2010 to 2014 for their study. The maximum and average concentrations of the water quality parameters of all the samples in each district were reported by Kawakami et al. (2014) [31], and used and reanalyzed in this study (Table S1).

Lithology, land use, and climatic conditions were quantitatively described for each district to examine the effects of environmental factors on groundwater quality (Table S2). The elevation data was obtained from EarthEnv-DEM90, that is a digital elevation model derived from CGIAR-CSI SRTM v4.1 and ASTER GDEM v2 data products [32]. To quantify the lithological distribution, a global lithological map was used and analyzed for Sri Lanka [29]. Land use dataset with a 10 m resolution was also used to obtain information for this study [33]. Land use types include cropland, forest, grassland, shrubland, wetland, water body, impervious surface, and bare land, among which cropland and impervious surface were mainly focused. Meteorological data including air temperature and precipitation

data used in this study were obtained from WorldClim version 2.1 at a spatial resolution of 30 arc-seconds and are the averages for years 1970–2000 [30].

#### *2.3. Statistical Methods*

Correlation Analysis (CA) and Principal Component Analysis (PCA) were used to examine the relationships between the groundwater quality parameters for exploring their possible sources. The suitability of data for PCA was checked by Kaiser–Meyer–Olkin (KMO) and Bartlett's sphericity tests. The average concentrations of groundwater quality parameters were used in the CA and PCA. Redundancy analysis (RDA) was selected to assess the explanation of variations of groundwater quality parameters by environmental factors. In the ordination diagram resulted from RDA, the arrows pointing in the same direction represent positive correlations, or vice versa. The data shown in Figure 2 were used in the RDA. The CA and PCA were performed through the statistical software package SPSS 22.0, and the RDA was carried out using CANOCO 5.0 software (Microcomputer Power Company, Ithaca, NY, USA).

**Figure 2.** Maps showing the spatial patterns of average concentrations of groundwater quality parameters in each district (maximum concentration (max) was showed for fluoride, As, and Cd). The ranges of data values were categorized into five classes based on the natural break method.

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

*3.1. Spatial Characteristics of Groundwater Geochemistry*

3.1.1. Nitrate, Hardness, and Fluoride

Groundwater quality parameters varied considerably in different regions (Figure 2). The average concentration of NO3 − in the district ranged from 1 mg/L to 42 mg/L. The highest concentration (366 mg/L) and highest average concentration (42 ± 90 mg/L) of

NO3 − were observed in the Puttalam district, where 49.7% of the land area is used for agriculture. The World Health Organization (WHO) guideline for NO3 − in drinking water is 50 mg/L [34]. The maximum NO3 − concentration in 13 out of 25 districts exceeded 50 mg/L limit. Higher concentrations of NO3 − were also found in other districts where a high percentage of land is used for agricultural purposes, for instance, the Jaffna (62.9%), Mannar (35.3%), and Vavuniya (28%) districts. Agricultural regions in the coastal zone are mainly underlain by recent unconsolidated sediments in which farmers excessively use a large amount of nitrate fertilizers for their cultivations, leading to nitrate pollution of groundwater [35]. Manure may be a cause for the high NO3 − concentration found in the Vavuniya district [36].

The average hardness ranged from 20 to 354 mg/L. The highest hardness (1734 mg/L) was found in the Hambantota district, while the highest average hardness (354 ± 190 mg/L) was noted in the Jaffna district, where the percentage of carbonate sedimentary rocks (37.1%) is the highest. The Puttalam and Vavuniya districts also showed a higher hardness. Since water hardness poses no apparent health concerns, there are no guidelines or regulations for optimum hardness in drinking water [34]. Considering the impact of hardness on water palatability, a non-health-based standard of 250 mg/L (as CaCO3) was defined by the Sri Lankan Standard for drinking water [37]. The geographic distribution of CKDu cases in the endemic region showed a strong association with the consumption of hard water [27].

The average concentration of fluoride ranged from <0.02 mg/L to 1.4 mg/L in Sri Lanka. The highest average concentration of fluoride (1.4 ± 1.3 mg/L) was recorded in the Moneragala district. High concentrations of fluoride (max: 7 mg/L; ave: 1.1 ± 0.9 mg/L) were also found in the Anuradhapura district. Spatially, fluoride showed a high concentration in the dry zone and a low concentration in the wet zone and also in the northern coastal zone. The spatial distribution of fluoride is closely resembled with that of the hardness. Although the maximum guideline value for fluoride in drinking groundwater as suggested by the WHO is 1.5 mg/L [34], and the limit suggested by the Sri Lankan Standard is 1.0 mg/L [37], both values cannot be applicable to the dry zone regions of Sri Lanka due to a higher water consumption for drinking under prevailing hot and dry ambient conditions [19]. The average concentrations of fluoride recorded in all the districts were lower than 1.5 mg/L, and only 9.9% of the samples had a concentration higher than 1.5 mg/L [36]. The fluoride concentration in Sri Lanka is much lower than that in many other regions in the world, however dental and skeletal fluorosis cases are common in these areas [10,16,19].

#### 3.1.2. Heavy Metals and Metalloids

In some studies of potential nephrotoxic effects of environmental exposure to heavy metals and metalloids, mainly As, Cd, and Pb have been identified as causative factors for CKDu [18,25,27,38,39], although some studies reject this hypothesis [7,8,15,20,40]. These metals in groundwater may be derived from rock weathering or by human activities, such as the use of agricultural chemicals [18]. Arsenic in groundwater has been proposed as a causal factor for CKDu in Sri Lanka [18,27]. The WHO guideline value and Sri Lankan standard limit for As are 10 μg/L [34,37]. The highest concentration (66 μg/L) and highest average concentration (7 ± 11.7 μg/L) of As were recorded in the Mannar district. High As concentrations were also found in the Batticaloa, Mullaitivu, and Puttalam districts that had the maximum As concentration higher than 10 μg/L. Some recent detail studies also noted higher As concentrations in groundwater, particularly in groundwater extracted from sedimentary aquifer systems [35,41].

Cd is known to be a toxic metal and has been listed as one of the causes of CKDu due to the well-known nephrotoxicity [12,25]. The highest Cd concentration in this study was 0.5 μg/L and recorded in the Moneragala district. All of the well water samples had a Cd concentration much lower than the WHO guideline value of 3 μg/L [34]. The highest Pb concentration (288 μg/L) was noted in the southern district of Galle, however, it may be an outlier. If this data was eliminated, the highest Pb concentration of 4.6 μg/L was recorded in the Anuradhapura district, but it is lower than the WHO guideline value of 10 μg/L [34].

Fe is one of the most abundant metals in the Earth's crust and is an essential human micronutrient [34]. Mn usually coexists with Fe and is also one of the most abundant metals in the Earth's crust. Cu is an essential nutrient and also a contaminant of drinking water. The highest concentrations and highest average concentrations of Fe (max: 828 μg/L; ave: 283 μg/L), Mn (max: 9772 μg/L; ave: 440 μg/L), and Cu (max: 443 μg/L; ave: 21 μg/L) were recorded in the Nuwara Eliya, Polonnaruwa, and Mullaitivu districts, respectively. The WHO guideline value for Cu in drinking water is 2 mg/L, but no WHO guideline values for Fe and Mn in drinking water are proposed [34]. The highest Cr concentration (14 μg/L) was found in the Ampara district, yet the level is lower than the WHO guideline value for Cr in drinking water of 50 μg/L [34]. The highest concentration (1457 μg/L) and highest average concentration (58 μg/L) of Al were found in the Moneragala and Nuwara Eliya districts, respectively. The maximum guideline value for Al in drinking water suggested by the WHO is 200 μg/L [34], and 2.3% of the samples exceeded the guideline value [36].

#### *3.2. Natural and Anthropogenic Controls on Groundwater Geochemistry*

The correlation matrix demonstrated that there were significant positive correlations among NO3 −, As, and Cr (*p* < 0.05) (Figure 3), indicating the higher contamination of As and Cr of groundwater in agricultural areas. Water hardness was positively correlated with fluoride (R = 0.57, *p* < 0.01), which suggests that they may be affected by similar geogenic factors and processes, such as the weathering of fluoride-bearing minerals under the hot climatic condition and excessive evaporation in the dry zone [19,42]. Cd and Cu showed a significant positive correlation (R = 0.54, *p* < 0.01), implying that they may be derived from the same source.

**Figure 3.** Correlation coefficients for Pearson correlation analysis among groundwater quality parameters. \* Correlation is significant at *p* < 0.05; \*\* Correlation is significant at *p* < 0.01.

PCA was employed to identify the associations between groundwater quality parameters and their potential sources [43]. Four principal components (PCs) with eigenvalues exceeding 1 were identified, and the PCA results are shown in Table 1. The first PC was dominated by NO3 −, water hardness, As, and Cr, with high loadings of 0.86, 0.67, 0.7, and 0.71, respectively, explaining 26.5% of the total data variance. The second PC explained 24.6% of the total data variance and was mainly contributed by Fe, Pb, and Al, with higher loadings of 0.81, 0.77 and 0.78, respectively. The third PC was mainly contributed by Cd and Cu, with high loadings of 0.88 and 0.76, respectively, explaining 12.6% of the total data variance. The fourth PC was mainly contributed by fluoride and Mn, with high loadings of 0.75 and 0.72, respectively, and explained 9.5% of the total data variance. The four PCs totally accounted for 73.2% of the total data variance. Agrochemicals especially phosphate fertilizer used in agriculture have been identified as the main source of As in areas affected with CKDu in Sri Lanka [18]. The groundwater in the dry zone has high fluoride concentrations due to the water–rock interaction with the aquifer [15,19]. In addition, the occurrence of CKDu is strongly associated with the consumption of hard water in the rural dry zone [22,27]. Therefore, the first PC may be related to the agricultural activities in the rural dry zone, suggesting the contributions of anthropogenic sources to NO3 −, As, and Cr. In contrast, other PCs may be dominated by natural geogenic sources and processes, especially bedrock weathering.

**Table 1.** Component matrix and explained variance (principal component analysis (PCA)) for groundwater quality parameters.

The groundwater quality in Sri Lanka is mainly influenced by environmental factors such as lithology, land use, and climatic condition, and anthropogenic activities including agricultural practices, disposal of domestic sewage, and industrial effluents [16–18]. To explore the potential influencing factors of groundwater quality in Sri Lanka, the relationships between the environmental factors and the groundwater quality parameters were examined. As shown in Figure 4, water hardness showed a negative correlation with the precipitation (for average hardness: R2 = 0.78, *p* < 0.01; for maximum hardness: R2 = 0.53, *p* < 0.01), following a power law relationship. Similar relationships were also found between fluoride and precipitation, and between As and precipitation. This suggests that in the dry zone with low precipitation, the concentrations of fluoride and As in groundwater were higher than those in the wet zone where excessive precipitation is characterized. Since lower precipitation is generally accompanied by high ambient temperature in Sri Lanka, the higher concentrations of hardness, fluoride, and As were also found in areas with high air temperature.

**Figure 4.** The relationships between hardness, fluoride, and As concentrations, and precipitation (**a**–**c**) and air temperature (**d**–**f**).

The positive correlations between NO3 − and hardness and cropland suggest the effects of agricultural activities on groundwater quality and that most high nitrate regions are underlain by sedimentary limestone sequences (Figure 5). In Sri Lanka, all wells with a high As concentration were located on a specific soil type of "sandy regosols on recent beach and dune sands" [36]. The positive correlation between As concentration and unconsolidated sediments and the spatial characteristics of As concentration indicate that As had a geological source. Amarathunga, et al. (2019) [44] noted that As is released from coatings of sand grains due to reductive dissolution under near natural pH condition in the Mannar region where a majority of wells exceeded the WHO recommended level of As.

**Figure 5.** The relationships between NO3 − concentration (**a**) and hardness (**b**) and cropland, and between As concentration and unconsolidated sediments (**c**).

RDA was used to assess the explanation of variations of the groundwater quality parameters by the environmental factors. The RDA results provided overall descriptions of the influences of environmental factors on groundwater geochemistry [45]. The results showed that the environmental factors can account for 58% of the spatial variation in the overall groundwater geochemistry (pseudo-F = 2.5; *p* = 0.002). The first two axes explained most of the spatial variation of groundwater geochemistry, with the first axis explaining 26% of the variation and the second axis explaining 16% of the variation. Ordination diagrams resulted from RDA showed the relationships between the groundwater quality parameters and the potential influencing factors (Figure 6). The results indicated that hardness was positively correlated with agricultural lands while fluoride concentration was positively correlated with ambient temperature. By contrast, precipitation and forest were negatively correlated with most of the groundwater quality parameters. These results

indicate the synergetic controls of lithology, land use, and climate on the groundwater quality in Sri Lanka.

**Figure 6.** Ordination diagrams showing the relationships between the groundwater quality parameters (represented by blue lines) and the potential influencing factors (represented by red lines) according to the redundancy analysis (RDA). The analyzed potential influencing factors include elevation (Ele), cropland (Cro), forest (For), impervious surface (IS), metamorphic rocks (MR), acid plutonic rocks (APR), unconsolidated sediments (US), carbonate sedimentary rocks (CSR), precipitation (Pre), and air temperature (AT).

#### *3.3. Relationships between Groundwater Quality and Occurrence of CKDu*

In the past three decades, possible etiology of CKDu has been extensively explored through multi-disciplinary studies. The CKDu distribution related to geography and socio-economy showed that the possible etiology is related to geogenic environment and occupational factors [6,9,39,46,47]. The long-term exposure to various nephrotoxic elements through drinking groundwater is widely suspected due to the geographical distribution of CKDu [7,8,15,16,20,21,28,40]. CKDu in Sri Lanka was first reported from the north central province of the rural dry zone, and the main victims of CKDu are young male farmers with substandard socio-economic background [6,10,39]. In the worst affected areas of CKDu, more than 98% of the population had completely relied on groundwater as their primary source for drinking or cooking at least five consecutive years between 1999 and 2018 [10]. This indicates that groundwater quality may interfere with human health in CKDu endemic region.

#### 3.3.1. Relationship between Fluoride and CKDu

Fluoride has been proposed as a causal factor for CKDu through drinking water. It is suggested that the source of F is bedrock weathering [19,21,22]. Climatic and hydrological conditions were found to be related to the anomalous concentrations of fluoride in the groundwater in Sri Lanka [19]. Fluoride may be an essential element for humans, for instance, fluoride is a beneficial element in the prevention of dental caries [34]. However, elevated fluoride intakes can have more detrimental effects such as dental and skeletal tissues [34]. Previous studies have revealed the dose–effect relationships between drinking water fluoride levels and damage to kidney functions [11,24]. High fluoride concentrations in drinking water elevated the levels of renal and liver function enzymes in serum and caused severe histological changes of the liver and kidneys [48]. Due to its high rank in the Hofmeister Series for denaturing proteins of the kidney membrane, fluoride may contribute to CKDu [20]. In this study, high fluoride concentrations were found in the dry zone where CKDu has its higher prevalence (Figure 1). Interestingly, no CKDu was observed in the areas in the southeast where the fluoride concentration in groundwater was greater than 1.0 mg/L [31]. This suggests that the high concentration of fluoride in groundwater is not only the factor affecting the occurrence of CKDu. The interaction of fluoride with aluminum to form nephrotoxic aluminum fluoride complexes has been proposed as a cause of CKDu [26]. In addition, when fluoride interacts with major cations

(Ca2+, Mg2+, and Na+) and metals such as cadmium, it is suspected that it will aggravate kidney failure [12,21,22,28,49]. Therefore, the synergistic effect between fluoride and major ions, and metals on kidney functions should be paid more attention.

#### 3.3.2. Water Hardness and Major Ions

Hardness in water is caused by a variety of dissolved polyvalent metallic ions, predominantly the cations of Ca and Mg. As classification of hardness, the ranges of 0–60 mg/L, 61–120 mg/L, 121–180 mg/L, and >181 mg/L are classified as soft, moderately hard, hard, and very hard, respectively [34]. Previous studies revealed the highly statistically positive correlation between the occurrence of CKDu and the hard water consumption in Sri Lanka, and found that 96% of CKDu patients consumed hard or very hard water from wells in shallow regolith aquifers for at least 5 years [12,27]. However, many studies found that hard to very hard groundwater is not exclusive to all CKDu endemic areas [14,16]. For instance, only 14% of the samples from Mahiyanganaya, a high CKDu endemic dry zone area of Sri Lanka, had hard to very hard groundwater [16]. Moreover, the groundwaters from some areas of Sri Lanka that are not seriously affected by CKDu had very hard waters, such as the districts of Puttalam (329 ± 174 mg/L) and Jaffna (354 ± 190 mg/L). Nevertheless, hardness might interfere with other chemical constituents in water to form various complexes harmful to human health [8,15].

#### 3.3.3. Nephrotoxic Heavy Metals and Metalloids

Heavy metals and metalloids such as As and Cd have been considered as important risk factors for CKDu [18,25,38]. The application of pesticides and fertilizers in rice paddy cultivation in Sri Lanka is a possible source of high Cd [25]. In this study, the Cd concentrations of all the well water samples were significantly below the WHO guideline value. Similar results were also reported by many other studies [16,40,46]. Moreover, previous studies noted that the levels of Cd and As of drinking water and rice samples collected from the affected areas of CKDu were within the levels recommended by WHO and Sri Lankan standard and did not indicate contamination in any form [40]. A number of studies also found that in the CKDu endemic areas of Sri Lanka, the As concentration in groundwater is significantly low [8,15,16,46]. Therefore, it is suggested that As and Cd in groundwater are unlikely to be risk factors for CKDu in Sri Lanka. In addition, Jayasumana et al. (2014) [27] have demonstrated the link between hardness and As. They proposed that As mainly derived from chemical fertilizers and pesticides can eventually damage kidney tissue when combined with Ca and/or Mg in groundwater. In addition, there is considerable evidence that agricultural workers in the CKDu endemic areas are exposed to As, but the exact source and entry mode of As are still controversial [27].

#### **4. Conclusions**

This study investigated the spatial characteristics of groundwater geochemistry in Sri Lanka. Groundwater quality data of 1304 water samples collected from 2010 to 2014 were statistically analyzed. Land use and lithology were quantitatively described to examine their relationships with the groundwater quality parameters. The data of ambient temperature and precipitation were used to explore the effects of climatic conditions on groundwater quality. The results showed that the high concentrations of NO3 − were found in the districts with a high percentage of cropland, especially in the districts in the coastal zone. The Puttalam district had the highest concentration and highest average concentration of NO3 −. The samples from the dry zone had a higher hardness, and in the Jaffna district, where the percentage (37.1%) of carbonate sedimentary rocks is the highest, the average hardness (354 ± 190 mg/L) of samples was the highest. Similar to the spatial distribution of hardness, fluoride also showed a high concentration in the dry zone, and 9.9% of the samples had a concentration higher than the WHO guideline value of 1.5 mg/L. Heavy metals and metalloids such as As, Cd and Pb have been suggested as causal factors for CKDu in Sri Lanka. The maximum As concentrations of the samples in the Mannar, Batticaloa, Mullaitivu, and Puttalam districts were higher than the WHO guideline and Sri Lankan standard value (10 μg/L). The Cd and Pb concentrations of all the samples were lower than the WHO guideline values of 3 μg/L and 10 μg/L, respectively. In addition, all of the samples had Cu and Cr concentrations below the WHO guideline values, and 2.3% of the samples had a concentration of Al above the WHO guideline value (200 μg/L).

Significant positive correlations (*p* < 0.05) among NO3 −, As, and Cr were observed, suggesting the effects of agricultural activities on As and Cr concentrations. The significant positive correlation between water hardness and fluoride indicated that they might be affected by rock weathering in the dry zone. According to the results of PCA, four PCs were identified and explained 73.2% of the total data variance. Considering the concentrations and spatial characteristics of groundwater quality parameters, PC1 with high loadings of NO3 −, hardness, As, and Cr was considered to be affected by agricultural activities, and other PCs was primarily attributed to natural sources and processes. Hardness, fluoride, and As concentration were positively correlated with precipitation and negatively correlated with air temperature, showing a clear difference in groundwater quality between different climatic regions, and suggesting that groundwater geochemistry in Sri Lanka is closely related to climatic conditions. In addition, NO3 − concentration and water hardness had positive correlations with cropland, and As concentration had a positive correlation with unconsolidated sediments, implying the effects of land use and lithology on groundwater geochemistry. The RDA results showed that 58% of the spatial variation in the overall groundwater geochemistry can be accounted by the environmental factors.

To our knowledge, there are few previous studies that quantitatively describe land use and lithology, and identify the natural and anthropogenic controlling factors of groundwater quality in Sri Lanka, especially at a national scale. This study provided key insights of the effects of environmental factors on groundwater quality in Sri Lanka, which considerably helps to investigating the relationships between groundwater quality and the occurrence of CKDu. In the future, high spatial-temporal resolution sampling is needed to unravel controlling factors of groundwater quality, and multidisciplinary studies are required to identify risk factors of CKDu in Sri Lanka.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/w13192724/s1, Table S1: Statistics of the groundwater quality parameters for each district in Sri Lanka, Table S2: Characteristics of the 25 administrative districts in Sri Lanka.

**Author Contributions:** Conceptualization, S.X. and S.-L.L.; methodology, S.X. and F.Y.; resources, S.-L.L. and F.Y.; data curation, S.X. and S.-L.L.; writing—original draft preparation, S.X., S.-L.L., F.Y., C.U. and R.C.; writing—review and editing, S.X., S.-L.L., F.Y., C.U. and R.C.; visualization, S.X. and S.-L.L.; supervision, C.U. and R.C.; project administration, S.-L.L. and R.C.; funding acquisition, S.-L.L. and R.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was carried out as collaborative research, funded by the National Natural Science Foundation, China (Grant No.<PHONE_NUMBER>6) and National Science Foundation, Sri Lanka (Grant No. ICRP/NSF-NSFC/2019/BS/02).

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

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

#### **References**

## *Article* **Effects of Seasonal Thermal Stratification on Nitrogen Transformation and Diffusion at the Sediment-Water Interface in a Deep Canyon Artificial Reservoir of Wujiang River Basin**

**Yongmei Hou 1,2, Xiaolong Liu 1,\*, Sainan Chen 3, Jie Ren 1,2, Li Bai 1, Jun Li 1, Yongbo Gu 1,2 and Lai Wei 1,2**

**Citation:** Hou, Y.; Liu, X.; Chen, S.; Ren, J.; Bai, L.; Li, J.; Gu, Y.; Wei, L. Effects of Seasonal Thermal Stratification on Nitrogen Transformation and Diffusion at the Sediment-Water Interface in a Deep Canyon Artificial Reservoir of Wujiang River Basin. *Water* **2021**, *13*, 3194. https://doi.org/10.3390/ w13223194

Academic Editor: Bommanna Krishnappan

Received: 20 September 2021 Accepted: 9 November 2021 Published: 11 November 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/).

**Abstract:** Watershed-scale nitrogen pollution in aquatic systems has become a worldwide concern due to its continuous impact on water quality deterioration, while the knowledge of key influencing factors dominating nitrogen transportation and transformation at the sediment-water interface (SWI) remains limited, especially in impounded rivers with an artificial reservoir. Hence, for a better understanding of the effects of thermal stratification on nitrogen transformation, we investigated the nitrogen species and isotopes in the sediment of a deep reservoir in Southwest China. Our results confirmed a significant difference in nitrogen species and isotopic composition in sediment between those in the thermal stratification period and non-thermal stratification period and indicated that the sediment biogeochemical process and transportation were clearly linked to the variations in water temperature and dissolved oxygen dominated by the process of thermal stratification. Significant seasonal differences in NH4 +-N and NO3 −-N in pore water of the upper layer (0–19 cm) revealed that nitrification exhausted NH4 + in the non-stratified period (NSP), and a potential low mineralization rate appeared when compared with those in the stratified period (SP). Seasonal differences in nitrogen species and isotope fractionation of δ15N-PON (about 2.3‰ in SP) in the upper layer sediment indicated a higher anaerobic mineralization rate of organic matter in SP than that in NSP. The diffusion fluxes of NH4 +-N at SWI were 9.48 and 15.66 mg·m−2·d−1 in NSP and SP, respectively, and annual NH4 +-N diffusion accounted for 21.8% of total storage in the reservoir. This study demonstrated that the nitrogen cycling processes, especially nitrification, denitrification, and mineralization, have been largely altered along with the changes in dissolved oxygen and that the diffusion of nitrogen species varied with the presence of the oxygen. The results contribute to the future study of watershed nitrogen budget evaluation and suggest that the endogenous nitrogen released from the sediment-water interface should be emphasized when aiming to fulfil water management policies in deep reservoirs.

**Keywords:** thermal stratification; nitrogen transformation; sediment-water interface; reservoir; NH4 + diffusion

#### **1. Introduction**

Excessive nitrogen loading at watershed-scale in aquatic systems, with the constant input of anthropogenically sourced active nitrogen, has posed serious environmental problems, such as eutrophication, greenhouse gas emission, and water quality deterioration [1,2]. Watershed-scale nitrogen pollution persistently attracts worldwide concern attributed to the enormous challenge of nitrogen assessment and management in global ecosystems. By intercepting runoff and sediment, massive dams and reservoirs worldwide would deposit the dissolved and adsorbed nutrients with the sediment and conserved

within the reservoir, which tends to reduce levels of nutrient transportation to downstream rivers and instead accumulate the nutrients in the reservoir [3,4]. Moreover, recent studies have emphasized that the accumulation of nitrogen in aquatic systems was affected by external nitrogen input as well as internal nitrogen release, which is dominated by nitrogen biogeochemical processes [5,6]. Evaluating the transformation and transportation of nitrogen species among the sediment-water interface (SWI) is essential if we aim to reduce nitrogen pollution in a reservoir and help to establish relevant policies for reservoir water quality protection. However, few studies have assessed the contribution of nitrogen releasing from the SWI in the reservoir. Currently, an increasing number of studies have raised the controversies surrounding the role of sediment in nitrogen transportation and transformation at the SWI, as well as its influence on the fate and distribution of nitrogen afterwards, which has become a hot topic in academic community [7,8].

Nitrogen experiences complex biogeochemical processes under different environmental conditions, e.g., DO has been recognized as the crucial factor that dominates nitrification and denitrification, and dissimilatory nitrate reduction to anammox (DNRA) generally occurs in sediment with a high C/N ratio and C loading [9]. Generally, ammonification, denitrification, and DNRA were regarded as being the dominant processes in a hypoxic environment in the stratified period (SP), while ammonification and nitrification could be favored by oxygenated conditions in the non-stratified period (NSP). In recent decades, intensive demands for hydropower resources led to the emergence of special landscape features caused by the cascading development in most rivers [10,11], which significantly posed a challenge for the synchronously high-efficiency management of water resources and the water environment. Over the past 60 years, approximately 58,000 large reservoirs have been built worldwide [12], which altered the original ecosystem conditions and attracted widespread concern [13].

It was well known that there was seasonal thermal stratification in deep canyon reservoirs, and water temperatures in the hypolimnion of deep lakes/reservoirs were generally low (4–6 ◦C) and remained nearly constant [13,14]. However, temperaturedriven hypoxia in the stratified period (SP) and oxygenated hypolimnion in the nonstratified period (NSP) may dominate the nitrogen exchange between sediment and water. Muller et al. [15] indicated that, in a lake, anaerobic environmental conditions were more conducive to internal nitrogen release than those in aerobic conditions. In addition, some studies reported that T and DO were the key factors controlling endogenous nitrogen diffusion of sediment in a shallow lake and that the diffusion fluxes of NH4 + were positively correlated with T, while there was a negative correlation with DO [16]. Cai et al. [17] investigated the relationship between nitrogen species and biological activities in shadow lake sediment, and reported that T and DO significantly alter the abundance of bacterial genes and community structures, subsequently affecting the transformation processes of sediment nitrogen. Hence, in this study, we hypothesized that the seasonal variation in DO and T would dominate the nitrogen transformation at SWI in a deep reservoir. However, as noted above, most current studies focused on nitrogen diffusion and influencing factors in thermal stratified lakes, while few studies have conducted targeted research in deep canyon reservoirs, especially the deep reservoirs in Southwest China. Because an artificial reservoir has a special hydrological regulation and water management regime for the purpose of power generation, the particulate matter sources and accumulation rates into the sediment are different from those in lakes. In order to optimize water protection and reduce nitrogen accumulation in reservoirs, sediment nitrogen transformation and diffusion to overlying water—which tends to play a role in endogenous pollution—urgently needed to be understood.

Water-soluble inorganic nitrogen and adsorbed inorganic nitrogen were the dominant nitrogen forms at the SWI; a growing number of studies have documented that the water-soluble and adsorbed nitrogen varied significantly along with the water chemistry changes [18]. Generally, there are two methods that have been employed universally to study nitrogen cycling at the SWI: incubation of intact sediment core experiments and

profiling of pore waters coupled to deterministic models [19,20]. The former methodology has revealed that the nitrogen species fluxes at the SWI varied among different aquatic environments, but it is difficult to infer the detailed relevant transformation processes involved at different depths of sediment. For the purpose of both identifying the nitrogen cycling processes and estimating the effluxes and influxes of nitrogen species at the SWI, the later methodology could provide more information on the nitrogen fate at the SWI. Additionally, stable isotopes approaches have been applied extensively in recent studies on watershed-scale nitrogen cycling for better understanding of the nitrogen transformation processes [21–23] and for identifying the source of particulate matter.

Most current studies on nitrogen transformation and diffusion of sediment paid close attention to aquatic systems, such as rivers, lakes, and oceans, while the knowledge on the fate of nitrogen and related key influencing factors at the sediment and water interface in deep reservoirs is still limited. In order to understand the effects of seasonal thermal stratification on nitrogen transformation and to quantitatively evaluate the contamination of endogenous nitrogen in a deep reservoir, the concentrations of nitrogen species and stable isotopes (δ15N-PON) of particulate organic matter were analyzed in sediments, which were collected during SP and NSP in Wujiangdu (WJD) reservoir, a typical deep canyon eutrophic reservoir in Southwest China. The study aimed to (1) clarify the seasonal variation characteristics of sediment nitrogen species in a deep reservoir; (2) to understand how the DO affected the transformation and diffusion of sediment nitrogen by comparative study between SP and NSP; (3) to quantitatively estimate the nitrogen exchanges at the SWI and assess the influence of seasonal thermal stratification. The results of this study provide implications useful for assessing the effects of endogenous nitrogen on nitrogen contamination in a deep reservoir that has characteristics of seasonal thermal stratification and nitrogen transformation and transportation that are similar to other deep reservoirs worldwide, and helping to fulfill water management policies in watersheds in consideration of the nitrogen contamination originating from the SWI.

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

#### *2.1. Study Area*

Wujiang River Basin is one of the largest tributaries on the south bank of the Yangtze River, and the WJD reservoir is the earliest (1979) seasonal regulating hydropower station built and put into use in China. The WJD reservoir is located in the middle reaches of the Wujiang River Basin (106◦47 E, 27◦18 N, 27,790 km2), which drains the Southeast Asian Karst Region [24,25]. The average annual flow rate of the WJD reservoir is 502 m3·s−1, and the storage area is 47.8 km2. There are abundant coal resources and industrial production bases in the area, such as a Wujiang manganese steel plant, phosphate fertilizer plant, and chemical fertilizer plant. Historically, intensive human activity has caused a serious eutrophication in the WJD reservoir, mainly on account of perennial cage fish aquaculture; hence, the hydrochemistry of the WJD reservoir is characterized by its high nitrogen loading [5].

The WJD reservoir, a deep canyon artificial reservoir, had significant seasonal variation in the profiles of water chemical parameters (especially the water T and DO). The WJD reservoir becomes thermally non-stratified from November to March, while remaining thermally stratified from April to October [5]. The water T difference between upper and lower layers reaches to 12.4 ◦C during the SP, while the water profile is homogeneous in the NSP, with small temperature differences (0.11 ◦C) among the vertical water mass. In Spring (typically in April) and Autumn (in October), the water layers are in a transitional stage of thermal stratification, and the most significant thermally stratification and nonstratification occurs in July and January, respectively. Thus, it was ideal to study the influences of seasonal thermal stratification on nitrogen transformation and transportation in the sediment-water interface by comparing those in July and January. The sediment type in the WJD reservoir is fine-grained, with a high proportion of clay [26]. The water chemical parameters in the WJD reservoir are shown in Table 1.

**Table 1.** Water chemical parameters in the WJD reservoir from our previous study [27].

#### *2.2. Sampling and Analysis*

A sampling site was chosen at the front of the dam (<1 km to the dam) in Wujiangdu Reservoir, where the depth of the water reached to 113 m and 95 m in the NSP and SP, respectively. Based on previous studies that indicated that sediment cores in lakes and reservoirs are generally consistent, single or two sediment cores would thus be enough for studying the transformation and spatial distribution of nutrients and heavy metals in a lake or reservoir [28–31]. Moreover, this study focused on nitrogen transformation and transportation mainly in the deepest parts of the reservoir (close to the dam), where the area was relatively small and the sediment was consistent. Referring to previous studies [28–31], two sampling campaigns were conducted on the NSP (January 2018) and SP (July 2019), respectively, and two sediment cores were collected by a sediment corer (Hydro-Bios, Inc., Altenholz, Germany) at the front of dam in the reservoir (<1 km from the dam), which is the deepest site in the reservoir (Figure 1). The sediment core collected in the NSP was 58 cm long and that in SP was 52 cm; the diameter of the sampling tube was 7 cm. Sample cores were sliced and sub-sampled immediately after sampling in the field (the first 20 cm was sliced at a 1 cm depth interval and at 2 cm intervals afterwards), then stored in anaerobic sample tubes, as in the studies of Hogarh et al. and Copetti et al. [28,31]. These samples were centrifuged by a centrifugal machine placed in a N2-filled container to separate pore water and were then filtered with a 0.45 μm cellulose acetate membrane. After these pretreatment process, they were stored at 4 ◦C away from light. In the laboratory, after the centrifuged sediment samples were freeze-dried at −30 ◦C, they were ground with an agate mortar, passed through a 100-mesh sieve, and then stored in polyethylene bags under 4 ◦C till analyzed. The overlying waters were sampled by two methods. First, for the overlying water close to the sediment (depth < 50 cm), we collected the water samples using a silicone tube (diameter 6 inches) by siphon method in the sediment corer before sediment subsampling; second, the overlying water samples above the sediment (>1 m) were sampled with a water sampler (Hydro-Bios, Altenholz, Germany).

**Figure 1.** Map showing the study area and the sampling sites at the WJD reservoir.

Water-soluble nitrogen was extracted from 2 g of the ground sediment samples using 10 mL deionized water, then 10 mL 2 mol·L−1 KCl was added and processed with the same method to obtain the absorbed inorganic nitrogen—the potential nitrogen of which could be replaced by special iron. Then, the concentrations of nitrogen species (NH4 +-N, NO3 −-N, and total nitrogen, TN) of pore waters and water-soluble and absorbed nitrogen species were determined by continuous-flow analysis (AA3 Auto Analyzer, SEAL, Norderstedt, Germany). The nitrogen detection limit was 0.001, 0.003, 0.003, and 0.01 mg·L−1 for NO2 −-N, NO3 −-N, NH4 +-N, and TN, respectively. NO2 −-N was under the detection limit. Laboratory standards and replicated samples were employed to keep the precision of NH4 +, NO3 −, and TN concentration analysis higher than ±5%. The dissolved organic nitrogen (DON) had a differential value to TN and total inorganic nitrogen (TIN, including NH4 +-N and NO3 −-N).

Processed sediment samples were weighed at 0.5 g and 15 mL 0.5 mol·L−1 HCl and 2 mol·L−1 KCl solution was added to remove inorganic carbon and nitrogen, then samples were washed repeatedly with deionized water until becoming neutral. The samples were freeze-dried and ground afterwards for nitrogen-stable isotopes (δ15N-PON) and C/N molar ratio analysis by stable isotopic mass spectrometer (MAT-253) and elemental analyzer (Elementar, Rhine main, Germany) at Tianjin Normal University.

#### *2.3. Calculation of Nitrogen Diffusion at SWI*

The diffusion of water-soluble nitrogen at the SWI was driven by the nitrogen gradients and calculated by Fick's first law, following equation [19,20]

$$\mathbf{F} = \mathbf{D}\_{\mathbb{S}} \times \boldsymbol{\varrho} \times (\mathbf{\hat{\partial}c/\partial z}) \tag{1}$$

where F (mg·m−2·d−1) is the diffusion flux of nitrogen, and *∂*c/*∂*z (mg·L−1·cm−1) is the gradient of dissolved substance across the SWI, which could be fitted by exponential function using the nutrient concentration of overlying water and the depth of SWI of about 5 cm. Then, the derivative of exponential function when Z = 0 was taken as the concentration gradient. C (mg·L−1) was the concentration of nitrogen, and Z (cm) was the vertical distance, which started from the upper boundary of the SWI and increased with depth. ϕ is the surface sediment porosity and can be approximately estimated by following equation [32]:

$$
\varphi = (\mathbf{r} - \mathbf{s}) / \mathbf{r} \tag{2}
$$

where r is the wet weight of sediment, while s is the dry weight. DS (cm2·s−1) is the actual diffusion coefficient under different T conditions, which can be expressed as follows [32,33]:

$$\mathbf{D}\_{\rm S} = \boldsymbol{\varphi} \cdot \mathbf{D}\_0 \left( \boldsymbol{\varphi} < 0.7 \right) \tag{3}$$

$$\mathrm{D}\_{\mathrm{S}} = \boldsymbol{\varphi}^{2} \,\mathrm{D}\_{0} \left( \boldsymbol{\varphi} \ge \mathrm{0.7} \right) \tag{4}$$

where D0 is the diffusion coefficient in infinite dilution, and the values for NO3 −-N and NH4 +-N at 25 ◦C are 19.0 × 10−6 and 19.8 × 10−6 cm2·s−1, respectively. Thus, it could be corrected according to the following formula due to the different T [34]:

$$\text{D}\_0\text{ (NO}\_3\text{ - N)} = 19.0 \times 10^{-6} + 0.4 \times (\text{T} - 25 \, ^\circ \text{C}) \tag{5}$$

$$\text{D}\_0\text{ (NH}\_4\text{ }^+\text{-N)} = 19.8 \times 10^{-6} + 0.4 \times \text{(T} - 25 \text{ }^\circ\text{C)}\tag{6}$$

#### *2.4. Data Processing*

IBM SPSS Statistics 19 was used for data analysis. Analysis of variance (ANOVA) was performed to examine the differences in the nitrogen concentrations in the SP and NSP. Additionally, *t*-tests were used to identify the variables of the C/N ratio and δ15N-PON in sediments of before and after artificial dam construction. Significance levels were reported to be *p* < 0.05 and *p* < 0.001. Microsoft Excel was used to analyze the linear correlation

of variables and obtain the value of the correlation coefficient (R2). All the data were completed by Origin 2018, Grapher 15, and Microsoft Office 2010. The concentrations of nitrogen species in the sediment of the WJD reservoir are listed in Table 2.

**Table 2.** The concentrations of nitrogen species in sediment of WJD reservoir.

#### **3. Results**

#### *3.1. Seasonal Variation of Inorganic Nitrogen in Sediment Pore Water*

The concentrations of NH4 +-N in pore water ranged from 5.7 to 27.6 mg·L−1 and from 15.3 to 31.4 mg·L−1 (averaged 14.9 and 21.0 mg·L−1) in the NSP and SP (Table 2), respectively. There were significant seasonal differences in vertical distribution of NH4 + and NO3 − in the upper layer (0–19 cm) of the sediment between the two periods, but the distribution nearly synchronously varied at deep layers (>19 cm) (Figure 2). The concentrations of NH4 +-N were lower in the upper layer of the sediment, then increased with depth in the NSP, while reversing in the SP. Inorganic nitrogen in pore water was mainly composed of NH4 +-N, accounting for 90.8~99.7% in the NSP. On the contrary, the concentration of NO3 −-N in sediment pore water was only determined at 0–19 cm depth (the concentration ranged from 0.056 to 0.58 mg·L−1) in the NSP, while it was under the detection limit (0.001 mg·L−1) in the SP and the deep layer of the sediment in the NSP. NO3 −-N was detected only in the NSP at the depth of 0–19 cm with a low concentration in pore water, while the concentration of NH4 +-N notably exceeded NO3 −-N (*p* < 0.001) in this period.

**Figure 2.** Characteristics of the vertical distribution of inorganic nitrogen in sediment pore water.

#### *3.2. Seasonal Changes in Water-Soluble and Adsorbed Nitrogen in Sediment*

Concentrations of adsorbed total nitrogen (TN) ranged from 120.7 to 398.56 mg·kg−1 in the NSP and from 91 to 255.7 mg·kg−1 in the SP, with an average value of 233.9 and 159.2 mg·kg−1, respectively. The water-soluble TN ranged from 57.4 to 339.4 mg·kg−1 and from 63.5 to 121.36 mg·kg−1 (averaged 144.5 and 86.8 mg·kg−1) in the NSP and SP (Table 2), respectively. The adsorbed TN was 1.4 to 2.7 times higher than that of the water-soluble TN. Meanwhile, the concentrations of nitrogen species displayed different seasonal variation; adsorbed TN in the sediment at 0–10 cm in the SP increased with depth, and there was no significant change in water-soluble TN (Figure 3).

**Figure 3.** Characteristics of the vertical distribution of absorbed nitrogen species in sediment.

The vertical distribution of adsorbed TN in two seasons was roughly the same as those of DON, and there were obvious seasonal differences in the upper layer (0–10 cm) of sediment (Figure 3). In addition, DON concentrations were higher and decreased with depth in the NSP, while displaying the opposite in the SP (Figure 4). Figure 4 shows the vertical variation of different nitrogen species in both the SP and NSP, with the statistics of each nitrogen species including the whole data of the profile. Moreover, the concentration of adsorbed inorganic nitrogen varied in a smaller range than that in DON, and NH4 +-N was the main component, accounting for 80.5% to 93.3% of the TN. The concentrations of absorbed NH4 +-N ranged from 48.0 to 90.1 mg·kg−1 in the NSP and from 58.1 to 85.7 mg·kg−1 in the SP, while the concentrations of absorbed NO3 −-N ranged from 5.8 to 24.75 mg·kg−1 in the NSP and from 4.8 to 12.7 mg·kg−1 in the SP. DON was the dominant component of TN (Figure 4), ranging from 52.8 to 300.1 mg·kg−1 in the NSP and from 13.0 to 176.4 mg·kg−1 in the SP, with a mean proportion ranging from 47% to 60%.

**Figure 4.** Concentrations of absorbed nitrogen species in the WJD reservoir.

*3.3. Characteristics of the Vertical Distribution of Stable Nitrogen Isotope (δ15N-PON) and the C/N Ratio in Sediment Particulate Organic Matter*

The profiles of the C/N ratio displayed no obvious seasonal variation when increasing with depth (Figure 5), and the values varied from 9.1 to 16.2 (mean value was 11.8) in the NSP and from 10.1 to 15.3 (mean value was 12.2) in the SP. The δ15N-PON values ranged from +5.7‰ to +6.8‰ and from +4.9‰ to +7.2‰, with an average value of +6.3‰ and +5.7‰ in the NSP and SP, respectively. Additionally, the δ15N-PON values decreased slightly when compared with those from 10 years ago (from +6.84‰ to +13.64‰), and there was inconspicuous difference in the vertical distribution of δ15N-PON in the NSP. Contrarily, it was roughly divided into two stages in the SP: the first stage was 0–15 cm (the mean value was +6.0‰, with the highest value of +7.2‰ and the lowest value of +4.9‰) and generally decreased with depth; the second stage was 16–52 cm (mean value was +5.4‰), where the δ15N-PON varied in a small range.

**Figure 5.** Characteristics of the vertical distribution of δ15N-PON and C/N ratio in particulate organic matter.