Background

The quality and availability of global renewable freshwater resources have been exposed to different anthropogenic burdens and non-sustainable practices, notably over-abstraction, land use/management change, and pollution [1,2,3,4]. Consequently, water has become limited, and as it is one of the most critical global resources, this problem expected to further escalate, owing to population growth and climate change [5, 6]. In the European Union (EU), ~ 1/3 of the territory already experiences a water imbalance in terms of availability/quality; this represents unbridged issues regarding the attainment of sustainable development in most Mediterranean countries [7] and ensuring increasing crop food demands (e.g., [8, 9]). The EU represents a large-scale example in which polices and approaches can be integrated within and among member states to protect and manage water resources more sustainably and aligned with modern climate change challenges. For example, one of the pillars of the EU water policy, the Water Framework Directive [10], is a comprehensive and holistic approach aimed to protect and foster the quality and quantity of hydro-resources, promote their more sustainable use, and protect them from ‘general’ pollution. However, diffuse pollution markedly impacts 90% of the EU river catchment districts, with agriculture as the primary source [11]; meanwhile, the agri-sector exerts the greatest pressure on renewable freshwaters, responsible for almost 60% of the total water use in the EU [12]. In addition, nitrate-N has been shown to be one of the most influential nutrients in the stability of the crop yield [13]; thus, N-fertilisation is often a favourable tool for ensuring food security. However, nitrates are very mobile in terrestrial ecosystems and, in the case of non-sustainable management practices (e.g., over-fertilisation) or unfavourable environmental conditions (e.g., farmlands on hillslope terrains after intensive precipitation) [8, 14], can rapidly reach and pollute particular ecological niches and compromise the health of living beings [15]. Some of the most recent reports have indicated that in many EU member states, the agri-sector (mostly intensive livestock and crop farming) represents the predominant source of nitrogen discharged into the (aquatic) environment [12, 16, 17]. In such contexts, very complex and dependent relationships in the agricultural-environmental nexus have been recognised at the EU level through the development of more ‘specific’ or ‘targeted’ agro-environmental (inter)national directives and incentives. Many of these EU policies have been applied to control and reduce the potentially negative effects of agriculture and agrochemicals on hydro-resources. Among them, the Nitrates Directive (ND) [18] is one of the oldest (almost 30 years in force) examples of EU anti-pollution targeted legislation, concerning the protection of hydro-resources against mineral-/organic-originating nitrate pollution from agroecosystems (crop/livestock production mostly) and striving to prevent the further environmental exacerbation of this issue [19]. As a directive document (i.e., a legislative act), it establishes a goal that all EU-28 member states must achieve; thus, it might be assumed that the situation regarding addressing the pollution of hydro-resources from agricultural nitrates has been controlled after almost three decades of effort. However, the situation at EU level is not yet that optimistic and sustainable, despite some well-defined milestones required of member states, including: (i) designation of nitrate-vulnerable zones (i.e., areas in which water resources are or could be nitrate-affected and eutrophicated), (ii) spatiotemporal analysis of nitrates in groundwater and surface waters, and (iii) implementation of a comprehensive water quality monitoring programme [18].

Most of the above-mentioned directives have recently resulted in the generation of many innovative methods and bioecologically/environmentally based models (without the use of data from national monitoring programmes) for the quality assessment of water resources in the EU [20]. Furthermore, the synthetisation of input variables (e.g., physicochemical and hydrological properties of river basins, agroecological conditions, land use, fertiliser application, crop/livestock production) by means of Geographical Information Systems (GIS) has been successful in predicting nitrogen emissions from different (non)point pollution sources (agriculture, urban areas, wastewater treatment plants, industrial discharges) at the catchment-scale level by using: INCA [21], SWAT [22], MONERIS [23], Generalised Linear and Boosted Regression Trees [24], Linear Mixed Effects [25], and other multi-criterial models or approaches [26, 27]. In general, methods based on (i) topology are recommended to help in pre-processing stations or sub-basins, (ii) geostatistics is recommended for specific monitoring with high spatial correlation, (iii) hydrology combined with optimization can facilitate statistical analysis to propose new candidate stations, while (iv) multivariate statistical approaches are most commonly used in designing water quality monitoring, but they fail when a new network is proposed for catchment characteristics [28]. For example, recently [26] developed a new approach with integrated geostatistics and multi-criteria evaluation to prioritise 15 water quality monitoring stations for a data-scarce region of ~ 3166 km2 in the mountainous (938–3670 m above sea level) catchment (W. Azerbaijan, Iran) with semi-arid climate and mean annual precipitation of 282 mm. Most recently, [27], using a similar methodology, identified six (out of 21) of the most optimal sampling stations for water quality monitoring in a ~ 2100 km2 tropical catchment with a mean annual precipitation of 2750 mm (Selangor, Malaysia). Nevertheless, new approaches to water quality monitoring are needed, especially for water bodies involved in complex systems [28] such as agroecosystems, as some methods that incorporate water quality variables and non-point sources are limited in predicting the effects of land use change on network sites [27, 29].

In this study, we propose a novel robust LUMNAqSoP approach for prioritising 151 candidate groundwater monitoring stations for implementation of the ND at national (not only catchment or regional) level in Croatia, covering an area of ~ 56,600 km2, from highly diverse Mediterranean over mountainous to continental (e.g., mean annual precipitation 600–3500 mm, 0–1831 m above sea level) ecosystems [30]. The geostatistical/multicriteria LUMNAqSoP methodology specifically responds to land use change and assesses input variables (Land Use and Management, net Nitrogen fertilisation, Aquifer vulnerability and Soil sensitivity to Pollution) that are considered to have the greatest agro-environmental impacts (i.e., weights, loads or risks) on the quality of national water resources in Croatia with respect to nitrates from agriculture. Croatia, as the youngest EU member state, has designated approximately 10% of its territory as a nitrate vulnerable zone [25]; however, the implementation of the ND is not yet complete, as a groundwater quality monitoring programme is not yet fully established, which motivated us to present this study. In Croatia, 111 stations have been identified within nitrate vulnerable areas as sites for groundwater quality monitoring, and potential new candidate stations are continuously being considered in order to optimise the criteria of ND [13], i.e., to determine the influence of loads from agriculture rather than from other sources of pollution. In addition, conducting water monitoring at groundwater vs. surface stations is much more technically, temporally, and financially challenging [26]; therefore, prioritising stations could expedite and improve the quality (i.e., conduct water testing at the most exposed sites) of the monitoring programme(s) (e.g., [27,28,29,30]).

Materials and methods

The study overview with description of input variables

Since Croatia is administratively organised as 21 counties and 556 municipalities, the LUMNAqSoP methodology was developed at the more spatially detailed municipality level, with 2017 as the reference year. In the approach, five input variables were recognised as dominant impacts (loads) from the agroecosystem and environmental vulnerabilities to groundwater quality with respect to nitrate pollution, as follows: (1) total agricultural land use and management (conventional/organic), (2) total amount of applied mineral nitrogen (N-min/ha) on agricultural land (excluding areas under organic farming management), (3) total amount of applied organic nitrogen (N-org/ha) on agricultural land, (4) natural vulnerability of aquifers to pollution (exclusively within the agroecosystem of a particular municipality), and (5) soil potential to pollutant transport (exclusively within the agroecosystem of a particular municipality) (Fig. 1). In the next, we show how raw data from different databases/sources were (pre-)processed for certain inputs.

Fig. 1
figure 1

Overall schematic flowchart of the study presented

Full size image

Land use and management (LUM)

To map the agricultural land use and management in 2017, two principal national databases were used: (i) ARKOD (www.arkod.hr), a geometric database consisting of 1,642,065 polygons with an average plot size of 1.52 ha and (ii) AGRONET (www.agronet.hr), a descriptive database with 224 different crops organised in 24 categories. Both databases serve as e-systems for the identification and classification of land parcels, i.e., records of agricultural land use/management and farm registry to the Agency for Payments in Agriculture, Fisheries, and Rural Development to award agricultural subsidies. However, for the land parcels that were not registered in 2017 in ARKOD/AGRONET, a classification of land use was performed using GIS tools and additional datasets, such as satellite images from Sentinel 2, Google Earth, a digital orthophoto map, and CORINE Land Cover 2018 (e.g., [24, 26]), providing a visual interpretation of Croatia. The overall mapping accuracy for each (of 25 total) category of agricultural land use was estimated, such that 500 points were randomly selected within the polygons of a specific crop category for examination (true/false) of each point to determine whether it was within a mapped polygon (true) or not (false) (more in [31]).

Such processed data generated a map of total agricultural land use and management in Continental and Mediterranean regions (Additional file 1: Figure S1), with a description of particular land and crop categories (Additional file 1: Table S1). Land use and management (LUM) was introduced into the model as input variable 1, in which LUM represents a ratio (%) of total agricultural land area (TALD) to total municipality area (TMA) as follows,

$${\text{LUM}} = \frac{{{\text{TALD}}}}{{{\text{TMA}}}}\;(\% ).$$
(1)

Generated data were distributed among eight classes using the statistical quantile method [32], where the 1st class consisted of municipalities for which LUM varied from 1 to 16% (the lowest share of agricultural land), and the 8th class consisted of municipalities for which LUM varied from 63 to 89% (the highest share of agricultural land) (Additional file 1: Figure S2). The same quantile approach was used for distribution of all other input variables presented in next sections.

Load from mineral N fertilisation (N-min)

This section describes the application (fertilisation) of the total amount of mineral nitrogen (N-min/ha) to the total agricultural land in 2017, excluding the area under organic farming management (Additional file 1: Figure S1b) given that N from mineral fertilisers (N-min) is not allowed in organic agricultural management [33]. The methodology to determine the N-min load at the municipality level was derived from the national total, like in some recent studies [17], and was explained in detail by [34], briefly considering the next: (i) the total amount of mineral fertiliser sold (domestically produced and imported) in Croatia during 2017 [30]; (ii) average N consumption by crops [35]; (iii) consultation regarding the N-fertiliser dosage with experts from the National Agricultural Advisory Service, and (iv) reference yields for the most represented national crops in 2017 [36]. Given the markedly different (agro)ecological conditions among Continental and Mediterranean regions (discussed later), the processing of the data resulted in an estimation of the total N-min fertilisation (kg/ha) for each land use (crop) category (Additional file 1: Table S2). The net mineral N fertilisation (N-min) was introduced into the model as variable 2, following the next calculation at the municipal level. The total area of each crop grown in a certain municipality was multiplied by the estimated N-min for the particular crop, to be than summarised and divided by TALD as follows;

$${\text{N-min}} = \frac{{\sum_{i = 1}^{i = 24} {\text{CROP}}_i \cdot {\text{N-min}}_i }}{{{\text{TALD}}}}\left( {{\text{kg}}/{\text{ha}}} \right).$$
(2)

The generated data were categorised at the municipality level into six classes, in which the 1st class included the municipalities with the lowest N-min rate (0.1–54.4 kg N/ha) and the 6th class represented municipalities with the highest N-min dosage (108–129 kg N/ha) (Additional file 1: Figure S3).

Load from organic N fertilisation (N-org)

The organic N fertilisation rate (N-org) on agricultural land as a result of livestock activity/animal farms was determined using the methodology explained in more detail by [37], which considered the following: (i) the total number of animal farms and their spatial (GIS) location provided by the National Farm Registry database; (ii) the number of animals kept on the farm during 2017 by species and production category determined using the Unified Register of Animals database; (iii) the amount of generated manure estimated according to common standards for a particular livestock species/category (e.g., [38, 39]), considering the technological specifications, production level, and specific management conditions for a particular municipality; and (iv) the number and power capacity of biogas plant facilities operating in 2017 [40] with the amount of livestock manure (N-org) required for their operation, estimating the amount of N-org that does not pose a direct impact (load) on the water resources. Considering these sub-criteria and using the same approach as for the previous input variable (N-min), the processed data resulted in an estimated application of net organic N (N-org) for each land use category, including the areas under organic farming management (Additional file 1: Figure S1b) given that organic N fertilisation is allowed in conventional/organic agriculture at a rate up to 170 kg N/ha [33] as follows;

$${\text{N-org}} = \frac{{\sum_{i = 1}^{i = 24} {\text{CROP}}_i \cdot {\text{N-org}}_i }}{{{\text{TALD}}}}\left( {{\text{kg}}/{\text{ha}}} \right).$$
(3)

The generated data were categorised into six classes, in which the 1st class represented the municipalities with the lowest N-org fertilisation rates (0–6.2 kg N/ha) and the 6th class municipalities exhibited the highest rates (38.1–154 kg N/ha) (Additional file 1: Figure S4).

Natural aquifer vulnerability to pollution (Aq)

Hydrogeologically, Croatia is divided into two markedly different basins (regions): (i) Pannonian (dominates over the continental areas) and (ii) Dinaric karst (dominates over the Mediterranean and mountainous areas) (Additional file 1: Figure S5). To map the natural aquifer vulnerability in the Pannonian region we used the generated outputs of the study by Brkić et al. [41]. Briefly, by using the SINTACS method [42,43,44] seven thematic layers were defined to represent the assessment of the natural properties of the hydrogeological system, as follows: (1) depth to groundwater, (2) effective rainfall infiltration, (3) characteristics of unsaturated aquifer zone, (4) characteristics of saturated aquifer zone, (5) soil properties, (6) hydraulic conductivity of aquifer, and (7) terrain slope (more in [41]). These parameters were multiplied by the appropriate weighting factors to calculate the natural vulnerability index [42,43,44]. Because of the significantly different hydrogeological conditions of the karstic aquifers and environment, the vulnerability of a natural aquifer for the entire Dinaric karst region was mapped from the study outputs by Biondić et al. [45] who considered: (1) hydrogeological characteristics of the aquifer (primarily based on lithologic composition), (2) degree of karstification (e.g., concentration of sinkholes, pits with water, sinks), (3) terrain slope, and (4) total precipitation.

Finally, given the well-known importance of hydrogeology and land use in transporting pollutants and their pollution potential [26, 46], in the last procedural step, the natural aquifer vulnerability (Aq) was overlapped with the agricultural land use map (Additional file 1: Figure S1), excluding the non-agricultural land within the municipality. By this operation, we maximised the inputs from agroecosystems to groundwater sources, striving to prevent as many other types of loads as possible (e.g., [20]). Consequently, a summarised map of aquifer vulnerability to N pollution on agricultural land areas, comprising the six classes was generated, in which the percentage of the dominant aquifer vulnerability class was assigned to the entire municipality (Additional file 1: Figure S5).

Soil potential to pollutant transport (SoP)

The last considered input variable represents the natural soil potential to transport (SoP) of a certain pollutant inside the relatively shallow solum profile (i.e., from the pedological perspective up to the parent material) in the descending/lateral direction, considering the soil hydraulic conductivity, particle size distribution, and soil organic matter content, as presented in detail by [47], 48. Briefly, soils with a high content of clay particles (> 40%) and very low (≤ 0.01 m/day) hydraulic conductivity pose the lowest risk (very poor) potential to pollution transport, whereas soils with high hydraulic permeability (> 1.0 m/day) and a lower share of clay particles (< 33%) pose the greatest risk (very high potential for pollution transport). In addition, important pedo-variables that affect the (de)sorption of pollutants beside of mentioned mineral clay fraction, represents and soil organic (complex humics) colloids [8]. Accordingly, texture-heavier clayey and organically-enriched (humic) soils, enable a higher sorption potential (i.e., very poor potential for pollution transport) vs. organically depleted and texture-lighter soils (i.e., very strong potential for pollution transport). Synthesising such pedo-data (related to hydraulic, texture, and organic properties) with the national soil map of Croatia (scale 1:300,000; [47, 48]), a new generalised map of soil potential for pollutant transport with six classes was produced (Additional file 1: Figure S6). The soil class that occupies the largest area within the total agricultural land (excluding non-agricultural lands, such as forests, watercourses, and urban areas) of a particular municipality was defined as the dominant class (load) for the entire municipality (Additional file 1: Figure S6).

Syntheses of input variables into the LUMNAqSoP approach and scoring of candidate stations and sensitivity analyses

The LUMNAqSoP methodology described in the previous section yielded a classification of five input variables, calculated at the municipal level and expressed as a numerical value (Table 1). Accordingly, the first input variable generated eight classes, whereas the rest of the inputs comprised six classes, wherein the first three variables were numerical and the last two were descriptive (Table 1). The maximal potential total loads (scoring) for a particular candidate groundwater station at the municipality level was calculated as presented in Eq. 2 and Table 1:

$${\text{TL}} = \sum {{\text{LUM}} + {\text{N-min}} + {\text{N-org}} + {\text{Aq}} + {\text{SoP}}.}$$
(4)
Table 1 Classification of input variables and scoring of candidate groundwater station according to the total agro-environmental loads at the municipality level
Full size table

Considering that some of the input parameters may be subject to significant uncertainties, it is common practise to perform the sensitivity analyses with some relevant and realistic scenarios [49, 50] for these parameters in order to assess their impact on the LUMNAqSoP outcome, i.e., the total agri-environmental loads at the municipality level (Table 1). Therefore, we considered the values generated from Table 1 as a baseline scenario. To assess the sensitivity, we altered selected input parameters (LUM, N-min and N-org due to their realistic annual variations) by a relative amount, i.e., an increase (+) or decrease (−) by 10% compared to the baseline scenario, and recalculated the total agri-environmental loads and finally compared these outputs with the baseline scenario. Briefly, in the second scenario we increased by 10% LUM per municipality assuming that conventional farming management increases, and conversely, in the third scenario we reduced by 10% LUM assuming that organic farming management increases, as well as the associated implications on N-min/N-org management (explained in “Load from mineral N fertilisation (N-min)” and “Load from organic N fertilisation (N-org)” sections) in both scenarios.

Data processing and methodology for the LUMNAqSoP testing and evaluation

A generated map of total loads obtained using the LUMNAqSoP methodology at the municipality level (previous section) was overlapped with a map of 151 new candidate stations (GPS-coordinated) for the groundwater quality monitoring network. Analysed nitrate concentrations from two existing national water quality monitoring programmes were used as benchmarks for the LUMNAqSoP testing and evaluation. In total, nitrate concentrations data from the national groundwater network at 391 stations (Water Framework Directive monitoring programme) and from the surface water network at 148 stations (ND monitoring programme) from the 2014 to 2017 period were examined. The sampling frequency for nitrate concertation detection in surface/groundwater samples was 4–12 times per year, whereas all collected samples were processed and analysed in accredited national laboratories according to standard procedures and methods [51]. Nitrate time data series were analysed using descriptive statistics in Statistica 10 software. Datasets with surface/groundwater nitrate concentration from particular existing monitoring locations, as well 151 new candidate stations for future monitoring, were obtained from a national authority of water management in Croatia (Hrvatske vode).

Finally, all maps were created using the ArcMap 10.5.1 platform, applying incorporated various Arc tools and operations (e.g., [26, 29, 43, 44]).

Results and discussion

Spatial variability and distribution of total agro-environmental loads to groundwater

The total agricultural land in Croatia in 2017 occupied 2,503,169 ha with 1,908,405 ha predominantly under conventional management and a smaller portion (104,253 ha) under organic farming, which has recorded constant annual growth (~ 9700 ha) since 2007 (Additional file 1: Figure S1, Table S1; [52]). For comparison, at the EU level in the 2012–2018 period, the portion constituting organic farming area in the total utilised agricultural area increased from 5.6 to 7.5%, covering approximately 13.4 M ha in 2018 [53]. The most dominant crop in the examined continental region (occupying 312,750 ha) was maize, followed by other cereals (247,318 ha), fodders (66,270 ha and oilseeds (57,613 ha) (Additional file 1: Table S1). In the Mediterranean region dominated permanent crops (grapes, olives) with organic pastures and meadows in the mountainous area (Additional file 1: Figure S1, Table S1). However, of the total agricultural land in 2017, a significant share (~ 25%, 594,764 ha) was characterised as unutilised (uncultivated; mostly overgrown with native vegetation), in both the continental and Mediterranean regions (Additional file 1: Figure S1 and Table S1). Such a large share of unutilised agricultural land is one the indicators of pronounced de-agrarianisation, which is very likely related to depopulation and de-ruralisation that Croatia has experienced over recent decades, albeit inadequate agricultural/demographic strategies and politics cannot be disregarded. The overall mapping accuracy was 92.3% (Additional file 1: Table S1) and ranged from 87% (vegetables) to 100% (sugar beets) (data not shown).

Generally, municipalities with the highest agricultural land use intensity (> 85% of total municipality area) were detected in the continental easternmost part of Croatia (Vukovarsko-Srijemska County) and in the upper north (Medimurska County and Koprivnicko-Krizevacka County) (Additional file 1: Figure S2). Next, municipalities with an agricultural land use score in the lowest class (up to 16% of total municipality area) were mostly located in the Mediterranean region, owing to well-known natural constrains and principally hill-sloped and karstified stony terrain with shallow low-fertility soils (e.g., [25]). However, municipalities in classes 6 and 7, mostly located in lower terrains (known as karst fields or poljes with fertile deep soils) near the seaside in the Istrian and central Dalmatian regions, were confirmed as well (Additional file 1: Figure S2; [54]).

High classes of total N-min load were mostly in accordance with high classes of land use, owing to the fact that dominant and widely grown crops (cereals, oilseeds) were located in traditionally intensive agricultural areas and fertilised with relatively higher amounts of N-min (Additional file 1: Figure S3, Table S2; [34]). Municipalities with the lowest N-min load were mostly located in the mountainous Mediterranean, with spatial hotspots of total N-min loads in the Istrian region that are associated with intensive horticultural production/fertilisation. Assessing the accurate amount of N-min fertilisation is a critical point in municipality N balancing, and the most sensitive factor in estimating regional N surplus [17]. Namely, In Croatia as well in some other EU Member States (e.g., Germany), mineral fertiliser sales statistics are available at the national level; thus derivation of the regional (i.e., at municipality/county level) mineral fertilising amounts remains among the biggest methodological challenges [17].

An analysis of the total N-org load indicates that the three municipalities with the highest score (range of 130–154 kg/ha) are located in the continental north and are characterised by intensive livestock production (Additional file 1: Figure S4). Other municipalities with a high score are small municipalities with < 2000 ha of utilised agricultural land (data not shown) but with significant livestock production (mainly poultry, pig, or dairy cattle farms) operating within large production systems [37]. In the Mediterranean region, as expected, most municipalities scored the lowest loads in N-org, whereas detected municipalities with higher loads are generally related to sheep farming [37].

Areas with very high and high natural aquifer vulnerabilities, covering approximately 3000 km2 in the continental portion, are located in the north-west part of Croatia (in Varaždinska and Međimurska County) and in the central, i.e., Zagreb area (Additional file 1: Figure S5). It is characteristic of alluvial aquifers of very good hydraulic properties, with relatively low depth to groundwater (< 5 m) and a poor protective function of the non-saturated zone and soil [41]. In the southern karstic part of Croatia, a very high natural vulnerability is defined near and within the immediate vicinity of swallow-holes (ponors) and their zones. Areas of high natural vulnerability are those developed on well-permeable carbonate rocks with a lower depth to groundwater, as in central Istria and the Kupa River basin ([45]; Additional file 1: Figure S5).

In summing up the highest classes of the input variables presented above (LUM with eight classes, N-min/N-org, and aquifer vulnerability and soil sensitivity with six classes each), the maximum possible total loads for a particular groundwater station candidate at the municipality level was set to 32 (Table 1). However, in the bassline scenario, none of 556 examined municipalities accumulated maximally possible loads, only two municipalities (in the northern territory) accounted 31 total loads, four municipalities accumulated 30 loads, while the largest number of municipalities (53) totally scored 17 loads (Fig. 2). In addition, by evaluating the sensitivity of total agri-environmental loads for all candidate stations to changes in inputs (LUM, N-min and N-org) in two opposite scenarios, we confirmed the responsivity of the LUMNAqSoP methodology as well (Fig. 2). For example, by increasing conventional agriculture by + 10% LUM (and N-org/-min management) on the national territory, nine municipalities scored 30 loads (5 more vs. baseline scenario), while the largest number of municipalities (55) totally accounted 18 loads (Fig. 2). In the third scenario, where conventional agriculture was reduced by − 10% LUM (at the expense of organic agriculture), only one municipality achieved 30 loads (4 less vs. baseline scenario), while the largest number of municipalities (60) achieved a total of 16 loads (Fig. 2).

Fig. 2
figure 2

Distribution of municipalities according to totally scored agro-environmental loads on groundwater resources with respect to nitrate pollution in the baseline and ± 10% Land Use and Management (LUM) change scenarios

Full size image

Accordingly, a final map with classified total loads at the municipal level was overlapped with GPS locations of 151 new candidate stations for the national ND groundwater quality monitoring programme, showing three priority levels (Fig. 2). Namely, since introducing all examined 151 stations into the national groundwater monitoring network is very likely going to be time-consuming and technically challenging (notably the installation of piezometers in deeper-positioned karstic aquifers), the prioritisation of locations is one of recommended approaches to national water management authorities (e.g., [26, 29, 55]). Thus, the groundwater monitoring process should be first introduced at locations in which targeted agricultural impact and environmental vulnerabilities were confirmed as the most intensive at the municipality level, i.e., candidate stations that scored: (i) the highest number of total loads (53 stations with the highest priority status), then (ii) those with less intense loads (59 with the moderate priority status) and (iii) finally stations with the least intensity (39 with the lowest priority status) in vulnerability to nitrate contamination (Fig. 3).

Fig. 3
figure 3

Priority classification map of 151 new candidate stations for the Nitrates Directive groundwater quality monitoring program in Croatia based on totally scored agro-environmental loads at the municipality level in the baseline scenario

Full size image

Moreover, the priority classification of candidate stations in the − 10% LUM change scenario responded different compared to the baseline scenario, i.e., 50 (11 less), 62 (3 more) and 39 (14 less) stations had the lowest, moderate and highest priority, respectively (Additional file 1: Figure S7), while in the + 10% LUM scenario even 64 new candidates (11 more vs. baseline scenario) were ranked with the highest priority (Additional file 1: Figure S8). A well-designed and conducted water monitoring programme of nitrates (NO3−) in water resources has proven to be a successful measure for controlling the threshold of 50 mg NO3−/L in the surface or ground water owing to agricultural activities [56]. However, it is clearly shown that a majority of proposed new candidate stations distributed inside the most vulnerable and intensive agroecosystems (44/53 with the highest priority status) are actually outside of currently designated nitrate-vulnerable zones (Fig. 2; [57]), i.e., areas in which water resources are already exceeding or is likely to exceed a threshold of 50 mg NO3−/L (e.g., [56]). The current situation, related to the nitrate-vulnerable zones, is not homogenous, and huge variations can be found across the EU territory [16, 58,59,60]. For example, some member states have decided to provide the same level of protection throughout their entire national territory (e.g., Slovenia, Austria, Ireland, Germany), whereas others such as Croatia (e.g., Italy, Spain, Portugal, Greece) have sporadically designated nitrate-vulnerable zones (Fig. 2; [2, 16, 56]). For instance, in Croatia, current vulnerable areas according to the ND are distributed in 75/556 municipalities, covering a total area of approximately 5091 km2 or < 10% of the national land territory ([57]; Fig. 3). Irrespective of the diverse criteria for the proclamation of nitrate-vulnerable zones, if they are accurately defined, then it is possible in such areas to implement proactive targeted measures (e.g., reinforced national action programmes, more sustainable agricultural practices; maximal annual use of 170 kg N from livestock manure/ha) for the protection of hydro-resources from the agricultural nitrate origin and control of the contamination risk [18, 56, 58]. However, Croatian agriculture has been very dynamic, not only in the past decades (e.g., net mineral fertiliser consumption in the 2012–2015 period dropped by > 30% over the 2008–2011 period, areas under organic farming increased > tenfold in the 2007–2017 period) but also in the shorter 2012–2015 period (recorded the sharpest decline in dairy cattle farms of > 19% at the EU-28 level) [16, 54, 61]. Other EU-28 member states have also recently experienced significant changes in agroecosystems inside nitrate-vulnerable areas (e.g., [2, 56, 58,59,60]). Consequently, the total area under nitrate-vulnerable zones has increased in the 2012–2015 period by 11.5% (223,963 km2), reaching 2,175,861 km2 or ~ 61% of the agricultural area [16].

Therefore, such strong changes in the agricultural structure and intensity should be reconsidered, and a revision of Croatian nitrates vulnerable zones would be highly recommended (Fig. 3). Furthermore, actions such as the revision of nitrate-vulnerable areas, implementation of action programmes, and estimation of mid-/long-term nitrate trends in national water resources, are obligatory for member states at least every 4 years [18]. Unfortunately, Croatia is one of the few EU member states that does not report mid-/long-term forecasts/trends on water quality with respect to nitrates in water resources as the result of measures in action programmes combined with the implementation of several agro-environmental measures incorporated into the Rural Development Programmes [16]. Thus, results provided by this approach can provide to national or (sub)regional water/environmental authorities a strong and reliable argument for: (i) the reassessment (i.e., derogation) of currently designated nitrate-vulnerable zones or (ii) the designation and extension of new nitrate-vulnerable zones in Croatia (Fig. 3). Similar recent examples have been confirmed in the UK (i.e., Northern Ireland), Belgium (i.e., Flanders), and Denmark, for which the European Commission has approved requested derogations pursuant to the ND in all cases [58,59,60]. In addition, there are lots of other agroecosystems, where action programmes over the last two decades were not effective in protection of groundwater against pollution caused by high nitrates loads [17].

Comparison and validation of LUMNAqSoP results with data from the groundwater quality monitoring program

Nitrates are generally one of the largest problems with regard to the groundwater quality in many EU member states [2, 12, 16, 53, 56], including Croatia; however, it is often very difficult to determine their agricultural/industrial/urban origin (documented below). Comparing all stations with the highest priority status (53 in total, which accumulated the most agro-environmental loads in 2017) with monitoring stations, where mean nitrate (NO3−) concentrations in groundwater were elevated (> 30 mg NO3−/L) during the 2014–2017 period, a very good match, i.e., overlapping of LUMNAqSoP vs. groundwater monitoring stations is clearly confirmed (Fig. 4). For example, higher average nitrate concentrations were found in Panonian (10.2 mg NO3−/L; σ 3.7) vs. Dinaric (4.9 mg NO3−/L; σ 1.41) regions at 391 observations sites during the 2014–2017 period (Additional file 1: Table S3, Fig. 4). In the Panonian region, concentrations above a threshold of 50 mg NO3−/L were recorded in a permeable alluvial aquifer in the northwestern part of Croatia at four monitoring stations (Fig. 4, Additional file 1: Table S3). Above the aquifer at those locations, there are mostly no protective semipermeable deposits; thus, strong loads and risks from agricultural pollution (and other) activities [62] are more pronounced. Additionally, in Varaždinska County, besides the alluvial aquifer, important carbonate aquifers exist in the mountainous part of the county, which contain groundwaters of low nitrate concentrations (on average up to 5 mg/L; Additional file 1: Table S3) as there are no pollutant sources in the recharge area [63, 64]. In north-eastern Croatia, the highest nitrate concentrations (112 mg/L) were identified only at a single station in the immediate vicinity of the state border with Hungary in Osječko-baranjska County (Fig. 4, Additional file 1: Table S3). Furthermore, increased average nitrates concentrations (104.9 mg/L) were registered at two locations in the eastern part of the Vukovarsko-srijemska County (Fig. 4, Additional file 1: Table S3), where generally, the most important groundwater resources are accumulated in the aquifers of intergranular porosity, above which relatively thick (~ 50 m) covering aquitard lies [65, 66]. Consequently, the aquifers in this area are dominated by reductive conditions; thus, nitrate concentrations are generally low [67]; meanwhile, increased nitrate content in groundwater is mainly related to shallow permeable deposits, which can be found within the covering deposits (Fig. 4, Additional file 1: Table S3).

Fig. 4
figure 4

Comparison of stations with the highest priority status (most accumulated agro-environmental loads) vs. groundwater monitoring stations for which mean nitrate concentrations were increased in the 2014–2017 period, with examples of deviation shown for the areas of Pula and Zagreb city

Full size image

In the Dinaric karst region, mean nitrate concentrations in groundwater generally did not exceed 5 mg/L, with the exception of the southern part of the Istria peninsula, where at certain locations, nitrate concentrations occasionally exceeded 50 mg/L (Fig. 3, Additional file 1: Table S3) and likely originated from pollution sources other than agriculture (discussed below). Such results are in agreement with recently reported findings that the lowest nitrate concentrations in the 2012–2015 period were observed in captive and karstic groundwater across the EU-28 member states, for which only 5% of groundwater stations recorded ≥ 50 mg NO3−/L [16].

Certain deviations are evident with regard to the smaller number of stations in the southern part of Istria (Pula city area) and in the wider Zagreb city area (Fig. 4), where elevated NO3− concentrations in groundwater were observed but not predicted by the LUMNAqSoP approach. For instance, in the 2014–2017 period at four monitoring stations in the vicinity of Pula city, average nitrate concentrations were 25–49 mg/L (Fig. 4). In addition, relatively higher nitrate concentrations (approximately 30 mg/L) were recorded at four monitoring stations on the left bank of the Sava River within the Zagreb urban area, where there is actually no land used for agricultural production (Fig. 3). Such exceptions can be explained by the fact that at both urban locations (Pula and Zagreb city), elevated nitrate concentrations in groundwater were most likely not exclusively of agricultural origin, but rather from some other additional or dominant source(s) (e.g., discharge of untreated industrial/urban wastewaters). However, further investigation through the application of appropriate research approaches and designs is needed to better understand these results (e.g., [54, 68]). For instance, [68] recently demonstrated that so-called Before-After-Control-Impact designs are very applicable for the assessment of natural-/human-induced spatiotemporal changes, notably when the observed locations cannot be randomly selected owing to different objective constraints, such as the selection of monitoring stations (e.g., access requirements, land ownership, costs, logistical support).

Namely, adequately designed (e.g., number and spatial distribution of locations) and performed (e.g., sampling frequencies, timing) monitoring of water quality is the starting point for the successful implementation of the ND. However, complete data regarding the contribution of agriculture to N discharge into water resources have not been provided by all member states thus far [16]; one of the causes is still (fully) unestablished monitoring networks in which many new candidate stations need to be installed (Fig. 4). For instance, in the last reporting period 2012–2015, the total number of elaborated groundwater monitoring stations incorporated within the ND in EU-28 was approximately 34,900 stations [16] or ~ 1246 per member state. The average density of the ND network in the EU-28 was ~ 8 stations/1000 km2 of land area, with the highest densities in Malta (130/1000 km2) and Belgium (97/1000 km2) and the lowest densities in Finland and Sweden (< 1 station/1000 km2) [16]. In such context, Croatia is above (12 stations/1000 km2) the EU-28 average density in respect to the ND groundwater network (Additional file 1: Table S3). Next, the annual sampling frequency per station at the EU level was also inconsistent (on average ~ 2) and varied from < 1 in Denmark, Latvia, Poland, and Sweden to ~ 5 in Belgium [16] and ~ 16 Croatia (Additional file 1: Table S3). Such information confirms that the comparability and interpretation of nitrate data among particular countries should be considered with caution owing to discrepancies in (ground)water monitoring approaches and strategies. However, results presented here (Fig. 4, Additional file 1: Table S3) are consistent with a recent report that Croatia belongs to a small group of EU-28 member states (i) with the lowest number of groundwater stations for which average concentrations exceeded the threshold of 50 mg NO3−/L and (ii) which mean nitrate concentrations are several fold below the EU-28 average [16].

Comparison and validation of LUMNAqSoP results with data from the surface water quality monitoring program

A comparison of stations that achieved the highest priority status and surface water locations with elevated nitrate concentrations also shows very good agreement, i.e., overlap (Fig. 5).

Fig. 5
figure 5

Comparison of stations with the highest priority status (most accumulated agro-environmental loads) vs. surface monitoring stations for which mean nitrate concentrations were increased in the 2014–2017 period, with an example of deviation provided by the area of Kutina city

Full size image

Monitoring of the nitrate concentration in surface waters within the ND program for this study was considered from 148 stations (104 in Panonian and 44 in Dinaric region) in the 2014–2017 period (Additional file 1: Table S4). As in the case of groundwater quality (previous section), mean nitrate concentrations in examined surface waters were generally higher in Pannonian (8 mg NO3−/L; σ 4.0) vs. Dinaric karst (0.71 mg NO3−/L; σ 0.30) region but were substantially below the threshold, rarely exceeding even 25 mg NO3−/L (Additional file 1: Table S4, Fig. 5). For instance, mean nitrate concentrations above 25 mg/L were the most concentrated in northern Croatia, Međimurska County (at five stations) and Varaždinska County (at four stations) than in the central part in Sisačko-moslavačka County (at one station; Kutina area, explained below) or eastern part of Osječko-baranjska County (at three stations) (Fig. 5; Additional file 1: Table S4). In the Dinaric karst region, the mean nitrate concentrations in surface water stations were 0.71 mg/L (Additional file 1: Table S4; Fig. 5), i.e., almost sevenfold lower than in groundwaters of this basin (Additional file 1: Table S3). Accordingly, lower nitrate concentrations in surface waters in Croatia are accompanied by less intensive surface water monitoring as well (Additional file 1: Tables S3, S4), similar to the EU level [16]. Namely, according to the last reporting period 2012–2015, the total number of examined stations in the surface waters in EU-28 reached 33,042 stations (increased by 23% compared with the 2008–2011 period; [16]), or 1180 stations per member state. The average density was 7.6 stations/1000 km2, with the highest densities in the Czech Republic, Belgium, and the UK and the lowest densities in Germany, Finland, and Croatia (3 stations/1000 km2); meanwhile, the sampling frequency over a year varied from once (Sweden) to approximately 20 in Ireland and 35 in Croatia (Additional file 1: Table S4, [16]).

However, at this level of comparison, some inconsistencies with presented LUMNAqSoP methodology and certain number of locations in the area of the eastern part of Istria and in the vicinity of Kutina city (i.e., Kutinica station; Fig. 5) were confirmed. Unlike most of the stations examined above, Kutinica station is definitely not a location suitable for monitoring nitrates of agricultural origin, given the strong evidence of the predominance of influence and loads from Kutina city and sources related to the petrochemical industry (e.g., [61]). Thus, such discrepancies can also be explained as very likely other more dominant (than agriculture) sources of nitrate pollution, such as industry/households; however, additional studies (as stated above for certain deviating groundwater stations) are needed for further elucidation.

Conclusions

The EU ND is one of the first large-scale initiatives for proactive protection from excessive nitrates from agricultural sources and their discharge into the aquatic environment. Over the last three decades, through the ND, many EU states and regions have succeeded in decreasing the nitrate concentration in surface and groundwater and remediating many eutrophicated surface waters. Such positive outcomes have hinged on well-structured and effectively executed national water quality monitoring programmes (e.g., spatial distribution of monitoring station, frequency and timing of water sampling) and strategies (e.g., methodology for detection of nitrate-vulnerable zones, proactive action plans/programmes, sustainable fertiliser/manure management) that are still hugely variable among the EU member states.

In this study, we developed a new very robust LUMNAqSoP approach for prioritising 151 groundwater stations for the comprehensive implementation of water quality monitoring programme in Croatia pursuant to ND. The integrative approach developed at the municipality level analyses and synthesises the most dominant agricultural loads (land use and management, total mineral and organic N fertilisation) and natural vulnerabilities (aquifer and pedosphere sensitivity to N pollution) to identify the monitoring stations with the highest risk of groundwater pollution by nitrate from the agroecosystems. Comparing the stations which scored the highest agro-environmental loads with actual water stations exhibiting elevated nitrate levels (from monitoring programmes) was showed a good fit across the entire national territory. However, several deviations in the vicinity of urban areas have also been detected, with elevated nitrates at groundwater stations (but not detected by the LUMNAqSoP), likely arising from other dominant (e.g., municipal/industrial) sources rather than from agroecosystems.

By implementing presented approach, it will be possible to strengthen the groundwater quality monitoring programme in Croatia for the purpose of the ND. Croatia is one of the few EU member states that still does not conduct forecast/establish trends of nitrates in groundwater and surface waters within the ND scope, as it lacks a quality and fully-implemented groundwater monitoring programme. The LUMNAqSoP approach can help improve the comparability of the data concerning the extent and trends of nitrates pollution from agricultural sources, as well provide a more realistic assessment of the quality of groundwater in Croatia, enabling the most critical locations (i.e., those with the highest priority status owing to the highest nitrate pollution risk from agricultural sources) to be included in the monitoring programme first. Additionally, generated results can facilitate and focus a dialogue between all the relevant stakeholders (agricultural and/or environmental authorities, farmers, water regulators) for the successful implementation of action programmes in the nitrate-vulnerable areas. Finally, the presented LUMNAgSoP approach can be a very useful tool and platform for national policy-makers and regulators to implement better management and a more realistic shaping (designation/derogation) of nitrate-vulnerable areas, notably in very diverse and dynamic (agro)environmental conditions that adorn Croatia and many other EU member states.