Tracing nitrate contamination sources and dynamics in an unconfined alluvial aquifer system (Velika Gorica well field, Croatia)

Abstract

Nitrate ions (NO₃⁻) are one of the most common contaminants in the groundwater of the Zagreb alluvial aquifer, which hosts strategic groundwater reserves of the Republic of Croatia and supplies drinking water to one million inhabitants of the capital city. To better understand the origin and the dynamics of NO₃⁻ in the unsaturated and saturated zones, the stable isotopes of nitrogen (δ¹⁵N) and oxygen (δ¹⁸O) in dissolved nitrate, combined with physico-chemical, hydrogeochemical and water stable isotope data, were used in the current work, together with statistical tools and mixing models. The study involved monthly sampling of groundwater, surface water, precipitation and soil water samples. Additionally, the isotopic composition of total nitrogen (δ¹⁵N_bulk) was determined in solid samples representing the local nitrate sources. The combination of a nitrous oxide isotopic analyzer and the titanium(III) reduction method provides reliable measurements of δ¹⁵N_NO3 and δ¹⁸O_NO3, with optimal stability achieved under specific conditions. Nitrate in the study area predominantly originates from organic sources, with nitrification as the main biogeochemical process, while denitrification was identified at sampling sites under specific anaerobic conditions. Although statistical analysis can be a valuable tool, it should be applied with caution if NO₃⁻ originates from multiple sources. The isotopic composition of water showed that groundwater is predominantly recharged by the Sava River but its contribution varied spatially. The results also show the existence of a different recharge source in the southern part of the aquifer. Our findings highlighted the importance of employing a diverse range of analytical methods to obtain reliable and comprehensive understanding of nitrate contamination. By integrating multi-method approaches, stakeholders can better understand the complexities of groundwater contamination and implement more targeted measures to safeguard the water supplies for future generations.

---

Environmental significance

Nitrate (NO₃⁻) is a highly stable and mobile form of nitrogen (N) and one of the more commonly found contaminants in the environment. Tracing NO₃⁻ transformations and sources is essential for better understanding of N cycling and gaining insights into water quality protection. Furthermore, it is critical that we identify and understand the processes affecting N forms in the soil zone and groundwater system. We demonstrate how integration of physico-chemical, hydrogeochemical and isotopic data with statistical tools and mixing models can enhance our understanding of N cycling in the environment and improve water resource protection in regions facing similar environmental challenges.

Introduction

Dissolved nitrate (NO₃⁻), which is one of the most prevalent forms of reactive nitrogen, is widespread in the environment and can originate from natural and/or anthropogenic sources.^1–5 Groundwater nitrate pollution is a global environmental problem^6–12 with increasing levels of NO₃⁻ in many freshwater systems.^2,13,14 The main anthropogenic sources of NO₃⁻ are sewage or industrial wastewaters, livestock manure and the intensive use of fertilizers.^15,16 Apart from anthropogenic NO₃⁻, nitrate is also naturally produced as a result of the decay of soil organic matter, production–fixation and atmospheric deposition.¹⁷ Tracing the sources and transformations of NO₃⁻ in groundwater, as well as in the soil zone, is important for an improved grasp of the nitrogen cycle, water quality protection and sustainable management of the aquifer.^18,19

Nitrogen and oxygen stable isotopes of NO₃⁻ (δ¹⁵N_NO3 and δ¹⁸O_NO3) are powerful tracers for identifying the source(s) of nitrogen (N) contamination, as well as N transformations and dynamics in aquatic and atmospheric systems.^16,20–25 Although nitrate isotopes are a valuable tool for tracking NO₃⁻ contaminants, the nitrate sources may show overlapping isotopic ranges and biogeochemical processes like nitrification and denitrification can complicate the identification of nitrate sources.²⁶ Biogeochemical processes predominantly occur in shallow aquifers, influencing nitrate content and δ¹⁵N_NO3 values. The nitrification process involves the oxidation of NH₄⁺ to NO₃⁻.²⁷ Denitrification generally occurs under anaerobic conditions and involves the reduction of NO₃⁻ to N₂, N₂O or NO.²⁸ Denitrification is a natural attenuation process that reduces NO₃⁻ concentrations and produces a linear relationship in the residual NO₃⁻ with a δ¹⁸O_NO3/δ¹⁵N_NO3 ratio which ranges from 1 : 2.1 to 1 : 1.3.²⁹ To improve the application of nitrate isotope techniques, the use of additional parameters, such as physico-chemical data and water stable isotopes is considered important.

The Zagreb aquifer is one of the most important aquifer systems in Croatia with high to very high vulnerability to contamination as evidenced by the presence of five major contaminants (toxic metals, nitrates, pesticides, pharmaceuticals and chlorinated aliphatics),³⁰ as well as elevated NO₃⁻ concentrations in the groundwater, particularly near industrial and agricultural activity zones.³¹ A preliminary assessment of nitrate contamination origin in the groundwater of the Zagreb alluvial aquifer showed that NO₃⁻ is predominantly of organic origin, particularly from wastewaters.³²

Previous studies in the groundwater of the Zagreb area did not focus on the investigation of nitrate ion distribution and migration through the soil zone (unsaturated zone), nitrate concentrations in the precipitation and Sava River, as well as the biogeochemical and hydrodynamic processes controlling NO₃⁻ content and dynamics in the aquifer. To determine both the origin and the dynamics of nitrate, it is necessary to measure stable isotopes δ¹⁵N and δ¹⁸O in dissolved nitrate in different types of water samples, which were evaluated within this research. However, previous study related to nitrate origin in groundwater of the Zagreb aquifer³² provided initial insights into the nitrate isotopic composition. This new research area was defined based on initial findings, which suggested an additional agricultural source beyond wastewater, primarily on the Sava River's right bank near the Velika Gorica well field. Furthermore, it is also essential to conduct an analysis of the isotopic composition of total nitrogen (δ¹⁵N_bulk) in the local solid samples of different potential nitrate sources if contribution from different nitrate sources has to be evaluated. Therefore, field and laboratory investigations included sampling and analysis of groundwater, surface water, precipitation and soil water. Here, we used physico-chemical, hydrogeochemical and isotopic data, together with statistical tools and mixing models to identify the origin of nitrate contamination and the conditions influencing the N dynamics. This was supported by the first determination of the nitrogen isotopic composition of the local sources around the Velika Gorica well field. There were four main objectives of this research: (1) to optimize laboratory procedures necessary to get reliable results using a nitrous oxide isotopic analyzer and the titanium(III) reduction method; (2) to define nitrate origin and associated biogeochemical processes as well as quantification of spatial and temporal proportional contributions from different N sources; (3) to define appropriate geochemical and statistical methods for interpretation in similar research; (4) to quantify aquifer recharge components. Our study provided new insights into the processes of nitrogen compound transformation in the soil zone, the identification of the dominant form of N entering the unsaturated zone of the aquifer from the soil, determination of the characteristic isotopic composition of N compounds appearing in the surface water, soil water and precipitation, and detailed quantification of nitrate sources which influence nitrate concentrations in the aquifer, with the aim of better understanding the nitrogen cycle to support sustainable management and protection of the Zagreb aquifer.

Methodology

Site description

The Zagreb aquifer is located in the north-western part of the Republic of Croatia within the Sava River catchment and covers an area of about 350 km² (Fig. 1A). It is an unconfined aquifer composed of quaternary sediments consisting mainly of sand, gravel and silt or silty clays. It comprises two hydraulically connected layers: (1) a shallow layer of Holocene alluvial deposits, primarily composed of a carbonate material, transported by the Sava River from the Alps and (2) a deep layer of Pleistocene lacustrine-marshy deposits, primarily composed of a siliciclastic material from the surrounding mountains.³³ The thickness of the aquifer varies from less than 10 meters in the western part up to 100 meters in the eastern part. The unsaturated zone comprises gravel in its lower section, with the upper part predominantly composed of silty to sandy materials, occasionally interspersed with clay layers. The thickness of the unsaturated zone varies from 2 to 11 meters.³⁴ In general, groundwater flows from west/northwest to east/southeast and coincides with the flow direction of the Sava River. The shallow layer is in direct hydraulic connection with the Sava River, which is the main source of groundwater recharge.^35,36 During high water levels, the Sava River recharges the aquifer along its entire course, while during medium and low water levels, it drains the aquifer in certain sections. The study area is characterized by three major pedological units, namely Fluvisols, Stagnosols and Eutric Cambisols on the Holocene deposits.^37,38 According to Ružičić et al.,³⁹ Fluvisols have higher permeability due to their lower clay content and higher sand content, which increases with depth. In contrast, Eutric Cambisols exhibit lower permeability and a slightly higher capacity for water retention. The areas where Fluvisols dominate allow for greater percolation of surface and soil water from precipitation. This is confirmed by the most recent study⁴⁰ in the zone of the Eutric Cambisols, which showed that only high average intensity precipitation events enable infiltration to the unsaturated part of the aquifer. The climate is classified as a humid continental climate, with an annual average precipitation of 967 mm and annual average air temperature of 11.9 °C, with a minimum in January and a maximum in July. The monthly average precipitation is around 80 mm and is evenly distributed throughout the year.⁴¹ For a sampling period from March 2021 to February 2022, the average monthly precipitation was 55.7 mm. The driest and the wettest months were September (29.8 mm) and May (102 mm), respectively. The average air temperature ranged from 1.2 (January) to 23.3 °C (July) with an average of 11.6 °C.⁴²

Fig. 1 (A) The location of the study area with the land use pattern; (B) the water table contour map of the study area during low and high groundwater levels, showing equipotential lines in meters above sea level (m a.s.l.); (C) groundwater flow lines during low and high water levels assuming a maximum pumping capacity (250 L s⁻¹ per well).

The Velika Gorica well field is one of the most important well fields in Zagreb County and is situated in the southern part of the Zagreb aquifer. It supplies drinking water to the capital city and the town of Velika Gorica. The well field consists of five pumping wells and observation wells generally used to monitor groundwater quality in the shallow aquifer. Within the first sanitary protection zone of the Velika Gorica well field, a pedological pit was constructed and equipped with devices for collecting soil water samples. The unsaturated zone thickness at the study site usually ranges from 5 to 8 meters and depends on the groundwater levels. At the top of the unsaturated zone, the following soil horizons are identified according to the World Reference Base classification: A (0–0.15 m), 2B (0.15–0.55 m), 3BC (0.55–0.9 m) and C (0.9–1.17 m). The upper 90 centimetres of the soil profile are predominantly composed of silty and sandy materials, with intermittent clay layers,³⁹ while gravels dominate in the lower part of the unsaturated zone.⁴³ According to Bogunović et al.³⁷ the pedological pit is located in Eutric Cambisols on Holocene deposits. The monitoring wells Vg-6, Lg-1 and Vg-5/2 are located in the area of Fluvisols, while all other monitoring wells are located on Eutric Cambisols on the Holocene deposits.

Based on previous research,^44–46 the groundwater in the wider area of the Velika Gorica well field belongs to the calcium–magnesium–hydrogen carbonate (CaMgHCO₃) water type. The land-uses, which constitute the potential sources of contamination, are presented in Fig. 1A.⁴⁷ The study area consists primarily of urban and agricultural land. The general direction of the groundwater flow in the study area is from north-west to south-east, with clear differences during low and high water level periods (Fig. 1B). During the high water level periods, the groundwater flows towards the south in the northern part of the study area, which is due to the increased influence from the Sava River. In addition, the pumping activities in the Velika Gorica well field also impact the groundwater flow direction near the pumping wells. By developing a groundwater flow model, Posavec⁴⁸ identified a change in the groundwater flow direction during low and high water level periods for a scenario with the maximum pumping capacity (Fig. 1C).

Sampling

Groundwater, surface water, precipitation and soil water were sampled each month from March 2021 to February 2022. Groundwater samples were collected from 10 monitoring wells, located in the catchment area of the Velika Gorica well field (Fig. 1A). Before collecting the groundwater samples, groundwater levels were manually measured while monitoring wells were pumped until pH and electrical conductivity (EC) stabilized to ensure collection of the representative groundwater sample. Surface water samples were collected from the Sava River at the western part of the Zagreb aquifer. The groundwater and Sava River samples were to be analysed for water chemistry, water stable isotopes and nitrate stable isotopes. In situ measurements included pH, EC, dissolved oxygen (DO) content and temperature (T) which were determined in situ using a WTW multi parameter 3630 IDS. Precipitation samples were collected within the first sanitary protection zone of the Velika Gorica well field using a Palmex Rain Sampler RS1 (Zagreb, Croatia) capable of preventing evaporation and appropriate for use in most hydrological studies.^49,50 For the precipitation samples, in situ parameters and nitrate stable isotopes were also determined. Soil water samples were collected from four soil horizons of a pedological pit using soil water samplers (suction cups; Eijkelkamp Soil & Water, Giesbeek, The Netherlands) and an integrated automatic vacuum pump unit AVP-100 (UGT GmbH, Müncheberg, Germany) to determine the isotopic composition of nitrate. The chemical and water isotopic data of soil water and precipitation were previously published.⁴² Each liquid sample was filtered in the field using a 0.22 μm nylon membrane filter to remove microbes, stored in a high-density polyethylene (HDPE) bottle and transferred to the lab for later analysis. The sample bottles were kept refrigerated at 4 °C for chemical and water isotope analysis or frozen at −20 °C for nitrate isotope analysis. During September 2021, ten (10) samples of nitrogen sources were collected in the study area for nitrate isotope analysis: six organic samples and four widely used commercial inorganic fertilizers. The organic samples included two manure samples, two samples from septic tanks and two sewage samples. The selected synthetic fertilizers were NP (nitrogen–phosphorus) 20–20, NPK (nitrogen–phosphorus–potassium) 15–15–15, NPK(S) (nitrogen–phosphorus–potassium–sulphur) 15–15–15(3) and NPK(MgO) (nitrogen–phosphorus–potassium–magnesium oxide) 13–10–12(4), in which most of the nitrogen is present as ammoniacal nitrogen. The isotopic signature of soil was defined based on the sampled soil water from the deepest soil horizon. A table summarizing the type of sample, location, number of sampling points, frequency, sampling period, type of treatment and analyses can be found in the ESI (Table S1†).

Chemical analysis

The major anion and cation concentrations (Cl⁻, NO₃⁻, SO₄²⁻, Na⁺, Mg²⁺, K⁺ and Ca²⁺) were determined with a Dionex ion chromatography system (ICS-90) at the LaGEMA laboratory of the Faculty of Mining, Geology and Petroleum Engineering, University of Zagreb. The ionic balance error was verified for all groundwater samples, with error values below 10% considered acceptable. The bicarbonate concentrations were determined using the titration method.

Water stable isotope analysis and mixing model

The water isotopic analyses (δ²H_H2O and δ¹⁸O_H2O) were carried out at the Laboratory for spectroscopy of the Faculty of Mining, Geology and Petroleum Engineering, University of Zagreb, using a liquid water isotope analyzer (LWIA-45-EP, Los Gatos Research, San Jose, California) by laser absorption spectroscopy. The analytical uncertainty of duplicate samples was ±0.9‰ for δ²H and ±0.19‰ for δ¹⁸O. The results were reported using the delta notation in per mil (‰) relative to Vienna Standard Mean Ocean Water (VSMOW). The data analysis and processing were performed using the Laboratory Information Management System (LIMS) for Lasers 2015.⁵¹

To quantify the contribution of precipitation and the Sava River to the aquifer recharge, a two-component mixing model was utilized which has been used in the study area in previous research,³⁶ but also in different hydrogeology applications.^52–55 In this research it was assumed that the aquifer recharge comes from two main sources, i.e., the Sava River and precipitation. The sum of the end member contributions is expressed as fractions (f) which are equal to 1. For the quantification of recharge, average values from the sampling period for each monitoring well, the Sava River and precipitation were used, based on the following equations:

f_river + f_{precipitation} = 1 (1)

f_river·δ¹⁸O_river + f_{precipitation}·δ¹⁸O_{precipitation} = δ¹⁸O_groundwater (2)

where f_river and f_{precipitation} present Sava River and precipitation fractions respectively, δ¹⁸O_groundwater presents isotopic composition of oxygen in the investigated monitoring well, while δ¹⁸O_river and δ¹⁸O_{precipitation} present isotopic compositions of oxygen in the Sava River and precipitation, respectively.

Nitrate stable isotope analysis

Instrumentation and data processing. The isotope analyses of δ¹⁵N_NO3 and δ¹⁸O_NO3 were performed using an ultrahigh precision advanced nitrous oxide isotopic EP analyzer (GLA451-N2OI3, ABB-LGR, Quebec, Canada) at the Laboratory for spectroscopy of the Faculty of Mining, Geology and Petroleum Engineering, University of Zagreb. The analyzer precision was 0.05 ppb for N₂O concentration, better than 1‰ for δ¹⁵N and better than 2‰ for δ¹⁸O. To convert dissolved NO₃⁻ to N₂O gas headspace the titanium(III) reduction method was used, first introduced and described in detail by Altabet et al.⁵⁶ After the analysis, the data was pre-processed by using an in-house MATLAB⁵⁷ script. This MATLAB script is designed to process multiple TXT files containing raw laser measurements by iterating throughout laser output TXT files and calculating mean values for essential data. This streamlined approach facilitates effective data processing and management, enabling further analysis. Subsequently, a new CSV file is created for each input file, containing the calculated mean values. A GRG nonlinear solving method in Excel Solver™⁵⁸ was employed to correct the raw instrumental data for isotopologues of interest (¹⁵N and ¹⁸O) to remove any concentration dependence in the δ values.

Testing and validation. To ensure laboratory setup functionality and reliability before research activities, optimization steps were conducted, i.e., target N₂O gas concentration inside the analyzer, reagent-to-sample ratio during sample preparation and reaction time. The testing of the optimal N₂O target that produced the most stable isotopic values was conducted using a control sample KNO₃ with different final NO₃–N sample concentrations injected into the analyzer at various injection volumes, ensuring that a broad range of N₂O concentrations were covered. N₂O concentrations ranging from 7 to 11 ppm N₂O exhibited the most stable isotopic values (δ¹⁵N and δ¹⁸O) with optimal precision (Fig. 2A). Instrumental δ¹⁵N and δ¹⁸O values were more precise using a 1 : 30 reagent-to-sample ratio with a mean s.d. of 1.9‰ and 1.7‰, respectively (Fig. 2B). The observed calibration slopes for δ¹⁵N were close to the theoretical value of 1 regardless of which reagent-to-sample ratio was used during sample preparation. However, the calibration plots generated with a reagent-to-sample ratio of 1 : 20 gave a higher slope for δ¹⁸O, i.e., 0.9 (Fig. 2C). The δ¹⁵N values were more accurate at a 1 : 40 reagent-to-sample ratio with a mean bias of −0.07‰, while a ratio of 1 : 20 resulted in a higher mean bias of 0.3‰. Conversely, the δ¹⁸O normalised values were most accurate at a 1 : 20 ratio with a mean autorun bias of 0.05‰ (Fig. 2D). Therefore, a 1 : 20 Ti(III) reagent-to-sample ratio is considered practical for all further measurements. Altabet et al.⁵⁶ noted that a 96-hour reaction yielded similar outcomes to a 21-hour reaction, suggesting that samples can be stored for at least 3 days before analysis without adverse effects. Our study yielded slightly different results (Fig. 2E), as both the 48-h and 72-h reaction yielded isotopic results comparable to those of a 24-h reaction, while the control sample left for longer than 72 hours did not yield acceptable results.

Fig. 2 Optimization steps taken to ensure laboratory setup functionality and reliability before initiating research activities include: (A) establishing the target N₂O gas concentration; (B)–(D) determining the reagent-to-sample ratio during the preparation method; (E) establishing the reaction time for optimal results. Note: USGS32 in (B) and the control sample in (E) are presented using the secondary axis.

Standards and samples. Groundwater and surface water nitrate samples were prepared with the Ti(III) reduction method to achieve a NO₃⁻ sample concentration of 0.55 mg_NO3–N L⁻¹ in each vial in order to get the target N₂O concentration inside the analyzer using a reagent-to-sample volume ratio of 1 : 20. Due to insufficient NO₃⁻ levels in two samples, the stable isotope composition of the nitrate was investigated for 10 precipitation samples. Additionally, due to very small sample volume amounts available for the Ti(III) preparation method and insufficient NO₃⁻ levels, the stable isotope composition of the NO₃⁻ was analysed for eight soil water samples. Due to the lower concentration of NO₃⁻, the precipitation and soil water samples were prepared to achieve a NO₃⁻ sample concentration of 0.1 mg_NO3–N L⁻¹.

In alignment with Altabet et al.⁵⁶ and in order to avoid non-representative results, an overnight reaction was chosen to be convenient and practical, allowing samples to be prepared on day 1 and analysed the following morning (day 2). Nitrate reference materials (USGS32, USGS34 and USGS35) were used for data normalisation to the AIR and VSMOW scales by applying least squares regression to the measured versus consensus δ values. Nitrate standards have consensus δ¹⁵N_AIR values of +180‰ (exactly), −1.8‰ (±0.1‰) and +2.7‰ (±0.1‰) and δ¹⁸O_VSMOW values of +25.7‰ (±0.2‰), −27.9‰ (±0.3‰) and +57.5‰ (±0.3‰), respectively.^59,60 Normalized results were least squares regression mean values ±1-sigma.

Solid nitrate sources. Ten solid samples of potential nitrate sources were dried and ground to homogenize them. The isotopic composition of total nitrogen (δ¹⁵N_bulk) was investigated by an external laboratory (The Faculty of Earth Sciences of the University of Barcelona) by using an Elemental Analyzer (EA) Flash IRMS coupled in continuous flow to an Isotope Ratio Mass Spectrometer (IRMS) Delta V Advantage with a universal interphase Conflo IV (CF) (EA IsoLink CN IRMS System, Thermo Scientific). Isotope ratios are reported in parts per mil (‰) relative to the international standard V-AIR for total nitrogen (δ¹⁵N_bulk). The precision was <±0.5%.

Statistical data analysis

The normality of 17 parameters (pH, EC, DO, T, Cl⁻, NO₃⁻, SO₄²⁻, HCO₃⁻, Na⁺, Mg²⁺, K⁺, Ca²⁺, δ¹⁵N_NO3, δ¹⁸O_NO3, δ²H_H2O, δ¹⁸O_H2O and d-excess) in the groundwater samples was assessed using the Shapiro–Wilk test. To evaluate temporal and spatial differences, normally distributed data were analysed using one-way ANOVA at a significance level of p-value <0.05. The test evaluated changes over time (from March 2021 to February 2022) and changes among sampling sites. If the normality assumption was not met, the non-parametric Kruskal–Wallis test was applied.

Correlation tests were conducted on the monitoring well data to examine the relationship between Cl⁻, NO₃⁻, SO₄²⁻, HCO₃⁻, Na⁺, Mg²⁺, K⁺ and Ca²⁺ in groundwater, using a significance level of p-value <0.05. All statistical analyses were performed using the TIBCO Software Inc. Statistica (Version 13.5.0.17).

Bayesian mixing model and uncertainly analysis

The proportional contribution of 3 nitrogen sources in the groundwater and the river was estimated via the application of a Bayesian isotope mixing model using the package simmr in R (Stable Isotope Mixing Model in R).^61,62 In brief, the model is expressed by using the following equations (eqn (1)–(4);⁶³):

S_ij ∼ N(μ_jk, ω_jk²) (4)

C_jk ∼ N(λ_jk, τ_jk²) (5)

ε_jk ∼ N(0, σ_j²) (6)

where X_ij is the isotope value j of the mixture i, in which i = 1, 2, 3,…, N and j = 1, 2, 3,…, J; p_k is the proportion of source k (k = 1, 2, 3,…, K) estimated by the model; S_jk is the source value k on isotope j normally distributed with mean μ_jk and standard deviation ω_jk; c_jk is the isotope fractionation factor for isotope j on source k normally distributed with mean λ_jk and standard deviation τ_jk; and ε_ij is the residual error of the additional unquantified variation between individual mixtures normally distributed with mean 0 and standard deviation σ_j.

The model was applied at each sampling site to identify spatial differences in proportional contributions, as well as across different seasons for each water type (groundwater and river). The isotopic average and standard deviation values of the 3 local nitrogen sources in the study area (soil, synthetic fertilizers and organic wastes) were used in the analysis and are shown in Table 1, while isotopic ranges are shown in Fig. 6 and were the following: soil (δ¹⁵N from −7.1 to −4.9‰, δ¹⁸O from −9.1 to −4.4‰); synthetic fertilizers (δ¹⁵N from 0.5 to 1.5‰, δ¹⁸O from −9.6 to −8.6‰); organic wastes from sewage, septic waste and manure (δ¹⁵N from 1.9 to 8.8‰, δ¹⁸O from −9.6 to −8.6‰); and precipitation (δ¹⁵N from −15.1 to 2.2‰, δ¹⁸O from 62.1 to 94.7‰). The δ¹⁵N values of the local solid sources were within expected ranges.⁶⁴

Table 1 Isotopic average values of the local nitrogen sources used in the mixing model

	δ ^15N (‰)	δ ^18O (‰)
Soil (N = 3)	−5.8 ± 0.9	−6.3 ± 2.0
Synthetic fertilizers (N = 4)	1.1 ± 0.4	−9.0 ± 0.3
Organic wastes (sewage, septic waste and manure) (N = 6)	3.7 ± 2.3	−9.0 ± 0.3

---

We excluded precipitation from the analysis, as the N isotope precipitation signal is quickly dampened.⁶⁴ For the sampling sites (i.e., Lg-1 and Lg-4) where denitrification was identified, since they fall within the optimal denitrification zone due to low DO levels³ and have characteristic slopes associated with denitrification,²⁹ an enrichment factor was calculated using the following equation:²⁷

δ ∼ δ_o + ε ln(f) (7)

where δ_o is the initial isotopic composition of the substrate, δ is the isotopic composition of the substrate, ε is the enrichment factor and f is the remaining fraction of the substrate. The average δ¹⁵N_NO3 value of the local sources (+1.5‰) was considered to be the initial isotopic composition.

Synthetic fertilizers were mainly NH₄⁺ based, which allowed determining only the δ¹⁵N value, whereas the δ¹⁸O value of this endmember was determined by considering the δ¹⁸O_H2O values of the local groundwater via nitrification. In general, the nitrification of NH₄⁺ in soil, precipitation or fertilizers utilizes two oxygens from water and one from atmospheric oxygen.^27,65 Thus, nitrification values may be estimated using the following experimental equation:^65,66

The soil signature was defined based on measured δ¹⁵N_NO3 and δ¹⁸O_NO3 data from soil water samples collected from the deepest soil horizon (C horizon) that was closest to the water table.

An uncertainty analysis of isotopic composition for each nitrogen source was conducted using the probability statistical method introduced by Ji et al.⁶⁷ In summary, the proportional contribution values from the simmr model for each potential nitrogen source were assigned frequency values based on the number of iterations. A plot of the cumulative frequency distribution versus the proportional contributions of the potential nitrogen sources was then created. The uncertainty index (UI90) for each nitrogen source was calculated by finding the difference between the proportional contributions at the 0.95 and 0.05 cumulative frequency distributions and then dividing this difference by 0.9.⁶⁸

Results and discussion

Evaluation of nitrate origin and the processes affecting nitrogen species concentrations is complex and requires comprehensive analysis. Effective data interpretation necessitates the integration of various analytical approaches and perspectives. This approach is particularly crucial in areas where nitrate concentrations arise from a combination of anthropogenic sources and natural environmental conditions.^1–4 Diverse data and methods enhance reliability, with isotopic signatures and the composition of main nitrate sources being critical. The results are discussed through hydrogeochemistry, statistical analysis, isotopic composition and Bayesian mixing models, contextualized for both urban and agricultural environments.

Hydrogeochemical aspects

The pH values ranged from 6.9 to 7.8 for groundwater, from 7.8 to 8.7 for Sava River water and from 6.5 to 8.1 for precipitation. The groundwater EC varies spatially from 675 to 938 μs cm⁻¹ with the lowest EC values observed at Lg-1 and Vg-10/2 and the highest at Vg-4. Surface water EC varies between 353 and 461 μs cm⁻¹, while precipitation EC ranged from 12.1 to 54.5 μs cm⁻¹. The DO content of the groundwater ranged from 0.6 to 9.0 mg_O2 L⁻¹; in the Sava River it varied between 7.1 and 14.1 mg_O2 L⁻¹, while in the precipitation it varied between 7.2 and 13.9 mg_O2 L⁻¹. The lowest groundwater DO values were observed at Lg-1 and Lg-4 sites with a mean value of 3.4 mg_O2 L⁻¹ and 3.6 mg_O2 L⁻¹, respectively. Groundwater temperatures varied between 12.4 and 15.0 °C, while river water temperatures ranged from 7.2 to 25.0 °C. Higher temperatures were recorded at Vg-4 and Vg-6 sites, while the lowest temperatures were recorded at the Vg-10/2 location.

The NO₃⁻ concentrations in groundwater ranged from 7.4 to 31.9 mg_NO3 L⁻¹ (Table S2†) with an average value of 17.1 mg_NO3 L⁻¹. None of the analysis exceeded the maximum permissible level for drinking water (50 mg_NO3 L⁻¹) as stipulated by EU and Croatian regulations. Elevated mean NO₃⁻ concentrations were found at Vg-10/2 (20.3 mg_NO3 L⁻¹) and Vg-5/2 monitoring wells (25.1 mg_NO3 L⁻¹), which was linked to fertilizers and wastewater, respectively. The NO₃⁻ concentrations in the Sava River ranged from 3.5 to 6.5 mg_NO3 L⁻¹. The lower NO₃⁻ concentrations in the Sava River compared to groundwater could be because of biogeochemical processes, such as assimilation (NO₃⁻ consumption by plants) and dilution with rain water. Given that the ambient background value of NO₃⁻ of the study area is around 7 mg_NO3 L⁻¹ (ref. 41) and the expected NO₃⁻ content in the soil zone is similar,⁴² the observed nitrate concentrations in the groundwater cannot be attributed only to natural conditions, i.e., sources. This is evident as the lowest recorded concentrations were around 9 mg_NO3 L⁻¹ in all groundwater samples apart from the Lg-4 site in August, where the NO₃⁻ concentration was 7.4 mg_NO3 L⁻¹. Most NO₃⁻ concentrations generally exceeded natural background levels.

Considering that monitoring wells Lg-1 and Lg-4 generally had DO levels lower than 4 mg_O2 L⁻¹, they were in the zone optimal for the denitrification process (Fig. 3A). Conversely, all other groundwater sampling sites had DO levels greater than 4 mg_O2 L⁻¹ and a pH range between 6.5 and 8, indicating that the measured NO₃⁻ concentrations could be a result of nitrification under aerobic conditions.³ The Sava River has a slightly higher pH values, with a maximum recorded value of 8.7 in January 2022.

Fig. 3 (A) Scatterplot of pH versus DO, with a dashed line separating the zone of nitrification and zone of denitrification; (B) TDS versus (NO₃⁻ + Cl⁻)/HCO₃⁻; (C) relationship between NO₃⁻/Na⁺ and the Cl⁻/Na⁺ molar ratio; (D) SO₄²⁻/Na⁺versus the NO₃⁻/Na⁺ molar ratio; (E) relationship of the NO₃⁻/Cl⁻ molar ratio and Cl⁻ molar concentration; (F) relationship between the NO₃⁻/K⁺ ratio and K⁺.

The lowest mean SO₄²⁻ concentrations were observed in Lg-1 (13.5 mg L⁻¹), P-7 (13.1 mg L⁻¹) and Vg-10/2 (14.0 mg L⁻¹) sites. Conversely, the Lg-4 and Vg-1 sites showed the highest content of sulphate with an average value of 23.1 and 23.3 mg L⁻¹, respectively. The content of SO₄²⁻ in the Sava River ranged from 6.3 to 17.0 mg L⁻¹.

The Cl⁻ concentrations in groundwater varied from 7.6 to 71.6 mg L⁻¹, with an average of 26.9 mg L⁻¹. The levels of Cl⁻ in the Sava River were significantly lower, from 5.9 to 12.4 mg L⁻¹. The highest Cl⁻ levels were observed at the monitoring wells Vg-4 and Vg-6, with average values of 55.6 mg L⁻¹ and 37.2 mg L⁻¹, respectively. In contrast, the lowest average Cl⁻ concentrations were found at Lg-1 (11.1 mg L⁻¹) and Vg-10/2 (10.1 mg L⁻¹) sites.

The concentration of K⁺ in groundwater ranged from 0.9 to 6.0 mg L⁻¹ with an average of 2.2 mg L⁻¹. The highest concentrations of K⁺ were detected at the monitoring wells Vg-5/2 and Vg-6 with average values of 4.8 and 3.7 mg L⁻¹, respectively. The K⁺ concentration in the Sava River varied between 1.1 and 3.7 mg L⁻¹.

The concentration of Na⁺ in groundwater ranged from 4.8 to 37.2 mg L⁻¹ with an average of 13.4 mg L⁻¹. The levels of Na⁺ in the Sava River were lower, from 4.6 to 10.1 mg L⁻¹. The highest concentrations of Na⁺ were observed at the sampling sites Vg-4 and Vg-6 with average values of 29.1 and 22.2 mg L⁻¹, respectively. The increase in Na⁺ concentration is due to sewage infiltration, as Vg-4 and Vg-6 are wells located in urban areas. In contrast, the lowest Na⁺ contents were recorded at Lg-1 and Vg-10/2 sites. Considering that these two locations had higher concentrations of NO₃⁻, the excess NO₃⁻ is related to agricultural or urban inputs (Fig. 3D).

As depicted in Fig. 3B the positive strong correlation between TDS and (NO₃⁻ + Cl⁻)/HCO₃⁻, R² = 0.64 and R² = 0.84 for the monitoring wells Vg-6 and P-7 respectively, indicated that the water chemistry at these locations is influenced by human activities, such as sewage infiltration or agricultural activities.⁶⁹ When considering the molar ratios of NO₃⁻/Na⁺ and Cl⁻/Na⁺ (Fig. 3C), the anthropogenic influence on the groundwater quality is evident in almost all groundwater samples.

The Cl⁻ content and NO₃⁻/Cl⁻ ratios have been used to evaluate the sources of NO₃⁻.⁷⁰ Here (Fig. 3E), the Vg-4 and Vg-6 sites exhibited high Cl⁻ content and low NO₃⁻/Cl⁻ ratios, indicating a contribution from sewage and/or manure, i.e., from organic contamination. Vg-6 is the only monitoring well that contains elevated Cl⁻ concentrations suggesting contamination from sewage⁷¹ and at the same time shows a correlation between Cl⁻ and NO₃⁻ suggesting contamination from domestic wastewater discharge.⁷² No monitoring well showed low Cl⁻ concentrations and a low NO₃⁻/Cl⁻ ratio. The samples from Lg-1 and Vg-10/2 sites, which showed low Cl⁻ content and high NO₃⁻/Cl⁻ ratios, could be associated with agricultural inputs.⁷³ The remaining six sampling sites did not fall within the range of any specific potential NO₃⁻ input, suggesting that NO₃⁻ was derived from a mixing between different sources.

As shown in Fig. 3F, a low NO₃⁻/K⁺ molar ratio for samples from sites Vg-6 and Vg-5/2, as well as for Sava River samples, confirms possible contamination by sewage effluent^32,74 from the urban area. Synthetic fertilizers and manure often have higher NO₃⁻ concentrations relative to K⁺. Thus, the high NO₃⁻/K⁺ molar ratios observed at the remaining sampling locations (Čp-23, Lg-4, P-7, Vg-10/2, Lg-1, P-3, Vg-1 and Vg-4) indicated multiple contamination inputs affecting NO₃⁻ and K⁺ contents.

The normality test showed significant differences between groundwater sites for all tested parameters (p < 0.05). However, no significant (p > 0.05) seasonal differences were observed for DO, EC, T, Cl⁻, NO₃⁻, Na⁺, δ²H_H2O and δ¹⁸O_H2O. The correlation coefficients provided a preliminary understanding of how different variables relate to nitrate sources. Detailed correlation matrices illustrating these relationships can be found in Tables S3 to S12.† For all monitoring wells, except for Vg-6, a significant and very strong positive correlation between Cl⁻ and SO₄²⁻ was observed indicating that these concentrations may be related to sewage disposal,^71,75 or agricultural activities, excretions of migrating livestock and domestic wastewater discharge.⁷² A significant and strong positive correlation between NO₃⁻ and Cl⁻ was observed for the sampling sites Vg-6 (r = 0.80) and Vg-5/2 (r = 0.69) but not between NO₃⁻ and other variables. This is probably due to contamination from septic tanks or wastewater discharge.⁷² Moreover, the sampling site P-7 showed a significant and strong positive correlation between NO₃⁻ and Cl⁻, K⁺, Na⁺ and SO₄²⁻ (r = 0.89, r = 0.89, r = 0.89 and r = 0.91, respectively) suggesting that most of these ions come from the same source of contamination, although this method does not allow the precise identification of the source. All other locations showed insignificant correlation between NO₃⁻ and K⁺ indicating that nitrate concentrations are related to the different non-point sources, for example manure, sewage or synthetic fertilizers.^76,77 Sampling site Vg-10/2 showed a significant and strong positive correlation between NO₃⁻ and Cl⁻ (r = 0.79) and SO₄²⁻ (r = 0.87) suggesting that these ions are likely influenced by similar sources. No significant correlation between NO₃⁻ and other ions was detected for the sampling sites Čp-23, Lg-1, Lg-4, Vg-1 and Vg-4. The monitoring well P-7 showed significant and strong positive correlation between K⁺ and Cl⁻ (r = 0.89), as well as between Na⁺ and Cl⁻ (r = 0.85), which could be related to salinization due to the application of KCl synthetic fertilizers.⁷⁸ Furthermore, the positive correlation between Na⁺ and Cl⁻ at P-7 is indicative of contamination from wastewater, while no correlation for all other sampling sites indicates the existence of natural sodium concentrations, probably due to rock weathering and cation exchange with Ca²⁺.⁷⁹ The strong positive correlation between NO₃⁻ and SO₄² for P-7 and Vg-10/2 sites agrees with the findings of Spoelstra et al.⁸⁰ who suggest that a significant proportion of the groundwater SO₄²⁻ originates from fertilizers. In general, SO₄²⁻ ions can originate from industrial products, fertilizers, and precipitation, but also from natural sulphide mineral dissolution.⁸¹ The correlation coefficients for the P-7 location are attributed to the intensified application of synthetic fertilizers.

Isotopic composition of water

The δ²H_H2O and δ¹⁸O_H2O values for the groundwater samples ranged from −65.0 to −57.5‰ and from −9.6 to −8.6‰, respectively. For the surface waters, the δ²H_H2O and δ¹⁸O_H2O values ranged from −61.8 to −59.2‰ and from −9.4 to −8.9‰, respectively. As documented earlier,⁴² the local meteoric water line (LMWL) of Velika Gorica is characterized by a linear evaporation trend described by δ²H_H2O = 7.7 δ¹⁸O_H2O + 9.3 and the average δ¹⁸O_H2O value of precipitation is −8.54‰ obtained from the Velika Gorica meteorological station. Additionally, the article also notes that the two shallowest soil horizons share isotopic signatures similar to precipitation, whereas the two deepest soil horizons have comparable isotopic compositions to each other but differ from those of precipitation. For both groundwater and surface water samples, the δ¹⁸O_H2Oversus δ²H_H2O average values were plotted close to the LMWL of Velika Gorica (Fig. 4) indicating their meteoric origin. The current work showed that the groundwater from seven sampling sites, namely, Čp-23, Lg-4, P-3, P-7, Vg-1, Vg-4 and Vg-5/2 and the Sava River exhibited similar isotopic compositions. Conversely, the monitoring well Vg-6 showed an isotopic composition similar to precipitation and the two shallowest soil horizons, while the Lg-1 and Vg-10/2 sites suggest the existence of a different recharge source.

Fig. 4 The relationship between δ¹⁸O_H2O and δ²H_H2O average values in groundwater, surface water, precipitation and soil water. LMWL is the local meteoric water line of Velika Gorica. The water isotopic composition of soil water and precipitation was published earlier in Buškulić et al.⁴²

Similar concentrations between sites Lg-1, Vg-10/2 and the Sava River do not imply a direct connection or common origin, as the δ¹⁸O_H2O and δ²H_H20 values of the Lg-1 and Vg-10/2 sites differed from those of the Sava River, indicating that these wells are primarily recharged by water inflow from the southwest part of the aquifer system (Fig. 1C). The average δ¹⁸O_H2O values for Lg-1 and Vg-10/2 were −9.3‰ and −9.5‰, respectively. For all other sampling sites, the values ranged from −9.0 to −8.7‰.

The results of the mixing model showed that the contribution of the Sava River to aquifer recharge spatially varied from 19.9% to 83.4%, with an average of 58.5%, while the recharge from precipitation varied from 16.6% to 80.1%, with an average of 41.5%. Additionally, the monitoring wells can be grouped as follows: Vg-6 showed the smallest fraction of the Sava River (19.9%); in Čp-23, P-3, Vg-4 and Vg-5/2, the Sava River fraction ranged from 47.1% to 54.2%; and in Lg-4, P-7 and Vg-1 group, the Sava River fraction ranged from 66.7% to 83.3%. The monitoring well Lg-4 showed the highest fraction of Sava River influence (83.4%). Conversely, for the sampling site Vg-6, precipitation represents the most important source of recharge, evidenced by the highest fraction of precipitation of 80.1‰. This was expected, given that Vg-6 is located in an area characterized by Fluvisols, which are characterized by a higher infiltration rate.³⁹

Isotopic composition of nitrate

The δ¹⁵N_NO3 values in the groundwater ranged from 2.6 to 18.3‰, while δ¹⁸O_NO3 values ranged between −6.7 and 9.7‰. The highest δ¹⁵N_NO3 and δ¹⁸O_NO3 values were recorded in the Lg-4 site in August 2021, with values of 38.9 and 17.6‰, respectively, which could be an outlier. Overall, the δ¹⁵N_NO3 values ranged from 4.4 to 13.1‰ for the Sava River, from −15.1 to 2.2‰ for precipitation and from −7.1 to 9.0‰ for soil water. The δ¹⁸O_NO3 values for the Sava River ranged from −0.8 to 3.9‰, between 62.1 and 94.7‰ for precipitation and between −12.7 and 10.1‰ for soil water.

The DO levels in eight monitoring wells (Čp-23, P-3, P-7, Vg-1, Vg-10/2, Vg-4, Vg-5/2 and Vg-6) ranged from 5.6 to 9.0 mg_O2 L⁻¹ indicating aerobic conditions that would favor nitrification, while denitrification is less likely to occur (Fig. 5A). In the Sava River, during the sampling period, the DO was measured in the range from 7.1 to 14.2 mg_O2 L⁻¹, which favors the nitrification process. Denitrification was detected in sampling site Lg-4 since it fell in the optimal zone for denitrification (Fig. 3A) and exhibited reductive conditions favourable for denitrification.³ Additionally, the site exhibited an increasing DO trend as δ¹⁵N_NO3 values decrease (Fig. 5B).

Fig. 5 (A) Scatterplot showing the δ¹⁵N_NO3 and DO values of groundwater and surface water samples; (B) the δ¹⁵N_NO3versus DO values and trendline for the Lg-4 sampling site; (C) the δ¹⁵N_NO3 and NO₃⁻ values of groundwater and surface water samples; (D) the δ¹⁵N_NO3versus NO₃⁻ values and trendlines for Lg-4, P-3 and soil zone samples; (E) the δ¹⁵N_NO3versus ln(NO₃⁻) values of groundwater and surface water samples; (F) the δ¹⁸O_NO3versus δ¹⁸O_H2O where the uppermost line indicates the limit of exchange with O₂, the middle line represents the limit of the nitrification process and the bottom line shows the limit of exchange with H₂O.

For almost all groundwater sampling locations and the Sava River, the δ¹⁵N_NO3 values were relatively constant with only minor variations as the NO₃⁻ concentration changes, which is characteristic of the dilution effect (Fig. 5C), suggesting mixing of nitrate from different sources with similar δ¹⁵N_NO3 signatures. Conversely, based on specific trends in δ¹⁵N_NO3 isotopes during NO₃⁻ attenuation, denitrification processes were identified in the monitoring wells Lg-4 and P-3 (Fig. 5D). This is indicated by an increase in δ¹⁵N_NO3 values along with a decrease in NO₃⁻ concentrations. Furthermore, when examining the relationship between δ¹⁵N_NO3 and the logarithmic concentration of NO₃⁻ (Fig. 5E), samples from Lg-4 and P-3 exhibit a negative slope with moderate correlations (R² = 0.66 and R² = 0.54, respectively) also suggesting natural attenuation, i.e., denitrification.⁷³ According to Böttcher et al.,²⁹ the slopes of δ¹⁵N_NO3versus δ¹⁸O_NO3 for Lg-4 and P-3 were characteristic of denitrification, i.e., 0.6 and 0.4, respectively. However, the sampling site P-3 showed high DO levels (7.6 to 8.3 mg_O2 L⁻¹) suggesting that denitrification is unlikely to occur and cannot be explained under such conditions. Conversely, the low DO levels in the monitoring well Lg-1 and a slope of 0.8 confirmed the occurrence of the denitrification process.²⁹ When observing the relationship between δ¹⁵N_NO3 and ln(NO₃⁻), sampling sites Čp-23, P-7, Vg-1, Vg-4 and Vg-6 exhibit a negative slope with low correlations (R² = 0.004, R² = 0.03, R² = 0.19, R² = 0.15 and R² = 0.29, respectively) (Fig. 5E) indicating a mixing process.⁷³

Nitrate consumption, i.e. the assimilation process, can be observed in the soil zone (Fig. 5D) through changes in the isotope signal. The deepest (C) horizon contains the highest NO₃⁻ concentrations and the lowest measured δ¹⁵N_NO3 values, whereas the shallowest (A) horizon shows the lowest NO₃⁻ concentrations and the highest δ¹⁵N_NO3 content. This suggests NO₃⁻ uptake by plants as the water infiltrates the deeper horizons.

The nitrification process was identified by using the local δ¹⁸O_H2O values to calculate the theoretical δ¹⁸O_NO3 values after applying eqn (8) and comparing them with the observed δ¹⁸O_NO3 values (Fig. 5F). The estimation of δ¹⁸O_NO3 values characteristic for the nitrification in the groundwater was performed using the δ¹⁸O_O2 value of atmospheric O₂ (23.5‰)²⁰ and the minimum and maximum δ¹⁸O_H2O values of groundwater (−9.6‰ and −8.6‰).

The calculated δ¹⁸O_NO3 values ranged from 1.5 to 2.1‰ (grey area shown in Fig. 6). Almost all groundwater samples fell between the line of the nitrification process and the water exchange line (Fig. 5F) confirming the occurrence of nitrification, except for monitoring well Lg-4, where all the isotopic signatures were higher than those theoretically calculated, confirming the denitrification process, i.e., nitrate attenuation.⁸² Similar to groundwater, the estimation of δ¹⁸O_NO3 values characteristic of nitrification in the river was conducted using the minimum and maximum δ¹⁸O_H2O values of −9.4‰ and −8.9‰, respectively. The resulting δ¹⁸O_NO3 values ranged from 1.6 to 1.9‰, which were quite similar to δ¹⁸O_NO3 content observed in the Sava River, where δ¹⁸O_NO3 values range from −0.8 to 3.9‰. As shown in Fig. 5F, the Sava River samples fell along the nitrification line. The estimation of δ¹⁸O_NO3 content in the soil zone was carried out using δ¹⁸O_H2O values of −10.7‰ and −4.2‰, i.e. minimum and maximum observed values. The calculated δ¹⁸O_NO3 values ranged from 0.7 to 5.0‰, while the observed δ¹⁸O_NO3 values ranged from −12.7 to 10.1‰. As depicted in Fig. 5F, soil samples from the first two horizons (A and 2B horizons) and the deepest horizon (C horizon) fell between the nitrification line and the water exchange line, thereby confirming the occurrence of nitrification. In contrast, the sample from the third soil horizon, i.e., the 3BC horizon, deviates from this range, precluding confirmation of nitrification for that particular sample. Considering that the investigated soil profile is situated within the first sanitary protection zone of the Velika Gorica well field, nitrification of NH₄⁺ from soil and/or precipitation could occur in the two upper soil horizons. Conversely, in the deepest horizon, nitrification is restricted to NH₄⁺ present in the soil due to the lack of isotopic signatures associated with precipitation.

Fig. 6 Bivariate plot of δ¹⁸O_NO3versus δ¹⁵N_NO3 values of groundwater, surface water, precipitation and soil water in the study area. The range of isotope compositions for five major N sources is derived from Kendall et al.⁶⁴ and indicated by boxes: (1) atmospheric NO₃⁻; (2) NO₃⁻ based fertilizers; (3) NH₄⁺ in fertilizers and precipitation; (4) soil NH₄⁺; and (5) manure and septic waste. Note: the brown, grey, green and blue boxes represent local N sources: soil, synthetic fertilizers, organic wastes and precipitation, respectively.

Specifically, the samples taken from sampling site Lg-4 indicated both anaerobic (March to September 2021) and aerobic conditions (from October 2021 to February 2022). During the anaerobic period, the average values of groundwater depth, DO, EC, NO₃⁻, δ¹⁵N_NO3, δ¹⁸O_NO3, δ²H_H2O and δ¹⁸O_H2O were 7.1 m, 1.9 mg_O2 L⁻¹, 777.6 μS cm⁻¹, 14.0 mg L⁻¹, 18.5‰, 8.0‰, −61.6‰ and −9.1‰, respectively. Nitrite (NO₂⁻) values during April to September ranged from 0.04 to 3.6 mg L⁻¹ but were below the detection limit in March. The phosphate ions ranged from 1.0 to 2.4 mg L⁻¹ and were below the detection limit in September. The ammonium ions were detected in May and June with an average value of 0.4 mg L⁻¹. Under aerobic conditions, the average values of groundwater depth, DO, EC, NO₃⁻, δ²H_H2O and δ¹⁸O_H2O were 7.8 m, 5.9 mg_O2 L⁻¹, 815.2 μS cm⁻¹, 17.9 mg L⁻¹, −60.3‰ and −8.9‰, respectively, significantly higher than those in the anaerobic period. The average values of δ¹⁵N_NO3 and δ¹⁸O_NO3 were lower, at 12.2‰ and 1.9‰, respectively. In summary, during high water periods, decreased concentrations of DO and NO₃⁻ were observed, whereas the concentrations of NO₂⁻, PO₄³⁻, NH₄⁺, δ¹⁵N_NO3 and δ¹⁸O_NO3 were higher. Conversely, during low water periods, elevated concentrations of DO and NO₃⁻ were recorded, alongside decreased levels of δ¹⁵N_NO3 and δ¹⁸O_NO3, and nearly undetectable levels of NO₂⁻, PO₄³⁻ and NH₄⁺ contents.

Based on the dual isotopic approach, typical δ¹⁵N_NO3 and δ¹⁸O_NO3 values from various N sources⁶⁴ and possible biogeochemical processes are shown in Fig. 6 together with the groundwater, surface water, precipitation and soil water samples. The δ¹⁸O_NO3 values observed in the groundwater and river samples were lower than the values reported for NO₃⁻ based fertilizers (ranging from 17 to 25‰) and atmospheric deposition of NO₃⁻ (ranging from 28 to 98‰) suggesting no contribution from these two sources. The lack of atmospheric NO₃⁻ deposition supports the statement that nitrate concentrations in the groundwater cannot be solely attributed to natural sources. Moreover, the δ¹⁸O_NO3 values were within the range for NH₄⁺ in fertilizers and precipitation, soil NH₄⁺ and manure/septic waste (ranging from −17 to 15‰), as well as around the estimated theoretical range of nitrification. Among the total number of groundwater and Sava River samples, 66% of groundwater samples and 42% of Sava River samples exhibited δ¹⁵N_NO3 isotope signatures higher than 9‰, indicating contamination predominantly originating from manure and/or septic waste. The remaining samples fell within the overlapping range reported for NH₄⁺ in fertilizers and precipitation and soil NH₄⁺.

Given that the investigated soil profile is located within the sanitary protection zone of the well field, where anthropogenic activities are forbidden and strictly controlled, the only possible sources of nitrates are soil NH₄⁺, NH₄⁺ in rain and atmospheric NO₃⁻. According to Kendall et al.,⁶⁴ the nitrate sources in the first two shallowest horizons (A and 2B) could be soil NH₄⁺ and/or NH₄⁺ in precipitation (Fig. 6). Since previous studies show that the two deepest soil horizons (3BC and C) are less influenced by precipitation,^40,42 the negative δ¹⁵N_NO3 isotopic signature is probably the result of nitrification, i.e., soil NH₄⁺ nitrification. As depicted in Fig. 5D, NO₃⁻ concentrations increase with depth while δ¹⁵N_NO3 decreases, which could be due to the nitrification process. The isotopic values in the deepest soil horizon confirmed the NH₄⁺ nitrification and supported a previous finding documented in Buškulić et al.⁴² that N primarily migrates to deeper parts of the unsaturated zone in the form of NO₃⁻.

Bayesian mixing model results and uncertainly analysis

A Bayesian isotope mixing model was employed to estimate spatial and seasonal variations in the proportional contributions of three local nitrogen sources to both the groundwater and river water.

Regarding spatial variations, the model outputs showed a similar pattern across all sampling sites (Fig. 7). The results for the Sava River are similar to those for groundwater, due to their continuous mixing. Across the entire study area, soil contributions ranged from 10 to 17%, making it the least significant source which is consistent with the results obtained using other approaches. In contrast, organic wastes showed the highest contribution, ranging from 56 to 68%, in agreement with Fig. 6 which illustrates that more than half of the groundwater samples indicate contamination due to leakage from the sewage system in urban areas, septic tank leakage in rural areas and/or the application of manure. The proportional contribution of synthetic fertilizers was obtained for each location and partly aligned with the results from the hydrogeochemical approach that detected the influence of synthetic fertilizers only on individual sampling sites (e.g., P-7). This is not a surprise given that nitrate concentrations represent the result of mixing of individual nitrate pollution sources, which makes it difficult to disentangle by using only hydrogeochemical parameters. The synthetic fertilizers contributed moderately, with a range of 21 to 27%. Thus, the contributions of NO₃⁻ sources can be ranked as follows: organic wastes > synthetic fertilizers > soil.

Fig. 7 Spatial variations in the proportional contributions of different local N sources per sampling site based on Bayesian isotope mixing model results. Note: the different local N sources are shown in different colours; the average proportional contributions are shown in red; boxplots illustrate the 25th, 50th and 75th percentiles.

When the model results are examined in the context of land use (urban and agricultural areas) and the average values from the corresponding sampling sites, a clear difference of the proportional contribution of organic waste between these two areas becomes apparent, but the presence of synthetic fertilizers in both areas is not very different. The average value of the organic waste in the urban area is 57%, while in the agricultural area it is 64%.

Regarding seasonal changes, there are noticeable variations in proportions for groundwater, as well as between groundwater and surface water (Fig. 8). The highest contributions of soil and synthetic fertilizers in groundwater were observed in the summer period (9% and 16%, respectively), while the lowest were in the spring (5% and 11%, respectively). Conversely, the proportional contribution of organic wastes is the lowest during summer (75%) and highest during spring (84%). For surface water, the model outputs showed a significant difference compared to groundwater, but a consistent pattern: soil contributions ranged from 26 to 30%, synthetic fertilizers ranged from 32 to 33% and organic wastes ranged from 38 to 41%. Such a significant difference compared to groundwater is attributed to the fact that only three values per season were available for Sava River water.

Fig. 8 Seasonal variations in the proportional contributions of different local N sources per water type based on Bayesian isotope mixing model results. Note: the different local N sources are depicted in distinct colours; the average proportional contributions are shown in red; boxplots illustrate the 25th, 50th and 75th percentile; gw stands for groundwater and sw stands for surface water (Sava River).

The results of the uncertainty analysis, i.e., the plots of the cumulative frequency distribution vs. the proportional contribution, of source apportionment per site are given in Fig. S1.† When all sites were considered, the average UI90 index per nitrate source followed the order: organic wastes (0.71) > NH₄⁺ in synthetic fertilizers and rain (0.63) > soil (0.36). Overall, the uncertainty analysis of the simmr model showed UI90 values comparable with those calculated in other case studies (e.g., ref. 24). However, the higher uncertainty in the first sources (organic wastes and NH₄⁺ in synthetic fertilizers and rain) is probably due to the fact that their isotopic signatures were relatively close in terms of δ¹⁵N_NO3 and δ¹⁸O_NO3. The average UI90 index per season (Fig. S2A†) was lower for the groundwaters (<0.5) showing the following order: summer (0.40) > winter (0.31) > autumn (0.28) > spring (0.18) but higher for the Sava River (>0.7) for all seasons (Fig. S2B†). This is because the amount of data for the Sava River used in the model was smaller (n = 12) than the groundwater (n = 119), which increased the uncertainty of the analysis for the surface waters.

Overall remarks

To summarize, the physico-chemical and hydrogeochemical data showed that groundwater and surface water quality is deteriorated by human activities, primarily due to urban or/and agricultural inputs. The concentrations of major anions and cations in the Sava River are lower compared to groundwater. However, certain sites (Lg-1 and Vg-10/2) contain Na⁺, SO₄²⁻ and Cl⁻ concentrations similar to those in the Sava River. Although these sites have similar concentrations, they are located furthest from the Sava River, suggesting that their primary recharge source may differ.

The statistical analysis revealed significant spatial differences in nitrate content, as well as strong correlations between NO₃⁻ and various ions, suggesting multiple contamination sources.

Furthermore, the nitrate isotopic composition indicated that NO₃⁻ contamination of groundwater and the Sava River during the study period mainly originated from three distinct nitrate sources: soil NH₄⁺, NH₄⁺ in fertilizers and manure/septic waste. This is especially important because previous research³² indicated that the nitrate isotopic composition suggested that nitrates originated solely from manure and wastewater sources. This was unusual due to the presence of a large proportion of agricultural areas, particularly on the right bank of the Sava River. The sites Čp-23, Lg-1, P-3, P-7, Vg-5/2 and the Sava River had δ¹⁵N_NO3 and δ¹⁸O_NO3 values which indicated contamination from NH₄⁺ based fertilizers, soil and manure/septic wastes, while the data from the other four locations (Vg-1, Vg-10/2, Vg-4 and Vg-6) had typical values that indicated contamination from soil and manure/septic waste. Furthermore, the δ¹⁵N_NO3 and δ¹⁸O_NO3 values from Lg-4 suggested contamination solely from manure and septic waste, as indicated by δ¹⁵N_NO3 values higher than 9%. The statistical analyses confirmed significant variations in δ¹⁵N_NO3 and δ¹⁸O_NO3 across groundwater samples from different locations (p < 0.05).

The results of the Bayesian model agree with findings from a previous study³² which stated that nitrate origin in the Zagreb aquifer area is predominantly organic. However, findings from this research for the first time cannot exclude synthetic fertilizers as the N source and possible cause of elevated NO₃⁻ concentrations. Given that synthetic fertilizers could not be identified as a potential nitrogen source at each sampling site using the hydrogeochemistry and δ¹⁵N_NO3vs. δ¹⁸O_NO3 bivariate plot, the mixing model results suggest that the biogeochemical process, most likely nitrification, masked the nitrogen isotopic signature typically associated with synthetic fertilizers and indicated the importance of employing multi-method approaches for more accurate N source identification.

Our results showed that separate usage of each method cannot give a clear understanding. In dynamic systems such as unconfined aquifers, where there are various sources of contamination, it is necessary to use and combine numerous types of data and methods, as well as perceive results from different aspects, if more reliable conclusions about nitrate origin and processes that affect them need to be achieved.

Additionally, it was not possible to consider other aspects, such as the inherent delay of N transport through the unsaturated and saturated zone,⁸³ in this work. This means that nitrate ions in groundwater may still be present even after years of the reduction or cease of fertilization rates.⁸⁴ This delay in the transport of water and nutrients in groundwater highly depends on the local hydraulic properties and heterogeneity of the unsaturated zone and the rates of biogeochemical processes.⁸⁵

Conclusions

This study presents a comprehensive integration of physico-chemical and hydrogeochemical data, isotopic analysis, statistical correlations and mixing models to elucidate the complex sources and processes affecting NO₃⁻ contamination in groundwater of an unconfined aquifer. The use of a combination of water stable isotopes alongside stable isotopes from NO₃⁻ dissolved in water, provides a unique way of identifying both natural and anthropogenic influences on groundwater quality, particularly in relation to agricultural and wastewater inputs. With respect to the main objectives of this research, the following can be concluded:

(1) A nitrous oxide isotopic analyser and the titanium(III) reduction method present a reliable combination for the determination of δ¹⁵N and δ¹⁸O from nitrate. It was shown that the results can vary based on the target N₂O gas concentration inside the analyser, reagent-to-sample ratio during sample preparation and reaction time. N₂O concentrations ranging from 7 to 11 ppm N₂O, 1 : 20 Ti(III) reagent-to-sample ratio and maximum 72-h reaction are found to generate the most stable isotopic values.

(2) The origin of nitrate in the research area is mainly related to the organic sources and the process of nitrification prevails in the most observed locations. Additionally, for the first time, the isotopic composition of local nitrate sources has been determined, which showed that the contribution from soil and synthetic fertilizer exists and varies from both a spatial and temporal point of view. From this perspective, a more detailed inspection of NO₃⁻ sources in future research must be performed in the study area with a focus on the evaluation of the isotopic composition of synthetic fertilizers, especially those which contain both ammonium and nitrate nitrogen in their formula. It must be noted that determination of the isotopic composition of NO₃⁻ sources is a crucial part from an interpretation point of view because it shows the range and variability of N which is introduced into the investigated system. The results of this research also showed that it is extremely important to monitor isotopic values in the soil horizons, especially in the deepest one, in order to get information about dominant nitrogen species which can enter the aquifer.

(3) Multitude of analysis is necessary when determination of nitrate origin and related processes presents the main subject of the research. Although correlation analysis was found to be useful in numerous previous research studies, in the areas where NO₃⁻ concentrations are the consequence of different sources it must be considered with great caution. Even though they can be indicative and helpful within the interpretation procedure, their usage without the isotopic composition should be avoided. Similar insights can be obtained from the hydrogeochemical data when observing ion ratios, but again not unambiguous. The initial insight into the nitrate isotopic composition clearly suggested that most of the nitrate concentrations were associated with manure and organic waste, which turned out to be correct, but at the same time did not recognize the influence of synthetic fertilizer. This is probably associated with the process of nitrification, which masked the nitrogen isotopic signature typically associated with synthetic fertilizers.

(4) The Sava River is the main source of recharge of the Zagreb aquifer, but the results also show that in the southern part of the aquifer groundwater flows from the south-western part of the aquifer and is recharged from the surrounding hills.

The significance of this work is underlined by its ability to distinguish between different sources of contamination and processes such as nitrification and denitrification in a specific hydrogeological context. This approach not only advances the understanding of nitrate dynamics in groundwater but also provides valuable insights for the management and protection of water resources in regions facing similar environmental challenges.

Data availability

Data will be made available on reasonable request. The MATLAB script used to preprocess data is available at the GITHUB repository: https://github.com/hrvojelukacic.

Author contributions

Patricia Buškulić: conceptualization, methodology, formal analysis, investigation, data curation, writing – original draft, writing – review & editing, visualization. Zoran Kovač: formal analysis, investigation, writing – review & editing. Ioannis Matiatos: formal analysis, investigation, writing – review & editing, visualization. Jelena Parlov: writing – review & editing, visualization, supervision. All authors have read and agreed to the published version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The work was part of the Young Researches' Career Development Project – Training New Doctoral Students (DOK-2020-01) supported by the Croatian Science Foundation (HRZZ) and IAEA TC project CRO7002 “The use of nitrogen and oxygen stable isotopes in the determination of nitrate origin in the unsaturated and saturated zone of the Velika Gorica well field”. The authors would like to thank Cedric Douence from IAEA, Travis B. Medor from Biology Centre CAS, and Francois Brisson and Robert Provencal from ABB. Thanks also go to Borna-Ivan Balaž, Saša Šipek, Hrvoje Čiček, Hrvoje Lukačić, Michaela Hruškova, Vinko Baranašić and Branka Prša from the Faculty of Mining, Geology and Petroleum Engineering, University of Zagreb.

References

1. D. Xue, J. Botte, B. De Baets, F. Accoe, A. Nestler, P. Taylor, O. Van Cleemput, M. Berglund and P. Boeckx, Present limitations and future prospects of stable isotope methods for nitrate source identification in surface- and groundwater, Water Res., 2009, 43(5), 1159–1170.
2. C. Fenech, L. Rock, K. Nolan, J. Tobin and A. Morrissey, The potential for a suite of isotope and chemical markers to differentiate sources of nitrate contamination: a review, Water Res., 2012, 46(7), 2023–2041.
3. O. Nikolenko, A. Jurado, A. V. Borges, K. Knöller and S. Brouyère, Isotopic composition of nitrogen species in groundwater under agricultural areas: a review, Sci. Total Environ., 2018, 621, 1415–1432.
4. B. Linhoff, Deciphering natural and anthropogenic nitrate and recharge sources in arid region groundwater, Sci. Total Environ., 2022, 848, 157345.
5. B. Mayer and I. Matiatos, Nutrient Dynamics in Rivers and Lakes, in Treatise on Geochemistry, ed. A. Anbar, D. Weis, Elsevier, ISBN 9780323997638, 3rd edn, 2024, pp 155–178, https://www.sciencedirect.com/science/article/abs/pii/B9780323997621000772?via%3Dihub.
6. J. Rockström, W. Steffen, K. Noone, Å. Persson, F. S. Chapin, E. F. Lambin, T. M. Lenton, M. Scheffer, C. Folke and H. J. Schellnhuber, et al., A Safe Operation Space for Humanity, Nature, 2009, 461, 472–475.
7. J. W. Erisman, J. N. Galloway, S. Seitzinger, A. Bleeker, N. B. Dise, A. M. Roxana Petrescu, A. M. Leach and W. de Vries, Consequences of Human Modification of the Global Nitrogen Cycle, Philos. Trans. R. Soc., B, 2013, 368,  DOI:10.1098/rstb.2013.0116.
8. X. Xia, S. Zhang, S. Li, L. Zhang, G. Wang, L. Zhang, J. Wang and Z. Li, The cycle of nitrogen in river systems: sources, transformation, and flux, Environ. Sci.:Processes Impacts, 2018, 20(6), 863–891.
9. A. Romanelli, D. X. Soto, I. Matiatos, D. E. Martínez and S. Esquius, A Biological and Nitrate Isotopic Assessment Framework to Understand Eutrophication in Aquatic Ecosystems, Sci. Total Environ., 2020, 715, 136909.
10. I. Matiatos, L. I. Wassenaar, L. R. Monteiro, J. J. Venkiteswaran, D. C. Gooddy, P. Boeckx, E. Sacchi, F. J. Yue, G. Michalski, C. Alonso-Hernández and C. Biasi, et al., Global patterns of nitrate isotope composition in rivers and adjacent aquifers reveal reactive nitrogen cascading, Commun. Earth Environ., 2021, 2(1), 52.
11. G. Balacco, G. D. Fiorese and M. R. Alfio, Assessment of groundwater nitrate pollution using the Indicator Kriging approach, Groundwater Sustainable Dev., 2023, 21, 100920.
12. R. N. Tabi, A. Gibrilla, P. Boakye, F. O. Agyemang, A. A. Foaah and S. Oduro-Kwarteng, Appraisal of groundwater quality and hydrochemistry in three regions of Ghana: implications for drinking purposes, Groundwater Sustainable Dev., 2024, 26, 101193.
13. E. Abascal, L. Gómez-Coma, I. Ortiz and A. Ortiz, Global Diagnosis of Nitrate Pollution in Groundwater and Review of Removal Technologies, Sci. Total Environ., 2022, 810, 152233.
14. M. J. Ascott, D. C. Gooddy, B. Marchant, N. Kieboom, H. Bray and S. Gomes, Regional scale evaluation of nitrate fluctuations in groundwater using cluster analysis and standardised hydrometeorological indices, J. Hydrol., 2024, 634, 131052.
15. G. Tokazhanov, E. Ramazanova, S. Hamid, S. Bae and W. Lee, Advances in the catalytic reduction of nitrate by metallic catalysts for high efficiency and N₂ selectivity: a review, Chem. Eng. J., 2020, 384, 123252.
16. R. Hussien, M. Ahmed and A. I. Aly, Tracking anthropogenic nitrogen-compound sources of surface and groundwater in southwestern Nile Delta: hydrochemical, environmental isotopes, and modeling approach, Environ. Sci. Pollut. Res., 2023, 30, 22115–22136.
17. J. A. Torres-Martínez, A. Mora, P. S. K. Knappett, N. Ornelas-Soto and J. Mahlknecht, Tracking nitrate and sulfate sources in groundwater of an urbanized valley using a multi-tracer approach combined with a Bayesian isotope mixing model, Water Res., 2020, 182, 115962.
18. H. Zhao, Q. Xiao, Y. Miao, Z. Wang and Q. Wang, Sources and transformations of nitrate constrained by nitrate isotopes and Bayesian model in karst surface water, Guilin, Southwest China, Environ. Sci. Pollut. Res., 2020, 27, 21299–21310.
19. I. Matiatos, L. R. Monteiro, M. Sebilo, D. X. Soto, D. C. Gooddy and L. I. Wassenaar, Isotopes reveal the moderating role of ammonium on global riverine water nitrogen cycling, ACS ES&T Water, 2024, 4(4), 1451–1459.
20. R. Aravena, M. L. Evans and J. A. Cherry, Stable isotopes of oxygen and nitrogen in source identification of nitrate from septic systems, Ground Water, 1993, 31(2), 180–186.
21. L. I. Wassenaar, Evaluation of the origin and fate of nitrate in the Abbotsford aquifer using the isotopes of ¹⁵N and ¹⁸O in NO₃, Appl. Geochem., 1995, 10(4), 391–405.
22. M. Zhang, Y. Zhi, J. Shi and L. Wu, Apportionment and uncertainty analysis of nitrate sources based on the dual isotope approach and a Bayesian isotope mixing model at the watershed scale, Sci. Total Environ., 2018, 639, 1175–1187.
23. G. Esquivel-Hernández, I. Matiatos, R. Sánchez-Murillo, Y. Vystavna, R. Balestrini, N. S. Wells, L. R. Monteiro, S. Chantara, W. Walters and L. I. Wassenaar, Nitrate isotopes (δ¹⁵N, δ¹⁸O) in precipitation: best practices from an international coordinated research project, Isot. Environ. Health Stud., 2023, 59(2), 127–141.
24. I. Matiatos, K. Lazogiannis, A. Papadopoulos, N. T. Skoulikidis, P. Boeckx and E. Dimitriou, Stable isotopes reveal organic nitrogen pollution and cycling from point and non-point sources in a heavily cultivated (agricultural) Mediterranean river basin, Sci. Total Environ., 2023, 901, 166455.
25. M. Rotiroti, E. Sacchi, M. Caschetto, C. Zanotti, L. Fumagalli, M. Biasibetti, T. Bonomi and B. Leoni, Groundwater and surface water nitrate pollution in an intensively irrigated system: sources, dynamics and adaptation to climate change, J. Hydrol., 2023, 623, 129868.
26. H. Jung, D. C. Koh, Y. S. Kim, S. W. Jeen and J. Lee, Stable Isotopes of Water and Nitrate for the Identification of Groundwater Flowpaths: A Review, Water, 2020, 12, 138.
27. C. Kendall, in Isotope Tracers in Catchment Hydrology, Elsevier, 1998, vol. 16, pp. 519–576.
28. J. A. Brandes and A. H. Devol, Isotopic fractionation of oxygen and nitrogen in coastal marine sediments, Geochim. Cosmochim. Acta, 1997, 61(9), 1793–1801.
29. J. Böttcher, O. Strebel, S. Voerkelius and H.-L. Schmidt, Using isotope fractionation of nitrate-nitrogen and nitrate-oxygen for evaluation of microbial denitrification in a sandy aquifer, J. Hydrol., 1990, 114(3), 413–424.
30. Z. Nakić, S. Ružičić, K. Posavec, M. Mileusnić, J. Parlov, A. Bačani and G. Durn, Conceptual model for groundwater status and risk assessment-case study of the Zagreb aquifer system, Geol. Croat., 2013, 66, 55–77.
31. Z. Kovač, M. Cvetković and J. Parlov, Gaussian simulation of nitrate concentration distribution in the Zagreb aquifer, J. Maps, 2017, 13, 727–732.
32. Z. Kovač, Z. Nakić, J. Barešić and J. Parlov, Nitrate Origin in the Zagreb Aquifer System, Geofluids, 2018, 2018, 2789691.
33. J. Velić and G. Durn, Alternating Lacustrine-Marsh Sedimentation and Subaerial Exposure Phases during Quaternary: Prečko, Zagreb, Croatia, Geol. Croat., 1993, 46, 71–90.
34. S. Ružičić, M. Mileusnić and K. Posavec, Building Conceptual and Mathematical Model for Water Flow and Solute Transport in the Unsaturated Zone at Kosnica Site, Rud.-Geol.-Naftni Zb., 2012, 25, 21–31.
35. K. Posavec, P. Vukojević, M. Ratkaj and T. Bedeniković, Cross-correlation Modelling of Surface Water – Groundwater Interaction Using the Excel Spreadsheet Application, Rud.-Geol.-Naftni Zb., 2017, 32, 25–32.
36. J. Parlov, Z. Kovač, Z. Nakić and J. Barešić, Using Water Stable Isotopes for Identifying Groundwater Recharge Sources of the Unconfined Alluvial Zagreb Aquifer (Croatia), Water, 2019, 11(10), 2177.
37. M. Bogunović, Ž. Vidaček, S. Husnjak and M. Sraka, Inventory of Soils in Croatia, Agric. Conspec. Sci., 1998, 63, 105–112.
38. S. Ružičić, PhD thesis, University of Zagreb, 2013.
39. S. Ružičić, Z. Kovač, D. Perković, L. Bačani and Lj. Majhen, The Relationship between the Physicochemical Properties and Permeability of the Fluvisols and Eutric Cambisols in the Zagreb Aquifer, Croatia, Geosciences, 2019, 9, 416.
40. Z. Kovač, L. Bačani, S. Ružičić, J. Parlov, K. Posavec and P. Buškulić, Using Water Stable Isotopes and Cross-Correlation Analysis to Characterize Infiltration of Precipitation through Unsaturated Zone at the Velika Gorica Site of Zagreb Aquifer, J. Hydrol. Eng., 2023, 28, 04023002.
41. P. Buškulić, J. Parlov, Z. Kovač and Z. Nakić, Estimation of Nitrate Background Value in Groundwater under the Long-Term Human Impact, Hydrology, 2023, 10(3), 63.
42. P. Buškulić, J. Parlov, Z. Kovač, T. Brenko and M. Pejić, Determination of Nitrate Migration and Distribution through Eutric Cambisols in an Area without Anthropogenic Sources of Nitrate (Velika Gorica Well Field, Croatia), Sustainability, 2023, 15, 16529.
43. A. Kukolja, Master's thesis, University of Zagreb, 2018.
44. Z. Nakić, S. Horvat and A. Bačani, Statistical Indicators of Groundwater Geochemical Characteristics in a Quaternary Aquifer from the Mala Mlaka Well Field Catchment Area (Zagreb, Croatia), Geol. Croat., 2005, 58, 87–99.
45. T. Vlahović, A. Bačani and K. Posavec, Hydrogeochemical Stratification of the Unconfined Samobor Aquifer (Zagreb, Croatia), Environ. Geol., 2009, 57, 1707–1722.
46. P. Buškulić and J. Parlov, Mathematical Methods and Terminology in Geology 2022, University of Zagreb, Zagreb, 2022, pp. 179–186.
47. CLC 2012, https://www.haop.hr/hr/baze-i-portali/pokrov-i-namjena-koristenja-zemljista-corine-land-cover, accessed August 2024.
48. Posavec 2016, https://www.vgvodoopskrba.hr/media/txadqokj/elaborat_zone_zastite_izvorista_vg_2016-1.pdf, accessed July 2024.
49. M. Gröning, H. O. Lutz, Z. Roller-Lutz, M. Kralik, L. Gourcy and L. Pöltenstein, A simple rain collector preventing water re-evaporation dedicated for δ¹⁸O and δ²H analysis of cumulative precipitation samples, J. Hydrol., 2012, 448–449, 195–200.
50. N. Michelsen, R. van Geldern, Y. Roßmann, I. Bauer, S. Schulz, J. A. C. Barth and C. Schüth, Comparison of Precipitation Collectors Used in Isotope Hydrology, Chem. Geol., 2018, 488, 171–179.
51. T. B. Coplen and L. I. Wassenaar, LIMS for Lasers 2015 for Achieving Long-Term Accuracy and Precision of δ²H, δ¹⁷O, and δ¹⁸O of Waters Using Laser Absorption Spectrometry, Rapid Commun. Mass Spectrom., 2015, 29, 2122–2130.
52. J. Jean-Baptiste, C. Le Gal La Salle and P. Verdoux, Use of mixing models to explain groundwater quality time and space variation in a narrowed fluctuating alluvial aquifer, Appl. Geochem., 2020, 121, 104700.
53. A. Behrouj-Peelya, Z. Mohammadia, L. Scheiberb and E. Vázquez-Suñé, An integrated approach to estimate the mixing ratios in a karst system under different hydrogeological conditions, J. Hydrol., 2020, 30, 100693.
54. L. Tian, Y. Gao, G. Yang, B. Schwartz, B. Cai, C. Ray, Y. Li and H. Wu, Isotopic tracers of sources of water for springs from the Edwards Aquifer, Central Texas, USA, Hydrol. Res., 2021, 52(3), 787–803.
55. X. Jian, S. Zhang, Q. Lao, F. Chen, P. Huang, C. Chen and Q. Zhu, Using dual water isotopes to quantify the mixing of water masses in the Pearl River Estuary and the adjacent northern South China Sea, Front. Mar. Sci., 2022, 9, 987685.
56. M. A. Altabet, L. I. Wassenaar, C. Douence and R. Roy, A Ti(III) reduction method for one-step conversion of seawater and freshwater nitrate into N₂O for stable isotopic analysis of ¹⁵N/¹⁴N, ¹⁸O/¹⁶O and ¹⁷O/¹⁶O, Rapid Commun. Mass Spectrom., 2019, 33(15), 1227–1239.
57. MathWorks Inc., MATLAB, Natick, Massachusetts, United States, 2023.
58. L. I. Wassenaar, C. Douence, M. A. Altabet and P. K. Aggarwal, N and O isotope (δ¹⁵N^α, δ¹⁵N^β, δ¹⁸O, δ¹⁷O) analyses of dissolved NO₃⁻ and NO₂⁻ by the Cd-azide reduction method and N₂O laser spectrometry, Rapid Commun. Mass Spectrom., 2018, 32(3), 184–194.
59. J. K. Böhlke and T. B. Coplen, in Reference and Intercomparison Materials for Stable Isotopes of Light Elements, International Atomic Energy Agency, Vienna, 1995, pp. 51–66.
60. J. K. Böhlke, S. J. Mroczkowski and T. B. Coplen, Oxygen isotopes in nitrate: new reference materials for ¹⁸O:¹⁷O:¹⁶O measurements and observations on nitrate-water equilibration, Rapid Commun. Mass Spectrom., 2003, 17, 1835–1846.
61. A. C. Parnell, D. L. Phillips, S. Bearhop, B. X. Semmens, E. J. Ward, J. W. Moore, A. L. Jackson, J. Grey, D. J. Kelly and R. Inger, Bayesian stable isotope mixing models, Environmetrics, 2013, 24(6), 387–399.
62. A. Parnell and M. A. Parnell, Package ‘simmr’, 2019, in Aquaculture Big Numbers, Food And Agriculture Organization Of The United Nations, Rome, Italy, 2016.
63. A. C. Parnell, R. Inger, S. Bearhop and A. L. Jackson, Source partitioning using stable isotopes: coping with too much variation, PLoS One, 2010, 5(3), e9672.
64. C. Kendall, E. M. Elliott and S. D. Wankel, in Stable Isotopes in Ecology and Environmental Science, Blackwell Science Publications, Oxford, 2007, vol. 12, pp. 375–449.
65. B. Mayer, S. M. Bollwerk, T. Mansfeldt, B. Hütter and J. Veizer, The oxygen isotope composition of nitrate generated by nitrification in acid forest floors, Geochim. Cosmochim. Acta, 2001, 65, 2743–2756.
66. K. K. Andersson and A. B. Hooper, O₂ and H₂O are each the source of one O in NO₂⁻ produced from NH₃ by Nitrosomonas: ¹⁵N-NMR evidence, FEBS Lett., 1983, 164(2), 236–240.
67. W. Ji, L. Shu, W. Chen, Z. Chen, X. Shang, Y. Yang, R. A. Dahlgren and M. Zhang, Nitrate pollution source apportionment, uncertainty and sensitivity analysis across a rural-urban river network based on δ¹⁵N/δ¹⁸O-NO₃⁻ isotopes and SIAR modeling, J. Hazard. Mater., 2022, 438, 129480.
68. X. Ji, R. Xie, Y. Hao and J. Lu, Quantitative identification of nitrate pollution sources and uncertainty analysis based on dual isotope approach in an agricultural watershed, Environ. Pollut., 2017, 229, 586–594.
69. Y. Xiao, X. Gu, S. Yin, X. Pan, J. Shao and Y. Cui, Investigation of geochemical characteristics and controlling processes of groundwater in a typical long-term reclaimed water use area, Water, 2017, 9, 800.
70. D. Widory, E. Petelet-Giraud, P. Negrel and B. Ladouche, Tracking the sources of nitrate in groundwater using coupled nitrogen and boron isotopes: a synthesis, Environ. Sci. Technol., 2005, 39, 539e548.
71. J. A. Tindall, K. J. Lull and N. G. Gaggiani, Effects of land disposal of municipal sewage sludge on fate of nitrates in soil, streambed sediment, and water quality, J. Hydrol., 1994, 163(1–2), 147–185.
72. T. L. L. J. Nenkam, R. Kringel, W. Y. Fantong, P. Nbendah, A. F. Takoundjou, Z. Elisabeth and B. T. Kamtchueng, Hydrochemistry of nutrients in groundwater under farmland in the Benue River Basin, North-Cameroon, Environ. Earth Sci., 2022, 81, 209.
73. J. A. Torres-Martínez, A. Mora, J. Mahlknecht, L. W. Daessle, P. A. Cervantes-Aviles and R. Ledesma-Ruiz, Estimation of nitrate pollution sources and transformations in groundwater of an intensive livestock-agricultural area (Comarca Lagunera), combining major ions, stable isotopes and MixSIAR model, Environ. Pollut., 2021, 269, 115445.
74. S.-L. Li, C.-Q. Liu, Y.-C. Lang, Z.-Q. Zhao and Z.-H. Zhou, Tracing the sources of nitrate in karstic groundwater in Zunyi, Southwest China: a combined nitrogen isotope and water chemistry approach, Environ. Earth Sci., 2010, 60, 1415–1423.
75. G. Venkatesan and G. Swaminathan, Review of chloride and sulphate attenuation in ground water nearby solid-waste landfill sites, J. Environ. Eng. Landscape Manage., 2009, 17(1), 1–7.
76. M. Adimalla and H. Qian, Groundwater chemistry, distribution and potential health risk appraisal of nitrate enriched groundwater: a case study from the semi-urban region of South India, Ecotoxicol. Environ. Saf., 2021, 207, 111277.
77. A. Rahman, N. C. Mondal and K. K. Tiwari, Anthropogenic Nitrate in Groundwater and Its Health Risks in the View of Background Concentration in a Semi Arid Area of Rajasthan, India, Sci. Rep., 2021, 11, 9279.
78. S. Buvaneshwari, J. Riotte, M. Sekhar, A. K. Sharma, R. Helliwell, M. S. M. Kumar, J. J. Braun and L. Ruiz, Potash fertilizer promotes incipient salinization in groundwater irrigated semi-arid agriculture, Sci. Rep., 2020, 10, 3691.
79. P. J. Sajil Kumar and K. Lemoon, Exposure and health risk assessment of nitrate contamination in groundwater in Coimbatore and Tirupur districts in Tamil Nadu, South India, Environ. Sci. Pollut. Res., 2021, 28, 10248–10261.
80. J. Spoelstra, K. A. Leal, N. D. Senger, S. L. Schiff and R. Post, Isotopic Characterization of Sulfate in a Shallow Aquifer Impacted by Agricultural Fertilizer, Groundwater, 2021, 59, 658–670.
81. M. Zendehbad, P. Cepuder, W. Loiskandl and C. Stumpp, Source identification of nitrate contamination in the urban aquifer of Mashhad, Iran, J. Hydrol., 2019, 25, 100618.
82. M. O. Rivett, S. R. Buss, P. Morgan, J. W. N. Smith and C. D. Bemment, Nitrate attenuation in groundwater: a review of biogeochemical controlling processes, Water Res., 2008, 42(16), 4215–4232.
83. A. Baillieux, C. Moeck, P. Perrochet and D. Hunkeler, Assessing groundwater quality trends in pumping wells using spatially varying transfer functions, Hydrogeol. J., 2015, 23(7), 1449–1463.
84. L. Wang, A. S. Butcher, M. E. Stuart, D. C. Gooddy and J. P. Bloomfield, The nitrate time bomb: a numerical way to investigate nitrate storage and lag time in the unsaturated zone, Environ. Geochem. Health, 2013, 35(5), 667–681.
85. I. Matiatos, C. Moeck, Y. Vystavna, H. Marttila, N. Orlowski, S. Jessen, J. Evaristo, M. Sebilo, G. Koren, E. Dimitriou, S. Müller, Y. Panagopoulos and M. P. Stockinger, Nitrate isotopes in catchment hydrology: insights, ideas and implications for models, J. Hydrol., 2023, 626, 130326.

---

Footnote

† Electronic supplementary information (ESI) available: Tables S1–S12 and Fig. S1 and S2. See DOI: https://doi.org/10.1039/d4em00527a

---