Nitrogen Mass Balance and Pressure Impact Model Applied to an Urban Aquifer

: The assurance of drinking water supply is one of the biggest emerging global challenges, especially in urban areas. In this respect, groundwater and its management in the urban environment are gaining importance. This paper presents the modeling of nitrogen load from the leaky sewer system and from agriculture and the impact of this pressure on the groundwater quality (nitrate concentration) in the urban aquifer located beneath the City of Ljubljana. The estimated total nitrogen load in the model area of 58 km 2 is 334 ton / year, 38% arising from the leaky sewer system and 62% from agriculture. This load was used as input into the groundwater solute transport model to simulate the distribution of nitrate concentration in the aquifer. The modeled nitrate concentrations at the observation locations were found to be on average slightly lower (2.7 mg / L) than observed, and in general reﬂected the observed contamination pattern. The ability of the presented model to relate and quantify the impact of pressures from di ﬀ erent contamination sources on groundwater quality can be beneﬁcially used for the planning and optimization of groundwater management measures for the improvement of groundwater quality.


Introduction
In 2018, 55% of the world's population lived in urban areas, and this proportion is expected to increase to 68% by 2050 [1]. Europe, with 74% of its population living in urban areas in 2018, is expected to reach an urban population of 80% in 2040 and approaching 85% by 2050 [1]. Considering this population growth in light of climate change projections, the estimates show that with current practices, the world will face a 40% shortfall between the forecasted demand and available supply of water by 2030. Water security is thus emerging as one of the biggest global risks in terms of developmental impact [2].
Groundwater will have an important role, as it represents over 95% of the world's available freshwater reserves [3]. Compared to surface water, it is better protected and a more stable water resource. Better quality, proximity of the resource, and the relatively low cost of pumping wells make groundwater an increasingly important source for city populations in the future.
Aside from socio-economic benefits, urbanization has also created environmental problems, manifested in the deterioration of air and water quality as well as a significant change in water balance [4][5][6][7]. Urban subsurface infrastructures strongly influence the chemical as well as the quantitative status of groundwater. The subsurface infrastructure of Bucharest (sewer system and water supply network) contributes more than 71% of the groundwater recharge. In addition, more than 77% of the groundwater discharge is associated with non-natural sinks (mainly the sewer system Water 2020, 12, 1171 2 of 15 and the drain of the main sewer collector) [4]. In the Baltimore metropolitan area, the infiltration of groundwater into wastewater pipelines was found to be the most important factor of urban development affecting the groundwater storage, which can cause its decrease by 11.1% [8].
The assurance of good quality groundwater in the urban environment is a complex task, due to numerous human activities and the number of potential sources of contamination [9,10]. Additionally, the limited ability for observation and the delay between the contamination and its detection make it complicated to define the causal relationship between human activities that exert pressures and their impacts on groundwater quality.
Wastewater exfiltration from the sewer system is an important factor influencing groundwater quality, with its impact being strongly variable both in space and time [11][12][13][14]. A review of exfiltration modeling approaches and measurement techniques in terms of both volume and pollutant load showed a wide range of exfiltration rates estimated on the catchment scale, ranging from 1% to 56% of dry-weather flow (DWF) [15].
Nitrate is one of the key groundwater contaminants, and measurements of its concentration are regularly included in groundwater quality monitoring. In Slovenia, about 37% of the groundwater in alluvial aquifers has a poor chemical status, mainly due to high concentrations of nitrate (above 50 mg/L) [16]. Background concentrations of nitrate largely depend on the geological properties of aquifers; in Slovenia, they are below 10 mg/L [17]. Increased concentrations are related to human activities, in particular to agriculture. The main source of diffuse nitrogen input on agricultural land in Slovenia are mineral fertilizers and livestock manure, which account for 84% of the total input of nitrogen [18]. Among the important non-agricultural sources of nitrate in groundwater, wastewater has been identified as one of the major sources of nitrogen in urban aquifers [19].
The impact of nitrogen sources on groundwater have been investigated in a number of studies. Different sources of nitrogen and their contribution to the nitrate content in the groundwater have been analyzed [20][21][22][23][24]. The research on the role of land use [25], vulnerability of groundwater [26,27], processes of leaching [28], and the migration and transformation of nitrogen in groundwater [29][30][31] have furthered the knowledge on factors controlling groundwater nitrate contamination and its mitigation. The spatial and temporal variation of nitrate in groundwater has been commonly simulated with numerical groundwater flow and transport models [24,29,32,33]. A fundamental step in the development of such models is the estimation of nitrogen loading. The broad range of nitrogen sources in urban areas makes reliable estimation of nitrogen load a difficult task in these areas [19].
The drinking water supply of the City of Ljubljana relies on groundwater, primarily on its abstraction from the Ljubljansko polje aquifer, situated beneath the city. The well fields Kleče (55%), Hrastje (12%), and Šentvid (9%) (Figure 1) assure around three-quarters of all the water needs in Ljubljana's drinking water supply system [34]. To determine the urban contamination in the Ljubljansko polje alluvial aquifer, an urban water balance modeling approach was applied [35]. The study focused on an area (approximately 0.76 km 2 ) located within the City of Ljubljana, and indicated that residential land use in an urban area with a high unsaturated thickness may have a significantly smaller impact on groundwater quality than agriculture or industry.
The first catchment scale pressure impact model for the Ljubljansko polje aquifer was developed in the frame of the INCOME project [36]. The aim of the model was to simulate the concentrations in groundwater for a number of pollutants. The model was able to explain practically the entire mass of nitrogen in the sewer system by the known sources. On the basis of these results and an encouraging decrease of observed nitrate concentration in groundwater, more optimistic groundwater management objectives for the year 2017 were set. The targeted average nitrate concentration in the aquifer was set below 15 mg/L and in all observation wells below 25 mg/L [37]. Nitrogen input from agriculture was identified as the most important source of uncertainty in the model. The explanation of high nitrate concentration at the local level, that is, individual observation wells, was pointed out as an important groundwater management activity.
In the present study, a refined and updated groundwater flow and transport model [38] was applied in the catchment scale approach. The aim of the modeling was to analyze the causal relationship between pressure (nitrogen load) arising from the leaky sewer system and from agriculture, and its impact on the quality of groundwater (nitrate concentration) in the Ljubljansko polje aquifer [39]. Nitrogen load from the leaky sewer system was estimated with the empirical model of spatially distributed sewer exfiltration [40] and together with estimated nitrogen surpluses from agriculture [41] used as input into the groundwater transport model. Modeled and observed nitrate concentrations in groundwater were analyzed to assess the uncertainty of the model and to identify priority areas for remediation. table is on average 25 m below the surface. The unconfined intergranular aquifer with high transmissivity (Figure 2) is mainly recharged from the Sava River in the northwestern part of the study area. It drains into the river in the northwestern part, which largely predisposes the groundwater flow direction (Figure 1). The main water balance components of the regional model area (Figure 1) for the 3 year modeling period (2013-2015) are presented in Table 1.

Study Area
The study area (57.91 km 2 ) is located in the flat part of the catchment of the alluvial Ljubljansko polje aquifer forming the subsurface of the City of Ljubljana, the capital of Slovenia. It corresponds to the area of the local groundwater model (Figure 1). The mean altitude of the area is 292 m above sea level (m a.s.l.) and prevailing land cover types are artificial surfaces (59%) and agricultural areas (39%) [42].
The Ljubljansko polje alluvial aquifer is developed in a tectonic basin, filled with highly permeable Quaternary gravel and sand beds, which are partly conglomerated (cemented with calcite cement) [43]. The thickness of the sedimentary fill is up to 100 m down to its very low permeable pre-Quaternary basement, composed of Carboniferous-Permian siliciclastic rocks. The groundwater table is on average 25 m below the surface. The unconfined intergranular aquifer with high transmissivity (Figure 2) is mainly recharged from the Sava River in the northwestern part of the study area. It drains into the river in the northwestern part, which largely predisposes the groundwater flow direction (Figure 1). The main water balance components of the regional model area (Figure 1 Table 1.

Nitrogen Mass Balance and Pressure Impact Model
The applied methodology [36,37,40] comprises the following main steps: (1) estimation of the spatial distribution of nitrogen load from the leaky sewer system and agriculture (nitrogen mass balance model), (2) use of the nitrogen load as input into the groundwater transport model, and (3) analysis of the modeling results with observed nitrate concentrations in the observation wells.

Nitrogen Load from the Leaky Sewer System
Assessment of nitrogen load from the leaky sewer system was based on the model of spatially distributed sewer exfiltration developed by Prestor et al. (2011) [36] for the sewer system of the City of Ljubljana. The main steps of workflow for the estimation of nitrogen load from the sewer system are presented in Figure 3. Firstly, the exfiltration rate (losses) from the sewer system was determined by the mass balance of wastewater released to sewer system and outflow from the sewer system. Then, the total mass of pollutant in the sewer system was calculated. In the next step, the exfiltration rate from the tested section of the sewer system was assessed and then upscaled to the spatial distribution of exfiltration of the entire sewer system.

Nitrogen Mass Balance and Pressure Impact Model
The applied methodology [36,37,40] comprises the following main steps: (1) estimation of the spatial distribution of nitrogen load from the leaky sewer system and agriculture (nitrogen mass balance model), (2) use of the nitrogen load as input into the groundwater transport model, and (3) analysis of the modeling results with observed nitrate concentrations in the observation wells.

Nitrogen Load from the Leaky Sewer System
Assessment of nitrogen load from the leaky sewer system was based on the model of spatially distributed sewer exfiltration developed by Prestor et al. (2011) [36] for the sewer system of the City of Ljubljana. The main steps of workflow for the estimation of nitrogen load from the sewer system are presented in Figure 3. Firstly, the exfiltration rate (losses) from the sewer system was determined by the mass balance of wastewater released to sewer system and outflow from the sewer system. Then, the total mass of pollutant in the sewer system was calculated. In the next step, the exfiltration rate from the tested section of the sewer system was assessed and then upscaled to the spatial distribution of exfiltration of the entire sewer system. of Ljubljana. The main steps of workflow for the estimation of nitrogen load from the sewer system are presented in Figure 3. Firstly, the exfiltration rate (losses) from the sewer system was determined by the mass balance of wastewater released to sewer system and outflow from the sewer system. Then, the total mass of pollutant in the sewer system was calculated. In the next step, the exfiltration rate from the tested section of the sewer system was assessed and then upscaled to the spatial distribution of exfiltration of the entire sewer system.

Exfiltration Rate (Losses) from the Sewer System
Daily flow rate measurements at the wastewater treatment plant (WWTP) for the period 2000-2010 were compared with in-house water consumption for that period, provided by the drinking water supply and the wastewater system operator. In the analysis, 13 DWF measurements were used, mainly in the winter period, when rainfall infiltration and outdoor water use have a low influence on the sewer system water balance. Nevertheless, the measured flow rate at the WWTP exceeded the quantity of in-house water consumption, which indicated the presence of extraneous water. Estimation of the extraneous water was based on the assumption that 0.01 km 2 of surface in the catchment of WWTP contributes 0.05 L/s of inflow of water [44], which accounted for the entire catchment 276.6 L/s. On the basis of this, the exfiltration rate from the entire sewer system in the catchment of WWTP was estimated to be 168 L/s (on average 0.21 L/s/km) [40].

Mass of Nitrogen in Wastewater
The estimation of the mass of nitrogen in water is based on the mass balance model [36,37]. The aim of this model was to investigate the shares of known and unknown sources of different pollutants (e.g., nitrogen, phosphorus, chloride, sulphate, copper, chromium, lead, organic carbon) in the sewer system. The mass of pollutants released into sewer system was calculated in the model as a sum of inputs from households, from roads and other traffic surfaces, and from industry. The model results were compared with two extensive analyses of chemical composition of wastewater at WWTP, made in 2008. More frequent measurements were performed only for nitrogen. The concentration of total nitrogen was correlated with discharge and both varied over time depending on the amount of rainfall (Figure 4). The yearly mass of other pollutants discharged at WWTP was a product of discharge and concentrations defined in the two mentioned analyses. Data from the first analysis was used as representative for the dry-weather periods (121 days), and from the second measurement as representative for the wet periods (121 days).
Good agreement between the modeled and observed content of phosphorous (−2.1%), nitrogen (+5.9%), and lead (13.1%) was found [36,37,40]. The content of other pollutants could not be well explained by the known sources. discharge and concentrations defined in the two mentioned analyses. Data from the first analysis was used as representative for the dry-weather periods (121 days), and from the second measurement as representative for the wet periods (121 days).
Good agreement between the modeled and observed content of phosphorous (−2.1%), nitrogen (+5.9%), and lead (13.1%) was found [36,37,40]. The content of other pollutants could not be well explained by the known sources.

Spatial Distribution of Sewer Exfiltration
The model of spatial distribution of sewer exfiltration [36] was based on a detailed inspection of a 14.82 km section of the sewer system, which represents 2% of the entire sewer system. Closed-circuit

Spatial Distribution of Sewer Exfiltration
The model of spatial distribution of sewer exfiltration [36] was based on a detailed inspection of a 14.82 km section of the sewer system, which represents 2% of the entire sewer system. Closed-circuit television (CCTV) provided by the sewer system operator was used to detect the pattern of sewer damages in this section. Among 1246 detected damages, 104 were selected as significant and characterized in more detail to quantify exfiltration. This pattern and information on type of material, age, and size of pipes [45] were used to estimate the spatial distribution of the exfiltration and nitrogen load from the entire sewer system ( Table 2). Table 2. Estimated exfiltration and nitrogen load from pipelines of the entire sewer system [40]. For the calculation of the exfiltration flow rate from the leaky areas, an approach based on Darcy's law was used. It assumes linear dependency between the hydraulic gradient (I), the hydraulic conductivity of the clogging layer (k f ), the leakage area of the pipe (A leak ), and the exfiltration flow rate (Q ex ) [46]:

Time of Construction
where the hydraulic gradient (I) is the ratio between the sum of the water level inside the sewer (h s ) and the pipe wall thickness (h w ), and the pipe wall thickness (h w ): In the calculation, three different hydraulic conditions in the pipelines were considered: (1) dry-weather conditions, days with mainly waste water discharge from water consumers (DWF); (2) days with significant rainfall inflow; and (3) extreme storms events with very high rainfall inflow. An analysis of rainfall data of a 10 year period (2000-2010) was performed. It showed that on average, in 1 year, DWF conditions last for 146 days, significant rainfall inflow occurs in 218 days, and once a year extreme storms are expected. It was assumed that sewer pipes are 10% full in DWF period, 50% in significant inflow period, and 100% in storm events ( Figure 5). These data were used together with data on the diameter of pipes and a representative sample of damages for different materials and ages of pipes to calculate the spatial distribution of yearly exfiltration (Figure 6). and the pipe wall thickness (hw), and the pipe wall thickness (hw): I = (hs + hw) / hw (2) In the calculation, three different hydraulic conditions in the pipelines were considered: (1) dryweather conditions, days with mainly waste water discharge from water consumers (DWF); (2) days with significant rainfall inflow; and (3) extreme storms events with very high rainfall inflow. An analysis of rainfall data of a 10 year period (2000-2010) was performed. It showed that on average, in 1 year, DWF conditions last for 146 days, significant rainfall inflow occurs in 218 days, and once a year extreme storms are expected. It was assumed that sewer pipes are 10% full in DWF period, 50% in significant inflow period, and 100% in storm events ( Figure 5). These data were used together with data on the diameter of pipes and a representative sample of damages for different materials and ages of pipes to calculate the spatial distribution of yearly exfiltration ( Figure 6).

Nitrogen Load from Agriculture
Nitrogen load from agriculture used in the model was the nitrogen surplus calculated as the difference between nitrogen added to an agricultural land (mineral fertilizers, livestock manure,

Nitrogen Load from Agriculture
Nitrogen load from agriculture used in the model was the nitrogen surplus calculated as the difference between nitrogen added to an agricultural land (mineral fertilizers, livestock manure, atmospheric deposition, and mineralization of organic matter) and nitrogen removed from the agricultural land (removed field crops and denitrification in the soil) [41]. On the basis of literature and statistical data, the nitrogen mass balance for different agricultural land use areas (arable fields, vineyards, orchards, and meadows and pastures) was calculated, and the results were provided in the polygon information layer. These data were used to generate a model input grid (Figure 7).

Groundwater Transport Model
Groundwater flow and transport simulations were carried out using the MIKE SHE/MIKE 11 modeling framework [47], which enabled transient simulation of the most important processes of the hydrological cycle in the study area. The solute transport model used in the study [38] is a refined and local version of the calibrated regional hydrological model [48]. The regional model (Figure 1) is a spatially distributed (200 m grid cells) and transient model with implemented daily values of time-dependent variables (temperature, precipitation, river levels and discharges at boundaries, groundwater level observations, abstractions, and leakage from the water supply system). The calibration and validation of the model was based on a comparison of the observed and simulated groundwater heads at the locations of 12 observation wells distributed over the entire study area, and the Sava River levels at one gauging station [48]. Additionally, the solute transport model was validated against the observed concentrations of the trichloroethylene pollution plume discovered in 2004 [49].
The local model (Figure 1) has one calculation layer and at boundaries uses hydraulic conditions from the regional model. Its model domain is discretized into horizontal grid cells of the size 50 × 50 m. The estimated nitrogen load was used as a spatially distributed, subsurface and constant mass input into the saturated zone of the aquifer. Neither retardation nor degradation processes were considered in the simulations. For the nitrate background concentration, 6 mg/L was used. This concentration is comparable to the lowest observed nitrate concentration (7.2 mg/L in observation location 18, Table S1 and is slightly higher than the observed concentration (4.9 mg/L) in the Sava River [50].
Scenarios for a 6 year simulation period (2010-2015) were performed. The last 3 years (2013-2015) of simulations were considered for the analysis. The first half of the simulation period was used for the model warm-up, which assured the avoidance of biased results affected by initial conditions.

Results
The estimated sources of nitrogen from the sewer system and from agriculture are distributed over an area of 19.84 km 2 (34%) and 27.81 km 2 (48% of model area), respectively. The mean nitrogen load in these model cells from the sewer system is 45 g/day and from agriculture is 52 g/day (Table 3). In the local model area of 58 km 2 , this accounts for a total nitrogen of 334.22 ton/year, 38% arising from the leaky sewer system ( Figure 6) and 62% from agriculture (Figure 7). Nitrate concentration in groundwater was available from 58 observation wells and 2 springs. The observed locations were numbered according to their geographical positions, in accelerating order from the west towards the east (Figure 1). The modeled total nitrate concentration (Figure 8) is the sum of the modeled impacts of nitrogen load from the sewer system and agriculture ( Figure S1, Figure S2 Table S1 and summary statistics in Table 4.  Table S1 and summary statistics in Table 4. The average difference between the observed and modeled concentrations at the observation locations was found to be 2.7 mg/L (Table 4, Figure 9). The distribution of differences was left-skewed (Sk = −0.73, Figure 10). The correlation coefficient between the observed and modeled concentrations was 0.51 ( Figure 11). * Sum of the modeled impacts from the sewer system and agriculture with an assumed natural background of 6 mg/L. ** Difference between observed and modeled total concentrations.   * Sum of the modeled impacts from the sewer system and agriculture with an assumed natural background of 6 mg/L. ** Difference between observed and modeled total concentrations.
The average difference between the observed and modeled concentrations at the observation locations was found to be 2.7 mg/L (Table 4, Figure 9). The distribution of differences was left-skewed (Sk = −0.73, Figure 10). The correlation coefficient between the observed and modeled concentrations was 0.51 ( Figure 11).

Discussion
The modeled spatial distribution of nitrate followed the groundwater flow pattern. It reflected the general characteristics of the observed concentrations, which are higher in the southeastern and central parts of the model area ( Figure 8). In spite of higher total nitrogen load from agriculture, the highest nitrate concentrations (>25 mg/L) were observed in areas (obs 2, 3, 4, 6 and 8) with predominant impact from sewer system (Figure 8).
The model gives insight into the impacts of two contamination sources in the capture zones of abstraction wells that are used for drinking water supply. Nitrate concentrations in observation wells in the Šentvid (obs 5 and 7) and Kleče well fields (obs 11-17, 19, 20, 21, 23, 25, 28 and 29) were found to be relatively low (on average below 15 mg/L) and, according to the model, arise predominantly from agriculture (72% of anthropogenic impact). Nitrate concentrations in observation wells in the Hrastje well field (obs 50-56) were higher (on average above 20 mg/L), and the shares of nitrogen load from agriculture (47%) and the sewer system (53%) are comparable. These findings give a bigger picture of the impact of different nitrogen sources on groundwater quality in comparison to the local study, focused on a small part of the urban area [35].
On average, the modeled concentrations were found to be 2.7 mg/L lower than the observed concentrations (Table 4, Figure 9). The discrepancy is relatively low and gives a good estimate of the nitrogen mass balance in the study area. Overestimation of observed nitrate concentrations by a factor of 4.5 was reported in the case study solely by using the nitrogen mass balance estimate for infiltration water [22]. This emphasizes the important contribution of the integration of groundwater dynamics in modeling in order to reduce the uncertainty in the estimate of the nitrate concentration in groundwater.
The modeled concentrations were most evidently underestimated in the area of the highest transmissivity and groundwater discharge (Figure 2). This zone stretches from the northeastern part of the model area, in which the most intensive recharge from the Sava River occurs and further follows the groundwater flow path. This discrepancy could be attributed either to the overestimated

Discussion
The modeled spatial distribution of nitrate followed the groundwater flow pattern. It reflected the general characteristics of the observed concentrations, which are higher in the southeastern and central parts of the model area ( Figure 8). In spite of higher total nitrogen load from agriculture, the highest nitrate concentrations (>25 mg/L) were observed in areas (obs 2, 3, 4, 6 and 8) with predominant impact from sewer system ( Figure 8).
The model gives insight into the impacts of two contamination sources in the capture zones of abstraction wells that are used for drinking water supply. Nitrate concentrations in observation wells in the Šentvid (obs 5 and 7) and Kleče well fields (obs 11-17, 19, 20, 21, 23, 25, 28 and 29) were found to be relatively low (on average below 15 mg/L) and, according to the model, arise predominantly from agriculture (72% of anthropogenic impact). Nitrate concentrations in observation wells in the Hrastje well field (obs 50-56) were higher (on average above 20 mg/L), and the shares of nitrogen load from agriculture (47%) and the sewer system (53%) are comparable. These findings give a bigger picture of the impact of different nitrogen sources on groundwater quality in comparison to the local study, focused on a small part of the urban area [35].
On average, the modeled concentrations were found to be 2.7 mg/L lower than the observed concentrations (Table 4, Figure 9). The discrepancy is relatively low and gives a good estimate of the nitrogen mass balance in the study area. Overestimation of observed nitrate concentrations by a factor of 4.5 was reported in the case study solely by using the nitrogen mass balance estimate for infiltration water [22]. This emphasizes the important contribution of the integration of groundwater dynamics in modeling in order to reduce the uncertainty in the estimate of the nitrate concentration in groundwater.
The modeled concentrations were most evidently underestimated in the area of the highest transmissivity and groundwater discharge (Figure 2). This zone stretches from the northeastern part of the model area, in which the most intensive recharge from the Sava River occurs and further follows the groundwater flow path. This discrepancy could be attributed either to the overestimated groundwater discharge rates or to underestimated nitrogen load, which in this area arises mainly from agriculture. In the urbanized part of the study area, underestimated concentrations could be related to unknown non-agricultural sources of nitrogen [19]. On the contrary, the overestimated nitrate concentrations are mainly in the southern part of the model area. This zone stretches from the southeastern part of the model area, in which the highest concentrations were observed (obs 2, 3 and 4) and follows the groundwater flow path towards the southeast.
The differences between the modeled and observed nitrate concentrations are related to the uncertainties arising from different sources. The level of complexity of used models was adjusted to available data and knowledge. The selection of the models followed the principle as simple as possible, as complex as necessary [51,52]. The used estimates of nitrogen load were based on the limited available data and a simplified nitrogen mass balance model. The wastewater exfiltration model used characteristics of 2% of the sewer system and scales up these findings to the entire sewer system. Due to this generalization and the lack of observations and the complexity of exfiltration process [14], the quantification of impact of sewer leaking on groundwater quality is uncertain. The available quantification of nitrogen load (surpluses) from agriculture was based on literature and statistical data [41], and thus it did not consider site-specific agricultural practices and their changes over time.
The groundwater model is a mathematical representation of an aquifer that is highly productive, characterized by a saturated zone up to 70 m deep and a fast groundwater flow with up to 20 m/day. The simulated pollution plumes are narrow and fast progressing, thus difficult to observe [49]. The modeling of the nitrate concentration is based on simulated groundwater discharge and flow direction. These two processes in the aquifer are difficult to measure and the lack of available data limits the possibility of validating the groundwater transport model. Especially in the areas close to the boundaries of the model, in which the transmissivity of the aquifer (Figure 2) is low, modeled concentrations were found to be very high ( Figure 8) and sensitive to the interpreted shape of the aquifer basement. It is assumed that in these areas, the uncertainty of the model is the highest.
Data on nitrate concentration in groundwater was acquired from three monitoring networks managed by the Slovenian Environmental Agency, the City of Ljubljana, and the water supply company JP VOKA SNAGA. The samples were taken from observation wells using withdrawal method. No multilevel sampling technique for collecting depth-specific groundwater samples was applied. Representative depth-specific groundwater sampling in unconsolidated and coarse-grained aquifers is difficult to obtain. In such aquifers, even the application of packers within a screened well does not assure representative depth-specific sampling, due to the high risk of vertical circulations in the annular space around the well screen or casing [53]. The available data on groundwater sampling were insufficient for characterization of vertical distribution of nitrate in the aquifer. Due to this, the median value of measurements (2013)(2014)(2015) in each well was considered as a representative concentration of the whole groundwater column at the observation location. This fact led to the simplification of the model and the use of a single calculation layer, which limited the simulation of vertical stratification and interpretation of the site-specific vertical concentration distribution.

Conclusions
The presented modeling approach enabled relation and quantification of the impact of pressure from contamination sources on groundwater quality. Modeled spatial distribution of nitrate content in groundwater pointed out the leaky sewer system as the prevailing source of nitrogen load inducing locally high nitrate concentration (above 25 mg/L). Reconstruction of the part of the sewer system in the catchments of observation wells with high nitrate content was recognized as a priority measure for the improvement of groundwater quality. The developed model can be used as a supportive tool in the optimization of mitigation measures, providing the simulation of effects of different remediation scenarios in a cost-benefit analysis. This approach could be applied also to assess the efficiency of groundwater management plans and to analyze the contribution of the stakeholders' activities to the improvement of groundwater quality.
Nitrate concentration in groundwater could be well explained by the known sources of nitrogen that approved the applied integrated mass balance and groundwater modeling concept. Application of this concept to other pollutants revealed a high share of unknown sources for several pollutants; among them, chromium (CrIV) was pointed out as the most threatening for the water supply [40]. Identification and remediation of unknown sources of pollutants whose content in groundwater is high, considering drinking water quality standards, or increases in time, is a priority for groundwater management and the assurance of safe drinking water supply for the City of Ljubljana.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/12/4/1171/s1, Table S1: Median values of modeled and observed nitrate concentrations (2013-2015), Figure S1: The spatial distribution of the modeled nitrate concentration arising from the sewer system; Figure S2: The spatial distribution of the modeled nitrate concentration arising from agriculture.