Introduction

The fresh waters of the world are collectively experiencing accelerating rates of degradation, affecting their physical and chemical properties and biodiversity (Wetzel 2001). Well established water monitoring programs provide important information on the effects of changing environmental, socio-economic, land-use, and climatic drivers on aquatic ecosystems (e.g., Magnuson 1990; Stow et al. 1998). Long-term records of water composition thus represent a valuable chronicle of environmental changes and a framework for understanding major terrestrial and aquatic processes affecting water quality in catchment-water systems (Lovett et al. 2007; Chapra et al. 2012; Kopáček et al. 2017). Water composition in rivers is fundamentally a reflection of hydrology, natural and anthropogenic catchment sources of water constituents, and their internal physical–chemical and microbial transformations, retention, removal, and production. Analyses of long-term data on water composition and the development of anthropogenic pollution enable estimates of natural background chemistry (Smith et al. 2003), as well as providing a basis for the most effective strategies for catchment management necessary for pollution control (Zessner et al. 2016; Dupas et al. 2018).

In this study, we evaluated trends of physical–chemical properties of water in the Slapy reservoir (Vltava River) during the last six decades (1960–2019). The Slapy catchment is in a central European region that underwent pronounced changes in land use, legislation, and inputs of sulphur (S), nitrogen (N), phosphorus (P), and base cations (BCs = sum of Ca2+, Mg2+, Na+, and K+) to the atmosphere, forests, agricultural soils, and wastewaters due to changes in socio-economic relationships from a market to a planned economy in the late 1940s and back to a market economy in the early 1990s (Kusková et al. 2008; Kopáček et al. 2017). The related changes in atmospheric pollution, agricultural practices, soil properties, sanitary systems, industrial and domestic wastewater production and treatment, legislation, and environmental protection have resulted in pronounced upward and downward trends in the concentrations and fluxes of all major water constituents in the Vltava River since the mid twentieth century (Kopáček et al. 2017). Previous studies on the Vltava water composition (Kopáček et al. 2013a, b, 2014a, b, 2017; Vystavna et al. 2017) mostly focused on anthropogenic effects on terrestrial exports of ions and P, but omitted the effects of climate change, seasonal variations in water properties, interactions of individual water constituents, and effects of changing eutrophication on water characteristics. The major aim of this study was thus to summarize and evaluate effects of all known historical changes in the intensity of point and diffuse pollution sources and climate in the catchment on the long-term trends and seasonality of physical–chemical properties of the Vltava water. Special attention was paid to (1) still non-evaluated water properties like water temperature, discharge, and concentrations of dissolved organic carbon (DOC), individual nitrogen forms, and chlorophyll a, and (2) joined effects of external and internal biogeochemical processes on the river water pH.

Materials and methods

Site description

The upper Vltava catchment (12,968 km2, elevation 271–1378 m) stretches from the mountain range along the borders of the Czech Republic with Austria and Germany to the Slapy reservoir (14.4°E, 49.8°N; Supplementary Information, SI, Fig. SI-1). The Slapy reservoir, filled in 1957, is a narrow, canyon-type reservoir ~ 43 km long, with a maximum depth of 58 m, volume of 270 million m3, and average water residence time of 24–84 days (Straškraba et al. 1973). The bedrock of the Slapy catchment is dominantly metamorphic and plutonic rocks (gneiss, schist, granulite, amphibolite, granite, and granodiorite). Currently, agricultural land (52%), forests (42%, mostly plantations of Norway spruce; Picea abies), surface waters (3%), and urban areas (3%) cover the catchment area (Kopáček et al. 2017). The area of the upper Vltava catchment is to a large extent identical to the area of the administrative South Bohemian Region (11,347 km2; Fig. SI-1), with available (1) annual statistics on population, sanitary systems, wastewater treatment, agricultural activities (drainage, application of synthetic fertilizers, livestock numbers, yields of major products), forestry and industry (Yearbooks by the Czech Statistical Office), and (2) regional mean precipitation and air temperature (data by the Czech Hydrometeorological Institute). For details on trends in these variables see Kopáček et al. (2013a, b, 2014a, b) and Hejzlar et al. (2016). Daily data on precipitation and snow cover in the mountain area of the study catchment are from Churáňov station (13.69°E, 49.08°N, elevation of 1118 m) operated by the Czech Hydrometeorological Institute.

The volume of surface waters in the Slapy catchment almost tripled from 0.57 to 1.66 km3 in 1959–1991 due to the construction of eight reservoirs (Table SI-1). The average residence time of surface water in the catchment thus increased from ~ 0.2 to ~ 0.7 year (Kopáček et al. 2013a). The biggest reservoir (Orlík, max. depth 74 m, mean water residence time of 89 days) was built ~ 9 km upstream of the Slapy inlet and filled in 1962. The outflow of water from the Orlík reservoir was permanently from its hypolimnion, except for the period from August 2002 to June 2003, when water was discharged by the upper overflow.

Water sampling and analysis

Water samples have been taken at 3-week intervals since March 1959 from the 0.5 m deep surface layer above the maximum depth, 8.8 km upstream from the Slapy dam. Sampling was usually performed in the morning (between 8 and 10 am). Water transparency (Secchi depth), pH (colorimetric determination), and vertical profiles of temperature and concentration of dissolved oxygen (O2) were measured during sampling. For titrimetric O2 determination, water was collected in a separate sample bottle. For other chemical analyses, water was immediately filtered through a 40-µm (or 200-µm since 1993 and 1995 for P and all water constituents, respectively) polyamide sieve to remove zooplankton and/or other coarse particles, transported to the laboratory, and analysed within 2–24 h after sampling. The effect of the new field filtration method was small (e.g., < 5% for total P concentration; see SI-Part-2) and neglected. For chlorophyll a determination (spectrophotometrically after acetone extraction according to Lorenzen 1967), samples were taken from the uppermost (0–4 m deep) water layer with a 4 m long tube and were not filtered (this sampling was irregular till 1975, then regular). The 21-day sampling interval was lower than the water residence time in the Slapy reservoir, and sample-to-sample variability in the concentrations of water constituents was usually small, except during floods. Previous studies have shown that the samples taken from the reservoir surface reasonably represented concentrations of most solutes (except for P) along the whole water column (Procházková and Blažka 1986; Vystavna et al. 2017).

In the laboratory, water was analysed for concentrations of acid neutralizing capacity (ANC; titration to the pH end-point of 4.5 till 1995, then Gran titration), ions (Ca2+, Mg2+, Na+, K+, NH4+, NO3−, NO2−, SO42−, Cl−), DOC, total phosphorus (TP), dissolved P (DP, since 1995), total organic nitrogen (TON), dissolved organic N (DON, since 1995), dissolved reactive silica (Si, since 1997), and chlorophyll a. For the determination of ions and dissolved Si, samples were filtered with membrane filters (pore size of 0.45 μm). DOC, DP, and DON were analysed after sample filtration with glass-fibre filters (pore size of 0.4 µm). Particulate C (PC) was determined in the seston retained on glass-fibre filters. Samples for pH (see SI-Part-2 and Fig. SI-2 for more details on its determination), electric conductivity (at 25 °C), ANC, TP and TON were not filtered beyond the field pre-filtration. Due to low concentrations of other buffering systems in the Vltava water, we assumed that ANC was equal to the HCO3– concentration. Particulate P (PP) and N (PN) were differences between their respective total and dissolved forms. Ratios of PN:PP and PC:PN were used to estimate the respective molar N:P and C:N ratios in seston.

Some water constituents (e.g., NO3−, NO2−, TON) were determined using the same method through 1960–2019. Other methods were changed during the monitoring period, but the old and new methods were always used in parallel for at least one year to be either sure that they provided identical data, or to obtain a significant relationship between the results, enabling homogenization of the whole long-term trend (see SI-Part-2 for more details on individual methods used in this study). The reliability of the analytical results was checked by means of an ionic balance approach (Kopáček et al. 2000), a comparison between measured and calculated conductivities, and a standard sample (SI-Part-2).

Daily data on the flow of the upper Vltava, surface water temperature, and occurrence of ice cover in the Slapy reservoir are from the databases of the Czech Hydrometeorological Institute and the Vltava River Board. Monthly data on mean precipitation and air temperature in the South Bohemian Region are from the databases of the Czech Hydrometeorological Institute.

Data evaluation and statistical methods

In this study, we used monthly means for all physical–chemical properties through the 1960–2019 period. The monthly data were arithmetical means of the 3-week data in the respective month. When solute concentrations were below their detection limits (SI-Part-2), half of these limits were used in the data evaluations. This occurred only exceptionally for NH4+ and NO2−. The per cent saturation of water with dissolved O2 was computed from its measured and equilibrium concentrations; the latter values were based on water temperature and elevation according to Mortimer (1981). Concentrations of H2CO3*, i.e. the sum of dissolved (and hydrated) carbon dioxide and carbonic acid, were computed from ANC and in-situ pH and water temperature according to Stumm and Morgan (1981). Annual mean NO3− and NO2− concentrations were discharge-weighted means computed from their respective monthly means and discharges.

Physical–chemical parameters were aggregated with mean monthly values within each decade to visualize seasonal patterns and to qualitatively assess changes in seasonal patterns across decades. These decadal means were calculated as geometric means for each unique month within each of the six consecutive decades, from 1960 to 2019. The geometric mean was used to eliminate the effect of non-normal distribution of values. The mean pH values were based on the geometric mean of H+ concentrations.

The trend analyses were done for monthly averages in the R environment for statistical computing (R Core Team 2019), using the function "decompose()". The Seasonal Mann–Kendall Trend Test (Hirsch–Slack Test) (Pohlert 2018) was used to detect seasonal trends in the series. Magnitudes of trends were computed by Sen's slope test. A non-parametric test after Pettitt was used to test for a breakpoint in the central tendency of a time series (Pohlert 2018). When the breakpoint was detected, the time series was split, and trend analysis was done for each part separately. Linear regressions against time were used to evaluate significance (p < 0.05) and slopes of trends in selected water characteristics for individual months over the study period (n = 60 for DOC and pH; n = 59 for precipitation and air temperature).

To evaluate the timing of spring snowmelt between January 1 and April 30, we determined the day of the year until when 75% of the cumulative inflow into the Slapy reservoir occurred, following the method by Hodgkins et al. (2003).

Thermal stratification of water column was assessed from vertical temperature profiles using the R package "rLakeAnalyzer" (Winslow et al. 2015). The water column was considered stratified when the maximum vertical temperature gradient in the water column was dtw/dz > 0.5 °C m−1 (tw = water temperature in °C, and z = depth in m) and the Schmidt stability index in the upper 30 m of the water column was at least 120 kg m−2.

The complete original 3-week data on chlorophyll a concentrations were used to analyse long-term changes in reservoir phenology from 1975 to 2019. We followed the method proposed by Mackas and Beaugrand (2010), which allows the determination of cardinal dates in seasonal succession. It is based on interval estimates and uses thresholds, the timing of cumulative percentiles, and peak maxima to index the start-of-season, peak maxima, and end-of-season. We calculated the length of the vegetation season from the start to the end of the season defined as the interval from 15 to 85 of cumulative percentiles of annual chlorophyll a concentration. We further analysed the position of spring and summer phytoplankton peaks, defined as the maximum of chlorophyll a concentrations for the periods from March 1 to June 20 and from June 22 to September 30, respectively.

Results

Basic physical–chemical variables

The median value of the Vltava discharge was 68 m3 s−1 in the Slapy reservoir from 1960 to 2019. The flow exhibited a pluvial regime, characterized by high discharge in winter and spring and a minimum in the late summer. In the 1960s, the trend of seasonal discharge was affected by the filling of the Orlík reservoir (mostly in the spring 1962 and 1963), which reduced the ten-year mean discharge by ~ 20% in March and April. During following two decades, the spring maximum flow fluctuated around April 1 (± 9 days), and then significantly (Mann–Kendall test; p = 0.04) shifted by 13 days to March 19 (± 6 days) between 1995 and 2019, while there were no changes in water management and operation of both Orlík and Slapy reservoirs. The shift in the spring maximum flow (Fig. 1a) occurred despite relatively stable seasonal patterns of precipitation from 1960 to 2019, with maximum values typically occurring from June to August (Fig. 1b). Trend analysis for the whole study period showed that significant (p < 0.05) changes in discharge occurred only in May and December (decreases), and precipitation significantly increased only in January (Fig. SI-3).

Fig. 1
figure 1

Seasonal variations of mean physical–chemical variables of water in the surface layer of the Slapy reservoir computed as geometrical means for individual months over six decades from 1960 to 2019: a discharge; b average precipitation in the Slapy catchment; c water temperature in the surface layer; d transparency (Secchi depth); e the concentration of dissolved oxygen; f per cent water saturation with dissolved oxygen; g conductivity at 25 °C; h the sum of dissolved carbon dioxide and carbonic acid (H2CO3*); and i pH

Full size image

Water temperature of the Slapy epilimnion reached maximum values exceeding 20 °C (up to 27 °C) in July–August, with the seasonal values higher in the 2000–2019 period than in the four preceding decades (Fig. 1c). Water temperature significantly increased in all months (except for October and November) and by 2.1 °C on an annual basis during the study period, reflecting a similar increase in air temperature in the South Bohemian Region (by 1.9 °C on an annual basis) (Fig. SI-4). The most pronounced increase in air temperature (by 2.6–2.7 °C) occurred in winter (December–January) and summer (July–August), but was insignificant (< 1.0 °C) during autumn (September–October).

The elevated air temperature in winter was the most likely reason for the shift in maximum discharge due to a decreasing proportion of winter precipitation in the form of snow. The annual cumulative amount of fresh snow significantly (p < 0.05) decreased from 416 to 296 cm year−1 (on average by 2 cm year−1) in the mountain area of the catchment (Churáňov station) from 1961 to 2019, representing a difference of ~ 114 mm year−1 of precipitation at the 1961–2019 average water equivalent of 1.09 mm cm−1 of snow. This means that ~ 10% of the average annual precipitation amount at Churáňov station of 1092 mm year−1 changed from snow to rain during the study period and resulted in lower snow cover (Fig. SI-5), more frequent thaws, earlier snowmelt, and elevated discharge already during winter, especially after 1995.

The elevated air temperature also significantly (p < 0.01) shortened ice-cover periods and advanced spring onset of thermal stratification (p < 0.001) of the Slapy reservoir by 19 days on average from 1960 to 2019 (Fig. 2).

Fig. 2
figure 2

Duration of ice cover (black bars) and stable thermal water column stratification (grey bars) in the Slapy reservoir during the 1960–2019 period. Nonparametric Mann–Kendall test indicated significantly decreasing the duration of ice cover (p < 0.01) and advancing spring onset of thermal stratification (p < 0.001)

Full size image

Seasonal variations in transparency were the least pronounced and exhibited the lowest values during the 1960s (i.e., the first decade after the reservoir filling) compared to during the rest of study (Fig. 1d). On a seasonal basis, the highest transparency usually occurred in winter, then sharply decreased, and reached their minima in April–May. The transparency increase corresponding to the "clear-water" phase of the seasonal succession (Sommer et al. 2012) usually occurred in June, and then decreased to the less pronounced minima in July–August, after which it continuously increased till winter.

Saturation of the Slapy surface layer with dissolved O2 exhibited seasonal variation similar to O2 concentrations, with the lowest values (35–80%; 4–6 mg L−1) during the autumn turnover of the water column in October–November and the highest values (oversaturation of 101–140%; 9–13 mg L−1) typically occurring from April to August (Fig. 1e, f).

Ionic composition

The long-term trends and seasonal variations in concentrations of ions were reflected by pronounced changes in electric conductivity (Fig. 1g). Trends in concentrations of individual ions, however, exhibited patterns differing from each other both in seasonal variations and long-term trends (Fig. SI-6). The lowest concentrations were observed for most ions in the 1960s (except for SO42− with minima in the 2010s), while the highest values typically occurred in the 1980s (except for Na+ and HCO3− with maxima in the 2010s). Seasonal variations in concentrations of SO42−, Cl−, and divalent base cations (Ca2+ and Mg2+) were similar and relatively small, with slightly elevated values in summer. In contrast, concentrations of monovalent base cations did not exhibit seasonal variations (except that for Na+ in the 2010s) and differed in their long-term trends. Concentrations of Na+ were lowest at the beginning and highest at the end of our study, while K+ concentrations were highest in the 1980s and lowest during the 2000–2019 period.

In contrast to the more conservative ions, concentrations of NH4+, NO3−, and NO2− exhibited more pronounced seasonal variations (Fig. 3). Concentrations of NH4+ were usually higher in the first half of a year, with maxima between March and June, and then decreased to averages < 2 µmol L−1 between July and December. This seasonal variation was most pronounced in the 1960s, while almost disappeared from 1990–2019 (Fig. 3a).

Fig. 3
figure 3

Seasonal variations of mean concentrations of nutrients in the Slapy surface layer computed as geometrical means for individual months over six decades from 1960 to 2019: a ammonium (NH4+); b nitrate (NO3−); c nitrite (NO2−); d dissolved organic carbon (DOC); e chlorophyll a; f total organic nitrogen (TON); g total phosphorus (TP); h dissolved reactive silica (Si); and i total nitrogen (TN = sum of NO3−, NO2−, NH4+, and TON)

Full size image

Concentrations of NO3− (the dominant N form throughout the study) were highest in April–May (Fig. 3b), i.e., about a month after the maximum discharges (Fig. 1a). The most pronounced (and consistent throughout the study) seasonal variations among ions occurred for NO2−, with negligible values in winter and maxima (up to 4 µmol L−1) in June–July (Fig. 3c). Despite different seasonal variations, NO2− concentrations exhibited a similar pattern as NO3− on a long-term basis, and NO3− explained 63% of the variability in annual mean NO2− concentrations (Fig. SI-7).

Organic matter, nutrients, and chlorophyll

Concentrations of DOC were highest in the Slapy reservoir in the 1960s (Fig. 4), with an average of 733 µmol L−1 and maxima in winter (Fig. 3d). The DOC concentrations decreased until the early 1990s, then with a breakpoint in 1992 began to increase, and both these trends were significant (Fig. 4). Despite the recent increase, DOC concentrations were still significantly lower in the 2010s than in the 1960s, especially in winter and spring (November–May). These months also exhibited a significant increase in DOC concentrations after 1992, while this increase was insignificant in summer (Fig. SI-8).

Fig. 4
figure 4

Schematic description of changes in major point and diffuse sources of DOC and TP in the upper Vltava catchment (a), and monthly mean values of DOC (b), TP (c), and pH (d) in the Slapy surface layer from January 1960 through December 2019. Red lines show significant trends of seasonal Mann–Kendall test (R Core Team 2019) and the green vertical line shows the breakpoint. Arrows in the panel (a) show: (1) decreasing DOC production with municipal wastewaters; (2) decreasing DOC production from wood pulp and paper productions, (2a) shut down of one paper mill in 1966 and beginning of partial wastewater treatment in 1968, (2b) beginning of recycling and treatment of all wastewaters in 1990, and (2c) end of pulp production in the catchment in 2002; (3) increasing DOC leaching from forests due to decreasing SO42− deposition and climate change since the 1980s, and tree dieback since the 2000s; (4) increasing P fertilization rate of farmland till 1989 and (4a) its sharp decrease since the 1990; and (5) P production with municipal wastewaters increasing till 1990, then decreasing due to (5a) partial reduction of P use in detergents and (5b) P control in effluents from large wastewater treatment plants (2005) and a ban of P-detergents (2006)

Full size image

Seasonal variations in TP, Si, and TON concentrations reflected the primary production of phytoplankton in the reservoir. Concentrations of chlorophyll a and TON had two maxima in April–May and July–September, while TP concentrations exhibited flat seasonal minima from May through September, and Si reached a sharp minimum in July (Fig. 3). On a long-term basis, TP concentrations significantly increased until the mid-2000s and then started to significantly decrease, with a breakpoint in 2004 (Fig. 4c). Concentrations of chlorophyll a exhibited a similar trend with increasing values until 2004, then decreasing and levelling off (Fig. SI-9A). In contrast, the seasonal Mann–Kendall trend test showed that TON concentrations increased significantly throughout 1960–2019 (Fig. SI-9B).

The timing of the seasonal phytoplankton succession showed remarkable temporal trends (Fig. 5). In 1975–2019, the start of phytoplankton succession advanced by 14 days, and its end delayed by ~ 22 days, which resulted in a substantial increase of the whole vegetation season by more than a month. Changes in phytoplankton phenology were the most pronounced in spring when the phytoplankton peak significantly advanced by 21 days (p < 0.01). A similarly large delay in the summer phytoplankton maximum was not, however, statistically significant.

Fig. 5
figure 5

Long-term changes in phytoplankton phenology in the Slapy surface layer from 1975 to 2019. Panel (a) depicts temporal trends in the start and end of the vegetation season, and panel (b) shows Julian days of spring and summer maxima in concentrations of chlorophyll a. Solid lines represent linear regressions and p values their significance. The dashed line shows a non-significant relationship. Slopes (day year−1) show average changes in the variables between 1975 and 2019

Full size image

The average (± standard deviation) molar DOC:TON ratios significantly (p < 0.001) decreased from 23 ± 8 in the 1960–1965 period to stable values of 12 ± 2 during the 1992–2019 period. The 1995–2019 average molar N:P and C:N ratios of the seston were ~ 20 ± 10 and ~ 10 ± 5, respectively, and did not exhibit any significant long-term trend.

Water pH and H2CO3*

The pH values increased in the Slapy surface layer from 1960 to the early 1990s, and then started to decrease, with a breakpoint identified by the seasonal Mann–Kendall test in 1992 (Fig. 4d). Similar breaks in long-term trends also occurred for the annual pH maxima, while the pH minima increased significantly through the whole 1960–2019 period (Fig. 6). Despite the slow pH decrease since 1992, the pH values still were significantly higher in the Slapy water in 2019 than at the beginning of the study in all months except for May (Fig. SI-10).

Fig. 6
figure 6

Annual minimum (a) and maximum (b) pH values in the Slapy surface layer from 1960 to 2019. Solid lines represent linear regressions and p values their significance

Full size image

Values of water pH were hightest (7.5–9.2) in summer months and lowest (6.7–7.1) from October to March (Fig. 4d). Due to parallel (or inverse) seasonal variations, pH positively correlated with concentrations of chlorophyll a and water saturation with dissolved O2 (both p values of < 0.001), while pH negatively correlated (p < 0.001) with Si concentrations.

Concentrations of H2CO3* showed a marked annual cycle (inverse to pH; Fig. 1h, i), with high values during winter and a significant decrease during the months with stable stratification (Fig. 2). On a long-term basis, H2CO3* concentrations were higher from 1960 to 1989 than in the 1990–2019 period.

Discussion

Major anions

The dominant anthropogenic sources of major ions for surface waters in the Slapy catchment were the application of synthetic fertilizers (for SO42−, NO3−, Cl−, and K+), liming of farmland (for Ca2+, Mg2+, and HCO3−) and its draining (for NO3−), atmospheric pollution/deposition (for SO42−), wastewaters (for NH4+, NO3−, and Cl−), and road de-icing (for Cl− and Na+) as discussed in detail by Kopáček et al. (2013a, b, 2014a, b, 2017). In short: the uneconomical use of synthetic fertilizers, excessive farmland drainage, and lack of environmental protection contributed to water pollution during the planned economy period until the late 1980s (Kopáček et al. 2017). As a result, concentrations of strong acid anions and BCs increased during this period (Fig. SI-6). In contrast, an economic recession immediately following the re-establishment of a market economy and continually improved legislation (regulation of emissions of S and N compounds and agricultural practice; Bičík et al. 2001; Kopáček et al. 2013b, 2014a, b) have resulted in dramatically reduced inputs of elements to farmland from fertilizers and to forests from atmospheric deposition; Figs. SI-11 and SI-12), and consequently, in decreasing (or levelling-off) trends in water pollution from diffuse sources since the 1990s. The long-term trends of major ions thus highlighted general relationships between the intensity of agriculture and land use on riverine element exports, well known from other large European and American catchments (e.g., Kaushal et al. 2005; Raymond et al. 2008; Aquilina et al. 2012). The exceptions from this pattern was Na+ due to continuously increasing use of NaCl for road de-icing from the early 1970s (Fig. SI-11F) that increased Na+ concentrations in surface waters, especially in the 2010s (Fig. SI-6G).

Seasonal variations in concentrations of major ions were relatively small. The seasonal cycles of SO42−, Cl−, Ca2+, and Mg2+ were similar, with maximum values in summer (Fig. SI-6), due probably to lower dilution by discharge (Fig. 1a) and higher evaporation from the epilimnion. Seasonal cycles of HCO3− and NO3− exhibited inverse patterns (Figs. 3b, SI-6A) due to the HCO3− dilution by elevated discharges in spring, while NO3− concentrations increased during this period (see “Nitrate” section for details).

Dissolved organic carbon

The highest DOC concentrations in the 1960s and their continuous decrease till 1992 (Fig. 4b) reflected decreasing surface water pollution with domestic and industrial wastewaters (Fig. SI-13). In the 1960s, ~ 36% of the Slapy catchment inhabitants were already attached to sanitary systems, but only 1–10% (3% on average) of these wastewaters were treated (Fig. SI-14A). The rest of the inhabitants used latrines, and biologically stabilized waste was hauled away from cesspools and used for soil fertilization in agriculture. By the mid-1990s, more than 80% of the catchment inhabitants were attached to sanitary systems, and the proportion of treated municipal wastewaters had increased to ~ 90% (Kopáček et al. 2013b). Consequently, the pollution of surface waters in the catchment with organic matter from municipal wastewaters steeply decreased in the 1990s (Fig. SI-14B).

Wood pulp and paper productions were other important DOC sources in the Slapy catchment. Most of the wood pulp production was based on the sulphite process and reached 72–171 Gg year−1 in the 1960–1981 period (Kopáček et al. 2014a), with two major paper mills. Their wastewaters were directly discharged into the Vltava. Sulphite pulp effluents are heavily polluted with various types of wood extracts, including readily biodegradable compounds, as well as non-degradable substances such as lignosulfonates, and have concentrations of biological and chemical oxygen demands in the hundreds and thousands of mg L−1, respectively (Pitter 1990; Suriyanarayanan et al. 2005). Organic pollution associated with sulphite pulp production in the Slapy catchment can be estimated at 13 to 31 Gg year−1 of biochemical oxygen demand (BOD5) in 1960–1981, assuming an average BOD5 production of ~ 180 g kg–1 of pulp (Pitter 1990). These wastewaters were not treated until the late 1960s. Their production, however, decreased in 1966 when one of the paper mills was shut down. Another decrease occurred in 1968 when wastewaters became partially treated using evaporators at the second paper mill. These measures reduced concentrations of brownish lignosulfonates in surface water and the water transparency, which was low in the Slapy reservoir in the 1960s, increased in the following decades (Fig. 1d). In 1986, the sulphite process was partly replaced by the sulphate process. In 1990, pulp effluents began to be efficiently recycled and treated, which led to a breakpoint in organic pollution of the Vltava above the reservoir cascade (Fig. SI-13) and a pronounced minimum in DOC concentrations in the Slapy reservoir over the next three years (Fig. 4b). Since this minimum, however, DOC concentrations have been steadily increasing.

The continuously increasing DOC leaching from diffuse sources in forested parts of the Slapy catchment (Hejzlar et al. 2003; Kopáček et al. 2013c) has reversed the decreasing DOC trend, associated with changes in surface water pollution by wastewaters from point sources, since the 1990s. The most likely reasons for the increasing DOC leaching were (1) climate, causing changes in temperature, precipitation, forest health, and soil moisture (Hejzlar et al. 2003; Kopáček et al. 2018), and (2) decreasing atmospheric deposition of SO42− since the 1980s (Kopáček et al. 2013c), resulting in elevated DOC mobility due to increasing pH and decreasing ionic strength of soil water (e.g., Hruška et al. 2009; Evans et al. 2012).

Climate-related terrestrial exports of DOC from soils are usually associated with elevated intensities and frequencies of lateral water flows through the uppermost soil horizons (Bearup et al. 2014). Their probability increases during high precipitation and runoff events and is currently magnified by climate change due to more frequent torrential rains (Raymond et al. 2016). The climate-related, more frequent thaws, snowmelt, and elevated discharge during winter probably contributed to the significantly increasing DOC concentrations in the Vltava water during the winter and spring months of the 1993–2019 period (Fig. SI-8B). DOC leaching from forested parts of the Slapy catchment was further amplified after bark beetle-induced spruce mortality in forests, which has increased since the late 1990s due to increasing air temperature. The tree mortality changes soil chemical and microbial characteristics, resulting in elevated DOC leaching (Kopáček et al. 2018). In parallel, the reduced water transpiration after tree death leads to increased soil moisture, shallower groundwater, and greater water fluxes through shallow, organic-rich soils, further accelerating terrestrial DOC exports (Hubbard et al. 2013; Mikkelson et al. 2013). The DOC concentrations thus continued their increase even after 2002, when pulp production in the Slapy catchment completely ceased.

Phosphorus, silica, and chlorophyll

Concentrations of TP and chlorophyll a increased from 1960 to the mid 2000s (Fig. 3), and then started to decrease, with breakpoints in 2004 (Figs. 4c, SI-9A). Primary production in the Slapy reservoir is P-limited (Komárková and Vyhnálek 1998). Hence, the increase in TP concentrations was paralleled with enhanced primary production, here proxied by the concentrations of chlorophyll a. The elevated TP availability resulted from increasing P exports from diffuse agricultural sources and point sources, especially domestic wastewaters (Fig. SI-14C; Hejzlar et al. 2016).

Increased exports of P from diffuse sources are usually associated with high flows after intense precipitation, which intensify the erosion of soil particles with adsorbed phosphate, as well as leaching of dissolved P from the uppermost soil layers during periods of lateral water flows (Bowes et al. 2008). The risk of P export associated with erosion increased with the increasing proportion of arable land in the catchment (Kopáček et al. 2017) and intensified P-fertilization until the late 1980s (Fig. SI-11D). Since the early 1990s, the risk has been decreasing due to the reduced use of both synthetic fertilizers as well as manure and slurry (due to a sharp reduction in cattle and pig production after 1989, Kusková et al. 2008), the increasing conversion of arable fields to pastures (Kvítek et al. 2009), and an agri-environment scheme implemented in the Czech law after the accession of the Czech Republic to the European Union in 2004 (Doležal and Kvítek 2004; Vystavna et al. 2017).

The P export from point sources increased in the Slapy catchment from the 1960s to the 1990s due to the increasing number of inhabitants attached to sanitary systems, the increasing use of P-rich detergents and the low P-retaining efficiency of wastewater treatment plants (Hejzlar et al. 2016). Then, the P export from point sources began to decrease (Fig. SI-14C) due to legislation changes. These included (1) the recommended and voluntary reduced use of P in detergents since the mid-1990s, followed by a complete legislative ban in 2006; (2) regulations of wastewater treatment efficiency, limits for P concentration in effluents, and fees for the amount of P discharged being set for large wastewater treatment plants (> 3 Mg year−1 of P) since 2005; and (3) by the reduced use of P in dishwasher detergents since 2014 (Hejzlar et al. 2016; Vystavna et al. 2017). These measures resulted in a continuing decrease in the Vltava TP concentrations, similarly to freshwaters in other European countries (Hukari et al. 2016), and consequently, in lower primary production in the Slapy reservoir compared to the 1990s (Fig. 3e).

The summer decrease in the epilimnetic TP concentrations (Fig. 3g) was caused by the uptake of this limiting nutrient by primary producers, and was accompanied by seasonally elevated chlorophyll a concentrations (Fig. 3e), low water transparency (Fig. 1d), elevated water saturation with dissolved O2 (Fig. 1f), decreased H2CO3* concentrations (Fig. 1h), depleted NH4+ (Fig. 3a), and increased pH (Fig. 1i). The elevated P uptake by phytoplankton resulted in a PP contribution to TP of 49% from May to September 1995–2019, while it only formed 29% of TP on an annual basis.

A substantial seasonal drop in concentrations of dissolved Si (Fig. 3h) was associated with the increasing dominance of diatoms in the reservoir phytoplankton. Diatoms, a group of algae with siliceous cell walls, can form dense populations under favourable conditions in early summer and deplete nutrients and Si to limiting levels (Wetzel 2001; Znachor et al. 2008). Due to the presence of siliceous frustules, diatom cells are substantially heavier than water and sink down despite a stable thermal water stratification in summer. Consequently, they are continuously replaced by other phytoplankton groups with no obligate Si requirement (Sommer et al. 2012), and therefore epilimnetic Si concentrations increase again later in the season (Fig. 3h).

Ongoing analyses of long-term data show widespread, systematic changes in the seasonal timing (phenology) of ecological events over recent decades (Winder and Schindler 2004). In the Slapy reservoir, phenological changes were most pronounced in spring (Fig. 5), which is in line with the commonly accepted notion that a shift in the spring biomass peak is the most sensitive indicator of phytoplankton response to climatic forcing (Winder and Schindler 2004).

Nitrogen forms

Concentrations of total N (TN) were dominated by NO3− and TON (76 and 21% on average, respectively), while NH4+ and NO2− contributed only negligibly (2% and 1%, respectively). Long-term trends and seasonal variations of TN thus predominantly reflected the NO3– dynamic (Fig. 3). The individual N forms, however, exhibited different time-series, reflecting different physical–chemical and microbial processes governing N transformations and uptake in the whole catchment-water system, as follows.

Ammonium

The NH4+ concentrations in the surface layer of the Slapy reservoir were highest in the early 1960s and probably originated from wastewaters and mineralization of the flooded organic matter. The untreated domestic wastewaters were a significant NH4+ source for surface waters, especially in the 1960s, with minimum sanitary systems attached to wastewater treatment plants (Fig. SI-14A). During the following five decades, the proportion of biologically treated domestic wastewaters in the catchment increased and intensive water aeration during biological treatment nitrified NH4+ and decreased its export from wastewater treatment plants to surface waters (Behrendt et al. 2000; Kopáček et al. 2013b).

Another factor contributing to more effective elimination of NH4+ from waste water treatment plants was the increasing volume of surface waters due to the construction of new reservoirs (Table SI-1), resulting in a tripled water residence time in the catchment between 1960 and 1991. More NH4+ thus was assimilated and nitrified prior to entering the Slapy reservoir.

Some of the NH4+ in the surface layer of the Slapy reservoir also temporally originated from mineralized vegetation and soil organic matter on the flooded bottom, especially in the early 1960s, entering the surface layer during spring and autumn turnovers (Procházková 1966). In addition, the water outlet from the Orlík reservoir always originated from the cold (~ 4 °C), dark, and anoxic hypolimnion. Under such conditions, the liberated NH4+ could not be effectively removed by assimilation and nitrification (Wetzel 2001) and thus added to the Slapy reservoir pool of NH4+.

Concentrations of NH4+ exhibited pronounced seasonal variations in the Slapy surface layer, with elevated winter values and pronounced spring maxima during the first three decades of this study (Fig. 3a). The major epilimnetic NH4+ sinks during the growing season were probably its assimilation (NH4+ is usually the primary source of inorganic N for freshwater phytoplankton) and nitrification, as in other circum-neutral waters (Procházková et al. 1970; Wetzel 2001).

Nitrate

Similarly to other large European catchments (e.g., Kronvang et al. 2008), NO3− leaching from agricultural land was the major source for the Slapy reservoir (Procházková 1966; Kopáček et al. 2013b). The increasing trend in NO3− concentrations from the 1960s to the 1980s (Fig. SI-7) predominantly reflected high net anthropogenic N inputs to farmland from synthetic fertilizers (Fig. SI-11C) and mineralization of soil organic N due to increasingly intensive drainage and deep tillage (Kopáček et al. 2013a). Point sources of N exhibited a similar increasing trend during this period due to the increasing population attached to sanitary systems (Fig. SI-14D), but this contribution was about an order of magnitude lower than that of diffuse sources (Kopáček et al. 2013b). Similarly, the subsequent decline in NO3− concentrations between the 1990s and the 2000s primarily reflected the de-intensification of agricultural production, reduced N fertilization rates, aging of drainage systems and the partial conversion of arable land to grassland (Bičík et al. 2001; Kopáček et al. 2013a), while increasing denitrification in wastewater treatment plants played a relatively minor role (Kopáček et al. 2013b). As observed elsewhere (e.g., Randall and Goss 2008; Kvítek et al. 2009), the reduced N inputs and decreasing soil aeration were the most effective measures rapidly reducing NO3− leaching from agricultural land also in the Slapy catchment.

The high seasonal variations in NO3− concentrations, with maximum values typically occurring in spring (Fig. 3b), were associated with elevated discharge and the preceding dormant season (Procházková and Brink 1991). Nitrate accumulated in soils during dry periods (especially in winter, with negligible N uptake by vegetation) because of low hydraulic conductivity in the soil profile was mobilized and leached in the following wet period (Meisinger and Delgado 2002; Randall and Goss 2008).

Nitrite

The NO2− concentrations in the Slapy surface layer exhibited pronounced seasonal variations (Fig. 3c). They started to increase in spring together with water warming, and this process accelerated in April–May when temperatures exceeded ~ 10 °C (Fig. 1c). The maximum NO2− concentrations corresponded with the highest transparency in the clear-water periods (i.e., during the lowest chlorophyll a concentrations; Figs. 1d and 3e) and the longest and most intensive sunshine in June. In the rest of summer, the NO2− concentrations continuously decreased with decreasing water transparency and finally reached their minima after water temperatures decreased to < 10 °C in November. In contrast to the terrestrial drivers dominating the NO3− dynamic, the NO2− variation was predominantly associated with the seasonal development of in-lake processes and probably resulted from nitrification and denitrification. Nitrification could have played an important role during elevated NH4+ concentrations at the beginning of the study. However, we expect that the role of denitrification in the NO2− production increased in subsequent years with increasing NO3− availability, as suggested by the parallel NO2– and NO3– long-term trends (Fig. SI-7). The NO2− production by both nitrification and denitrification can be expected to exhibit maxima in summer for the following reasons.

The nitrification rate is low at < 10 °C and very sensitive to small changes in temperature over a range from 10 to 17 °C, with the inhibitory effect of low temperature greater for NO3− formers, such as Nitrobacter, than for NO2− formers, such as Nitrosomonas (Randall and Buth 1984). The inhibitory effect of low temperature thus temporarily results in a NO2− build-up during temperature increases above or decreases back to ~ 10 °C. Moreover, the rate of microbial NO2– oxidation is more inhibited by light then its production from NH4+ (Olson 1981; Vrba 1990).

The NO2− production associated with NO3− reduction occurs during photochemical and microbial denitrification (Spokes and Liss 1999; Stief et al. 2016). Photochemical denitrification is a first-order reaction related to the NO3− concentration and also depends on the concentration of organic matter present in the system and the intensity of solar radiation (Spokes and Liss 1999; Porcal et al. 2014). At high (and compared to NO2− relatively stable) NO3− and DOC concentrations during the year (Fig. 3), the intensity of photochemical NO3− reduction thus primarily reflects seasonal variations in transparency and sunshine intensity, both of which are high in June. We also cannot exclude microbial denitrification, despite the high summer concentrations of dissolved O2 (Fig. 1e), due to the formation of local anoxic microsites that may occur inside algal aggregates, especially diatoms (Stief et al. 2016). This process would also result in the highest NO2− production during the summer months with the highest primary production and algal concentrations in the epilimnion.

Organic nitrogen

Concentrations of TON exhibited a small but significant (p < 0.001) increase between 1960 and 2019 (Fig. SI-9B). Data on PN and DON, available from 1995 to 2019, showed that DON increased (p < 0.001), while PN concentrations significantly (p < 0.05) decreased during this period. This decrease in PN concentrations was similar to PC, PP, and chlorophyll a from 2004 to 2019. The long-term increase in TON was thus associated with the DON representing 92% of TON during the 1995–2019 period.

The decrease in molar DOC:TON ratios from ~ 23 in the 1960s to ~ 12 during the 1992–2019 period was apparently caused by reduced inputs of DOC-rich wastewaters. The stable DOC:TON ratios since 1992 indicated that the long-term increase in concentrations of organic N was probably associated with increasing terrestrial DOC export from diffuse forest sources. In-lake N assimilation by algae, however, contributed to the elevated TON (and also PN and DON) concentrations during the growing period, resulting in seasonal variations similar to chlorophyll a (Fig. 3).

pH

The long-term trend in pH and its seasonal variations, with winter minima and maxima in July and August (Fig. 4d), integrated several effects. They include changing external acidity sources, in-lake H+ producing/consuming processes, the P inputs (affecting the intensity of P-limited primary production), and changes in water buffering capacity. These processes had different effects on trends in annual pH minima and maxima (Fig. 6).

There were at least three reasons contributing to lower pH values during the 1960s–1980s than in the three following decades: (1) acidic deposition and the intensive use of synthetic fertilizers [including (NH4)2SO4)] acidified soils and the receiving waters (e.g., Reuss and Johnson 1986). (2) Decomposing organic matter on the bottom of the flooded valleys of the Slapy and Orlík reservoirs and microbial mineralization of organic pollution from domestic and industrial wastewaters were more important CO2 sources (resulting in higher H2CO3* concentrations in the Slapy water) until the 1990s than in the following years, with a breakpoint in the H2CO3* trend in 1993 (Fig. SI-15). The elevated H2CO3* concentrations lowered water pH (Stumm and Morgan 1981), especially in winter, when CO2 assimilation by algae was limited. (3) The H+ release associated with NH4+ assimilation and nitrification (Reuss and Johnson 1986) was higher during the 1960s–1980s than in the following decades due to higher NH4+ concentrations (Fig. 3a).

The increasing pH trend between the 1960s and 1990s thus probably resulted from the recovery of headwaters from atmospheric acidification since the mid-1980s (Kopáček et al. 2013c); the decreasing use of (NH4)2SO4 for soil fertilizing and intensifying liming of farmland since the 1970s (Kopáček et al. 2017); the decreasing pollution of surface waters with DOC and NH4+ from point sources; and the decreasing mineralization of flooded organic matter during aging of the Slapy and the Orlík reservoirs since the 1960s. Another process contributing to the pH increase during this period was an increasing proportion of NO3− assimilation. The average molar N:P ratio of the seston was higher (~ 20) than the water NH4+:TP ratio (from < 1 to 8) and stable. A substantial contribution of N2 fixing cyanobacteria to the PN budget was rather unlikely due to their low abundance (Komárková and Vyhnálek 1998). Hence, the stable sestonic N content was mostly maintained by NO3− assimilation as in other ammonium-poor waters (Wetzel 2001). The microbial NO3–N reduction consumes an equivalent amount of H+ (Reuss and Johnson 1986). With the decreasing NH4+ availability (Fig. 3a) and increasing primary production (Fig. SI-9A), the role of NO3− as an additional N source for algae (Wetzel 2001) continuously increased and contributed to the in-lake H+ consumption and increasing pH, especially in summer.

The breakpoint in the trends of Slapy water pH occurred in 1992 (Fig. 4d), despite the continuing increase in the annual pH minima (Fig. 6a). The decreasing trend in pH since 1992 thus predominately resulted from the decreasing summer pH maxima (Fig. 6b). Processes contributing to decreasing pH between the 1990s and 2019 include at least four processes: (1) The intensive liming of farmland almost ceased in the early 1990s (Fig. SI-11G), and agricultural soils in the Slapy catchment have started to re-acidify (Kopáček et al. 2017). (2) The increasing leaching of DOC from forested areas in the Slapy catchment (Hejzlar et al. 2003; Kopáček et al. 2018) has supplied more organic acids to surface waters since the early 1990s. The surface pH is usually higher than that of the original soil water, leading to higher dissociation of organic acids and H+ liberation when more acidic soil water enters less acidic surface water (Stumm and Morgan 1981). (3) The decreasing P loading of the Slapy reservoir since the mid-2000s (Fig. 4c) has reduced the P-limited primary production, resulting in a lower depletion of epilimnetic concentrations of dissolved CO2 by growing algae, and hence in lower summer pH maxima (e.g., Wetzel 2001). (4) The decreased primary production also reduced NO3− assimilation and the associated H+ consumption, which further contributed to the decreasing summer pH maxima.

A shift in the water buffering capacity was another important factor contributing to higher annual pH maxima in the 1990s, compared to the 1960s and 2010s (Fig. 6b). The lower pHs in the 1960s were closer to the value 6.35 at which the buffering capacity of the carbonate system in water reaches a maximum, while in the 1990s, pH values were closer to the 8.3 (Fig. 4c), at which the buffering capacity is minimum (Stumm and Morgan 1981). The ability of water to mitigate pH changes caused by H+ production/consumption thus decreased from the 1960s to the 1990s, allowing a more sensitive pH response to H+ sources and sinks. Then, the increasing water buffering capacity due to decreasing pH values contributed to the less pronounced summer pH maxima between the 1990s and 2010s.

These results show that long-term trends and changes in the seasonal variations of pH may reflect anthropogenic and climatic effects, reservoir aging, and changes in water eutrophication, not only in acid-sensitive regions (Psenner and Catalan 1994) but also in circum-neutral systems.

Conclusions

Our six-decade study on the composition of Vltava water in the Slapy reservoir showed a complex influence of various stressors on the water quality. Effects of agricultural activities and water and air pollution deteriorated the water quality during the 1960s–1980s period, followed by improvements over the following decades. The ionic water composition predominantly reflected changes in fertilizing, liming, draining, land-use, and road de-icing. The TP, NH4+, and DOC concentrations were primarily affected by municipal and industrial wastewater production and treatment, i.e., the development of industrial and treatment technologies and legislation. Reservoir aging further affected concentrations of DOC, NH4+ and H2CO3*. Climate change resulted in earlier spring snowmelt and related discharge maxima, increasing water temperatures, earlier thermal stratification of water column, and a longer period for primary production in the Slapy epilimnion. The combined effect of climate change and decreasing acidic deposition resulted in increasing DOC concentrations during the last three decades. The integrated effects of all major stressors (water pollution, reservoir aging, climate change) significantly affected long-term trends and seasonal variations of water pH even in the circum-neutral Vltava River.