Gaussian simulation of nitrate concentration distribution in the Zagreb aquifer

ABSTRACT Nitrates present one of the most common groundwater contaminants in the world and one of the five major groups of contaminants in the study area. Gaussian simulation (GS) algorithm was used for determining the spatial distribution of average nitrate concentrations from 2010 to 2015 on 95 sampling points. Results indicate two main focus areas of nitrate contamination, located on the left and right bank of the Sava River. Those areas generally extend according to groundwater flow, while areas near Sava River have much smaller concentrations. GS showed that they can be useful for this kind of mapping because they favor abrupt changes in data values which are in this case a result of heterogeneous lithological composition of the aquifer.


Introduction
Nitrate contamination presents one of the most common groundwater quality problems in the world (Almasri, 2007;Peña-Haro, Pulido-Velazquez, & Sahuquillo, 2009) and in the Croatia (Larva, Marković, & Brkić, 2010;Nakić et al., 2013). Nitrates do not have the ability to bond on soil by adsorption and they are also very soluble and have high potential for leaching from surface through unsaturated zone to groundwater (Chowdary, Rao, & Sarma, 2005;Mkandawire, 2008). Numerous studies showed a strong relationship between agriculture and increased nitrate concentration in groundwater (Hosono et al., 2013;Li, Liu, Lang, Zhao, & Zhou, 2010;Liu, Ming, & Ankumah, 2005;Peña-Haro et al., 2009). However, nitrate concentrations can be a consequence of different simultaneous sources, e.g. manure, fertilizers, atmospheric deposition, irrigation, septic tanks, sewage systems, lawns, gardens, dairy lagoons and landfills (Almasri, 2003(Almasri, , 2007Peña-Haro et al., 2009). Different natural factors can also have influence on their stability, mobility and transformation in unsaturated and saturated zones. Lake et al. (2003) determined surface leaching, soil characteristics, drift cover and aquifer type as four main groups of parameters that have influence on groundwater vulnerability to nitrate pollution. Almasri (2003Almasri ( , 2007 noted that nitrate leaching from unsaturated zone presents consequence of complex interaction between on-ground nitrogen loading, land use, nitrogen dynamics, soil characteristics, recharge and depth to groundwater table. In groundwater, they mostly depend on transport (e.g. advection, dispersion and diffusion) and geochemical (e.g. denitrification) processes, but also on aerobic and anaerobic conditions in unsaturated and saturated zones. Kovač, Pavlić, and Nakić (2016) have showed that there is good to very good statistically significant positive correlation between nitrate and dissolved oxygen concentrations in Zagreb aquifer. Kovač, Nakić, and Pavlić (2017) have concluded, based on the results of multivariate statistical analysis, that nitrate concentrations in the Zagreb aquifer are more related to dissolved oxygen concentrations, pH and electrical conductivity values than oxidation-reduction potential and groundwater temperature. Also, Kovač et al. (2017) have pointed out that lower nitrate and dissolved oxygen concentrations in areas near Sava River, and in the eastern part of the aquifer, are probably the consequence of Sava River influence and increase of aquifer depth toward the east where mixing of aerobic and anaerobic water occurs. Nakić et al. (2013) determined nitrates as one of the five main groups of contaminants in groundwater of Zagreb aquifer. Generally, spatial distribution of nitrate concentration in groundwater is a function of numerous physical, chemical and biological processes that take place at the surface, unsaturated and saturated zones. The main objective of this paper was to construct the most representative nitrate concentrations map using Gaussian simulations (GS) algorithm and identify the most contaminated areas with increased nitrate concentrations. The other objective was to test whether the GS could be used for the construction of maps that show spatial distribution of chemical substances in groundwater when the spatial distribution of observation wells is not ideal and when lithological and hydraulic properties vary widely.

Research area
Zagreb aquifer is located in NW part of the Republic of Croatia (Main Map) between mountain Medvednica in the north and Vukomeričke Gorice hills in the south ( Figure 1). It presents the only source of potable water for about 850,000 inhabitants of the City of Zagreb and part of Zagreb County, and it is designated as part of country's strategic water reserves. It covers an area of approximately 350 km 2 . In general, the region is characterized by large variability in lithology, hydraulic properties of the aquifer, pedological features and land use.
Zagreb aquifer is made of two Quaternary age aquifers deposited during the Middle and Upper Pleistocene (lacustrine-marshy deposits) and Holocene (alluvial deposits) (Velić & Durn, 1993;Velić, Saftić, & Malvić, 1999). They present one hydrogeological unit, but indicate geochemical stratification along the depth (Marković, Brkić, & Larva, 2013). Lower Pleistocene deposits are mostly composed of clayey silts and silty clays with sporadic interbeds and lenses of gravelly sands. Lower and middle parts of Middle Pleistocene deposits are dominantly made of sands, while the upper part is made of silts and clays. In the Late Pleistocene, frequent lateral changes of gravels, sands, silts and clays can be found. The Holocene part is mostly made of gravels and sands (Velić & Durn, 1993).
Main stratigraphic determination was performed based on microfaunal and microfloral analysis (Hernitz, Kovačević, Velić, & Urli, 1981;Sokač, 1978). The boundary between Pleistocene and Holocene was set based on the change in the sedimentary environment from lacustrine-marsh deposits to fluvial, and the change in origin in the petrography of the clasts in gravels and sands. With the onset of Holocene, it is presumed that Sava River started to flow and with it came the material from the Alps which was dominantly carbonate in contrast to material which was deposited during Pleistocene, which was mostly siliciclastic (Velić & Saftić, 1991).
In hydrogeological terms, Quaternary deposits contain three units. First unit is the overburden generally made of silt and clay. The thickness of the first unit varies from 2 to 8 m (Ružičić, Mileusnić, & Posavec, 2012). In some areas, overburden is destroyed by anthropogenic influence. Second unit is shallow aquifer made of gravel and sand, while third unit is the deeper aquifer characterized with lateral and vertical alterations of sand, gravel and clay. Hydraulic conductivities of the shallow aquifer go up to 3000 m/day (Nakić et al., 2013).
Zagreb aquifer presents an unconfined aquifer which is in the direct contact with Sava River, the main source of recharge. General groundwater flow is from W/NW to E/SE. The thickness of the shallower aquifer is from 5 to 40 m, while deeper aquifer extends up to 60 m in thickness in the eastern area (Nakić, Posavec, Parlov, & Bačani, 2011). Posavec (2006) stated that during high water levels Sava River generally gives water to the aquifer, while during medium and low water levels Sava River is draining water from some parts of the aquifer. Also, Posavec (2006) determined no flow boundary in the north of Zagreb aquifer, inflow boundary in the south and west, and outflow boundary in the east. Determination of those boundaries was based on the analyses of head contour maps for high, medium and low groundwater levels. However, determination of areas where Sava River is recharging or draining the aquifer strongly depends on hydrologic conditions and its duration and intensity, as well as anthropogenic influence, for example groundwater abstraction. That means that local groundwater flow direction can vary depending on the hydrological conditions and Sava River fluctuations, what is also very important for studying the movement of potential contamination. There are six major well fields which generally pump groundwater from shallower aquifer for water supply purposes. Bačani, Posavec, and Parlov (2010) noted that groundwater levels are declining on average for 1-2 m every 10 years, while the decrease in permanent groundwater reserves is about 4% based on the data from 1976 to 2006. Reasons for groundwater decrease are numerous. The main reasons are associated with deepening the Sava riverbed caused mainly by construction of the power plants and its reservoirs upstream from Zagreb, increase in groundwater abstraction and construction of dikes along the Sava River to prevent occasional flooding of the floodplain. The last reason can consequently lower the potential infiltration of water from the flooded areas into the aquifer. It can be seen that groundwater levels are lower than before what can suggest that some part of potential contamination is trapped in the unsaturated zone and can be slowly released due to intensive rainfall and high Sava River water levels. The thickness of unsaturated zone varies from 8 meters in NW part to 2 meters in SE part. The upper part of this zone is composed mainly of silty to sandy material while the lower part consists of gravels. Both parts are in places intersected with clay layers (Ružičić et al., 2012). Predominantly, three pedological units are developed on these sediments: Fluvisols, Stagnic Podzoluvisols (Pseudogley) and Eutric Cambisols (Sollitto, Romić, Castrignano, Romić, & Bakić, 2010). Except from nitrates, Nakić et al. (2013) have identified pesticides, potentially toxic metals, pharmaceuticals and chlorinated aliphatics as main groups of contaminants. Even though there are different possible potential contamination sources, agricultural activity and leakage from septic tanks and sewage system present the main potential sources of nitrate contamination. On the left bank of the Sava River, there is dominantly urban area what suggests waste water as the main source of nitrate contamination. On the right bank, there are more agricultural areas, where nitrate contamination is probably the consequence of several different and simultaneous sources. Regardless of different anthropogenic sources at the study area, it has been shown that nitrate concentrations are strongly related to dissolved oxygen concentrations in the study area. This suggests that they probably have influence on the stability and mobility of nitrates in Zagreb aquifer (Kovač et al., 2016(Kovač et al., , 2017.

Data and methods
Data used as input for nitrate concentration mapping were spread out across more than 200 km 2 which covers the most of the Zagreb aquifer. All data were obtained from National monitoring program of Croatian Waters. Nitrate concentrations were derived from water samples obtained from 93 observation wells and 2 sampling sites at Sava River, in period from 2010 to 2015. Map was produced based on the average values derived from 10 to 67 measured values, dependent on the sampling interval of certain observation wells. Only observation wells that had data in at least five of six years were used for map production. Average nitrate concentrations varied from 3.52 to 35.59 mg/l of NO − 3 . It is important to note that values indicate that groundwater is safe for drinking according to nitrates concentration and its maximum allowed concentrations (50 mg/l NO − 3 ) defined in Croatian regulations.
As previously stated, nitrate concentrations rely on many factors, one of which is lithological distribution with its respective petrophysical parameters. Therefore, GS were selected rather than some other conventional and deterministic algorithm, which favors gradual transition of parameters trough space. GS are often used in mapping of petrophysical properties and lithofacies in hydrocarbon reservoir modelling and facies definition (Cvetković, 2016;Novak Zelenika, Cvetković, Malvić, Velić, & Sremac, 2013;Novak-Zelenika, Vidaček, Ilijaš, & Pavić, 2017). The basics of the methodology is presented in Hand, Yang, and Moritz (1994), Damayanti and Hicks (1996), Deutsch and Journel (1998) and Sahin and Al-Salem (2001). The first step for building a map using any kind of kriging or GS is building a variogram, which determines the influence of the data point on the map values. More specifically, how much and how far away from the data point. For the successful building of a variogram, a normal distribution of the input data is needed. As it did not follow the normal distribution (Figure 2), the data had to be transformed. Box-Cox transformation was firstly employed, as with it a near normal transformation can be obtained (Horváth, Borka, & Geiger, 2017), but for any kind of γ value (Box & Cox, 1964), the results were unsatisfactory. Box-Cox transformation of the data with most favorable γ value of 0.4 is shown in Figure 2. Normal score transformation proved successful for data to be fitted to the normal distribution ( Figure 2). Experimental variogram was built based on the all of the analyzed data within Petrel software by testing various values of computation of experimental variogram to achieve optimum variogram parameters (e.g. nugget value, major and minor range). An anisotropic variogram was to be expected due to dynamic and complex interrelation between the Sava River and Zagreb aquifer. Most favorable was spherical, anisotropic, with W-E orientation.
Major direction axis reach was 7100 (W-E) and 4600 m on the smaller axis (N-S) with a nugget effect of 0.017 (Figure 3). In other words, sampling points have an influence on assigning cell values in maximum extent of 7100 eastward or westward while the maximum spatial relevance of the sampled point was 4600 to the South or North. Orientation of the variogram was defined based on the Sava River influence on the distribution of groundwater flow, and respectively nitrates concentration in the aquifer. GS, unlike traditional approach, are a stochastic algorithm which means that for the same input data, there can be several output maps/realizations of which all are equally  probable (Deutsch & Journel, 1998). Number of realizations which is needed to adequately define the distribution of a parameter in space varies for different applications (Jakab, in press). For the purpose of nitrate distribution visualization 50 realizations were performed out of which one was selected. The criteria for selecting one realization out of 50 performed was the comparison of the distribution of nitrate data similarity between the input data set and the result map and similarity of mapped nitrate values in Sava River proximity to the ones measured in the Sava River. For better representation, topographic overlay was used to represent urban areas. Grid resolution of the output surface was 50 × 50 meters.

Conclusions
Nitrates present one of the five main groups of contaminants in the study area. By determining the areas with increased nitrate concentrations in groundwater, researchers could be able to map high-and low-risk areas. Areas with the highest risk probably have combination of high nitrogen input and thin unsaturated zone. Therefore, local variations in characteristics of unsaturated and saturated zones have a large impact on the nitrate concentration in groundwater. GS algorithm was used for the construction of nitrates concentration map. Two main focus areas of nitrates contamination in Zagreb aquifer can be observed (Main Map), one in the dominantly urban area located on the left bank of the Sava River, and the other on the right bank, in a dominantly agricultural area. Those contaminated areas mainly extend according to general groundwater flow, i.e. from W/NW to E/SE. Due to very little or no agricultural activity on the left bank of the Sava River, it can be assumed that nitrate concentrations are probably the outcome of leakage from sewage or septic tanks. The highest nitrate concentrations are observed in the central and eastern part of the City of Zagreb where sewage is very old. On the right bank, agricultural activity is much more prevalent so it can present additional possible source of nitrate contamination. In this area, the highest nitrate concentrations are observed south of the Sava River and the urban part of the City of Zagreb, what suggests increased use of organic manure and/or fertilizer, as well as the leakage from septic tanks in rural settlements. Areas near Sava River have smaller nitrates concentrations, probably due to bigger influence of the Sava River, and associated geochemical environment and hydrodynamic mixing in its hyporheic zone, which is consistent with the results of previous research. Generally, it can be concluded that GS can be used for this kind of mapping because they favor abrupt changes in data values much in relation to changes in lithological composition that can be observed at the investigated location.

Software
Prior to mapping, input data were firstly averaged in Microsoft Excel to be statistically viable. Nitrate concentration mapping was performed in Schlumberger Petrel with GS algorithm. The end map was constructed in Corel Draw X6 while the geocoded terrain image was obtained from geoportal of Croatian Geodetic Administration.