The Italian DSLs Lake Como, Lake Garda, Lake Idro, Lake Iseo and Lake Maggiore are a group of lakes of glacial origin at the southern slopes of the Alpine chain (Figure 1; Table 1). They are vital for the Alpine hydrological balance and for the water uses of the Po River plain in Northern Italy, one of the most anthropized areas in Europe. They collect extensive water runoff through large catchments, which include one or 2 MTs, whose watersheds cover most of the overall lake ones. The catchments of the MTs originate from small streams, mainly fed by nival and glacial fusion, which flow rapidly from high altitudes to low-sloping plains, where they combine into larger rivers which flow into the DSLs (Salmaso and Mosello, 2010; Terzi et al., 2021). The high-mountain areas where the MTs originate are pristine environments, whereas the downstream ones are made up of urban, industrial and rural areas.

The MTs considered in this study are the Adda River for Lake Como, the Sarca River for Lake Garda, the Toce River and the Ticino River for Lake Maggiore, the Chiese River for Lake Idro and the Oglio River for Lake Iseo (Figure 1; Table 1). Lake Lugano was not considered in the present study, despite being usually considered part of the DSLs (Salmaso and Mosello, 2010), because all its tributaries are small-sized, the largest one Creek Vedeggio having a mean discharge of 3.6 m³ s⁻¹ in the 2004–2023 period (FOEN ² ), hence incomparable to the flowrate of the other MTs (Table 1).

As anticipated before, despite the high number of studies performed on the DSLs, direct experimental assessments of the input nutrient loads from their catchments have been scarce. A quantification of the total nutrient loads flowing into Lake Como in the 2000–2019 period was carried out by Fenocchi et al. (2023), adopting a regionalization procedure to estimate loads for minor unmonitored tributaries as well. For Lake Garda, the only study on nutrient load was undertaken by Decet and Salmaso (1997), encompassing only 14 months of sampling. Some considerations on nutrient mismanagement over the Chiese River catchment were published for chemically meromictic Lake Idro in Nizzoli et al. (2018) and Severini et al. (2023). Critical nutrient loading conditions also characterize Lake Iseo, the other chemically meromictic DSL, despite the recent efforts to reduce nutrient load inputs (Pilotti et al., 2021; Scibona et al., 2022; Valerio et al., 2022). In fact, nutrient inputs from the Oglio River over 2017–2018 were experimentally assessed by Valerio et al. (2022) determining that current loading levels are still too high. Last, studies on Lake Maggiore (Rogora et al., 2001; Rogora et al., 2003) determined that N concentrations in the lake and its tributaries are strongly associated with atmospheric deposition, this being related to the comparatively low anthropization of the catchment.

European national and regional authorities have focused on actions to reach higher water quality standards and reduce nutrient releases, conforming to the WFD ³ . Due to the monitoring and assessment requirements of the WFD itself, data series with a minimum duration of around 15 years and monthly resolution are available for nutrient concentrations in the MTs near their outlets into the DSLs, obtained from the analysis of water samples. Together with available daily discharge series in the same locations from gauging stations of the hydrological monitoring network (Supplementary Material S1), they allow studies on nutrient load quantification and evolution for the MTs.

Nutrient concentration and discharge series are available for Adda River and Oglio River from the local environmental protection agency ARPA Lombardia. The Oglio River is partially diverted 18 km upstream of its mouth to the Italsider industrial channel, which also flows into Lake Iseo with a comparable discharge to the residual one of the natural river, up to a 50 m³ s⁻¹ hydraulic capacity threshold. Both the diverted discharge and the nutrient concentrations measured in such channel were considered, performing a weighted average of the resulting loads with those of the natural Oglio River at its mouth, for an overall assessment of the Oglio catchment. Daily discharges for the industrial channel Italsider were provided by Consorzio dell’Oglio ⁴ , whereas nutrient concentrations by ARPA Lombardia (yet sampled with quarterly frequency).

For the Sarca River and Chiese River, nutrient concentrations were provided by the local environmental protection agency APPA Trento ⁵ , whereas discharges are available from the Civil Protection Department of the Autonomous Province of Trento ⁶ .

For Lake Maggiore, monthly nutrient concentration sampling on the MTs Ticino River and Toce River started in 1981 within the framework of the International Commission for the Protection of Italian-Swiss Waters (CIPAIS ⁷ ), work being performed by the Water Research Institute of the Italian National Research Council (CNR-IRSA) of Verbania, so that the available concentration series span more than 40 years. Daily discharges for the entire period were available for the Toce River from ARPA Piemonte, whereas for the Ticino River from the Swiss Federal Office for the Environment (FOEN). As both discharge and water quality of the Toce River are measured upstream of the junction with the Strona River, its largest tributary, which occurs 4 km upstream of Lake Maggiore, we also considered Strona River data for a weighted-average overall assessment of nutrient loads from the Toce catchment. Discharges and nutrient concentrations of the Strona River were obtained from ARPA Piemonte and CNR-IRSA, respectively.

The nutrient substances considered in this study are nitrate (N-NO₃), ammonium (N-NH₄) and total nitrogen (TN) for N, orthophosphate (P-PO₄) and total phosphorus (TP) for P.

For this study, the LOAD ESTimator (LOADEST ⁸ ) freeware software by the United States Geological Survey (USGS ⁹ ) was employed to estimate daily nutrient load series from daily discharges and monthly concentration samples. LOADEST has been widely applied and validated over a multitude of case studies (e.g., Chen et al., 2015; Park and Engel, 2016; Yu et al., 2023; Song et al., 2024). This software infers daily load series through a regression analysis on input data, generating a model that fits them, and employing it to turn the available daily discharge series into a load series. LOADEST reports statistical indicators of the goodness-of-fit of the regression, as well as the regression parameters themselves.

The expression for observed load computation considered by LOADEST is reported in Equation 1, where L _τ is the measured transported load, τ is the considered integration period, equal to 1 day if daily discrete loads are produced, Q and C are the measured flowrate and the nutrient concentration, respectively. L τ = ∫ 0 τ Q C d t (1)

LOADEST then applies the regression models on the series of discrete loads obtained from Equation 1, expressing them in the logarithmic form and considering the inferences of discharge, of seasonality and of temporal trends. LOADEST includes ten different fitting models, which include a variable portion of the terms on the right-hand-side of Equation 2 with increasing complexity. These span from the simplest model, which includes only the first two terms, leading to the well-known power-law dependence on discharge (Godsey et al., 2009), up to the most complicated possible regression model expressed by the full Equation 2. Ln L τ = a 0 + a 1   L n Q + a 2   L n Q 2 + a 3 ⁡ sin   2 π   d t i m e + a 4 ⁡ cos   2 π   d t i m e + a 5   d t i m e + a 6   d t i m e 2 (2)

As per Equation 2, through LOADEST dependence on discharge can be considered up to a quadratic regression, dependence on seasonality is expressed by the two terms including trigonometric functions, whereas the last two terms on the right-hand side of Equation 2 express possible temporal trends, up to a quadratic relationship. The variable dtime expresses time in a decimal form, the variable being centered on the midpoint of the considered time series.

LOADEST allows choosing among the ten different regression models, depending on the specific goal of the calibration, or it can select the best-performing model itself, according to the Akaike Information Criterion (AIC) (Mohammed et al., 2015). This last possibility was followed in this study, to obtain the most accurate possible load series. The software automatically computes the AIC value for all ten fitting models (which depends on likelihood, number of records and number of estimated parameters) and selects the one associated with the lowest AIC value. Depending on the features of the observed series, LOADEST selected different fitting models for different MTs and for different substances for the same MT, aiming at the best fit.

The Adjusted Maximum Likelihood Estimation (AMLE) method was herein selected to obtain the regression model parameters in LOADEST. The load estimation process was applied homogenously for all nutrient substances on all the MTs, obtaining daily series of nutrient loads starting from the daily discharge data series and adopting the obtained specific regression parameters.

Table 2 reports the mean and standard deviation, the maximum and the minimum annual discharges for the considered MTs, obtained from available time series. The annual flow variability which surfaces out of the statistics highlights the influence of interannual precipitation variation, resulting in an up to an order-of-magnitude difference between minimum and maximum annual values. This variability plays a fundamental role in nutrient delivery (Stackpoole et al., 2021; Yates et al., 2022).

As N and P concentrations can show either a dilution behavior (source-limited), a concentration behavior (transport-limited) or an uncertain behavior with discharge (Godsey et al., 2009; Fenocchi et al., 2023; Yu et al., 2023; Song et al., 2024), sampled concentration data series were analyzed in relation to the flowrates, to determine the specific Q-C (Flowrate–Concentration) behavior of each catchment for each considered nutrient substance (Table 3). The catchment of the Oglio River stands out as the only one in which the transport-limited concentration behavior is clearly observed from Q-C data. This is a signal of enduring nutrient pollution reaching the river waters during rainfall periods, being either excess fertilizers, resuspended sewer sediments, uncollected sewage or untreated sewage from the activation of CSOs, confirming the evaluations of Valerio et al. (2022).

An unclear behavior emerges for all nutrient substances for the Chiese River, potentially showing that nutrient pollution is still considerable on this catchment as well, as per the judgement of Nizzoli et al. (2018) and Severini et al. (2023). The other MTs instead display a clear dilution behavior, indicating a less dramatic anthropogenic nutrient pollution, with an essentially steady ordinary input being diluted by surface runoff in rainy periods.

As regards River Oglio and the industrial channel Italsider, nutrient concentration data sampled in the same day showed significant differences between the two streams (Supplementary Material S2). Such discrepancy should be attributed to different nutrient delivery to the two streams downstream of their fork, enforcing the weighted-average approach employed in this study for overall nutrient load estimation from the Oglio catchment.

The Pearson correlation coefficients (R) between the observed daily nutrient loads and those estimated with LOADEST for the MTs are reported in Table 4. All values are R ≥ 0.39, highest values with R > 0.98 resulting for N-NO₃ and TN, likely due to pollution sources and transport dynamics for nitrates being simpler (TN is primarily made up by N-NO₃ for the studied catchments). Lower correlation coefficients result for N-NH₄, to which analytical errors due to the low observed concentrations surely contribute, and for P-PO₄ and TP, which have multifaceted pollution sources and more complex transport dynamics.

Estimated loads of the considered nutrient substances for the studied MTs, cumulated over an annual basis, are reported in Figures 2–6, together with the linear trends of the series.

Dependence of nutrient loads on rainfall is evident from Figures 2–6, as revealed by the high and very low values obtained for almost all nutrient substances and MTs in 2014, which was a very rainy year (Corsini et al., 2017), and in 2022, which was an extremely dry year (Koehler et al., 2022; Avanzi et al., 2024), respectively.

Such dependence is fully disclosed in Figure 7, in which the Pearson correlation coefficients between the series of annual mean discharges and of cumulated annual nutrient loads are displayed for all considered nutrient substances and MTs. All correlation coefficients are in fact positive, nutrient loads largely depending on flood events (Moatar et al., 2017; Torres and Baronas, 2021), this being due to discharges varying by more orders of magnitudes than concentrations (Godsey et al., 2009; Moatar et al., 2017).

Correlation coefficients close to unity mostly result in Figure 7 for N-NO₃ and TN, echoing their more straightforward dependence on rainfall due to both atmospheric deposition and higher solubility from agricultural soil wash-off. A larger dispersion of correlation coefficients, with lower overall values, results for P-PO₄ and TP, likely due to a more evident source-limited behavior, which affects also N-NH₄. Similar results among nutrient substances were obtained for the correlation coefficients between mean annual discharge series and total nutrient loads delivered to Lake Como from its major and minor tributaries in Fenocchi et al. (2023).

The Toce River and the Oglio River show the highest correlation coefficients overall among the nutrient substances, highlighting a stronger dependence on discharge variations. This could be due to the higher side slope of the valley for the Toce River, favoring nutrient delivery to the stream with rainfall. For the Oglio River, instead, the cause could be a relevant amount of dispersed nutrients being stocked in the soil and sediments and being subsequently released with flood events, although this last hypothesis would need direct observations to go beyond the speculation stage.

The correlation coefficients for Adda River, Sarca River and Ticino River display a higher variability of the coefficients among the nutrient substances, still abiding to the general behavior with higher N-NO₃ and TN values and lower N-NH₄, P-PO₄ and TP values. This could be due to higher withdrawal of P by vegetation for these basins, due to the higher relevance of agriculture inside them. Finally, the R values are lower for the combined Oglio River and Italsider channel than for the individual streams, the same happening for separate Toce River and Strona River. This issue must be ascribed to the intuitively more straightforward dependence of nutrient loads on rainfall over simpler catchments.

Past and present nutrient reduction interventions have been performed in the catchments of the MTs. However, the decreasing trends visible in Figures 2–6 must be considered in light of trends in mean annual discharges, as given in Supplementary Material S1, also given the outcomes in Figure 7.

Decreasing loads result for the Adda River and the Ticino River for all nutrient substances except N-NO₃ and TN, even though with concurrent mean annual discharge decrease. For both basins, such increasing trends for N could be due to increasing agricultural pressures. Specifically, for the Adda River, the slight increase in estimated TN annual loads complies with the higher concentrations measured in recent years, which could be attributed to rising organic nitrogen from agricultural waste in the catchment.

For the Chiese River, estimated loads show a general increase for all nutrient substances in the first half of the time series, peaking on the rainy year 2014, followed by a decreasing trend. The cumulated annual loads for this MT markedly oscillate with discharge for all nutrient substances, with high correlation coefficients between them (Figure 7). This could be due to the steeper and smaller nature of the Chiese catchment compared to those of the other MTs (Table 1; Figure 1), which enforces a behavior which depends more on transport limitation (Torres et al., 2015; Fenocchi et al., 2023).

As regards the Oglio River, it should be noted first that the annual loads from the Italsider channel are in the same order of those from the natural river for all nutrient substances, being even higher overall for P-PO₄ (Figure 6). This highlights the relevance of the industrial channel in terms of nutrient delivery to Lake Iseo, due to both discharges being slightly higher than those of the natural river under ordinary conditions and to the usually higher nutrient concentrations (Supplementary Material S2). Broadly speaking, external nutrient loading to Lake Iseo is known to be still relevant, leading together with chemically induced meromixis (Ambrosetti and Barbanti, 2005) to critical water quality conditions with deep-water anoxia and bottom nutrient accumulation (Scibona et al., 2022). Efforts to tackle this situation have been undertaken only recently on the watershed (Valerio et al., 2022), and a decreasing trend can be indeed seen from Figures 2–6, even though simultaneous with the decrease of mean annual discharge (Supplementary Material S1). It should be last mentioned that the water quality sampling location for the Oglio River is placed upstream of the Pisogne wastewater treatment plant outlet, so that a minor load underestimation is furthermore expected.

The cumulated annual loads of the Sarca River are the only ones which seem to increase through the time series, especially for TP and P-PO₄ (Figures 5, 6), yet together with a rising trend in discharges (Supplementary Material S1). A longer time series and further assessments would be needed to fully state if such an increase in nutrient loading depends solely on discharges or if the pressures are indeed rising. Figure 7 reveals that the lowest correlation coefficient between annual loads and discharges results for N-NH₄ in the Sarca River. This may highlight the presence of a steady, possibly rising, pollution source, more independent than others from discharge. As a matter of fact, increasing N-NH₄ concentrations have been observed throughout the time series (Figure 3). This could also explain the comparatively low correlation coefficients obtained for TP and P-PO₄ with flowrate in Figure 7.

The annual nutrient loads of the Toce River show a slow steady reduction over the last 40 years, still echoing the behavior of mean annual discharges (Supplementary Material S1). The Toce catchment is characterized by a lower degree of anthropisation compared to those of other MTs, especially in the riparian areas, which are largely occupied by natural vegetation, which also works as a buffer to nutrient pollution. The Toce valley is also quite narrow and steep, therefore urban development is limited by morphology (Bach et al., 2003). Furthermore, industrial activities in the Toce catchment are mostly of chemical nature (Riva et al., 2010), producing little nutrient pollution to the surrounding environment.

To properly compare the nutrient loads and their evolution among the considered MTs, the statistical features and the linear trends of cumulated annual loads for all nutrient substances and of mean annual discharges were computed over the common period of data availability 2008–2022 for all the MTs (Supplementary Material S3, Figure 8). To make trends truly comparable among catchments of different sizes, they have been normalized by the mean annual values over the common period 2008–2022.

From Figure 8, it is clear that all MTs exhibited a decreasing nutrient load trend in the 2008–2022 period, except for the Sarca River, for which significantly increasing normalized trends were found for N-NH₄, P-PO₄ and TP. Mean annual discharges have been similarly decreasing over the same period for the MTs, the normalized reduction being negligible for the Sarca River, likely due to its basin being at the eastern boundary of the DSL region, hence being exposed to a slightly different meteorology. Despite this, such a result deserves further analyses and represents a reason of concern for catchment management, as the N-NH₄, P-PO₄ and TP nutrient substances for which the positive trend is present are those more closely associated to sewage pollution. On the opposite side, the most relevant decreasing normalized trends over 2008–2022 result for N-NH₄ and especially P-PO₄ and TP for the Oglio River, likely reflecting the awaited actions on wastewater collection and treatment improvement undertaken in the last decade on the catchment, as highlighted in Valerio et al. (2022). For the Oglio River, the relevance of interventions to reduce nutrient loading from the catchment is further proven by the reducing trend in normalized discharge not being the largest one among the MTs.