Analysis of Nitrate Pollution Pathways on a Vulnerable Agricultural Plain in Slovenia: Taking the Local Approach to Balance Ecosystem Services of Food and Water

: Groundwater pollution with nitrate of agricultural origin is a major problem in many countries. A great deal of e ﬀ ort is focused on ﬁnding ways to reduce leaching from agricultural land. In this study, di ﬀ erent land management scenarios were evaluated with the SWAT model in order to determine which are the most e ﬀ ective in reducing nitrate leaching on speciﬁc soil types in the Krška kotlina alluvial plain (Slovenia). The area is very important both for agriculture production and drinking water resources. The model was calibrated for three soil moisture ﬁeld trial sites, each representing one major soil type of the area. Simulated soil moisture values were in good agreement with the observed values (PBIAS (percent bias) ± 25%). Of the nine land management scenarios that were evaluated, vegetable rotation caused the most nitrate leaching on all soil types, but it fared better on Cambisol than on Fluvisol. Orchards on the other hand leached the least amount of nitrate, but also fared better on Cambisol. Presented studies should be considered as a preliminary stage in the study of nitrate pollution in the investigated area. Results show that nitrate leaching varies for di ﬀ erent land management scenarios on di ﬀ erent soil types. Further work should concentrate on ﬁeld trials to evaluate the impacts of reduced fertilization on nitrate leaching and both crop yield and quality on di ﬀ erent soil types.


Introduction
Nitrate pollution of groundwater, despite many regulations and measures applied in previous decades, is still a concern in many places around the globe. Areas with significant nitrate leaching oftentimes provide ecosystem services of both water and food. This is certainly the case for most of the vulnerable areas in Slovenia, where shallow aquifers that lie underneath the best agricultural lands are the main source of drinking water [1]. To combat that, the EU Nitrates Directive [2] requires a system of water protection zones to be implemented in order to protect the areas around pumping wells. With fertilization being the main source of diffuse nitrate pollution through leaching losses in agricultural areas, water protection strategies are usually based on reduction or modification of fertilizing. Despite the water protection regime that was introduced in Slovenia more than a decade ago, groundwater in several aquifers in Slovenia is still not in a 'good state' [3], according to the official monitoring reports [4].
Research in nitrate leaching field has been going on for decades and countless solutions were proposed [5][6][7][8][9], but all this just to prove that a one-size-fits-all solution does not exist. Many studies   Three field trials were set-up in the central area of the plain, each on a different soil type (Figure 1b-d)-Žadovinek (apple orchard, organic fertilizers, Calcaric Fluvisol-shallow-medium depth), Brege E (vegetable rotation, mixed fertilization, Calcaric Fluvisol-deep) and Brege SW (arable rotation, mixed fertilization, Eutric Cambisol). Soil properties for the three trial locations are laid out in Table 1. The three management regimes found on trial sites are quite typical for the area. With a large portion of farms being family owned, a mixed production structure is common (arable rotation in Brege SW is a good example). A bigger apple production company is present in the area, and the orchard in Žadovinek illustrates their typical management operation. Vegetable production field was also included in the analysis (Brege E) because vegetable production is present in the area, and in order for our country to be more self-sufficient with vegetables, it is important to increase production. Management operations for the three sites (crops grown and fertilization) were replicated based on the information we obtained from the owners, and are presented in Table 2. Irrigation is only used in the area to a limited extent (mostly in the orchards and vegetable production), but the government is evaluating possibilities to build a larger irrigation system in the future. In the case of the three sites only the orchard was irrigated (drip irrigation with shallow groundwater, average irrigation amount during summer droughts was approximately 6 mm every two days).

Modeling Approach
In order to quantify differences between scenarios, the widely used Soil and Water Assessment Tool (SWAT [26,27]) model was used. The model is considered very suitable for analyzing nitrate leaching and evaluating impacts of agricultural activities. It is a semi-distributed process-based model that operates on a spatial basis of sub basins and hydrologic response units (HRUs) in daily time interval. Input data required include topography, soil and land use characteristics, climate and land management. Detailed explanation of the SWAT model strengths and weaknesses is given in a previous publication [28]. The model is based on hundreds of parameters, many of which should resemble the modelled environment as closely as possible, to ensure a representative simulation. Local data needed for model functioning include information on land elevation, land use, soil properties, management operations and weather. Local data sources and their resolutions are presented in Table 3, along with some of the important parameters for each group. Soil hydraulic properties, weather and land management data have a great influence on leaching, thus, great care was taken to acquire quality local data (field sampling, measurements, on site observations). The great majority of parameters was not changed in the process of the model setup, as default values were proven suitable for similar areas in previous studies [15,28]. The SWAT model was utilized over field scale models because it is much more practical for evaluation of processes on a larger area, which we are planning for the future and for which a field scale model would not be optimal [15]. In Slovenia, most data are accessible publicly upon request from different governmental agencies, but to improve the resolution of the soil map and land management data, additional field work was executed. In the available soil map, the main area of the plain is only presented in three soil types. Additional sampling with soil auger was performed to determine the location of a representative soil samples, and then a soil profile was dug for each of them to collect the necessary data for SWAT model soil input data. Soil profiles were also dug for each of the field trial sites (Table 1). Land management scenarios were determined based on information of local farmer's management strategies. Irrigation was applied in the orchard for the model calibration and validation periods, but at this point, further evaluation and comparison of scenarios was done without irrigation, to exclude its effect on nitrate leaching and make the scenarios more comparable by excluding the influence of precipitation.
Based on current practices, a basic scenario was formulated. Water protection limitations based on the Nitrate directive (mostly based on limitations of fertilization) were taken into account where applicable. Of the three field trial locations, only the orchard in Žadovinek is protected as a zone 2 water protection area. Water protection regimes are described in more detail in a previous publication [29]. The model was run from beginning of the year 2010 to the end of 2018, with a 5-year warm up period to ensure a good representation of processes.

Calibration and Evaluation of the Model
Calibration of the SWAT model was performed on a daily time step. It was not performed in the way it is usually done, based on two factors. First, the study area is not a catchment itself; therefore, the streamflow would not give us any useful feedback on the model capabilities. Second, soil moisture data were the only measurements available to us for the three study locations. As the leaching of nitrate through soil is dependent on movement of soil water, a soil moisture calibration was implemented. Calibration of the model with satellite sensed soil moisture data has been getting more attention in recent years [30], but the measured data comes in quite rough resolutions; thus, we used data from time-domain reflectometry (TDR) probes (several probes in each trial location), similar to what was described before [15]. Probes were installed at 30 cm depth in the three field trial locations described earlier. Each trial was set up on different soil type, thus, all three main soil types of the area were characterized in the soil moisture data samples. Data from these locations was converted from volumetric % to mm of water in the designated soil layer in order to make measured data comparable to simulated values for corresponding layer depths from the output.swr SWAT output file. For conversion, the Equation (1) was applied. It was defined by the following reasoning: • SWAT model only outputs the soil water (SWR) amounts that exceed the wilting point (WP) value of soil water (hence we subtract it) [31], and does not account for capillary rise from deep layers (over 60 cm), • To determine the volumetric amount of water in a certain depth of soil, the percentage must be multiplied by the depth (h) of that soil layer.
The calibration itself was done on an HRU level (HRU in which the TDR measurements were performed) in SWAT-CUP software, with the SUFI-2 method [32]. To enable calibration of soil water parameters in SWAT-CUP, additional soil water package for reading output.swr file from SWAT was provided by the SWAT-CUP developer upon request. The program was then prepared and run via a command line interface, because the graphical user interface is not adapted for the use of the soil water package. When calibrating, major focus was given to soil and evaporation parameters, which are the most influential on soil moisture. Parameter ranges are listed in Table 4. Based on soil hydraulic properties being very influential on the value of soil moisture, we (since the trial locations were located on each main soil type) assumed that with using three "point" measurements of soil moisture, we got a reasonable approximation of processes in each trial site. Calibration was performed for the year 2018; for the orchard for the period from 7 June till 7 August; for the arable field from 21 June till 15 July; and for the vegetable field from 21 June till 22 August. Different periods were chosen for each study site because the TDR probes failed to provide a stable series of measurements throughout the period. With the use of delicate equipment such as TDR probes, measurement errors are definitely an important source of model uncertainty. For evaluation of model performance, we used the PBIAS (percent bias) coefficient, as recommended by Moriasi et al. [33]. The same data source was used for validation-the periods for each trial location were: 7 June till 7 August (orchard), 16 July till 31 July (arable) and 23 August till 15 October (vegetable). Moreover, a visual comparison of soil water results to measured daily rainfall for the area was also performed, to ensure the peaks in soil water corresponded to its cause-the precipitation. Additional model evaluation was done by comparing model results with local estimations for evapotranspiration, fertilization and nitrate leaching.

Scenario Analysis
To evaluate options for nitrogen leaching mitigation based on soil properties, a scenario analysis was performed. Each field trial location features a different land use (an orchard, vegetable field and an arable field), thus, for the basic scenarios the current situation was evaluated. To determine the influence of soil properties on leaching, the three original land uses from trial sites were interchanged, so that we tested each of them on each soil type, thus gaining six alternative scenarios (nine in total). This was done in order to see if any of the existing land management regimes is better suited for one specific soil type than the other Land management regimes (basic scenario) for each of the three trial sites, which is described in Table 2.

Description of Scenarios
Each of the management scenarios was tested on all three soil types, for example orchard land use was simulated on shallower Fluvisol in Žadovinek, on deeper Fluvisol in Brege E and on Cambisol in Brege SW; same for arable and vegetable regimes. All combinations are listed in Table 5, for clarity. Table 5. Scenario combinations (land use + soil type) used for further evaluation (land management is described in Table 2).

Soil Type
Land Use

Model Evaluation
Soil water simulation results from SWAT model are reasonable, both when comparing it to precipitation measurements ( Figure 2) and to measured soil water data from the sites (Figure 3). Comparison with precipitation is useful because we can visually see how soil water in the model responds to rain. Drought in the year 2017 is visible in all three sites, and it is also noticeable that soil water storage differs between soil types. The soil in Žadovinek has the best storage capacity, which is in line with measured soil properties such as texture and depth (Table 1). In Figure 3 we can observe an increase of soil water storage after rain events, and a decrease during a drought in the second part of September. August drought was successfully prevented by irrigation in Žadovinek, but in Brege E the soil water levels dropped significantly. The irrigation was not applied in further scenario analysis on purpose, to eliminate the differences in precipitation amounts, which could give a false comparison of nitrogen leaching between the sites. It would be interesting to study the use of irrigation as a nitrate leaching mitigation measure, as nitrate leaching usually occurs during strong rainfall after a period of drought [21]. As plants are not able to use as much nitrogen in periods of drought, it accumulates, only to be washed out of the soil profile with later rain events. Slow and uniform water flow from irrigation, though, might decrease such extremities, but more research should be done on this in the future to verify it.
The PBIAS coefficient was under ±25% for all sites both in calibration and validation periods: for Žadovinek 8.6% (calibration) and 8.2% (validation); for Brege SW 2.6% (cal.) and −0.3% (val.); and for Brege E −1.0% (cal.) and 12. (val.). From what we can observe, the model performance seems adequate for our needs, both by visual evaluation and by statistics [33]. Still, it has to be noted that the calibration and validation periods were quite short, thus not capturing all of the possible weather conditions of the area. Such an example might be the validation period for Brege E, which shows simulated soil water continuously higher than measured. This might be caused by a combination of factors. Precipitation distribution might not have been the same in the two sites, and since precipitation was measured in Žadovinek, actual weather might have differed somewhat in Brege E (local storms are known to vary in intensity, even in small areas). Crop development might be another source. Crop development stages could differ between model and real conditions, which influences evapotranspiration (and consequently soil water storage). It is unclear why the simulated soil water decreases rapidly during drought in August (following measurements), but during the drought in September it decreases much slower. We believe this could be connected with development stage of the crop, as described above. Based on this we figured the bias was caused by uncertainties connected with data, but also the rather simplified conversion of measured data from % to mm.
(local storms are known to vary in intensity, even in small areas). Crop development might be another source. Crop development stages could differ between model and real conditions, which influences evapotranspiration (and consequently soil water storage). It is unclear why the simulated soil water decreases rapidly during drought in August (following measurements), but during the drought in September it decreases much slower. We believe this could be connected with development stage of the crop, as described above. Based on this we figured the bias was caused by uncertainties connected with data, but also the rather simplified conversion of measured data from % to mm.

Water and Nitrogen Balance
Soil water and nitrogen balance results are daily time step variables from the SWAT output files for each HRU. They are presented in Table 6. Simulated evapotranspiration (ET) is expectedly somewhat lower than the estimated potential ET (650-750 mm) for the area and in the range of findings from a study in a nearby catchment [15]. Daily soil water content responds well to the seasonal distribution of precipitation ( Figure 2) in terms of it increasing with rain events and decreasing during droughts. During serious drought soil water in SWAT can be shown as 0 mm, which is impossible in nature, but is not an error, since SWAT only outputs water content above

Water and Nitrogen Balance
Soil water and nitrogen balance results are daily time step variables from the SWAT output files for each HRU. They are presented in Table 6. Simulated evapotranspiration (ET) is expectedly somewhat lower than the estimated potential ET (650-750 mm) for the area and in the range of findings from a study in a nearby catchment [15]. Daily soil water content responds well to the seasonal distribution of precipitation (Figure 2) in terms of it increasing with rain events and decreasing during droughts. During serious drought soil water in SWAT can be shown as 0 mm, which is impossible in nature, but is not an error, since SWAT only outputs water content above wilting point. The soil water simulation results for the three sites show that soils in the area definitely react differently even under the same weather conditions. The soil hydraulic conductivity characteristics influence the way water moves so this is expected. The profile depths are also different, which plays an important role as well. Nitrogen balance is more closely connected to different management regimes in use, and as we can see from the top section of Table 6, nitrogen fertilizer applications (N_APP) are quite diverse in the three trial sites. Applied fertilizer amounts slightly exceed the actual fertilizer inputs shown in Table 2, because by this the model accounts for the slow release of N from organic fertilizers. Nitrogen fixation (NFIX) only occurs in scenarios with the arable management regime, because of the use of red clover as winter cover crop. Mineralization (F-MN) is also the most significant in the arable regime, because of larger quantities of cattle manure applied almost every year. Despite the biggest applications of nitrogen in the arable rotation, average yearly nitrogen leaching (NO3L) is greater in the vegetable rotation by nearly a factor of 2, either because inputs of fertilizer are not done in several rations, which causes leaching and consequently a lower plant nitrogen uptake (NUP), or because no winter cover crop was grown in this scenario. No measured data was available for comparison of results with realistic leaching amounts for the area, but simulated results were in line with findings from another modeling study done in a similar area nearby [15] and with local expert estimations. Models, despite being quite reliable for evaluating possible trends, are still not able to predict future, and in order to get a better understanding of how fertilization influences yield, field trials are still a more definite way to go. They are also very important because they show us the quality of yield, which cannot be neglected if a farmer is planning to sell his vegetables to the costumer. Here is where the balancing of ecosystem services really shows its importance: we cannot influence water quality without impacting food production in such areas, and vice versa.

Scenario Analysis
In the scenario analysis, each of the management regimes was tested on each of the soil types in the area. Table 6 shows us the results. The highest difference in average yearly nitrate leaching was more than 60 kg N/ha, depending on scenario. The highest leaching (102.5 ± 16.7) was observed in the base scenario combination of vegetable production on deeper Fluvisol, while the lowest (27.9 ± 28.9) was simulated for combination orchard on Cambisol. Of the three management regimes, vegetable production was the most prone to leaching losses, while average results for orchard were the least prone to leaching. Of the soil types Cambisol, probably because of its shallow depth, was the most susceptible to leaching, while the Fluvisols fared quite equally. The differences are small though, thus, further research (possibly field trials) should be considered. Vegetable production showed similar results in a study in one of the other alluvial plains in Slovenia [15], thus, the current practices of vegetable producers in Slovenia do not seem to be very effective for groundwater protection. This does not mean, however, that vegetable production should not be considered in these areas, because many other factors (fertile soil, proximity of irrigation water sources and market etc.) make it a highly suitable choice. As Slovenia's self-sufficiency with vegetables (only 39%), but also fresh fruits (21%) [34], is very low, limiting such production is not sustainable. Considering the low self-sufficiency with fruit in light of this study's results, more orchards could be another way to deal with this problem as a whole. At least one field trial has previously evaluated the impact of vegetable production on groundwater quality in Slovenia [21], but to better evaluate the influence of reduced fertilizer inputs on yield and quality of the vegetables, more field trials should be established. Without knowing how reduced fertilization impacts quality of produce, we cannot know if such vegetables would be marketable at all, even if we managed to reduce the nitrate leaching.
Similar conclusions can also be observed in a monthly distribution (Figure 4). The highest leaching peaks can be observed in vegetable rotation, though the Arable rotation on shallower Fluvisol also features a significant peak in July. The nitrate loads have three common peaks-in February, probably because of warmer weather, in July (probably because of fertilization) and another one in the autumn after harvest. Overall, the results seem to prove our hypotheses that with careful adjustment of management operations to each soil type nitrate leaching could be diminished. Furthermore, it might be necessary to discuss the current WPZ delineation, was not done with focus on nitrate leaching, but rather on point source contamination. Current WPZs in the area are located only in the NW, thus, only the field trial in Žadovinek is protected (under zone 2). In the future, we are planning to study the spatial extent of WPZ in the area in order to also consider the soil factors, as previously suggested by other authors [15,34,35].
Soil conditions play a major role in leaching processes, thus, good understanding of local soil conditions is crucial for the adjustment of agricultural practices to achieve the best environmental effects. It is important to note that soil conditions also reflect partially the crops that grow there; for example, the topsoil in an orchard is better structured, while an arable field can be cracked, which influences the hydraulic conductivity. Actual leaching in different scenarios might therefore vary, not only because of unavoidable measurement errors (as in TDR soil water measurements, soil sampling etc.), but even more so because of the use of modeling.
Use of modeling for scenario analysis is a common approach [36], but results should never be looked upon as definite, as a degree of error is always present when simulating complex processes. In this regard we are using nitrate loads or other results for comparative purposes only. Another possible source of uncertainty is the lack of calibration for nitrate. This was not done in this study for two reasons, the first and more important being the lack of reliable measurements. Since nitrate leaching is closely connected to soil water dynamics, the risk of having nitrate cycle completely wrong was greatly diminished, and to represent the nitrate leaching loads as close as possible, the SWAT management files were edited with attention to real crop rotations and fertilization in the area. In the process of calibration, soil properties were changed slightly, to better represent measured soil moisture. It would be interesting to see how the model performance would change if different methods of determining soil hydraulic properties were used (measured/calculated by pedotransfer functions) as performing measurements is both time and cost-consuming. For this study, we used measured data, but to our knowledge, no studies concerning this question were published to date, thus, it is not clear whether this extra step improves the model performance in the end.

Conclusions
In conclusion, this is the first study to evaluate nitrate leaching from different soil-land management combinations in this specific area. The focus of research was on discovering which of the common management regimes works best for each of the common soil types in the area. The study results suggest that designation of "recommended" management regimes for each soil type could minimize leaching, while a diverse array of crops could still be grown in the area. That is good news, as monoculture cropping can have other unwanted consequences [37], such as pests and diseases, which in turn demand larger pesticide consumption.
An extensive monitoring network from three field trials was utilized for measurements of soil moisture used for calibration of the SWAT model, which was used to estimate the nitrate leaching loads. The model was validated by comparing soil moisture, evapotranspiration, fertilization and nitrate leaching with local estimates. Average yearly nitrate leaching amounted to 34.5 ± 6.1 kg N/ha (orchard in Žadovinek), 58.7 ± 53.5 kg N/ha (arable field in Brege SW) and 102.5 ± 16.7 kg N/ha (vegetable field in Brege E), for the three trial locations, respectively (business as usual scenario).

Conclusions
In conclusion, this is the first study to evaluate nitrate leaching from different soil-land management combinations in this specific area. The focus of research was on discovering which of the common management regimes works best for each of the common soil types in the area. The study results suggest that designation of "recommended" management regimes for each soil type could minimize leaching, while a diverse array of crops could still be grown in the area. That is good news, as monoculture cropping can have other unwanted consequences [37], such as pests and diseases, which in turn demand larger pesticide consumption.
An extensive monitoring network from three field trials was utilized for measurements of soil moisture used for calibration of the SWAT model, which was used to estimate the nitrate leaching loads. The model was validated by comparing soil moisture, evapotranspiration, fertilization and nitrate leaching with local estimates. Average yearly nitrate leaching amounted to 34.5 ± 6.1 kg N/ha (orchard in Žadovinek), 58.7 ± 53.5 kg N/ha (arable field in Brege SW) and 102.5 ± 16.7 kg N/ha (vegetable field in Brege E), for the three trial locations, respectively (business as usual scenario).
Generalizations and limitations connected to the modeling approach always present some uncertainty. Results should therefore not be viewed as a prediction of the future, but they offer us insight into trends connected with nitrate leaching in the area. The SWAT model was not calibrated in its most common way by using flow measurements, but with soil moisture measurements, which still need to be studied more in order to find all of its uncertainties. The most prominent uncertainties we discovered were connected with quality or lack of measurements. Even though we have concentrated on acquiring quality local input data, we found it difficult to acquire reliable data for calibration and evaluation. This goes to show that even in the age when modeling is an increasingly popular tool for assessing environmental impacts of agriculture, field trials and in situ measurements are of great importance to fill the gaps in knowledge of local areas.
This study has reinforced the findings of other authors that soil characteristics should be taken into account when applying groundwater protection measures [11,13,15]. It has revealed that certain soil types are more suitable for the production of certain crops, while other might respond differently.
Care should therefore be taken when proposing policy and regulations connected with water protection, as slight differences in natural conditions can greatly influence water quality, but also agricultural production economy. Therefore, a focused examination of local conditions proposed by other studies [13,15] is the key. Further research should focus on studying soil characteristics and fertilization reduction possibilities in specific vulnerable areas. Further education of agricultural advisors, farmers and other stakeholders, but also the general public should, also be considered, as too much emphasis is given solely on water protection, while it should be studied hand in hand with economy and sociology of the agricultural areas.