Estimation of nutrient loads with the use of mass-balance and modelling approaches on the Wełna River catchment example (central Poland)
Case study areaThe studied catchment (2 621 km2) is located in the central-western part of Poland, and constitutes a part of the Oder River basin. The Wełna River (118 km) discharges to the Warta River near the town of Oborniki18, with an average flow rate of 8.1 m3s−1 (1980–2019) in this profile19. The natural conditions in this catchment favour the development of intensive agriculture, which covers almost 72% of this area (1888 km2) and contributes to the high consumption of mineral fertilizers20. Forest areas cover another 22% of this catchment (589 km2), while urbanised ones only 4% (93 km2) (Fig. 1). The Wełna River catchment is inhabited by approx. 230,000 people, of which only approx. 74% is served by wastewater treatment facilities21.Figure 1Localisation of the Wełna River catchment with its land use forms and nutrient sources. This figure was created using ArcGIS 10.2.1 for Desktop available at https://www.esri.com/en-us/home. Licence granted to Institute of Meteorology and Water Management.Full size imageInput dataBoth the mass-balance method and the modelling method require a similar amount and type of input data (Supplementary Table S1). Basic information on the Wełna River daily flow rates and nutrient concentrations in the closing profile of the catchment (Oborniki) has been obtained from the state monitoring services (Institute of Meteorology and Water Management—National Research Institute—IMGW-PIB13 and State Environmental Monitoring22—SEM) (Supplementary Table S1). They have formed the basis for the estimation of the share of individual sources in the mass-balance method, as well as for the calibration of the Macromodel DNS/SWAT in the modelling method. Other data, such as maps of elevation, river network and soil maps, as well as meteorological data, necessary for the development of an accurate representation of the studied catchment area on the Macromodel DNS/SWAT digital platform, were also obtained from state repositories. Data on the land use comes from the Corine Land Cover8, while detailed information on nutrient sources has been obtained mostly from the Local Data Bank of statistical information. The utilisation of the collected database has been presented in Fig. 2, and described in the following text. The comparison of the results for nutrient loads from both method was based on the year 2017, which was characterised by the maximum amount of monitoring data for both flows (365 measurements) (IMGW-PIB) and total nitrogen (TN) and total phosphorus (TP) (12 measurements–SEM). The average air temperature in 2017 in Poland was 1.5 °C higher than the long-term average (1971–2000) and was over 10 °C which resulted from the warm autumn and the end of the year. The time of the snow cover presence was shorter than the long-term data, and the rest of the year was classified as thermally normal.Figure 2Methodology diagram with relevant chapters marked in grey ovals (green—steps and data used for Mass Balance method, blue—steps data used for Modelling method, green/blue—steps and data used for both methods).Full size imageIn terms of precipitation, 2017 was assessed as wet, similarly due to rainy autumn and summer. In the Wełna River catchment area, the annual rainfall was about 770 mm, however high variability of precipitation conditions in particular months, from extremely wet to very dry, should be noted23. Therefore hydrologically, 2017 was considered normal with the flows only slightly lower than the long-term average.Mass-balance methodThe first method used for the quantification of sources and loads in the studied area was the static mass-balance method. It is widely used by the Polish administration authorities responsible for water management17. This method is based on the assumption that the sum of the nutrient loads in the river’s closing profile (selected based on access to the monitoring data) and its retention in the catchment equals the emission of nutrients in a given time. Such assumption allows the apportioning of the river loads among identified sources and the estimation of their contribution to the total loads, based on known or assumed values of their retention.River load calculationThe total load of nutrients discharged from the catchment was calculated using the daily flow rate and nutrient concentrations in the closing profile of the catchment area (Oborniki, Fig. 1) from the SEM (Supplementary Table S1). The daily river load was calculated using the following Eq. 5:$${L}_{river}=0.0864sum_{t=1}^{n=365}{({Q}_{t}cdot {C}_{t})}_{t}$$
(1)
where: Lriver is the annual load [kga-1], n is the number of days, t is the consecutive day, Ct is the concentration [mg L-1], Qt is the mean daily flow rate [Ls-1], and 0.0864 is the unit conversion.Due to the fact that the flow rate is measured daily and nutrient concentrations only 12 times a year, the linear interpolation method5 was used to obtain the daily concentration values:$${x}_{k}={x}_{a}+kcdot frac{{x}_{b}-{x}_{a}}{n+1}$$
(2)
where: xk is the interpolated concentration value, xa is the first of the two measured concentration values between which the concentrations are interpolated, xb is the second of the two measured concentration values between which the concentrations are interpolated, k is the consecutive number of missing value and n is number of missing values.Source apportionmentFor the mass-balance method, data on nutrient loads for source apportionment (emission inventory) was collected in order to proceed with further calculations. The calculations were performed for 2017, due to the availability of river monitoring data and the nutrient sources were divided into 7 categories, based on the HELCOM guidelines5: municipal (MWS) and industrial (INS) point sources, municipal diffuse sources (SCS), forestry (MFS), agriculture (AGS), natural background (NBS) and atmospheric deposition (ATS). The category of “unknown sources” (UKS) was taken into account, in order to include possible discrepancies in nutrients load apportionment, and to cover eventual differences between calculated river load and inventoried emission.The MWS loads were calculated on the basis of the number of inhabitants served by the wastewater treatment plants (WWTPs)21. In the Wełna River catchment, 151 771 inhabitants were served by the 12 WWTPs covered by the National Wastewater Treatment Program (NWTP)24, which provides information on the total discharge volume from each facility. For 5 of these plants, information on influent/effluent nutrient concentrations was also available, allowing the direct calculation of discharged loads. For the remaining seven facilities, the loads were calculated on the basis of the mean influent concentration information, available for the WWTPs covered by the NWTP (80 mgL−1 and 12 mgL−1 for TN and TP, respectively), and approximated nutrient reduction level in non-biological WWTPs. This reduction level, based on data from the NWTP, was set at 65% for TN and 35% for TP24. Another 19 350 inhabitants of this catchment were connected to the small WWTPs, not included in the NWTP. This part of the MWS load was calculated using the mean daily wastewater outflow (0.12 m3day−1 per person), the same mean nutrient concentrations and reduction levels as presented above. Additionally, the remaining 25% of the catchment’s population (58,000) is not connected to any WWTP and uses septic tanks and other types of individual wastewater treatment systems. The load from this source was expressed as SCS, and calculated using unit loads set on 11 gday−1 per person and 1.6 gday-1 per person for TN and TP, respectively17. The industrial nutrient input information (INS) was gathered directly from the Statistics Poland office database21.The AGS load was calculated using nitrate and phosphate concentrations in shallow groundwater (90 cm below the ground surface), from 22 sampling points located on agricultural areas in the Wełna River catchment17. Concentrations were recalculated to TN and TP respectively and averaged. Thus, the obtained mean concentrations were 8.25 mgL−1 of TN and 1.92 mgL−1 of TP. Subsequently, load values were calculated by multiplying these concentrations by the outflow from agricultural areas, calculated as a share of the total catchment outflow respective to the agricultural use of the area. The calculated loads were multiplied by coefficients reflecting the share of monitored outflow (groundwaters and tile drainage) from the agricultural runoff (1.11 for TN and 4.17 for TP)25. Subsequently, the natural background (NBS) was subtracted from the AGS load.The load corresponding to NBS was calculated using the total catchment outflow and nutrient concentration values reflecting conditions in undisturbed areas of pre-human activity, set as 0.15 mgL−1 and 0.02 mgL−1, for TN and TP respectively17. The MFS load was also calculated in a similar way, using nutrient concentrations set to represent forest catchment as 0.31 mgL−1 and 0.038 mgL−117 and the outflow calculated as the share of the total catchment area, respective for the catchment part covered by forest. Also in this case, the NBS load has been subtracted. As for the ATS load, data on pollutant deposition into the ground from precipitation was taken from the SEM network26. This data was based on precipitation and its chemistry measurements taken from 22 monitoring stations covering the entire territory of Poland. The total load from the point and diffuse sources was calculated by adding the loads mentioned above. The eventual difference between the river load (“River load calculation” Section) and inventoried emission (“Mass-balance method” Section) accounted for the other sources (UKS).Load apportionmentThe contribution of each source to the calculated river load was calculated based on a simplified equation modified from HELCOM5:$${L}_{river}=DP+LOD-RP-RD$$
(3)
where: Lriver is the river load [kga−1], DP is the load from point sources (MWS and INS) [kga−1], LOD is the load from diffuse sources (SCS, ATS, MFS, AGS and, NBS) [kga−1], RP is the point source retention [kga−1] and RD is the diffuse source retention [kga−1].In the adopted mass-balance method, it is assumed that nutrient load from the point sources (DP) is introduced directly into the river bed phase, while load from the diffuse sources (LOD) is discharged into both phases of the catchment, land and river bed ones. In both phases, self-purification processes are taking place, resulting in the reduction of nutrient loads on the way from the source to the catchment closing profile. However, due to the limited amount of data, the self-purification processes in the river have been omitted, therefore the point source retention (RP) equalled 0 kga−1. Subsequently, the diffuse source retention (RD) has been estimated on the basis of the difference between each nutrient load of the river (Lriver) and the point sources (LOD). The remaining river load has been then attributed proportionally to the contribution of the particular diffuse sources to the total source apportionment (emission inventory).Modelling methodThe digital platform, the Macromodel DNS with the SWAT module27,28,29,30,31,32, was used for comparison for the nutrient balancing method described in “Mass-balance method” Section. This advanced dynamic tool tracks nitrogen and phosphorus migration paths in the river basin taking into account their spatial and temporal variability. For this purpose, it takes into account a very extensive input database, similar to that used in the mass balance method (Supplementary Table S1). Natural and anthropogenic processes that affect the transport and transformation of nutrients, are also part of this platform. The SWAT module (version 2012) is a tool which operates in the spatial information system (GIS) and is fully integrated with it. Using the digital elevation model (DEM), the SWAT module divided the entire analysed Wełna River catchment into a total of 225 sub-catchments of an average area of 11.5 km2. The subsequent use of data on land use (forests, agriculture and urbanised areas) and the types of soils (31 classes) allowed the authors to identify a total of 2824 hydrological response units (HRUs), homogeneous in terms of vegetation, soil and topography33. Afterwards, a simulation of soil water content, evapotranspiration, surface runoff, primary and underground flows was carried out in accordance with the water balance Eq. (4), which represents the basis for the quantitative component and the HRU.$${SW}_{t}={SW}_{0}+sum_{i=1}^{t}({R}_{day}-{Q}_{surf}-{E}_{a}-{W}_{seep}-{Q}_{gw})$$
(4)
where: SWt is the final soil water content (mm H2O), SW0 is the initial soil water content (mm H2O), Rday is the amount of precipitation (mm H2O), Qsurf is amount of surface runoff (mm H2O), Ea is the amount of evapotranspiration (mm H2O), Wseep is the amount of water entering the vadose zone from the soil profile (mm H2O), Qgw is the amount of return flow (mm H2O).The model directs all runoff flows generated by each HRU through the channel network, thus simulating a catchment. The water balance equation also represents a basis for the simulation, transport and transformation of nutrients required for the quantitative component of the model. This tool models forms of nitrogen, organic and inorganic , different forms of phosphorus in soil34, as well as organic nitrogen and phosphorus forms associated with plant residues, microbial biomass and soil humus35,36,37,38. Final results of simulations are produced by the SWAT model as all the forms of nitrogen and phosphorus (in kilograms of N and P per a time unit, respectively) are then summed up to give TN and TP values. To verify that the model properly predicts TN and TP values its results are calibrated with the TN and TP values resulting from SEM, as described in Sect. 2.4.1. Moreover, the particular forms of nitrogen and phosphorus have also been compared with the modelling results (Supplementary Table S4). A detailed overview of the migration and transformation pathways of nitrogen and phosphorus forms in the catchment has been presented in the Supplementary Information (Sect. S1), while mathematical description of these processes is included in the theoretical documentation of the SWAT model39.Similarly, as in the case of the mass-balance method, diffuse sources of nutrients from agriculture (AGS), forestry (MFS) or urban areas (URB) in SWAT were simulated in the land phase of the catchment. In the land phase, the model simulates both the infiltration of nutrients into the soil (fertilization, plant biomass, precipitation) and their removal from it (volatilization, denitrification, erosion, surface runoff). Additionally, changes in the distribution of nutrients in the soil (uptake by plants) and the low mobility of phosphorus itself are also taken into account39,40,41.Pollutants from municipal and industrial point sources (MWS, INS) are introduced directly into the river bed phase. The exception here is the nutrient load from municipal diffuse sources (SCS) which, reduced as a result of the self-purification processes taking place in the land phase, is also treated in the model as point sources. The SCS nutrient load mainly derives from leaking or illegally emptied septic tanks. For this purpose, septic tanks have been divided into three types: leaky, partially illegally emptied, and sealed septic tanks, legally emptied. The loads from the legally emptied tanks are included in the effluents from WWTPs reported in the catchment. While for the remaining types, their loads are calculated using factors depending on their effectiveness in nutrient removal (15 – 50%). The final nutrient load derived from these types of facilities is then calculated, taking into consideration the number of inhabitants using the different types of septic tanks and the average chemical composition of wastewater21.The load of nutrients from the atmospheric deposition (ATS) affects both land and river phases due to the presence of two deposition mechanisms in the SWAT module, i.e., wet and dry deposition. The model also allows for the determination of nutrient loads generated as a result of natural processes of nitrogen and phosphorus transformation and transport in the soil, with the omission of all anthropogenic pressure—natural background (NBS).Calibration, verification and validationThe SWAT module for the Wełna River has been calibrated, verified and validated using the SWAT-CUP software42. For the quantitative component (water circulation in the catchment), the implemented daily flow data (source: IMGW-PIB) for the period of 18 years (2001–2018) came from the water gauge stations on the Wełna River (Pruśce and Kowanówko) and its tributary (the Flinta River-Ryczywół) (Fig. 1). The qualitative component (nitrogen and phosphorus concentration in the catchment) was gathered from the SEM stations localised at the Wełna River (Oborniki and Rogoźno) (Fig. 1) and covered a period of 13 years (2005–2018). Three statistical measures, coefficient of determination (R2)43, percent bias (PBIAS)44, and Kling-Gupta efficiency (KGE)45, have been used to indicate the Wełna River model performance (Supplementary Tables S2 and S3). In terms of the quantitative component, the calibration and verification coefficients R2, KGE and PBIAS classified the model performance generally as good and very good for the main river (Wełna), and satisfactory and good for its tributary (Flinta). During the validation procedure, all coefficient values rated the model performance for daily flow simulations as very good. In terms of qualitative components, the model performance for TN simulations can be considered as very good or good, according to the all-applied coefficients. Lower model performance, mostly satisfactory, was observed for TP mainly due to the variability of phosphorus temporal distribution patterns (Supplementary Table S2). The entire process was described in detail in Orlińska-Woźniak et al46.Variant scenariosIn order to determine the contribution of individual sources to the total load of nutrients in the profile closing the analysed catchment, a final simulation of the model was used and subjected to calibration, verification and validation procedures, and called the baseline scenario (A0). Subsequent so-called variant scenarios (A1–A5), i.e. model simulations, were developed. In variant scenarios the values of selected parameters were changed in relation to the A0 scenario. This was used both in the river bed phase for point sources (A1) and for individual diffuse sources (A2–5), thus imitating surface changes for particular types of land use in the land phase of the catchment (Fig. 3).Figure 3Variant analysis diagram for assessment of nutrient loads (L) for particular modelling scenarios and sources described in the text: MWS, INS, SCS—point sources, URB—urban, AGS—agricultural, MFS—forest.Full size imageIn the A1 scenario, all parameterized and aggregated point sources (MWS, INS, SCS) for each relevant sub-basin (LMWS,INS,SCS), were removed from the model to calculate their contribution to the total nutrient loads in the closing profile of the studied catchment (LA1).In the next two scenarios (A2 and A3), concerning urban and agricultural land use, their surface areas (5 663 ha and 192 917 ha, respectively) were successively replaced by the forest land use. This procedure was based on the assumption that the forest is the primary type of land use for this catchment area47. In order to completely eliminate the influence of these areas, the nutrient loads from the relevant surface area occupied by forest land use were subtracted, in order to estimate the contribution of urban and agricultural land (LURB and LAGS, respectively).The change in land use from urbanised and agricultural areas to forest areas increased their percentage of the catchment area to almost 100%, thus the original image of the catchment area and the nutrient load at its mouth. On this basis, in scenario A4, the nutrient load from forests LMFS, which currently occupy only 20% of the catchment area (A0), flowing to the closing profile, was calculated from the proportion.The A5 scenario is the difference between the nutrient load from the A0 scenario and the sum of the remaining loads from the subsequent variant scenarios (A1–A4). In this way, both the natural background (NBS) and atmospheric deposition (ATS) were taken into account. More
