Towards a multi-basin SWAT model for the migration of nutrients and pesticides to Puck Bay (Southern Baltic Sea)

Study area

The study area is located in northern Poland on the southern coast of the Baltic Sea (Fig. 1). It covers the catchment area of Puck Bay within the Puck commune district with an average population of 107 people per km². The basin area is 176 km². The use of land for agriculture dominates (60%), in addition to forests (29%) and urban areas (11%). This is a young glacial area of varied relief, characterized by a large denivelation (from −0.5 to 113.5 m a.s.l.) and cut by smaller valleys with steep slopes (Potrykus et al., 2018). The bottoms of the valleys contain the beds of rivers, e.g., Gizdepka, Błądzikowski Stream, and Płutnica. The geological structure consists of fluvioglacial sediments, mainly loamy sands and sandy loams interlaced with clay. Peat soils dominate in the valleys. There are two nature reserves in the catchment area: Beka Nature Reserve and Darżlubskie Buki Nature Reserve. Waters from the studied area flow into the Bay of Puck. It is a brackish shallow bay separated from the Baltic Sea by the Hel Peninsula. It is an area of very high biodiversity included in the NATURA 2000 network as a Special Protection Area and Special Area of Conservation (Zima, 2019). Average flows in the analyzed streams are: Gizdepka 0.178 m³ s⁻¹, Błądzikowski Stream 0.035 m³ s⁻¹, Płutnica 0.718 m³ s⁻¹ (Bogdanowicz & Cysewski, 2008). The water quality of each river varies depending on the season of the year, the intensity of vegetation and the agricultural practices (Wojciechowska et al., 2019a). The specific climate of this region, with moderate winters and mild summers, is due to the proximity of the Baltic Sea. In the period 2011–2019 the annual average temperature was about 7.5 °C and the average annual precipitation was 712 mm.

SWAT model

The catchment model was created based on a digital terrain model, soil maps and meteorological data. The digital elevation model (DEM) was created from LIDAR data supplied by the Polish Head Office of Geodesy and Cartography. The model was set up using the QGIS SWAT interface. The resolution of the used DEM is 10 m (Kalinowska et al., 2018; Mustafa & Szydłowski, 2020). The basin boundaries determined on the basis of the DEM were compared with topographic and hydrological maps and were also verified by fieldwork (Fig. 2) (Kalinowska et al., 2018; Dzierzbicka-Głowacka et al., 2019a; Dzierzbicka-Glowacka et al., 2019b).

The soil maps were supplied by the Polish Geological Institute. Due to the very fragmented “mosaic” soil, an attempt to accurately reproduce the complicated soil system needed almost 50 different soil profiles. Such a large diversity resulted in the release of over 1000 HRUs, thus the decision was made to generalize. The first step of the generalization was to leave only 2 soil layers in the soil profile. Then 13 types of soil profiles were selected, covering over 95% of the analyzed area. The remaining area was assigned to the main groups based on similarities in the share of the ratio of fraction content (silt, sand and clay fractions) and the saturated hydraulic conductivity. The soil parameters of every soil class were taken to be the average values reported in the Polish literature (Table 1) (Zawadzki & Bednarek, 1999).

The land use map was based on the soil-agricultural map from the Institute of Soil Science and Plant Cultivation, and the digital surface model from the Polish Head Office of Geodesy and Cartography, as well as fieldwork.

Meteorological data from the period 2011–2019, used in the simulations, were supplied by the Interdisciplinary Centre for Mathematical and Computational Modelling, University of Warsaw (Fig. 3). In addition, data provided by the Institute of Meteorology and Water Management for the years 2000–2010 were used. This period was treated as the time needed to start the model (warm-up period). Potential evapotranspiration was estimated with the Hargreaves method (Arnold et al., 2013).

The possible full division into HRU units was adopted based on the soil-agricultural map. The slope criterion was not included in our research because crop areas, meadows and pastures located in the flat part of the basin were the most important terrain. The biggest variation in terrain occurs in the western part covered with forests. Finally, a simulation was made consisting of 17 sub-basins and 353 HRUs.

Crop calendars and fertilizing

Information was obtained from farmers on the type and method of cultivation, dates of agrotechnical operations and doses of fertilizers used in the catchment area (Dzierzbicka-Głowacka et al., 2019a; Dzierzbicka-Glowacka et al., 2019b). Based on the above data, agricultural calendars were created for the most popular crops, containing average agricultural treatment dates and fertilizer doses (Table 2 presents an example crop calendar for winter wheat). The SWAT model allows for defining the doses and dates of fertilization or using the autofertilization mode (Neitsch et al., 2011). The autofertilization mode supplies as many nutrients as the plants need. The auto mode requires only the specification of the type of fertilizer and the start date of autofertilization. In the autofertilization mode, the fertilizers were composed of elemental nitrogen and elemental phosphorus.

In order to investigate the effects of different land use on the time variability of loads of nutrients and pesticides in the surface water, 13 scenarios were defined, based on typical agricultural practices for the investigated area. Scenario S1: rotation of main crops in the area (winter canola, winter wheat, silage corn). Scenario S2 is S1 in autofertilization mode. Scenarios S3–S13 assume only one type of crop on agricultural land without any rotation. The crop types are as follows: S3 winter wheat (also representing winter triticale), S4 winter wheat in autofertilization mode, S5 silage corn, S6 silage corn in autofertilization mode, S7 winter canola, S8 winter canola in autofertilization mode, S9 mixture of spring cereals (represented by barley), S10 mixture of spring cereals (represented by barley) in autofertilization mode, S11 potatoes, S12 potatoes in autofertilization mode, S13 peas (Pisum). For silage, dry biomass was considered as the yield. SWAT gives the yields of potatoes as dry mass, thus we assumed a 20% content of dry mass in potatoes in our calculations (Wierzbicka, 2012). Peas occur in just one scenario, without autofertilization, because Polish law limits the fertilization of this plant.

Grassland areas were assigned a permanent plant cover (fescue species, hay cut twice a year). Forests were also modeled as a permanent plant cover. Pine was selected as the representative species.

Autofertilization and doses practiced by farmers did not exceed the standards contained in the Nitrates Directive. Fertilizers available on the local market were added to the implemented fertilizer database, by specifying the necessary parameters in the SWAT model. Farmers mainly used slower acting forms of nitrogen. For natural fertilizers, the SWAT database was used.

Pesticides

The SWAT model has been used many times to model the surface runoff of pesticides. However, many of the pesticides from the SWAT database have been withdrawn from use in Poland and the European Union, e.g., DDT, atrazine (EU Pesticides database—European Commission, 2019; Barchanska et al., 2017).

The researched area was examined for the presence of pesticides in rivers (Wojciechowska et al., 2019b). Among the detected compounds, pesticides found in the SWAT database (glyphosate amine), as well as new ones not defined in the SWAT database, were determined. Their necessary parameters were adopted based on the literature (Table 3) and then entered into the database used by the model (Pesticide Properties DataBase, 2019; Sawicki & Zima, 1997; Luo & Zhang, 2009; Neitsch et al., 2011; Thomas, Mineau & Juraske, 2011; Arnold et al., 2013; Fohrer et al., 2014).

Calibration and validation

Any mathematical model requires calibration in order to reflect reality. In the case of the SWAT model, with a sufficiently large data resource, it is possible to use ready-made programs (e.g., SWAT-CUP) to automate this process. The situation is complicated by the multitude of issues analyzed and the number of parameters included in the calculations, and the related amount of data needed to verify the model. For the needs of the WaterPUCK project, the model was required to deliver data about hydrology, water quality, infiltration, and crop simulation (Dzierzbicka-Glowacka et al., 2018; Dzierzbicka-Głowacka et al., 2019a; Dzierzbicka-Glowacka et al., 2019b). Although monitoring existed at the turn of 2018/2019, in conjunction with several field studies conducted in 2017–2019, the complete data needed for comprehensive model calibration cannot be collected. This is related to the specificity of the area. The collection of data for calibration coincided with a severe drought, during which, in the summer of 2018, almost half of the watercourses dried up. The water levels on the Błądzikowski Stream were so low that most of the time it was impossible to perform hydrometric measurements. Even after heavy rain, the surface runoff was barely noticeable because the non-urbanized drainage basin made up for shortages. Moreover, large slopes, especially in the initial forest fragments of watercourses, affect the rapid response of the drainage basin and the very quick discharge of excess water to Puck Bay.

Due to the lack of sufficient data on stream discharge, it was decided to perform calibration by a manual trial-and-error procedure (Szymkiewicz et al., 2020). For calibration, meteorological data from the years 2000–2009 were used. In the calibration process, the following parameters were adjusted: SCN curve numbers (increased), maximum canopy storage (CANMX, increased in forests), and water in the shallow aquifer returning to the root zone REVAP (increased in forests), as well as the parameters describing pine growth in the “plant.dat” file. The objectives for calibration were the annual groundwater recharge from each HRU (constrained to the range of between 3% and 30% of yearly precipitation), the average annual evapotranspiration and the average annual biomass production in pine forests. A similar approach was successfully used in other studies (Cao et al., 2006; Smarzyńska & Miatkowski, 2016).

Calibration and validation

In the present study, the model calibrated for the period 2000–2009 has been used to simulate another time period (2011–2019), which can be considered as a validation of the model. In Table 4 the model outputs for the period 2011–2019 are compared with the reference values obtained from the literature. It can be seen that the agreement is good, except the groundwater recharge, for which the maximum model value significantly exceeds the reference value. This is due to higher precipitation in the years 2011–2019 (from 506 to 979 mm, average 712 mm), compared to the calibration period (average 620 mm). The increased difference between the precipitation and evapotranspiration led to increased surface runoff and groundwater recharge.

The results obtained differ slightly in comparison with the reference data. This is caused by different weather conditions than those reported in the data used for calibration, e.g., higher precipitation.

Information from the monitoring network was used to validate the model. At the mouth of the watercourse an automatic level station was located, then the water level was converted to the flow rate, based on the rating curve. The impact of the cultivation method (type of cultivated plant, type of agrotechnical treatments as well as method of fertilization and plant protection) on the environment, particularly on the quality of the waters of Puck Bay, could be determined after the full agrotechnical cycle, i.e., after a year. Therefore, the results related to the main product based on the SWAT model—the agricultural calculator, are presented in the form of the annual nutrient load. ${\displaystyle \mathrm{BIAS}=\frac{\sum _{i=1}^n\left(Q_i^{\mathrm{sim}}-Q_i^{\mathrm{obs}}\right)}{n}.}$

where: Qi^sim—the ith observation; Qi^sim—the ith simulated value; Qmean—the mean of the observed data observation; n—the total numer of observations.

Both on a monthly and daily basis, the BIAS parameter = −0.012, which was considered sufficient to compensate for the model error in a longer time unit. This is confirmed by a similar result of the average annual flow: simulated $Q_{\mathrm{avg}}^{\mathrm{sim}}=0.171$ m³ s⁻¹ and observed $Q_{\mathrm{avg}}^{\mathrm{obs}}=0.172$ m³ s⁻¹.

Water balance

The highest precipitation (979 mm) was recorded in 2017. In this wet year, the surface runoff averaged 238 mm (24% of precipitation), infiltration 293 mm (30% of precipitation), and evapotranspiration 436 mm (45% of precipitation) (Table 5). The proportions of the water balance are less affected by the type of crop and the fertilization method used. For wheat (scenarios S3 and S4), evapotranspiration was the highest and amounted to 47% for calendar fertilization and 48% for autofertilization. For each scenario, the share of the water balance components differs by a maximum of 3% from the average values.

The lowest precipitation was recorded in 2018—506 mm. In this dry year, the surface runoff averaged 76 mm (15% of precipitation), and infiltration 101 mm (20% of precipitation). Plants evapotranspired an average of 373 mm (74% of precipitation) (Table  5). The proportions of the water balance were less affected by the method of fertilization. For wheat (scenarios S3 and S4), evapotranspiration was the highest and amounted to 80% for calendar fertilization and 81% for autofertilization. However, the share in the balance of evapotranspiration in the dry year for a particular crop ranges from 70% (S13 peas) to 81% (S4 winter wheat autofertilization). Depending on the crop, evapotranspiration varies considerably.

The share of individual components of the water balance changes each year depending on the weather conditions. There are slight differences between the water balance in the same year for cultivation in the autofertilization mode, and fertilization specified in the calendar of agricultural practices (for example, wheat in Table 6). For wheat in the autofertilization mode, the share of evapotranspiration is on average 2% (maximum 3%) higher than in the calendar mode. Other components of the balance sheet differ by a maximum of 2%.

Yields

Average yields are presented in Table 7. The SWAT model calculates the yield for each HRU unit taking into account soil and hydrological conditions. In the scenarios with odd numbers, fertilization was assumed according to the implemented agricultural calendar, while in the scenarios with even numbers, the autofertilization option was used. Biomass was used as the corn yield because only silage corn is grown in the catchment area.

Nutrient doses

Comparing the total nitrogen doses, the doses given according to the crop calendar are higher for crop rotation S1, silage corn S5, winter canola S7, and potatoes S11 (Appendix A). However, in autofertilization mode, they are higher for winter wheat and spring cereals. The fertilizers used by farmers in the crop calendar are dominated by ammonium. As fertilizers in autofertilization mode, elemental nitrogen and elemental phosphorus were used. The doses of phosphorus applied according to the crop calendar are several times higher than phosphorus fertilization in the autofertilization mode. Organic forms of nitrogen and phosphorus are applied only in crops fertilized with manure by farmers, i.e., S5 silage corn and S11 potatoes.

Nutrient losses

The largest total nitrogen losses are in silage corn (S5) and potatoes (S11) according to the crop calendar. The highest nitrogen infiltration into groundwater also occurs for these scenarios. It is worth emphasizing that organic nitrogen dominates in the runoff of nutrients. Only during the cultivation of silage corn and potatoes does autumn manure fertilization occur. Losses of nitrate-nitrogen to surface waters are on average three times higher in autofertilized crops that are fertilized only with the nitrate form of nitrogen. The amount of nitrogen uptake by plants is independent of the fertilization mode (Appendix B).

According to the model, the amount of phosphorus uptake by plants is similar for both fertilization modes (Appendix C). The losses of phosphorus and its uptake by plants are similar for different crops, regardless of the fertilization method. According to the model, mainly organic phosphorus is loaded into the stream. The largest phosphorus losses are in silage corn (scenarios S5 and S6).

Gizdepka River investigation

Despite the fact that with the adopted criteria for the calibration and validation of the model, the reliable time step of the analyzed results is one year, the influence of the type of autofertilization mode on the daily and monthly results was also checked. The paper analyzes the Gizdepka River in detail (Table 8). This is due to the fact that this stream has all the typical features of rivers in the studied area. In addition, it had the best hydrological conditions (no periodic drying).

Analyzing the average monthly nutrient concentrations in the outflow of the Gizdepka, for winter wheat (scenarios S3 and S4) in the autofertilization mode (mainly nitrate-nitrogen), N-NO₃ concentrations several times higher can be seen. Nitrate values are from 4.3 to 9.1 mg l⁻¹ in the autofertilization mode, while in fertilization according to the crop calendar, from 1.3 to 3.2 mg l⁻¹. Organic nitrogen concentrations are similar (Figs. 4 and 5). Organic phosphorus losses are comparable between scenarios.

Similar to the mean values from the 2011–2019 period, the trend of a change in the proportion of the forms of nitrogen observed in the daily outflow is confirmed in 2019, with a significant increase in the concentration of nitrates. A comparison of the monthly sum of nutrient loads clearly indicates the dominance of nitrates in the autofertilization mode with the simultaneous relative compatibility of phosphorus loads (Figs. 6 and 7).

Pesticides

Not all pesticides used in the analyzed area were modelable. Table 9 presents the results available for the most popular plant protection products: Diflufenican and Glyphosate. Dose refers to the dose of the substance applied to the field, according to the agricultural calendar. The proportion of the dose that was able to react is about 38%. Most of the dose, ranging from 60% to 75%, is destroyed (broken down) before it starts affecting the environment. These losses are caused, among others, by atmospheric conditions (e.g., wind) (Neitsch et al., 2011). There was no significant effect of autofertilization on the reaction of pesticides. In the model, no pesticides were found at the outflow from watercourses, which is consistent with the results of field studies (Pazikowska-Sapota et al., 2020). Only trace amounts of the substance remain on the plant surface (28–38 mg ha⁻¹). A large part of the pesticide remains on the ground, amounting to 3% of the substance dose in the case of glyphosate and 85% of diflufenican.

Water balance

The water balance and the proportions between its components are adequate to the amount of precipitation. In years where rainfall is lower than the average rainfall, the share of evapotranspiration in the balance sheet is much greater. The average rainfall in the studied period was 712 mm. In the analyzed period, an average of 420 mm evapotranspired (60% of the annual average precipitation). An average of 143 mm reached the surface runoff (which is 20% of the average precipitation) and an average of 154 mm infiltrated (21% of the average precipitation—Table 10). By comparison, for another basin in Poland for the years 1963–2010, the average rainfall is 611 mm, while the surface runoff is 106 mm (17% of rainfall) (Banasik & Hejduk, 2012).

Yields

The average yields for spring crops, represented by barley, are higher in the autofertilization mode scenario (scenario S10) by an average of 2300 kg than when fertilization was based on the agricultural calendar (scenario S9) (Table 7). This difference is most likely due to the type of fertilizer and rapidity of reaction. For the remaining crops, the impact of the fertilization mode is much smaller and without a clear advantage of one fertilization method. The obtained yields are close to the values given in the survey (Dzierzbicka-Głowacka et al., 2019a; Dzierzbicka-Glowacka et al., 2019b).

Nutrient doses

The greatest difference between fertilization in auto mode and the crop calendar is for spring cereals. In the S9 and S10 scenarios, doses in autofertilization with nitrogen are twice as high. None of the scenarios exceed the maximum doses specified in the Nitrates Directive (Appendix A) (“Council Directive 91/676/EEC of 12 1991 concerning the protection of waters against pollution caused by nitrates from agricultural sources”). In the analyzed area, the use of nitrogenous fertilizers was higher than the average consumption for the whole of Poland and for the Pomeranian Voivodeship. The mean phosphorus fertilizer consumption was higher than in the Pomeranian Voivodeship, but lower compared to the entire country (Dzierzbicka-Głowacka et al., 2019a; Dzierzbicka-Glowacka et al., 2019b).

Nutrient losses

The biggest difference in nitrogen uptake by spring cereals (S9 and S10) may be due to the fact that nitrate-nitrogen from autofertilization is more easily and faster available for spring cereals. Also, the yields of spring cereals (S9 and S10) in the autofertilization scenario are much higher (Table 7). This confirms that nitrate-nitrogen works faster and is better absorbed even during spring water shortages. The diagrams (Figs. 6 and 7) show a higher ratio of nitrates in the surface runoff in the autofertilization scenario, which indicates greater leaching of nitrates from the catchment.

Gizdepka River investigation

Analyzing nutrient loads, it is plain to see increased values of organic forms (Table 8). This may be related to the fact that the catchment has an agricultural character, arable lands are located on the slopes of moraine hills and in valleys characteristic of a young glacial landscape. These valleys were covered with peat bogs and drained by watercourses. Currently, they are transformed into agricultural areas, becoming a reservoir of soil with a high content of organic carbon (nearly 20% of the catchment area) (Table 1). For the SWAT model, organic nitrogen levels are assigned assuming that the C:N ratio for humic materials is 14:1. Organic phosphorus levels are assigned assuming that the N:P ratio for humic materials is 8:1 (Neitsch et al., 2011). Dependence on the reaction between organic and mineral forms of nutrients as a function of carbon content may be too simplistic (Kemanian et al., 2011). This may explain the high proportion of organic nutrient forms in the runoff.

For the needs of the model, a generalization was made, consisting of the unification of crops in the monoculture in the entire area, classified as farmland. In reality, agricultural land is a mosaic of many crops, including those not included in the analyzed scenarios. A large range of nutrient outflow variability, depending on the adopted scenario, from 872 to 6078 kg N km⁻² year⁻¹ and from 110 to 939 kg P km⁻² year⁻¹ (Table 8), leads to the assumption that the percentage distribution of crops in the catchment area may also have a large effect on the amount of nutrient outflow.

Pesticides

As part of the fieldwork related to the WaterPUCK project, in 2018 the catchment and adjacent bay were tested for the presence of pesticides. Samples were tested for the presence of 309 substances of different generations, both those currently used and discontinued. Samples of soil, groundwater, water from drainage ditches, watercourses flowing into the bay and water from the bay itself were taken. Only in one sample (Błądzikowski Stream, July 2018), was diflufenican detected with a concentration of 0.13 µg dm⁻³. Substances found in plant protection products were also detected in drainage ditches, in just 17 samples. The majority of them were glyphosate and its derivatives (AMPA), with a maximum concentration of 20 µg dm⁻³ (Gizdepka basin, August 2018), as well as metazachlor with a maximum concentration of 2 µg dm⁻³ (Błądzikowski Stream basin, August 2018). The concentration of other substances was below the limit of quantification (Pazikowska-Sapota et al., 2020).