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.
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,
9,
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,
14,
15,
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,
26,
27,
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 km
2, most of it (142,000 km
2) 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.
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:
an increase in the amount of metal in the wastewater of all localities with controlled wastewater discharges;
a salvo of a significant amount of HM into the watercourse as a result of an emergency discharge in one of the localities;
complete exclusion of anthropogenic impact on the catchment.
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].
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.