Nitrate contamination risk of the Zagreb aquifer

Nitrates present one of the main groundwater contaminants in the world and in the Zagreb aquifer. In order to reduce nitrate concentrations in groundwater, it is necessary to spatially define main nitrogen sources and areas which have the highest risk of nitrate contamination. This paper presents a map of nitrate contamination risk in the area of the Zagreb aquifer. It was constructed based on nine different layers that include natural characteristics and anthropogenic pressures. For the construction of the Main map, which has been defined as the most representative one, 15 different variations have been tested. The Main map has shown that the urban part of the City of Zagreb, especially central and eastern parts, together with the area of Jakuševec landfill and marshalling station, present areas with the highest risk of nitrate contamination, which is consistent with the results of previous research. ARTICLE HISTORY Received 4 March 2019 Revised 1 July 2019 Accepted 2 July 2019


Introduction
Nitrates present one of the most common groundwater contaminants in the world (Hosono et al., 2013;Pennino, Compton, & Leibowitz, 2017;Peña-Haro, Pulido-Velazquez, & Sahuquillo, 2009;Xue et al., 2009) and in the Zagreb aquifer (Nakić et al., 2013). Nitrates (NO − 3 ) are negatively charged and they do not have the ability to bond to the soil by adsorption. They are very soluble and have high leaching potential through unsaturated zone (Chowdary, Rao, & Sarma, 2005;Mkandawire, 2008). Almasri (2007) showed that nitrate leaching from unsaturated zone is very complex and it presents the interaction between different factors, mostly surface nitrogen loading, nitrogen dynamics, land use, soil characteristics, groundwater recharge and depth. When they reach groundwater, their movement generally depends on transport processes, geochemical processes and conditions that prevail in the aquifer. Nitrates can reach the aquifer from different sources, such as fertilizers, manure, waste water, atmospheric deposition, lawns and gardens. However, their stability in groundwater also depends on the existence of the aerobic/anaerobic environment, i.e. processes such as nitrification and denitrification.
Nitrates have been recognized as one of the main groups of contaminants in the Zagreb aquifer (Nakić et al., 2013), which is defined as part of country's strategic water reserves and presents the only source of potable water for the inhabitants of the City of Zagreb and a part of the Zagreb County. The existence of elevated nitrate concentrations in the Zagreb aquifer have been confirmed through recent research (Kovač, Cvetković, & Parlov, 2017;Kovač, Nakić, & Pavlić, 2017). It has been shown that nitrate trends are declining in most parts of the Zagreb aquifer (Kovač, Nakić, Špoljarić, Stanek, & Bačani, 2018), and that nitrate origin is mostly organic, generally associated with wastewater leaching from permeable sewage and septic tanks (Kovač, Nakić, Barešić, & Parlov, 2018). Furthermore, nitrate stable isotopes composition has shown that the denitrification process is negligible in most parts of the Zagreb aquifer. All results have shown that nitrates contamination exists in the Zagreb aquifer, especially in two areas. One area is located in the urban part of the City of Zagreb, on the left bank of the Sava River, while the other is located on the right bank of the Sava River, where agricultural areas are also present. It has been shown that aerobic conditions prevail in the Zagreb aquifer, except in the eastern part where the aquifer deepens and in the vicinity of the Sava River, and that nitrate concentrations are strongly related to dissolved oxygen concentrations Kovač, Pavlić, & Nakić, 2016).
The main goal of this paper was to construct a map which includes variables that are related to the vulnerability of the Zagreb aquifer and those related to the risk of nitrate contamination. For that purpose, the most important physical and chemical parameters of the Zagreb aquifer were taken into account together with spatial distribution of main nitrate sources, i.e. nitrate sources related with wastewater and agriculture. The Zagreb aquifer system is located in the northwest part of the Republic of Croatia (Map), between Mountain Medvednica in the north and Vukomeričke Gorice hills in the south. It covers approximately 350 km 2 of an area which is characterized by great variability in land use, lithology, pedology and hydraulic properties. The Zagreb aquifer is an unconfined aquifer and it consists of three main units: unsaturated zone, shallow aquifer layer and deep aquifer layer. Unsaturated zone thickness generally varies from 2 to 8 m (Ružičić, Mileusnić, & Posavec, 2012) and is disintegrated by anthropogenic influence. Nitrate concentrations are mostly higher in the upper part of the unsaturated zone, and they decrease with depth (Ružičić, Kovač, & Tumara, 2018). This is consistent with the research made in the Fluvisol area of the Zagreb aquifer where it was shown that the upper part of the unsaturated zone is generally permeable, except in dry periods (Ružičić, Kovač, Nakić, & Kireta, 2017), which allows the infiltration of contaminants into the aquifer. Apart from nitrates, pesticides, potentially toxic metals, pharmaceuticals and chlorinated aliphatics have been identified as main groups of contaminants in the Zagreb aquifer (Nakić et al., 2013). Shallow Holocene aquifer consists of alluvial deposits, while deeper Pleistocene aquifer consists of lacustrine-marshy deposits (Velić & Durn, 1993;Velić, Saftić, & Malvić, 1999). The whole system is up to 100 m thick, the thickness of the shallow aquifer varies from five to 40 m, while the deeper part extends up to 60 m (Nakić, Posavec, Parlov, & Bačani, 2011). Both aquifers are hydraulically connected, and the shallow aquifer is in direct contact with the Sava River, which presents the main source of recharge and the influence on the groundwater levels is more pronounced in its vicinity (Posavec, Vukojević, Ratkaj, & Bedeniković, 2017). General regional groundwater flow is from W/ NW to E/SE, but local groundwater flow is subject to the Sava River fluctuations, i.e. hydrologic conditions.

Data and methods
Data used to create the risk of nitrate contamination map in the area of the Zagreb aquifer were adjusted to the area of approximately 300 km 2 . Data were collected from previous research (Kovač, 2017;Kovač, Nakić, Barešić, et al., 2018;Ružičić, 2013), official documents of the Republic of Croatia and various registers of the Republic of Croatia: Hydrogeological map of the Republic of Croatia M 1:300,000 (Biondić, Brkić, Biondić, & Singer, 1996), Outlines of Water Management of the Republic of Croatiasection Groundwater (Biondić & Brkić, 2001), CORINE Land Cover 2012 (Croatian Environmental and Nature Agency, 2018a), WFS layer of The Register of industrial plants in which the presence of dangerous substances was established (Croatian Environmental and Nature Agency, 2018b) and Farm register used in Kovač (2017). Higher concentrations of nitrates in groundwater are thought to be a consequence of multiple different sources, as well as the aquifer properties. The properties of the Zagreb aquifer that were taken into further consideration were: hydraulic conductivity, unsaturated zone thickness, soil, groundwater temperature and dissolved O 2 concentrations. In this research, the main nitrogen sources that were taken into account are land cover, industry, sewage system, septic tanks and farms.
Hydraulic conductivity was the first property that was taken into deeper consideration. Hydraulic conductivity values and its distribution were taken from the Hydrogeologic map of the Republic of Croatia M 1:300,000 (Biondić et al., 1996) and the Outlines of Water Management of the Republic of Croatiasection Groundwater (Biondić & Brkić, 2001). Hydraulic conductivity values are reaching more than 1000 m/ day. It can be seen that most of the evaluated aquifer area has very high hydraulic conductivity that enables very fast movement of contaminants in the saturated zone (Main Map 1), which can be specifically addressed to nitrates in an aerobic environment. Hydraulic conductivities are smaller in the edge sections of the aquifer, where slower groundwater velocities can be expected.
At the Zagreb aquifer area, there are six types of soil (Bogunović, Vidaček, Racz, Husnjak, & Sraka, 1998). Certain parts of the Zagreb aquifer area are not covered with soil, but with urban surfaces and water bodies (Main Map 2). Soil and its properties represent natural protection for groundwater systems. In the Zagreb aquifer area, Stagnasols represent a type of soil through which the leakage is the lowest, and water with contaminants stagnates the longest. The reason for stagnation is the silty-clay-loamy soil structure, density and low permeability. Fluvisol and Gleyic Fluvisols represent very permeable soils, through which contaminants move easily. On the other hand, the properties of Eutric Cambisol, Gleyic Stagnosols and Vertic Pseudogleys are between two previously mentioned groups (Lončarić et al., 2014;Ružičić, 2013).
Unsaturated zone thickness, and its changes over time, is a very important property to evaluate in the analysis of its role in groundwater protection. All contaminants can be trapped here, as long as the groundwater level is below this zone. When groundwater level rises, it can rinse contaminants that are located in the unsaturated zone, which can result in percolation by which contaminants move to the groundwater systems. Taking into account the safety factor, while creating the map of the thickness of the unsaturated zone, the data used were showing the equipotential map for high water levels (Kovač, 2017). The high water levels occurred on 5 November 2012. On the map of the thickness of the unsaturated zone, values vary from −1.93 to 16.78 m. Values less than zero (0), show flooded areas on the above-mentioned date (Main Map 3).
Temperature is a key factor for every chemical reaction. Temperature is considered to be a key factor that affects the kinetics of nitrification and denitrification (Almasri, 2007). In general, reaction rates are assumed to double for every 10°C increase in temperature (Arrhenius rate law). It was determined that groundwater temperature is significantly changing due to changes in depth. In addition, it has been shown that in the first 10 m from the soil surface, groundwater temperature is affected by seasonal changes in temperature. In the interval of 10-20 m of depth, groundwater temperature is 1-2°C lower than the average annual temperature of the atmosphere. Below the depth of 20 m, groundwater temperature is changing depending on the geothermal gradient (Domenico & Schwartz, 1997). Also, it was proved that denitrification in subsoil temperatures was significantly slower than in laboratory conditions (Lind, 1983). Although it is difficult to observe temperature dependency of denitrification rates in temperature stable groundwater environment, Saunders and Kalff (2001) observed that 5°C increase can result in 10% increase of denitrification rate. Furthermore, it was shown that the correlation between water temperature and denitrification rate exists in a permeable reactive barrier system. Denitrification rates at lower temperatures (2-5°C) were approximately 5 mg-N/l/day, while at higher temperatures (10-20°C) they increased up to 30 mg-N/l/day (Robertson, Blowes, Ptacek, & Cherry, 2000). However, Rivett, Buss, Morgan, Smith, and Bemment (2008) have emphasized that one of directions of the future research should be aimed at the exploration of environmental conditions and their influence on denitrification rates. In the Zagreb aquifer area groundwater temperature is not measured at different depth intervals, but average values of groundwater temperature were taken into account for each piezometer. There were 153 values, some of which have the same location, due to nested piezometers. In nested piezometers, only the highest values were taken into further research, which resulted in 138 values. Using Kriging interpolation method the map of the distribution of groundwater temperature at the Zagreb aquifer was created. The values of groundwater temperature shown on the map vary from 11.92°C to14.92°C (Main Map 4).
The real values measured in the aquifer vary from 11.43°C to 17.02°C. The interval of the values on the map is smaller due to interpolation and extrapolation that were carried out on the data. Values obtained in the analysis process were divided into five categories, by 0.6°C. The highest temperatures were measured in the urban part of the City of Zagreb.
Dissolved O 2 concentrations are very important in the research of nitrogen species in groundwater because the stability of nitrogen forms depends on the existence of aerobic/anaerobic environment. Dissolved O 2 concentrations control nitrification and denitrification processes (Balderacchi et al., 2012). NO 3 are the most mobile nitrogen compound in groundwater. NO − 3 form during the nitrification and due to a negative charge, they are not adsorbed onto soil fragments and they can travel to long distances (Lee, Lee, Hyun, Clement, & Hamilton, 2006). Dissolved O 2 data were collected at the same piezometers as the temperature values, after filtrating the values there were 138 values of dissolved O 2 , with the highest values taken into account in the nested piezometers. These 138 values were used to create a map of the dissolved O 2 distribution in the Zagreb aquifer. This map was also created using Kriging interpolation method. Both dissolved oxygen concentrations map and groundwater temperature map were created in the same way. To get more realistic results, data were tested with two different Kriging methods: Simple and Ordinary, which are better explained in Malvić (2008). The results derived from the Simple method were used for the groundwater temperature map, while the results derived from Ordinary method were used for the dissolved oxygen concentrations map. The results have shown similar maps for the middle part of the aquifer, where the piezometers were located, and different, but expected relationships for the edges of the aquifer. According to previous research, the amount of dissolved O 2 determines the type of conditions that are dominant in the aquifer: < 1 mg/l O 2 represent dominantly anaerobic conditions, 1-2 mg/l O 2 represent transitional zone, and > 2 mg/l O 2 represent dominantly aerobic conditions (Dimkić, Brauch, & Kavanaugh, 2008). To obtain a clearer presentation for the Zagreb aquifer, data were divided into five categories: <1, 1-2, 2-4, 4-6 and >6 mg/l O 2 (Main Map 5).
The research area is partly covered with vegetation and partly with urban areas. Different types of cover represent different possible sources of nitrogen compounds. In the area covered with vegetation (agricultural land, grasslands, forests) the main sources of nitrogen compounds are fertilizers and animal waste. In the urban area, the main sources of nitrogen compounds are leaking sewage system, septic tanks, industry, construction pits, etc. To create the map of the land cover of the Zagreb aquifer system, the CORINE Land Cover 2012 maps were used. Data in these maps were divided into four categories (Main Map 6). The first category covers grassland and forests, the second category covers agricultural land, the third category covers urban area, and the fourth category covers locations where exploitation of mineral resources occurs, such as construction pits, waste landfills, and construction sites, i.e. places where natural land cover is seriously damaged and/or is lacking.
The industry is a significant source of contamination of groundwater with various compounds due to various purposes that the water in industry is used for cleaning, heating or cooling. According to the official data of the Republic of Croatia, WFS layer of The Register of industrial plants in which the presence of dangerous substances was established (Croatian Environmental and Nature Agency, 2018b), there are 12 potential industrial contaminants in the Zagreb aquifer area. The sources of industrial contaminants are located in different parts of the aquifer, with the biggest concentration in the industrial area of the City of Zagreb (Main Map 7). To illustrate their area of influence, general course of groundwater flow was taken into detailed observation. In the direction of regional groundwater flow, the influence of industrial objects is assumed to be bigger. The direction of the regional groundwater flow at the Zagreb aquifer area is NW -SE, which means that the area of influence should be bigger to the SE of the object. The area of influence was marked with two arcs. The first arc has a radius of 3000 m, the angle of the arc is 60°, and the orientation of the arc is 135°. The second arc has a radius of 500 m, the angle of the arc is 300°, while its orientation is 315° (Main Map 7). The values used for radius in the direction of groundwater flow were given by Bhadra, Pathak, and Sharma (2013) where the value of 3000 m was used to determine the buffer zone in the area around Pali City (India), mostly in quaternary alluvium, which corresponds to deposits in the Zagreb aquifer area.
In the observed area, there are many farms. According to the Register of farms of Croatian Agricultural Agency (Kovač, 2017), in the wider Zagreb aquifer area there are 1174 farms, mostly on the right bank of the Sava River (Main Map 9). Farms are considered to be a significant source of nitrogen compounds because nitrogen is excreted from animal waste in the form of urea. As some previous research state, if livestock density exceeds 2.2-4.3 cows/ha, then the amount of nitrogen excreted from animal waste can reach up to 800 kg/ha N per year (Roche et al., 2016). However, the data for farms in the Zagreb aquifer area are not precise. Therefore, neither the species bred at these farms nor the number of animals per hectare of land, are known. These findings led to determining the area of influence of each farm in the same way that the area of influence of industrial objects was determined, based on the research of Bhadra et al. (2013). The area of influence of each farm was illustrated as two arcs, one with a bigger radius (3000 m) in the direction of regional groundwater flow (135°), and another one with a smaller radius (500 m) in the direction opposite to the direction of the regional groundwater flow (315°).
Sewage system and septic tanks are considered to be the most significant source of nitrogen compounds in the Zagreb aquifer area. The agricultural activity follows sewage systems and septic tanks with the amount of nitrogen compounds produced (Balderacchi et al., 2012). Some authors consider septic tanks to be the main source of nitrates and nitrates contamination, others find problems in both sewage systems and septic tanks of the big cities, which is also the case with the Zagreb aquifer (Kovač, Nakić, Barešić, et al., 2018). Sewage systems and septic tanks tend to be overloaded with high amounts of wastewater, which leads to leaking through the systems and infiltration in groundwater systems (Almasri, 2007;Delleur, 2007). Many research has shown that sewage systems and septic tanks are the main sources of nitrogen compounds in groundwater systems (Gutiérrez, Biagioni, Alarcón-Herrera, & Rivas-Lucero, 2018;Kovač, Nakić, Barešić, et al., 2018;Lee et al., 2006). The polygon that shows the area at the Zagreb aquifer system where the sewage system was built, was used in previous research (Kovač, 2017;Kovač, Nakić, Barešić, et al., 2018). For this research, the area which is not covered with the sewage system polygon was assigned to the area where there are dominantly septic tanks (Main Map 8).

Map design
The collected data came in forms of point, line and polygon shapefiles. The goal for each data layer was to put it in the form in which it could be easily read. Data that were gathered in the form of polygon shapefile were the hydraulic conductivity, sewage system, land cover and soil. Data in form of line shapefiles were: equipotential lines on 5 November 2012data that were used to create maps of the thickness of the unsaturated zone at the Zagreb aquifer area.
Point data, such as readings from piezometers, were used to create maps with interpolation methods. Those maps were the map of the distribution of dissolved O 2 concentrations and the map of the distribution of groundwater temperature. Furthermore, the maps of farms and industry in the observed area were also made from point shapefiles.
After the data were adjusted, each layer and each category were given a weight factor. The values of weight factors were chosen in the range of 1-10, where value 1 presents the smallest load, and value 10 the biggest load. In order to determine the most representative risk of nitrate contamination map, 15 variations were made (Main Map 1 -Main Map 15, Table S1). In each variation, there were changes in weight factors for one or several layers. Weight factors were changed several times for temperature, dissolved O 2 concentrations, unsaturated zone thickness, farms and industry. In the soil layer, only the weight factor for urban area and water bodies was changed since the weight factors for urban area and water bodies were determined in the land cover layer. That is why urban area and water bodies have weight factor 0 in the soil layer. It was determined that variation 15 (Main Map 15) was the most representative and it is shown as the Main Map of this research.
In the soil layer, given weight factors are 3, 6, and 9. The difference between each group was made taking into consideration other soil types with soil properties between the explored types of soil. Weight factors for the unsaturated zone are 2, 4, 6, 8, and 10. Weight factor allocation started in the category <0 m, to which factor 10 was given. The absence of unsaturated zone represents a great danger for the aquifer from the perspective of risk of contamination, so the weight factors for other categories were given accordingly.
The same principle was applied for dissolved oxygen concentrations layer, except in this case the allocation started with the highest value of dissolved O 2 . Since the values of dissolved oxygen can be higher than 7.2 mg/l, the highest given value was 9. In this case, it is assumed that more stable aerobic conditions provide more stable conditions for nitrate transport. Weight factors for groundwater temperature were given continuously from 4 to 8. The lowest value of 4 was given to the category 11.92-12.52°C. Temperature values for this category show that groundwater temperature can be much lower, but also much higher than 14.92°C, which is why this category was given the weight factor 8.
In the farms and industry layers, it was very difficult to determine which weight factors should be given. In most variations, farms and industry were not used to create the final map. In the last four variations, weight factors for farms and industry significantly differ, because it was unclear how to make sure that their impact was shown realistically. In the end, it was concluded that their weight factors should be 6 for farms and 8 for the industry. Farms and their waste are not controlled and the amount of animal waste deposited on the land is unknown. Because of that, farms are considered a dangerous source of contamination, but not as dangerous as an industry. Previous research (Kovač, Nakić, Barešić, et al., 2018), defined wastewater as the main source of nitrogen, and thus the industry is considered to be a more dangerous source of nitrogen contamination. Consequently, the urban part of the City of Zagreb, especially the eastern part where most of the industry is located, is defined as the area with a very high risk of nitrate contamination. However, the type of industry that is present in the Zagreb aquifer area is not directly connected with nitrogen production or nitrogen compounds, that is why weight factor assigned to the industry is not 10 or 9, but smaller. Although agricultural activities have an impact on nitrates concentrations, they do not present the main source of nitrate contamination, especially in the urban part of the City of Zagreb.
The land cover layer was one of the layers for which weight factors were not changed. The factor for each category is increasing by 2 compared to the previous category. The factors are 2, 4, 6 and 8. Weight factor 2 was given to grassland and forests. Namely, grassland and forests are areas where human activity is lacking or is not significant but exists. These factors were not changed because they did not make a significant difference in the final map. Weight factors for the sewage system and septic tanks were also not changed. Their weight factors have to be very high because they present the main source of nitrogen in the Zagreb aquifer area. However, septic tanks are maintained by citizens and possible leaking should be controlled more often. Weight factors were not changed for hydraulic conductivity either. The values of hydraulic conductivity in the Zagreb aquifer area are very high, sometimes reaching more than 1000 m/day. Due to the existence of the aerobic environment which allows nitrates to be stable and transport to long distances, it was decided to give hydraulic conductivity the three highest values of weight factors.
After determining which layer should be assigned which weight factor, the formula for calculating the index of risk was created. The created formula has elements of DRASTIC and SINTACS methods (Aller, Bennett, Lhr, Petty, & Hackett, 1987;Cvita & De Maio, 1997): where VI stands for the index of risk, n for the weight factor, S for the type of soil, UZ represents unsaturated zone, K hydraulic conductivity, T stands for temperature, O marks the values of dissolved O 2 , C stands for the type of land cover according to categories, SS for the sewage system, ST for septic tanks, F for farms and I stands for industry. In the temperature layer, unsaturated zone layer and dissolved O 2 layer, the values for temperature degrees, the concentration of O 2 and unsaturated zone thickness were obtained for the entire research area of the Zagreb aquifer. After that, each of those values was multiplied by a weight factor. This process ensured that values of temperature, concentrations of dissolved O 2 and unsaturated zone thickness have a bigger impact on the reactions in the aquifer. On the other hand, soil layer, hydraulic conductivity layer, land cover layer, sewage system and septic tanks layer, farms and industry layers, did not carry any values in each pixel and there were not enough data to calculate their values, and in order to calculate their influence on the aquifer, only their weight factors were taken into account.
As mentioned before, variation 15 was chosen as a variation that shows the Zagreb aquifer's risk of nitrate contamination the best, in accordance with previous research. The Main Map shows that the values of risk indexes increase on the southern border of the aquifer. This increase can also be seen in Kovač, Cvetković, et al. (2017). Also, a significant increase in the risk of contamination indexes is visible in the area of the Jakuševec landfill, marshalling station and the City of Zagreb. Significant decrease in the risk of contamination and the lowest values of indexes can be found in the eastern part of the aquifer, which is also consistent with the map provided by Kovač, Cvetković, et al. (2017). The map shows a moderate risk of nitrate contamination, which is in accordance with a decrease in dissolved O 2 concentrations not allowing for the nitrification to occur. Using the software, all values were multiplied by weight factors, after which they were cumulated, and the map of the distribution of the risk of contamination index was created. After that, the values of the index of risk of contamination, that range from 100.03 to 220.4, were divided into five categories which represented a very low, low, moderate, high and very high risk of nitrate contamination in the area of the Zagreb aquifer. Very low risk of nitrate contamination has the values of the indexes below 44.08, low risk of contamination has the values of the indexes varying from 44.80 to 88.16, moderate risk of contamination shows the values of the indexes from 88.16 to 132.24, while high risk of contamination varies from 132.24 to 176.32 in the values of the index of risk of nitrate contamination. Finally, a very high risk of nitrate contamination marks places where the index of risk of contamination varies from 176.32 to 220.4 (Main Map).

Conclusions
Main map (Variation 15) presents the most representative map because it shows most realistically the places where the Zagreb aquifer has the greatest risk of nitrate contamination. This map shows that the urban part of the City of Zagreb, as well as the area near landfill Jakuševec and marshalling station present areas in the Zagreb aquifer where the risk of nitrate contamination is the highest. In addition, these results are consistent with all previous research. In the west, nitrogen load is not as high as in the central part, while in the eastern part the risk of nitrate contamination is in some places moderate due to deepening of the aquifer and the existence of a more anaerobic environment. In those areas, dissolved oxygen concentrations and groundwater temperature are lower, which results in denitrification and slower chemical reactions.

Software
All maps were made using ESRI ArcGIS 10.1 for Desktop. The final map was made in Adobe Illustrator CC 2019.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
The publication process is supported by the Development Fund of the Faculty of Mining, Geology and Petroleum Engineering, University of Zagreb.